* progress on MUGAST analysis for Commissionning 2019

NPLib/Detectors/Mugast/TMugastPhysics.cxx

@@ -246,16 +246,21 @@ void TMugastPhysics::BuildPhysicalEvent() {
         PosZ.push_back(GetPositionOfInteraction(i).z());
       }
 
-      if (m_Take_E_Y)
+      if (m_Take_E_Y){
         DSSD_E.push_back(DSSD_Y_E);
-      else
+        TotalEnergy.push_back(DSSD_Y_E);
+      }
+      else{
         DSSD_E.push_back(DSSD_X_E);
+        TotalEnergy.push_back(DSSD_X_E);
+      }
 
       if (m_Take_T_Y)
         DSSD_T.push_back(DSSD_Y_T);
       else
         DSSD_T.push_back(DSSD_X_T);
 
+
       for (unsigned int j = 0; j < SecondLayerEMult; ++j) {
         if (m_PreTreatedData->GetSecondLayerEDetectorNbr(j) == N) {
           if (Match_SecondLayer(X, Y, m_PreTreatedData->GetSecondLayerEStripNbr(j))) {
@@ -284,8 +289,8 @@ void TMugastPhysics::BuildPhysicalEvent() {
         SecondLayer_E.push_back(-1000);
         SecondLayer_T.push_back(-1000);
       }
-    }
-  }
+    } // loop on couples
+  } // if (CheckEvent)
   return;
 }
 
@@ -379,15 +384,17 @@ vector<TVector2> TMugastPhysics::Match_X_Y() {
 
 ////////////////////////////////////////////////////////////////////////////
 bool TMugastPhysics::IsValidChannel(const int& Type,
+// Uses raw channel number
     const int& telescope, const int& channel) {
   if (Type == 0){
-    return *(m_XChannelStatus[m_DetectorNumberIndex[telescope] - 1].begin() + channel - 1);
+    return *(m_XChannelStatus[telescope].begin() + channel - 1);
   }
-  else if (Type == 1)
-    return *(m_YChannelStatus[m_DetectorNumberIndex[telescope] - 1].begin() + channel - 1);
+  else if (Type == 1){
+    return *(m_YChannelStatus[telescope].begin() + channel - 1);
+    }
 
   else if (Type == 2)
-    return *(m_SecondLayerChannelStatus[m_DetectorNumberIndex[telescope] - 1].begin() + channel - 1);
+    return *(m_SecondLayerChannelStatus[telescope].begin() + channel - 1);
 
   else
     return false;
@@ -398,25 +405,35 @@ void TMugastPhysics::ReadAnalysisConfig() {
   NPL::InputParser parser("./configs/ConfigMugast.dat");
   vector<NPL::InputBlock*> blocks = parser.GetAllBlocksWithToken("ConfigMugast");
 
+  cout << endl << "//// Read MUGAST analysis configuration" <<endl;
 
   for (unsigned int i = 0; i < blocks.size(); i++) {
     if(blocks[i]->HasToken("MAX_STRIP_MULTIPLICITY"))
       m_MaximumStripMultiplicityAllowed = blocks[i]->GetInt("MAX_STRIP_MULTIPLICITY");
+    
     if(blocks[i]->HasToken("STRIP_ENERGY_MATCHING"))
       m_StripEnergyMatching = blocks[i]->GetDouble("STRIP_ENERGY_MATCHING","MeV");
-    if(blocks[i]->HasToken("DISABLE_CHANNEL")){
-      vector<int> v = blocks[i]->GetVectorInt("DISABLE_CHANNEL");
-      *(m_XChannelStatus[v[0] - 1].begin() + v[1] - 1) = false;
+    
+    if(blocks[i]->HasToken("DISABLE_CHANNEL_X")){
+      vector<int> v = blocks[i]->GetVectorInt("DISABLE_CHANNEL_X");
+      *(m_XChannelStatus[v[0]].begin() + v[1] - 1) = false;
+    }
+    
+    if(blocks[i]->HasToken("DISABLE_CHANNEL_Y")){
+      vector<int> v = blocks[i]->GetVectorInt("DISABLE_CHANNEL_Y");
+      *(m_YChannelStatus[v[0]].begin() + v[1] - 1) = false;
     }
+    
     if(blocks[i]->HasToken("DISABLE_ALL")){
       int telescope = blocks[i]->GetInt("DISABLE_ALL");
       vector<bool> ChannelStatus;
       ChannelStatus.resize(128, false);
-      m_XChannelStatus[telescope - 1] = ChannelStatus;
-      m_YChannelStatus[telescope - 1] = ChannelStatus;
+      m_XChannelStatus[telescope] = ChannelStatus;
+      m_YChannelStatus[telescope] = ChannelStatus;
       ChannelStatus.resize(16, false);
-      m_SecondLayerChannelStatus[telescope - 1]  = ChannelStatus;
+      m_SecondLayerChannelStatus[telescope]  = ChannelStatus;
     }
+
     if (blocks[i]->HasToken("TAKE_E_Y"))
       m_Take_E_Y = blocks[i]->GetInt("TAKE_E_Y");
 
@@ -732,9 +749,20 @@ void TMugastPhysics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1,
   TVector3 Strip_1_1;
 
   //   Geometry Parameter
-  double Face          = 98; // mm
+  double Base,Height;
+  if(DetectorType[m_NumberOfTelescope-1]==MG_TRAPEZE){
+    Base          = 91.48; // mm
+    Height        = 104.688; // mm
+  }
+    
+  if(DetectorType[m_NumberOfTelescope-1]==MG_SQUARE){
+    Base          = 91.716; // mm
+    Height        = 94.916; // mm
+  }
+    //double Face          = 98; // mm
   double NumberOfStrip = 128;
-  double StripPitch    = Face / NumberOfStrip; // mm
+  double StripPitchBase    = Base / NumberOfStrip; // mm
+  double StripPitchHeight  = Height / NumberOfStrip; // mm
   //   Buffer object to fill Position Array
   vector<double> lineX;
   vector<double> lineY;
@@ -745,8 +773,9 @@ void TMugastPhysics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1,
   vector<vector<double>> OneTelescopeStripPositionZ;
 
   //   Moving StripCenter to 1.1 corner:
-  Strip_1_1 = C_X1_Y1 + (U + V) * (StripPitch / 2.);
-  Strip_1_1 += U + V ;
+  //Strip_1_1 = C_X1_Y1 + (U + V) * (StripPitch / 2.);
+  Strip_1_1 = C_X1_Y1 + U  * (StripPitchBase / 2.) + V * (StripPitchHeight / 2.);
+  //Strip_1_1 += U + V ;
 
   for (int i = 0; i < 128; ++i) {
     lineX.clear();
@@ -754,7 +783,8 @@ void TMugastPhysics::AddTelescope(TVector3 C_X1_Y1, TVector3 C_X128_Y1,
     lineZ.clear();
 
     for (int j = 0; j < 128; ++j) {
-      StripCenter = Strip_1_1 + StripPitch * (i * U + j * V);
+      //StripCenter = Strip_1_1 + StripPitch * (i * U + j * V);
+      StripCenter = Strip_1_1 + i*U*StripPitchBase  + j*V*StripPitchHeight;
       lineX.push_back(StripCenter.X());
       lineY.push_back(StripCenter.Y());
       lineZ.push_back(StripCenter.Z());
@@ -778,13 +808,24 @@ void TMugastPhysics::InitializeStandardParameter() {
 
   ChannelStatus.resize(128, true);
   for (int i = 0; i < m_NumberOfTelescope; ++i) {
-    m_XChannelStatus[i] = ChannelStatus;
-    m_YChannelStatus[i] = ChannelStatus;
+    auto it=m_DetectorNumberIndex.begin();
+    for(unsigned int j = 0 ; j < i; j++){
+      it++;
+      }
+    int det= it->first;
+    m_XChannelStatus[det] = ChannelStatus;
+    m_YChannelStatus[det] = ChannelStatus;
   }
 
   ChannelStatus.resize(16, true);
   for (int i = 0; i < m_NumberOfTelescope; ++i) {
-    m_SecondLayerChannelStatus[i]  = ChannelStatus;
+    auto it=m_DetectorNumberIndex.begin();
+    for(unsigned int j = 0 ; j < i; j++){
+      it++;
+      }
+
+    int det= it->first;
+    m_SecondLayerChannelStatus[det]  = ChannelStatus;
   }
 
   m_MaximumStripMultiplicityAllowed = m_NumberOfTelescope;

Projects/MUGAST/56Nidp.reaction

-%%%%%%%%%%%%%%%%%%%%%% S1107 at Triumf %%%%%%%%%%%%%%%%%%%%%%%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-Beam
-  Particle= 56Ni 
-  ExcitationEnergy= 0 MeV
-  SigmaThetaX= 3 mrad
-  SigmaPhiY= 3 mrad
-  SigmaX= 25 mm
-  SigmaY= 25 mm
-  MeanThetaX= 1 mrad
-  MeanPhiY= 1 mrad
-  MeanX= 0 mm
-  MeanY= 0 mm
-  EnergyProfilePath= eLise.txt EL
-  %XThetaXProfilePath=
-  %YPhiYProfilePath=
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-TwoBodyReaction
- Beam= 56Ni
- Target= 2H
- Light= 1H
- Heavy= 57Ni
- ExcitationEnergyLight= 0.0 MeV
- ExcitationEnergyHeavy= 1.0 MeV
- DoubleDifferentialCrossSectionPath= 56Nidp.root h
- ShootLight= 1
- ShootHeavy= 0
- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- 
- 

Projects/MUGAST/Analysis.cxx

@@ -36,17 +36,20 @@ Analysis::~Analysis(){
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::Init() {
+
+  agata_zShift=51*mm;
+
   // initialize input and output branches
   InitOutputBranch();
   InitInputBranch();
-  myInit = new TInitialConditions();
   // get MUST2 and Gaspard objects
   M2 = (TMust2Physics*)  m_DetectorManager -> GetDetector("M2Telescope");
   MG = (TMugastPhysics*) m_DetectorManager -> GetDetector("Mugast");
+  CATS = (TCATSPhysics*) m_DetectorManager->GetDetector("CATSDetector");
 
   // get reaction information
-  myReaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
-  OriginalBeamEnergy = myReaction.GetBeamEnergy();
+  reaction.ReadConfigurationFile(NPOptionManager::getInstance()->GetReactionFile());
+  OriginalBeamEnergy = reaction.GetBeamEnergy();
   // target thickness
   TargetThickness = m_DetectorManager->GetTargetThickness();
   string TargetMaterial = m_DetectorManager->GetTargetMaterial();
@@ -54,9 +57,9 @@ void Analysis::Init() {
   WindowsThickness = 0;//m_DetectorManager->GetWindowsThickness(); 
   string WindowsMaterial = "";//m_DetectorManager->GetWindowsMaterial();
 
- // energy losses
-  string light=NPL::ChangeNameToG4Standard(myReaction.GetNucleus3()->GetName());
-  string beam=NPL::ChangeNameToG4Standard(myReaction.GetNucleus1()->GetName());
+  // energy losses
+  string light=NPL::ChangeNameToG4Standard(reaction.GetNucleus3()->GetName());
+  string beam=NPL::ChangeNameToG4Standard(reaction.GetNucleus1()->GetName());
   LightTarget = NPL::EnergyLoss(light+"_"+TargetMaterial+".G4table","G4Table",100 );
   LightAl = NPL::EnergyLoss(light+"_Al.G4table","G4Table",100);
   LightSi = NPL::EnergyLoss(light+"_Si.G4table","G4Table",100);
@@ -86,7 +89,6 @@ void Analysis::Init() {
   Y = 0;
   Z = 0;
   dE = 0;
-  dTheta=0;
   BeamDirection = TVector3(0,0,1);
 
 }
@@ -95,14 +97,18 @@ void Analysis::Init() {
 void Analysis::TreatEvent() {
   // Reinitiate calculated variable
   ReInitValue();
-  //double zImpact = Rand.Uniform(-TargetThickness*0.5,TargetThickness*0.5);
-  double zImpact = 0 ;
-  BeamImpact = TVector3(0,0,zImpact); 
+  double XTarget = CATS->GetPositionOnTarget().X();
+  double YTarget = CATS->GetPositionOnTarget().Y();
+  TVector3 BeamDirection = CATS->GetBeamDirection();
+
+  BeamImpact = TVector3(XTarget,YTarget,0); 
   // 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));
+  //reaction.SetBeamEnergy(Rand.Gaus(myInit->GetIncidentFinalKineticEnergy(),myInit->GetIncidentFinalKineticEnergy()/235));
 
-  //////////////////////////// LOOP on MUST2 //////////////////
+  ////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////// LOOP on MUST2  ////////////////////////////
+  ////////////////////////////////////////////////////////////////////////////
   for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
     /************************************************/
     //Part 0 : Get the usefull Data
@@ -145,7 +151,7 @@ void Analysis::TreatEvent() {
     // Evaluate energy using the thickness
     ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface);
     // Target Correction
-    ELab   = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2.-zImpact, ThetaNormalTarget);
+    ELab   = LightTarget.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
 
     if(LightWindow)
       ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
@@ -153,24 +159,22 @@ void Analysis::TreatEvent() {
 
     /************************************************/
     // Part 3 : Excitation Energy Calculation
-    Ex = myReaction.ReconstructRelativistic( ELab , ThetaLab );
+    Ex = reaction.ReconstructRelativistic( ELab , ThetaLab );
     ThetaLab=ThetaLab/deg;
 
     /************************************************/
 
     /************************************************/
     // Part 4 : Theta CM Calculation
-    ThetaCM  = myReaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+    ThetaCM  = reaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
     /************************************************/
   }//end loop MUST2
 
   ////////////////////////////////////////////////////////////////////////////
+  //////////////////////////////// LOOP on MUGAST ////////////////////////////
   ////////////////////////////////////////////////////////////////////////////
-  //////////////////////////// LOOP on MUGAST //////////////////
-  
-/*  
-for(unsigned int countMugast = 0 ; countMugast < MG->DSSD_E.size() ; countMugast++){
-  if(MG->GetEnergyDeposit(countMugast)>0){
+  for(unsigned int countMugast = 0 ; countMugast < MG->DSSD_E.size() ; countMugast++){
+
     // Part 1 : Impact Angle
     ThetaMGSurface = 0;
     ThetaNormalTarget = 0;
@@ -189,22 +193,44 @@ for(unsigned int countMugast = 0 ; countMugast < MG->DSSD_E.size() ; countMugast
     Energy = ELab = 0;
     Energy = MG->GetEnergyDeposit(countMugast);
     // Target Correction
-    ELab   = LightTarget.EvaluateInitialEnergy( Energy ,TargetThickness*0.5-zImpact, ThetaNormalTarget);
+    ELab   = LightTarget.EvaluateInitialEnergy( Energy ,TargetThickness*0.5, ThetaNormalTarget);
 
     if(LightWindow)
       ELab = LightWindow->EvaluateInitialEnergy( ELab ,WindowsThickness, ThetaNormalTarget);
-    
+
     // Part 3 : Excitation Energy Calculation
-    Ex = myReaction.ReconstructRelativistic( ELab , ThetaLab );
+    Ex = reaction.ReconstructRelativistic( ELab , ThetaLab );
 
 
     // Part 4 : Theta CM Calculation
-    ThetaCM  = myReaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
+    ThetaCM  = reaction.EnergyLabToThetaCM( ELab , ThetaLab)/deg;
     ThetaLab=ThetaLab/deg;
 
   }//end loop Mugast
-}
-*/
+  
+  ////////////////////////////////////////////////////////////////////////////
+  ///////////////////////////////// LOOP on AGATA ////////////////////////////
+  ////////////////////////////////////////////////////////////////////////////
+  if(nbTrack==1){ // keep only multiplicity one event
+    TLorentzVector GammaLV;
+    // Measured E
+    double Egamma=trackE[0];
+    // Gamma detection position
+    TVector3 GammaHit(trackX1[0],trackY1[0],trackZ1[0]); 
+    // Gamma Direction 
+    TVector3 GammaDirection = GammaHit-BeamImpact;
+    // Beta from Two body kinematic
+    TVector3 beta = reaction.GetEnergyImpulsionLab_4().BoostVector();
+    // Construct LV in lab frame
+    GammaLV.SetPx(Egamma*GammaDirection.X());
+    GammaLV.SetPy(Egamma*GammaDirection.Y());
+    GammaLV.SetPz(Egamma*GammaDirection.Z());
+    GammaLV.SetE(Egamma);
+    // Boost back in CM
+    GammaLV.Boost(-beta);
+    // Get EDC
+    EDC=GammaLV.Energy();
+    }
 
 }
 
@@ -214,6 +240,7 @@ void Analysis::End(){
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitOutputBranch() {
   RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+  RootOutput::getInstance()->GetTree()->Branch("EDC",&Ex,"Ex/D");
   RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
   RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
   RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
@@ -222,20 +249,54 @@ 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("dTheta",&dTheta,"dTheta/D");
+  // Vamos 
+  RootOutput::getInstance()->GetTree()->Branch("LTS",&LTS,"LTS/l");
+
+  // Agata
+  // Time stamp of the agata trigger
+  RootOutput::getInstance()->GetTree()->Branch("TStrack",&TStrack,"TStrack/l");
+
+  // Array of reconstructed tracks
+  RootOutput::getInstance()->GetTree()->Branch("nbTrack",&nbTrack,"nbTrack/I");
+  RootOutput::getInstance()->GetTree()->Branch("trackE",trackE,"trackE[nbTrack]/F");
+  RootOutput::getInstance()->GetTree()->Branch("trackX1",trackX1,"trackX1[nbTrack]/F");
+  RootOutput::getInstance()->GetTree()->Branch("trackY1",trackY1,"trackY1[nbTrack]/F");
+  RootOutput::getInstance()->GetTree()->Branch("trackZ1",trackZ1,"trackZ1[nbTrack]/F");
+  RootOutput::getInstance()->GetTree()->Branch("trackT",trackT,"trackT[nbTrack]/F");
+  RootOutput::getInstance()->GetTree()->Branch("trackCrystalID",trackCrystalID,"trackCrystalID[nbTrack]/I");
+
+  // Array of reconstructed core
+  RootOutput::getInstance()->GetTree()->Branch("nbCores",&nbCores,"nbCores/I");
+  RootOutput::getInstance()->GetTree()->Branch("coreId",coreId,"coreId[nbCores]/I");
+  RootOutput::getInstance()->GetTree()->Branch("coreTS",coreTS,"coreTS[nbCores]/l");
+  RootOutput::getInstance()->GetTree()->Branch("coreE0",coreE0,"coreE0[nbCores]/F");
+  //
 }
 
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitInputBranch(){
-  RootInput::getInstance()->GetChain()->SetBranchStatus("Run",true);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("Run",&Run);
-  RootInput::getInstance()->GetChain()->SetBranchStatus("InitialConditions",true);
-  RootInput::getInstance()->GetChain()->SetBranchAddress("InitialConditions",&myInit);
+  // RootInput:: getInstance()->GetChain()->SetBranchAddress("GATCONF",&vGATCONF);
+  // Vamos
+  RootInput::getInstance()->GetChain()->SetBranchAddress("LTS",&LTS);
+  // Agata
+  RootInput::getInstance()->GetChain()->SetBranchAddress("TStrack",&TStrack);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("nbTrack",&nbTrack);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("trackE",trackE);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("trackX1",trackX1);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("trackY1",trackY1);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("trackZ1",trackZ1);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("trackT",trackT);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("trackCrystalID",trackCrystalID);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("nbCores",&nbCores);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("coreId",coreId);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("coreTS",coreTS);
+  RootInput::getInstance()->GetChain()->SetBranchAddress("coreE0",coreE0);
 }
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::ReInitValue(){
   Ex = -1000 ;
+  EDC= -1000;
   ELab = -1000;
   ThetaLab = -1000;
   ThetaCM = -1000;
@@ -243,7 +304,6 @@ void Analysis::ReInitValue(){
   Y = -1000;
   Z = -1000;
   dE= -1000;
-  dTheta= -1000;
 }
 
 

Projects/MUGAST/Analysis.h

@@ -28,7 +28,7 @@
 #include"RootInput.h"
 #include "TMust2Physics.h"
 #include "TMugastPhysics.h"
-#include "TInitialConditions.h"
+#include "TCATSPhysics.h"
 #include <TRandom3.h>
 #include <TVector3.h>
 #include <TMath.h>
@@ -50,10 +50,11 @@ class Analysis: public NPL::VAnalysis{
  
   private:
   double Ex;
+  double EDC;
   double ELab;
   double ThetaLab;
   double ThetaCM;
-  NPL::Reaction myReaction;
+  NPL::Reaction reaction;
     //	Energy loss table: the G4Table are generated by the simulation
   NPL::EnergyLoss LightTarget;
   NPL::EnergyLoss LightAl;
@@ -85,13 +86,31 @@ class Analysis: public NPL::VAnalysis{
   double X ;
   double Y ;
   double Z ;
- 
+  // Vamos Branches
+  unsigned long long int LTS;
+  
+  // Agata branches
+   double agata_zShift;
+   unsigned long long int TStrack;
+   int nbHits;
+   int nbTrack;
+   float *trackE= new float(100);
+   float *trackX1= new float(100);
+   float *trackY1= new float(100);
+   float *trackZ1= new float(100);
+   float *trackT= new float(100);
+   int *trackCrystalID = new int(100);
+   int nbCores;
+   int *coreId= new int(100);
+   ULong64_t *coreTS= new ULong64_t(100);
+   float *coreE0= new float(100);
+   //
   double dE;
   double dTheta;
   // Branches and detectors
   TMust2Physics* M2;
   TMugastPhysics* MG;
-  TInitialConditions* myInit ;
+  TCATSPhysics* CATS;
 
 };
 #endif

Projects/MUGAST/RunToTreat.txt

 TTreeName 
-	AutoTree
+  TreeMaster	
 RootFileName 
-  /data/mugastX/mugast/acquisition/root/run_0040_0.root 
+  /data/mugastX/com2019/acquisition/MERGED_DATA/MergedTree_run_0107_0000.root

Projects/MUGAST/configs/ConfigMugast.dat

+ConfigMugast
+ TAKE_E_X= 1 
+ %raw channel 86 corresponds to stripY 46
+ DISABLE_CHANNEL_Y= 3 123
+ MAX_STRIP_MULTIPLICITY= 100
+ STRIP_ENERGY_MATCHING= 1 MeV
+ DSSD_X_E_RAW_THRESHOLD= 8250
+ DSSD_Y_E_RAW_THRESHOLD= 8100

Projects/MUGAST/configs/ConfigMust2.dat

 ConfigMust2
-	MAX_STRIP_MULTIPLICITY 10
-	STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA 5
-	STRIP_ENERGY_MATCHING_SIGMA 0.5
-	CSI_SIZE 256
-	SI_X_E_RAW_THRESHOLD 0
-	CSI_E_RAW_THRESHOLD 0
+ MAX_STRIP_MULTIPLICITY 100
+ STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA 50
+ STRIP_ENERGY_MATCHING_SIGMA 0.05
+ DISABLE_CHANNEL MM1STRX12
+ SI_X_E_RAW_THRESHOLD 8200
+ SI_Y_E_RAW_THRESHOLD 8180
+ CSI_E_RAW_THRESHOLD 8250
+ CSI_SIZE 34
+ TAKE_T_X
+ TAKE_E_X

Projects/MUGAST/eLise.txt

-0  0
-613.219 709.855
-614.514 1313.08
-615.809 2238.92
-617.104 3326.08
-618.399 5471.15
-619.694 7846.54
-620.988 10682.6
-622.283 14367.4
-623.578 18476.6
-624.873 24443.9
-626.168 30839.4
-627.463 38091.2
-628.757 46423.8
-630.052 55296.8
-631.347 65853.4
-632.642 76681.1
-633.937 88051
-635.232 99897.7
-636.527 111983
-637.821 123753
-639.116 135331
-640.411 146525
-641.706 157397
-643.001 168109
-644.296 177356
-645.59 186232
-646.885 194365
-648.18 202193
-649.475 209870
-650.77 216265
-652.065 222212
-653.36 227261
-654.654 232041
-655.949 236686
-657.244 241899
-658.539 247179
-659.834 252593
-661.128 257121
-662.423 261204
-663.718 265745
-665.013 270121
-666.308 274168
-667.603 278329
-668.898 282547
-670.192 286326
-671.487 289804
-672.782 292678
-674.077 295509
-675.372 298320
-676.667 300010
-677.961 301688
-679.256 303343
-680.551 303699
-681.846 303407
-683.141 304337
-684.436 304835
-685.73 304469
-687.025 302633
-688.32 300063
-689.615 298270
-690.91 295638
-692.205 291329
-693.5 287682
-694.794 284365
-696.089 279889
-697.384 274682
-698.679 268013
-699.974 261764
-701.269 255725
-702.563 247673
-703.858 239926
-705.153 232786
-706.448 225368
-707.743 217812
-709.038 208636
-710.333 199665
-711.627 191105
-712.922 181974
-714.217 172557
-715.512 164076
-716.807 155603
-718.101 147146
-719.396 138130
-720.691 128834
-721.986 120126
-723.281 111349
-724.576 102434
-725.871 94069.8
-727.165 85981.5
-728.46 78207.3
-729.755 70298.4
-731.05 62120.2
-732.345 54567
-733.64 47326.1
-734.934 40399.1
-736.229 33915.6
-737.524 28319.1
-738.819 23384.8
-740.114 18781.6
-741.409 14615.9
-742.704 10933.7
-743.998 8218.64
-745.293 5802.54
-746.588 3535.92
-747.883 2403.3
-749.178 1431.15
-750.473 779.938
-2000 0

Projects/MUGAST/macros/DrawFunction.C

+void DrawEfficiencyMugast(){
+  PhysicsTree->Draw("acos(PosZ/sqrt(PosX*PosX+PosY*PosY+PosZ*PosZ))*180./3.141592>>h(180,0,180)","",""); 
+  TH1F* h = new TH1F();
+  gDirectory->GetObject("h",h);
+TF1 *fa2 = new TF1("fa2","sin(x*TMath::Pi()/180.)",0,180);
+h->Divide(fa2,1);
+}
+
+void DrawImpactMugast_XY(){
+  PhysicsTree->Draw("Mugast.PosY:Mugast.PosX","","colz"); 
+}
+void DrawImpactMugast_XY(int tel){
+  PhysicsTree->Draw("Mugast.PosY:Mugast.PosX",Form("Mugast.TelescopeNumber==%d",tel),"colz"); 
+}
+void DrawImpactMugast_YZ(){
+  PhysicsTree->Draw("Mugast.PosY:Mugast.PosZ","",""); 
+}
+void DrawImpactMugast_XZ(){
+  PhysicsTree->Draw("Mugast.PosX:Mugast.PosZ","",""); 
+}
+void DrawComparisonThetaMG_RealTheta(){
+  // Theta rec - Mugast alone
+  PhysicsTree->Draw("acos(PosZ/sqrt(PosX*PosX+PosY*PosY+PosZ*PosZ))*180./3.141592>>h(180,0,180)","",""); 
+  // Theta rec - combine MG,M2 and CATS
+  PhysicsTree->Draw("ThetaLab>>h2(180,0,180)","","same"); 
+  TH1F* h2 = new TH1F();
+  gDirectory->GetObject("h2",h2);
+  h2->SetLineColor(kRed);
+
+}
+
+void DrawPhiAnnular(){
+  //PhysicsTree->Draw("atan(Mugast.PosY/Mugast.PosX)*180./3.141592>>h(180,0,180)","",""); 
+  PhysicsTree->Draw("Mugast.PosY/Mugast.PosX","Mugast.TelescopeNumber==11",""); 
+}

Projects/MUGAST/macros/Kine.cxx

+#include"NPReaction.h"
+
+TCutG* ETOF=NULL;
+TCutG* EDE=NULL;
+TChain* chain=NULL ;
+TCanvas* c1 = NULL;
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadCuts(){
+TFile* File_ETOF = new TFile("cuts/ETOF.root","READ");
+ETOF = (TCutG*) File_ETOF->FindObjectAny("ETOF");
+
+TFile* File_EDE = new TFile("cuts/EDE.root","READ");
+EDE= (TCutG*) File_EDE->FindObjectAny("EDE");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void LoadChain(){
+chain = new TChain("PhysicsTree");
+chain->Add("../../Outputs/Analysis/Example1.root");
+}
+////////////////////////////////////////////////////////////////////////////////
+void LoadEventList(){
+}
+////////////////////////////////////////////////////////////////////////////////
+void plot_kine(NPL::Reaction r, double Ex,Color_t c,int w, int s){
+  r.SetExcitation4(Ex);
+
+  TGraph* g= r.GetKinematicLine3();
+  g->SetLineColor(c) ;
+  g->SetLineStyle(s) ;
+  g->SetLineWidth(w) ;
+
+  g->Draw("c");
+  }
+////////////////////////////////////////////////////////////////////////////////
+void plot_state(double Ex,double max,Color_t c,int w, int s){
+  TLine* line = new TLine(Ex,0,Ex,max) ; 
+  line->SetLineColor(c) ;
+  line->SetLineStyle(s) ;
+  line->SetLineWidth(w) ;
+  line->Draw();
+
+  }
+////////////////////////////////////////////////////////////////////////////////
+void Kine(){
+// Load stuff
+LoadChain();
+LoadEventList();
+//LoadCuts();
+
+// for Kinematic calculation
+NPL::Reaction O17("16O(d,p)17O@96");
+NPL::Reaction O16dd("16O(d,d)16O@96");
+NPL::Reaction O16pp("16O(p,p)16O@96");
+NPL::Reaction O16C12("16O(12C,12C)16O@96");
+NPL::Reaction O15("16O(d,t)15O@96");
+Color_t color_dp=kOrange+7;
+Color_t color_dd=kRed;
+Color_t color_pp=kGreen;
+Color_t color_dt=kAzure+7;
+Color_t color_ca=kBlack;
+
+c1 = new TCanvas("Example1","Example1",0,0,600,600);
+c1->Divide(2,2);
+// Kinematic Plot
+c1->cd(1);
+TH2F* hKine = new TH2F("hKine","hKine",1000,0,180,1000,0,40);
+hKine->Draw("col");
+hKine->GetXaxis()->SetTitle("#Theta_{Lab} (deg)");
+hKine->GetYaxis()->SetTitle("E_{Lab} (MeV)");
+chain->Draw("ELab:ThetaLab","","col");
+// O17
+plot_kine(O17,0,color_dp,2,1);
+plot_kine(O17,0.87073,color_dp,2,1);
+plot_kine(O17,3.05536,color_dp,1,1);
+// Sn
+plot_kine(O17,4.143,color_dp,1,2);
+
+// O15
+plot_kine(O15,0,color_dt,2,1);
+//plot_kine(O15,5.183,color_dt,2,1);
+//plot_kine(O15,5.241,color_dt,1,1);
+//plot_kine(O15,6.176,color_dt,1,1);
+// Sp
+//plot_kine(O15,7.297,color_dt,1,2);
+// Elastic
+plot_kine(O16dd,0,color_dd,2,1);
+plot_kine(O16pp,0,color_pp,1,1);
+plot_kine(O16C12,0,color_ca,1,2);
+
+// Kinematic Plot
+c1->cd(2);
+int bin = 1000;
+double start=-10;
+double end =  10;
+double binning = (end-start)/bin;
+TH1F* hEx = new TH1F("hEx","hEx",bin,start,end);
+hEx->FillRandom("gaus", 10000);
+hEx->GetXaxis()->SetTitle("E_{17O} (MeV)");
+hEx->GetYaxis()->SetTitle(Form("counts / %.0f keV",binning*1000));
+
+hEx->Draw();
+double max = hEx->GetMaximum();
+plot_state(0,max,color_dp,2,1);
+plot_state(0.87073,max,color_dp,2,1);
+plot_state(3.05536,max,color_dp,1,1);
+plot_state(4.143,max,color_dp,1,2);
+
+
+
+}

Projects/MUGAST/mugast.detector

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Target
+ THICKNESS= 3.2 micrometer
+ ANGLE= 20 deg
+ RADIUS=	24 mm
+ MATERIAL= CD2
+ X= 0 mm
+ Y= 0 mm
+ Z= 0 mm
+
+%%%%%%% Telescope 1 %%%%%%%
+M2Telescope
+ X1_Y1=	     10.85	105.03	162.16 mm
+ X1_Y128=	 22.80	  9.84	191.94 mm
+ X128_Y1=	104.09	105.03	124.76 mm
+ X128_Y128=	116.04	  9.84	154.54 mm
+ SI=	1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 2 %%%%%%%
+M2Telescope
+ X1_Y1=	    -116.04   9.84	154.54 mm
+ X1_Y128=	 -22.80	  9.84	191.94 mm
+ X128_Y1=	-104.09	105.03	124.76 mm
+ X128_Y128=	 -10.85	105.03	162.16 mm
+ SI=	1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 3 %%%%%%%
+M2Telescope
+ X1_Y1=	    -10.85	-105.03	162.16 mm
+ X1_Y128=	-22.80	 -9.84	191.94 mm
+ X128_Y1=	-104.09	-105.03	124.76 mm
+ X128_Y128=	-116.04	 -9.84	154.54 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%% Telescope 4 %%%%%%%
+M2Telescope
+ X1_Y1=  	116.04	  -9.84	154.54 mm
+ X1_Y128=	22.8	  -9.84	191.94 mm
+ X128_Y1=	104.09	-105.03	124.76 mm
+ X128_Y128=	10.85	-105.03	162.16 mm
+ SI= 1.00
+ SILI= 0.00
+ CSI= 1.00
+ VIS= all
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid 
+ DetectorNumber= 1
+ X128_Y128=       42.4    22.21   -98.74 mm
+ X1_Y128=         24.58   40.02   -98.74 mm
+ X128_Y1=	  122.38  55.31   -31.63 mm
+ X1_Y1=		  57.69  120.01   -31.63 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid 
+ DetectorNumber= 3
+ X128_Y128=       -42.4  -22.21   -98.74 mm
+ X1_Y128=        -24.58  -40.02   -98.74 mm
+ X128_Y1=	-122.38  -55.31   -31.63 mm
+ X1_Y1=		 -57.69 -120.01   -31.63 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid      
+ DetectorNumber= 4
+ X128_Y128=     -22.21     42.4   -98.74 mm
+ X1_Y128=       -40.02    24.58   -98.74 mm
+ X128_Y1=	-55.31   122.38   -31.63 mm
+ X1_Y1=	    	-120.01   57.69   -31.63 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid      
+ DetectorNumber= 5
+ X128_Y128=      22.21     -42.4  -98.74 mm
+ X1_Y128=        40.02    -24.58  -98.74 mm
+ X128_Y1=	 53.91   -121.57  -29.29 mm
+ X1_Y1=	        120.01    -57.69  -31.63 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Trapezoid      
+ DetectorNumber= 7
+ X128_Y128=    -14.28    -45.68   -98.74 mm
+ X1_Y128=       10.92    -45.68   -98.74 mm
+ X128_Y1=      -47.42   -125.65   -31.63 mm
+ X1_Y1=	        44.06   -125.65   -31.63 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Square      
+ DetectorNumber= 9 
+ X1_Y1=      143.81   8.21   91.70 mm
+ X128_Y1=    107.49  95.88   91.70 mm
+ X1_Y128=    143.81   8.21     0.0 mm
+ X128_Y128=  107.49  95.88     0.0 mm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Square 
+ DetectorNumber= 10
+ X1_Y1=     143.81   -8.21  91.7 mm
+ X128_Y1=   107.49  -95.88  91.7 mm
+ X1_Y128=   143.81   -8.21   0.0 mm
+ X128_Y128= 107.49  -95.88   0.0 mm
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+Mugast Annular  
+ DetectorNumber= 11
+ Center=	0 0 -125 mm
+
+ModularLeaf 
+ DefaultValue= -5000
+ Leafs= GATCONFMASTER GATCONFSLAVE GATCONFSLAVE2 ADC1_9 CORREL_INVAMOS ADC_CATS_0
+ X= ADC1_9              ADC1_9               ADC_CATS_0
+ Y= CORREL_INVAMOS    ADC_CATS_0           CORREL_INVAMOS
+