* Progress on Fixing Gaspard trapezoid by angle definition

NPLib/Detectors/GASPARD/GaspardTrackerTrapezoid.cxx

@@ -30,10 +30,13 @@
 #include <cmath>
 #include <stdlib.h>
 
+// NPLib
+#include "NPSystemOfUnits.h"
 // Gaspard
 #include "TGaspardTrackerPhysics.h"
 
-
+//Root
+#include"TRotation.h"
 GaspardTrackerTrapezoid::GaspardTrackerTrapezoid(map<int, GaspardTrackerModule*> &Module,
                                                  TGaspardTrackerPhysics* EventPhysics) 
    : m_ModuleTest(Module),
@@ -87,7 +90,7 @@ void GaspardTrackerTrapezoid::ReadConfiguration(NPL::InputParser parser){
         double Phi = blocks[i]->GetDouble("PHI","deg");
         vector<double> beta = blocks[i]->GetVectorDouble("BETA","deg");
       
-        AddModule(R,Theta,Phi,beta[0],beta[1],beta[2]);
+        AddModule(Theta,Phi,R,beta[0],beta[1],beta[2]);
         m_ModuleTest[m_index["Trapezoid"] + m_NumberOfModule] = this;
       }
     }
@@ -289,10 +292,6 @@ void GaspardTrackerTrapezoid::AddModule(double theta,
 {
    m_NumberOfModule++;
 
-   // convert from degree to radian:
-   theta *= M_PI/180.;
-   phi   *= M_PI/180.;
-
    // Vector U on Module Face (paralelle to Y Strip) (NB: remember that Y strip are allong X axis)
    TVector3 U ;
    // Vector V on Module Face (parallele to X Strip)
@@ -312,21 +311,21 @@ void GaspardTrackerTrapezoid::AddModule(double theta,
 
    W = C.Unit();
    V = W.Cross(YperpW);
-   U = W.Cross(U);
+   U = V.Cross(W);
 
    U = U.Unit();
    V = V.Unit();
 
-   U.Rotate( beta_u * M_PI/180. , U ) ;
-   V.Rotate( beta_u * M_PI/180. , U ) ;
-
-   U.Rotate( beta_v * M_PI/180. , V ) ;
-   V.Rotate( beta_v * M_PI/180. , V ) ;
-
-   U.Rotate( beta_w * M_PI/180. , W ) ;
-   V.Rotate( beta_w * M_PI/180. , W ) ;
+  C=W*10*NPUNITS::mm+C;
 
+   TRotation R;
+   R.Rotate(beta_u,U);
+   R.Rotate(beta_v,V);
+   R.Rotate(beta_w,W); 
 
+   //W=R*W;
+   V=R*V;
+   U=R*U;
    vector<double> lineX;
    vector<double> lineY;
    vector<double> lineZ;
@@ -338,19 +337,20 @@ void GaspardTrackerTrapezoid::AddModule(double theta,
    double X, Y, Z;
 
    // Moving C to the 1.1 corner:
-   C.SetX( C.X() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.X() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.X() )   ;
-   C.SetY( C.Y() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Y() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Y() )   ;
-   C.SetZ( C.Z() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Z() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Z() )   ;
-
+   TVector3 Strip_1_1;
+   Strip_1_1.SetX( C.X() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.X() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.X() )   ;
+   Strip_1_1.SetY( C.Y() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Y() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Y() )   ;
+   Strip_1_1.SetZ( C.Z() - ( m_FirstStageBaseLarge/2 - m_StripPitchX/2 ) * U.Z() - (m_FirstStageHeight/2 - m_StripPitchY/2 ) *V.Z() )   ;
+  
    for (int i = 0; i < m_NumberOfStripsX; i++) {
       lineX.clear();
       lineY.clear();
       lineZ.clear();
 
       for (int j = 0; j < m_NumberOfStripsY; j++) {
-         X = C.X() + m_StripPitchX*U.X()*i + m_StripPitchY*V.X()*j ;
-         Y = C.Y() + m_StripPitchX*U.Y()*i + m_StripPitchY*V.Y()*j ;
-         Z = C.Z() + m_StripPitchX*U.Z()*i + m_StripPitchY*V.Z()*j ;
+         X = Strip_1_1.X() + m_StripPitchX*U.X()*i + m_StripPitchY*V.X()*j ;
+         Y = Strip_1_1.Y() + m_StripPitchX*U.Y()*i + m_StripPitchY*V.Y()*j ;
+         Z = Strip_1_1.Z() + m_StripPitchX*U.Z()*i + m_StripPitchY*V.Z()*j ;
 
          lineX.push_back(X);
          lineY.push_back(Y);

NPSimulation/Detectors/GASPARD/GaspardTracker.cc

@@ -199,15 +199,9 @@ void GaspardTracker::InitializeScorers()
 // Called at in the EventAction::EndOfEventAction
 void GaspardTracker::ReadSensitive(const G4Event* event)
 {
-  // Before looping on each sub-detector, clear the static variable
-  // ms_InterCoord
-  // This is done on the first element of the m_Modules vector.
-  // This should be done here since this variable (of type TIneractionCoordinates)
-  // deals with multiplicity of events > 1.
-  m_Modules[0]->GetInterCoordPointer()->Clear();
-
-  // We do the same for the static variable ms_Event
-  m_Modules[0]->GetEventPointer()->Clear();
+  // Before looping on each sub-detector, clear the static variable ms_Event
+  if(m_Modules.size()>0)
+    m_Modules[0]->GetEventPointer()->Clear();
 
   // loop on sub-detectors belonging to GaspardTracker
   int nbDetectors = m_Modules.size();

Projects/BeLise/Analysis.cxx

@@ -78,7 +78,8 @@ void Analysis::TreatEvent(){
   // Beam energy is measured using F3 and F2 plastic TOF
   double BeamEnergy = Rand.Gaus(Initial->GetIncidentInitialKineticEnergy(),4.5);
   //BeamEnergy = Li11CD2.Slow(BeamEnergy,TargetThickness/2.,0);
-
+  OriginalThetaLab = Initial->GetThetaCM(0);
+  OriginalELab =  Initial->GetKineticEnergy(0);
   myReaction->SetBeamEnergy(BeamEnergy);
   //////////////////////////// LOOP on MUST2 Hit //////////////////
     for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
@@ -174,6 +175,10 @@ void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
   RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
   RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
+  RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
+ RootOutput::getInstance()->GetTree()->Branch("OriginalELab",&OriginalELab,"OriginalELab/D");
+
+
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -189,6 +194,8 @@ void Analysis::ReInitValue(){
   ELab = -1000;
   ThetaLab = -1000;
   ThetaCM = -1000;
+  OriginalThetaLab = -1000;
+  OriginalELab = -1000;
 }
 
 ////////////////////////////////////////////////////////////////////////////////

Projects/BeLise/Analysis.h

@@ -64,7 +64,8 @@ class Analysis: public NPL::VAnalysis{
     double Si_Y_M2;
     double ZTarget;
     double TargetThickness;
-
+    double OriginalThetaLab;
+    double OriginalELab;
     NPL::EnergyLoss He3CD2  ;
     NPL::EnergyLoss He3Al   ;
     NPL::EnergyLoss He3Si   ;

Projects/MUGAST_LISE/Analysis.cxx

@@ -102,7 +102,8 @@ void Analysis::TreatEvent() {
   // determine beam energy for a randomized interaction point in target
   // 1% FWHM randominastion (E/100)/2.35
   myReaction.SetBeamEnergy(Rand.Gaus(myInit->GetIncidentFinalKineticEnergy(),myInit->GetIncidentFinalKineticEnergy()/235));
-
+  OriginalThetaLab = myInit->GetThetaLab_WorldFrame(0);
+  OriginalELab = myInit->GetKineticEnergy(0);
 
   //////////////////////////// LOOP on MUST2 //////////////////
   for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
@@ -230,6 +231,10 @@ void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("Y",&Y,"Y/D");
   RootOutput::getInstance()->GetTree()->Branch("Z",&Z,"Z/D");
   RootOutput::getInstance()->GetTree()->Branch("dE",&dE,"dE/D");
+  RootOutput::getInstance()->GetTree()->Branch("OriginalThetaLab",&OriginalThetaLab,"OriginalThetaLab/D");
+  RootOutput::getInstance()->GetTree()->Branch("OriginalELab",&OriginalELab,"OriginalELab/D");
+
+
   RootOutput::getInstance()->GetTree()->Branch("dTheta",&dTheta,"dTheta/D");
 }
 
@@ -252,6 +257,8 @@ void Analysis::ReInitValue(){
   Z = -1000;
   dE= -1000;
   dTheta= -1000;
+  OriginalThetaLab = -1000;
+  OriginalELab = -1000;
 }
 
 

Projects/MUGAST_LISE/Analysis.h

@@ -79,7 +79,8 @@ class Analysis: public NPL::VAnalysis{
   double Si_E_M2 ;
   double CsI_E_M2  ;
   double Energy ;
-  
+  double OriginalThetaLab;
+  double OriginalELab;  
   double ThetaGDSurface ;
   double X ;
   double Y ;

Projects/MUGAST_LISE/DetectorConfiguration/MUGAST_LISE.detector

@@ -22,7 +22,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 45 deg
@@ -32,7 +32,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 90 deg
@@ -42,7 +42,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 135 deg
@@ -52,7 +52,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 180 deg
@@ -62,7 +62,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 225 deg
@@ -72,7 +72,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 270 deg
@@ -82,7 +82,7 @@ GaspardTracker Trapezoid
  THIRDSTAGE= 0
  VIS= all
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%GaspardTracker Trapezoid
+GaspardTracker Trapezoid
  R= -90 mm
  THETA= 60 deg
  PHI= 315 deg
@@ -99,5 +99,18 @@ Tiara Barrel
  X= 0 mm
  Y= 0 mm
  Z= 0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X1_Y1=   115.12	-11.35	1540.94 mm
+ X1_Y128= 24.33	-11.60	1900.81 mm
+ X128_Y1= 103.64	-103.42	1206.51 mm
+ X128_Y128= 12.89	-103.90	1602.29 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
 
 

Projects/MUGAST_LISE/show.cxx

+void show(){
+
+  TFile* file = new TFile("../../Outputs/Analysis/PhysicsTree.root","READ");
+  TTree* tree = (TTree*) file->FindObjectAny("PhysicsTree");
+  
+  tree->Draw("ELab:ThetaLab>>h(360,0,180,360,0,20)","ELab>0","colz");
+  NPL::Reaction r("68Ni(d,p)69Ni@680");
+  r.GetKinematicLine3()->Draw("c");
+  TCanvas*c =new TCanvas();
+  c->Divide(2,1);
+  c->cd(1);
+  tree->Draw("ELab:OriginalELab>>hE(500,0,16,500,0,16)","ELab>0","colz");
+  TLine* lineE = new TLine(0,0,16,16);
+  lineE->Draw();
+  c->cd(2);
+  tree->Draw("ThetaLab:OriginalThetaLab>>hT(500,100,180,500,100,180)","ELab>0","colz");
+  TLine* lineT = new TLine(100,100,180,180);
+  lineT->Draw();
+}