* Adding MUGAST analysis

NPAnalysis/MUGAST/Analysis.cxx

+#include "Analysis.h"
+
+int main(int argc, char** argv){
+  // command line parsing
+  NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv);
+
+  // Instantiate RootInput
+  string runToReadfileName = myOptionManager->GetRunToReadFile();
+  RootInput:: getInstance("RunToTreat.txt");
+  TChain* Chain = RootInput:: getInstance()->GetChain();
+  // if input files are not given, use those from TAsciiFile
+  if (myOptionManager->IsDefault("DetectorConfiguration")) {
+    string name = RootInput::getInstance()->DumpAsciiFile("DetectorConfiguration");
+    myOptionManager->SetDetectorFile(name);
+  }
+
+  if (myOptionManager->IsDefault("EventGenerator")) {
+    string name = RootInput::getInstance()->DumpAsciiFile("EventGenerator");
+    myOptionManager->SetReactionFile(name);
+  }
+
+  // get input files from NPOptionManager
+  string detectorfileName    = myOptionManager->GetDetectorFile();
+  string OutputfileName      = myOptionManager->GetOutputFile();
+
+  // Instantiate RootOutput
+  RootOutput::getInstance("Analysis/"+OutputfileName, "ResultTree");
+  // RootOutput::getInstance()->GetFile()->SetCompressionLevel(0);
+  // Instantiate the detector using a file
+  NPA::DetectorManager* myDetector = new DetectorManager();
+  myDetector->ReadConfigurationFile(detectorfileName);
+  // Attach new branch
+  InitOutputBranch();
+  InitInputBranch();
+
+  //	Instantiate the Reaction
+  NPL::Reaction*  myReaction = new Reaction ;
+  myReaction	->	ReadConfigurationFile(myOptionManager->GetReactionFile());
+
+  ////////////////////////////////////////////////////////
+  // Get pointer to the different detector
+
+  TMust2Physics* M2  = (TMust2Physics*) myDetector -> GetDetector("MUST2");
+  GaspardTracker* GD = (GaspardTracker*)  myDetector -> GetDetector("GASPARD");
+
+  // Beam Energy
+  double BeamEnergy = 4.0* 74; // AMEV
+
+  // intermediate variable
+  TRandom3 Rand = TRandom3();
+  int DetectorNumber = 0 ;
+  double ThetaNormalTarget = 0 ;
+  double ThetaM2Surface = 0; 
+  double X_M2 = 0 ;
+  double Y_M2 = 0 ;
+  double Z_M2 = 0 ;
+  double Si_E_M2 = 0 ;
+  double CsI_E_M2 = 0 ; 
+  double Energy = 0;
+  double E_M2 = 0;
+  double Si_X_M2 = 0;
+  double Si_Y_M2 = 0;
+  double ZTarget = 0;
+  double TargetThickness = 9*micrometer;
+
+  double ThetaGDSurface = 0; 
+  double X_GD = 0 ;
+  double Y_GD = 0 ;
+  double Z_GD = 0 ;
+  double Si_E_GD = 0 ;
+  double E_GD = 0;
+  double Si_X_GD = 0;
+  double Si_Y_GD = 0;
+ 
+  
+  // Get number of events to treat
+  cout << endl << "///////// Starting Analysis ///////// "<< endl;
+  int nentries = Chain->GetEntries();
+  cout << " Number of Event to be treated : " << nentries << endl;
+  clock_t begin = clock();
+  clock_t end = begin;
+  cout.precision(5);
+  //////////////////////////////////////////////////////////////////////////////
+  // main loop on entries //
+  for (int i = 0 ; i < nentries; i++) {
+    if (i%10000 == 0 && i!=0)  {
+      end = clock();
+      long double TimeElapsed = (long double) (end-begin) / CLOCKS_PER_SEC;
+      double percent = (double)i/nentries;
+      double TimeToWait = (TimeElapsed/percent) - TimeElapsed;       
+      cout  << "                                                      "<< flush;
+      cout  << "\r Progression:" 
+        << percent*100 << " % \t | \t Remaining time : ~" 
+        <<  TimeToWait <<"s |  Analysis Rate : "
+        << (double) i/TimeElapsed << flush;
+    }
+    else if (i == nentries-1)  cout << "\r Progression:" << " 100% " <<endl;
+    // Get the raw Data
+    Chain -> GetEntry(i);
+    // Clear previous Event
+    myDetector->ClearEventPhysics();
+    // Build the current event
+    myDetector->BuildPhysicalEvent();
+    // Reinitiate calculated variable
+    ReInitValue();
+    
+    // Get the Init information on beam position and energy
+    // and apply by hand the experimental resolution
+    // This is because the beam diagnosis are not simulated
+    // PPAC position resolution on target is assumed to be 1mm
+    double XTarget = 0;
+    double YTarget = 0;
+    TVector3 BeamDirection = TVector3(0,0,1);
+
+    // Beam energy is measured using F3 and F2 plastic TOF
+    BeamEnergy = BeamCD2.Slow(BeamEnergy,TargetThickness/2.,0);
+    myReaction->SetBeamEnergy(BeamEnergy);
+
+    //////////////////////////// LOOP on MUST2 //////////////////
+    for(unsigned int countMust2 = 0 ; countMust2 < M2->Si_E.size() ; countMust2++){
+      /************************************************/
+      //Part 0 : Get the usefull Data
+      // MUST2
+      int X = M2->Si_X[countMust2];
+      int Y = M2->Si_Y[countMust2];
+      int TelescopeNumber = M2->TelescopeNumber[countMust2];
+      Si_X_M2 = X ;
+      Si_Y_M2 = Y ;
+
+      /************************************************/
+      // Part 1 : Impact Angle
+      ThetaM2Surface = 0;
+      ThetaNormalTarget = 0;
+      if(XTarget>-1000 && YTarget>-1000){
+        TVector3 BeamImpact(XTarget,YTarget,0);
+        TVector3 HitDirection = M2 -> GetPositionOfInteraction(countMust2) - BeamImpact ;
+        ThetaLab = HitDirection.Angle( BeamDirection );
+
+        ThetaM2Surface = HitDirection.Angle(- M2 -> GetTelescopeNormal(countMust2) );
+        ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+        X_M2 = M2 -> GetPositionOfInteraction(countMust2).X() ;
+        Y_M2 = M2 -> GetPositionOfInteraction(countMust2).Y() ;
+        Z_M2 = M2 -> GetPositionOfInteraction(countMust2).Z() ;
+      }
+
+      else{
+        BeamDirection = TVector3(-1000,-1000,-1000);
+        ThetaM2Surface    = -1000  ;
+        ThetaNormalTarget = -1000  ;
+      }
+
+      /************************************************/
+
+      /************************************************/
+      // Part 2 : Impact Energy
+      Energy = ELab = 0;
+      Si_E_M2 = M2->Si_E[countMust2];
+      CsI_E_M2= M2->CsI_E[countMust2];
+
+      // if CsI
+      if(CsI_E_M2>0 ){
+        // The energy in CsI is calculate form dE/dx Table because 
+        // 20um resolution is poor
+        Energy = 
+          LightSi.EvaluateEnergyFromDeltaE(Si_E_M2,300*micrometer,
+              ThetaM2Surface, 0.01*MeV, 
+              450.*MeV,0.001*MeV ,1000);
+        E_M2=CsI_E_M2;
+      }
+
+      else
+        Energy = Si_E_M2;
+
+      E_M2 += Si_E_M2;
+
+      // Evaluate energy using the thickness 
+      ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaM2Surface); 
+      // Target Correction
+      ELab   = LightCD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+      /************************************************/
+
+      /************************************************/
+      // Part 3 : Excitation Energy Calculation
+      Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
+      ThetaLab=ThetaLab/deg;
+      /************************************************/
+
+      /************************************************/
+      // Part 4 : Theta CM Calculation
+      ThetaCM  = myReaction -> EnergyLabToThetaCM( ELab , 0)/deg;
+      /************************************************/
+    }//end loop MUST2
+    
+    ////////////////////////////////////////////////////////////////////////////
+    ////////////////////////////////////////////////////////////////////////////
+    //////////////////////////// LOOP on GASPARD //////////////////
+    if(GD->GetEnergyDeposit()>0){
+      /************************************************/
+      // Part 1 : Impact Angle
+      ThetaGDSurface = 0;
+      ThetaNormalTarget = 0;
+      if(XTarget>-1000 && YTarget>-1000){
+        TVector3 BeamImpact(XTarget,YTarget,0);
+        TVector3 HitDirection = GD -> GetPositionOfInteraction() - BeamImpact ;
+        ThetaLab = HitDirection.Angle( BeamDirection );
+
+        ThetaGDSurface = HitDirection.Angle( TVector3(0,0,1) ) ;
+        ThetaNormalTarget = HitDirection.Angle( TVector3(0,0,1) ) ;
+        X_GD = GD -> GetPositionOfInteraction().X() ;
+        Y_GD = GD -> GetPositionOfInteraction().Y() ;
+        Z_GD = GD -> GetPositionOfInteraction().Z() ;
+      }
+
+      else{
+        BeamDirection = TVector3(-1000,-1000,-1000);
+        ThetaGDSurface    = -1000  ;
+        ThetaNormalTarget = -1000  ;
+      }
+
+      /************************************************/
+
+      /************************************************/
+      // Part 2 : Impact Energy
+      Energy = ELab = 0;
+      Energy = GD->GetEnergyDeposit();
+      
+      ELab = LightAl.EvaluateInitialEnergy( Energy ,0.4*micrometer , ThetaGDSurface); 
+      // Target Correction
+      ELab   = LightCD2.EvaluateInitialEnergy( ELab ,TargetThickness/2., ThetaNormalTarget);
+      /************************************************/
+
+      /************************************************/
+      // Part 3 : Excitation Energy Calculation
+      Ex = myReaction -> ReconstructRelativistic( ELab , ThetaLab );
+      if(Ex!=Ex){
+        cout << ELab << " " << ThetaLab << " " << Ex<< endl;
+      }
+      ThetaLab=ThetaLab/deg;
+      /************************************************/
+
+      /************************************************/
+      // Part 4 : Theta CM Calculation
+      ThetaCM  = myReaction -> EnergyLabToThetaCM( ELab , 0)/deg;
+      /************************************************/
+    }//end loop GASPARD
+
+    if(ELab>0)
+      RootOutput::getInstance()->GetTree()->Fill();
+   }// loop over events
+
+  cout << "A total of " << nentries << " event has been annalysed " << endl ;
+
+  RootOutput::getInstance()->Destroy();
+  RootInput::getInstance()->Destroy();
+  NPOptionManager::getInstance()->Destroy();
+  /////////////////////////////////////////////////////////////////////////////
+  return 0 ;
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void InitOutputBranch() {
+  RootOutput::getInstance()->GetTree()->Branch("Ex",&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");
+}
+
+////////////////////////////////////////////////////////////////////////////////
+void InitInputBranch(){
+  RootInput:: getInstance()->GetChain()->SetBranchAddress("InitialConditions",&Init );
+  RootInput:: getInstance()->GetChain()->SetBranchStatus("InitialConditions",true );
+  RootInput:: getInstance()->GetChain()->SetBranchStatus("fIC_*",true );
+}
+////////////////////////////////////////////////////////////////////////////////     
+void ReInitValue(){
+  Ex = -1000 ;
+  ELab = -1000;
+  ThetaLab = -1000;
+  ThetaCM = -1000;
+}
+
+
+////////////////////////////////////////////////////////////////////////////////

NPAnalysis/MUGAST/Analysis.h

+// You can use this file to declare your spectra, file, energy loss , ... and whatever you want.
+// This way you can remove all unnecessary declaration in the main programm.
+// In order to help debugging and organizing we use Name Space.
+
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -------------------------------------- VARIOUS INCLUDE ---------------------------------------
+
+// NPL
+#include "DetectorManager.h"
+#include "NPOptionManager.h"
+#include "NPReaction.h"
+#include "RootInput.h"
+#include "RootOutput.h"
+#include "TMust2Physics.h"
+#include "GaspardTracker.h"
+
+#include "TInitialConditions.h"
+#include "NPEnergyLoss.h"
+using namespace NPL ;
+// STL C++
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include <string>
+#include <cmath>
+#include <cstdlib>
+using namespace std;
+// ROOT
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TVector3.h>
+#include <TRandom3.h>
+#include <TMath.h>
+#include <TObject.h>
+
+// ----------------------------------------------------------------------------------------------
+void InitOutputBranch() ;
+void InitInputBranch() ;
+void ReInitValue() ;
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// ----------------------------------- DOUBLE, INT, BOOL AND MORE -------------------------------
+namespace VARIABLE{
+  double Ex;
+  double ELab;
+  double ThetaLab;
+  double ThetaCM;
+  TInitialConditions* Init = new TInitialConditions();
+}
+
+using namespace VARIABLE ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+// -----------------------------------ENERGY LOSS----------------------------------------------
+namespace ENERGYLOSS{
+  //	Energy loss table: the G4Table are generated by the simulation
+  EnergyLoss LightCD2 = EnergyLoss("proton_CD2.G4table","G4Table",100 );
+  EnergyLoss LightAl = EnergyLoss("proton_Al.G4table","G4Table",10);
+  EnergyLoss LightSi = EnergyLoss("proton_Si.G4table","G4Table",10);
+  EnergyLoss BeamCD2 = EnergyLoss("Kr74[0.0]_CD2.G4table","G4Table",100);
+}
+
+using namespace ENERGYLOSS ;
+// ----------------------------------------------------------------------------------------------
+/////////////////////////////////////////////////////////////////////////////////////////////////
+
+

NPAnalysis/MUGAST/Makefile

+# include same architecture file than for NPLib
+# so that consistency is ensured
+include $(NPLIB)/Makefile.arch
+
+# additional libraries
+LIBRARY  = `$(NPLIB)/liblist`
+LIBRARY += -L$(CLHEP_LIB_DIR) -l$(CLHEP_LIB)
+
+PROGRAMS	= NPAnalysis
+
+all:	$(PROGRAMS)
+
+NPAnalysis:	Analysis.o
+	$(LD) $(LDFLAGS) $^ $(LIBS) $(LIBRARY) $(OutPutOpt) $@
+	@echo "$@ done"
+	
+
+# rule for creating .o from .cxx
+.SUFFIXES: .$(SrcSuf)
+.$(SrcSuf).$(ObjSuf):
+	$(CXX) $(CXXFLAGS) $(INCLUDE) -c $<
+
+# some cleaning
+clean:
+	rm -rf *.o
+
+distclean:
+	make clean; rm $(PROGRAMS)
+
+# dependences
+Analysis.o:	Analysis.cxx	Analysis.h

NPAnalysis/MUGAST/RunToTreat.txt

+TTreeName 
+	SimulatedTree
+RootFileName 
+	../../Outputs/Simulation/myResult.root
+

NPAnalysis/MUGAST/configs/ConfigMust2.dat

+ConfigMust2
+	MAX_STRIP_MULTIPLICITY 1
+	STRIP_ENERGY_MATCHING_NUMBER_OF_SIGMA 5
+  	STRIP_ENERGY_MATCHING_SIGMA 0.02
+	DISABLE_CHANNEL MM1STRY12
+	DISABLE_CHANNEL MM2STRY12
+	DISABLE_CHANNEL MM3STRY12
+	DISABLE_CHANNEL MM4STRY12
+	DISABLE_CHANNEL MM1STRX12
+   DISABLE_CHANNEL MM2STRX12
+   DISABLE_CHANNEL MM3STRX12
+   DISABLE_CHANNEL MM4STRX12
+	DISABLE_CHANNEL MM1STRY124
+	DISABLE_CHANNEL MM4STRX3 
+	DISABLE_CHANNEL MM4STRX4
+	DISABLE_CHANNEL MM4STRX5
+	DISABLE_CHANNEL MM4STRX6
+	DISABLE_CHANNEL MM4STRX8
+	DISABLE_CHANNEL MM4STRX31
+	DISABLE_ALL MM5
+	DISABLE_ALL MM6
+	DISABLE_ALL MM7
+	DISABLE_ALL MM8
+	SI_X_E_RAW_THRESHOLD 0
+  	CSI_E_RAW_THRESHOLD 0
+  	CSI_SIZE 256

NPAnalysis/MUGAST/configs/ConfigSSSD.dat

+ConfigSSSD
+  PEDESTAL_THRESHOLD 0
+  DISABLE_CHANNEL SI2STR1
+  DISABLE_CHANNEL SI4STR13