diff --git a/Inputs/EventGenerator/60Fe.reaction b/Inputs/EventGenerator/60Fe.reaction
index 467bd51d9509590ead464e0f56adc36d461b14b9..ea0874438d3a70b841bdbd8852d269b7e4d3e949 100644
--- a/Inputs/EventGenerator/60Fe.reaction
+++ b/Inputs/EventGenerator/60Fe.reaction
@@ -10,10 +10,10 @@ Transfert
ExcitationEnergy= 2.0
BeamEnergy= 1200
BeamEnergySpread= 0
- SigmaX= 1
- SigmaY= 1
- SigmaThetaX= 2
- SigmaPhiY= 2
+ SigmaX= 3
+ SigmaY= 3
+ SigmaThetaX= 0.5
+ SigmaPhiY= 0.5
CrossSectionPath= ni69_g7_01.n
ShootLight= 1
ShootHeavy= 0
diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh
index fd9d4b26a378eb1bf7cdd13be6564d2fe29c458a..47e09a3cb2df16b98579a2ef3f7a8f651ef9c268 100644
--- a/NPAnalysis/10He_Riken/include/ObjectManager.hh
+++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh
@@ -75,7 +75,7 @@ namespace GRAPH
{
// Declare your Spectra here:
- TH1F *myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 1000 , -5 , 5 ) ;
+ TH1F* myHist1D = new TH1F("Hist1D","Histogramm 1D ; x ; count", 100 , -40 , 40 ) ;
TH2F *myHist2D = new TH2F("Hist2D","Histogramm 2D ; x ; y ", 128 , 1 , 128 , 128 , 1 , 128 ) ;
@@ -111,15 +111,18 @@ namespace ENERGYLOSS
// 3He Energy Loss
EnergyLoss He3TargetWind = EnergyLoss ( "He3_Mylar.G4table" ,
"G4Table",
- 10000 );
+ 1000,
+ 3 );
EnergyLoss He3TargetGaz = EnergyLoss ( "He3_D2.G4table" ,
"G4Table",
- 10000 );
+ 1000,
+ 3 );
EnergyLoss He3StripAl = EnergyLoss ( "He3_Aluminium.G4table" ,
"G4Table",
- 10000 );
+ 1000,
+ 3 );
// EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" ,
// "LISE",
@@ -143,21 +146,21 @@ namespace ENERGYLOSS
EnergyLoss He3StripSi = EnergyLoss ( "3He_Si.txt" ,
"LISE",
- 10000 ,
- 1 ,
- 3 );
+ 100 ,
+ 3 ,
+ 1 );
// proton Energy Loss
EnergyLoss protonTargetWind = EnergyLoss ( "proton_Mylar.txt" ,
"LISE",
- 1000 ,
+ 100 ,
1 ,
1 );
EnergyLoss protonTargetGaz = EnergyLoss ( "proton_D2gaz_1b_26K.txt" ,
"LISE",
- 1000 ,
+ 100 ,
1 ,
1 );
diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc
index 0fe0d7c40e4efbc22e7d52b433013164740836da..8cbb0f434dcfe4c703a601b3a8d721e8e79af847 100644
--- a/NPAnalysis/10He_Riken/src/Analysis.cc
+++ b/NPAnalysis/10He_Riken/src/Analysis.cc
@@ -34,13 +34,14 @@ int main(int argc,char** argv)
myDetector -> ReadConfigurationFile(detectorfileName) ;
// Instantiate the Calibration Manger using a file
- CalibrationManager* myCalibration = new CalibrationManager(calibrationfileName) ;
+ CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ;
// Attach more branch to the output
double ELab[2], ExcitationEnergy[2] ;
double ThetaLab[2] , ThetaCM[2] ;
double X[2] , Y[2] ;
+ double EventWeight=0;
// Exitation Energy
RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",ExcitationEnergy,"ExcitationEnergy[2]/D") ;
@@ -52,11 +53,16 @@ int main(int argc,char** argv)
RootOutput::getInstance()->GetTree()->Branch("X",X,"X/D[2]") ;
RootOutput::getInstance()->GetTree()->Branch("Y",Y,"Y/D[2]") ;
+ // For phasespace event
+ RootOutput::getInstance()->GetTree()->Branch("EventWeight",&EventWeight,"EventWeight/D") ;
// Get the formed Chained Tree and Treat it
TChain* Chain = RootInput:: getInstance() -> GetChain() ;
// Open the ThinSi Branch
+ Chain->SetBranchStatus("EventWeight",true);
+ Chain->SetBranchAddress("EventWeight" ,&EventWeight );
+
Chain->SetBranchStatus("InitialConditions",true) ;
Chain->SetBranchStatus("fIC_*",true) ;
@@ -72,8 +78,38 @@ int main(int argc,char** argv)
TMust2Physics* M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ;
TPlasticPhysics* Pl = (TPlasticPhysics*) myDetector -> m_Detector["Plastic"] ;
TSSSDPhysics* ThinSi = (TSSSDPhysics*) myDetector -> m_Detector["SSSD"] ;
- cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
-
+
+
+RootOutput::getInstance()->GetList()->Add(myHist1D);
+
+ TCutG *cut3He_MUST2 = new TCutG("Cut3HeMUST2",11);
+ cut3He_MUST2->SetPoint(0,-1.49426,24.2781);
+ cut3He_MUST2->SetPoint(1,15.3161,14.3151);
+ cut3He_MUST2->SetPoint(2,47.069,7.6732);
+ cut3He_MUST2->SetPoint(3,110.575,4.35222);
+ cut3He_MUST2->SetPoint(4,308.563,2.25477);
+ cut3He_MUST2->SetPoint(5,310.431,1.20604);
+ cut3He_MUST2->SetPoint(6,232.917,1.11864);
+ cut3He_MUST2->SetPoint(7,89.0948,2.42955);
+ cut3He_MUST2->SetPoint(8,32.1264,6.71186);
+ cut3He_MUST2->SetPoint(9,1.30747,11.9555);
+ cut3He_MUST2->SetPoint(10,-1.49426,24.2781);
+
+ TCutG *cut3He_M2_SSSD = new TCutG("Cut3HeM2SSSD",11);
+ cut3He_M2_SSSD->SetPoint(0,7.44252,1.45432);
+ cut3He_M2_SSSD->SetPoint(1,40.6322,0.684984);
+ cut3He_M2_SSSD->SetPoint(2,95.9483,0.406912);
+ cut3He_M2_SSSD->SetPoint(3,218.649,0.216896);
+ cut3He_M2_SSSD->SetPoint(4,362.471,0.101033);
+ cut3He_M2_SSSD->SetPoint(5,355.431,0.0315148);
+ cut3He_M2_SSSD->SetPoint(6,100.977,0.0315148);
+ cut3He_M2_SSSD->SetPoint(7,25.546,0.314221);
+ cut3He_M2_SSSD->SetPoint(8,3.41954,0.995498);
+ cut3He_M2_SSSD->SetPoint(9,3.41954,1.40797);
+ cut3He_M2_SSSD->SetPoint(10,7.44252,1.45432);
+
+
+ cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
int i ,N=Chain -> GetEntries();
cout << " Number of Event to be treated : " << N << endl ;
@@ -84,8 +120,8 @@ int main(int argc,char** argv)
// Clear local branch
for(int hh = 0 ; hh <2 ; hh++)
{
- ELab[hh] = -1 ; ExcitationEnergy[hh] = -1 ; ThetaLab[hh] = -1 ;
- X[hh] = -1000 ; Y[hh] = -1000 ; ThetaCM[hh] = -1 ;
+ ELab[hh] = -10000 ; ExcitationEnergy[hh] = -10000 ; ThetaLab[hh] = -10000 ;
+ X[hh] = -10000 ; Y[hh] = -10000 ; ThetaCM[hh] = -10000 ;
}
// Minimum code
@@ -122,7 +158,6 @@ int main(int argc,char** argv)
// Must 2 And ThinSi //
for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
{
- ELab[hit] = -1 ; ThetaLab[hit] = -1;
// Get Hit Direction
TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(XTarget,YTarget,0);
// Angle between beam and particle
@@ -135,21 +170,24 @@ int main(int argc,char** argv)
if(M2 -> GetPositionOfInteraction(hit).Z() > 0)
{
- if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0)
+ if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0 )
{
ELab[hit] = M2 -> Si_E[hit] ;
ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 2*0.4*micrometer , // Target Thickness at 0 degree
+ 0.4*micrometer , // Target Thickness at 0 degree
ThetaMM2Surface );
if(ThinSi -> Energy.size() > 0)
{
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
ELab[hit] += ThinSi-> Energy[hit];
ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
0.4*micrometer , // Target Thickness at 0 degree
ThetaMM2Surface );
- }
+ }
ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
10*micrometer , // Target Thickness at 0 degree
@@ -161,6 +199,11 @@ int main(int argc,char** argv)
ThetaCM[hit] = He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ;
ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) ;
+
+ if(ThinSi -> Energy.size() > 0 )
+ if(cut3He_M2_SSSD->IsInside(ThinSi -> Energy[hit], M2 -> Si_E[hit]) )
+ myHist1D->Fill(ExcitationEnergy[hit],EventWeight);
+
X[hit] = HitDirection . X();
Y[hit] = HitDirection . Y();
ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -186,6 +229,9 @@ int main(int argc,char** argv)
if(ThinSi -> Energy.size() > 0)
{
+ ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
+ 0.4*micrometer , // Target Thickness at 0 degree
+ ThetaMM2Surface );
ELab[hit] += ThinSi-> Energy[hit];
ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
0.4*micrometer , // Target Thickness at 0 degree
@@ -199,9 +245,13 @@ int main(int argc,char** argv)
ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
1.5*mm , // Target Thickness at 0 degree
ThetaN );
-
+
ThetaCM[hit]= He10Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ;
ExcitationEnergy[hit] = He10Reaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit]) ;
+
+ if( cut3He_MUST2->IsInside(M2 -> Si_E[hit], M2 -> CsI_E[hit]) )
+ myHist1D->Fill(ExcitationEnergy[hit],EventWeight);
+
X[hit] = HitDirection . X();
Y[hit] = HitDirection . Y();
ThetaLab[hit] = ThetaLab[hit] / deg ;
@@ -215,12 +265,12 @@ int main(int argc,char** argv)
{} */
}
-
RootOutput::getInstance()->GetTree()->Fill() ;
}
cout << " A total of " << i << " event has been annalysed " << endl ;
cout << endl << "/////////////////////////////////"<< endl<< endl ;
+ myHist1D->Write();
RootOutput::getInstance()->Destroy();
return 0 ;
}
diff --git a/NPAnalysis/e530/RunToTreat.txt b/NPAnalysis/e530/RunToTreat.txt
index 1a9cc7e2870b0895210c65798057d40d0dcfe585..c56d1ce44cd471573db652b261e0f3be3adc4cd8 100644
--- a/NPAnalysis/e530/RunToTreat.txt
+++ b/NPAnalysis/e530/RunToTreat.txt
@@ -1,5 +1,7 @@
TTreeName
- AutoTree
+ SimulatedTree
RootFileName
- ../../../../root/run_0003.root
-
+ %../../Outputs/Simulation/3He_source.root
+ %../../Outputs/Simulation/alpha_source.root
+ %../../Outputs/Simulation/REALIST_DATA.root
+ ../../Outputs/Simulation/mySimul.root
diff --git a/NPAnalysis/e530/pipo.txt b/NPAnalysis/e530/pipo.txt
index 7b4859f4567953d6a47cd82929be6fa61bda9f41..a6ffecd1bfe9a754e329aa65d0f5fedaec9284f8 100644
--- a/NPAnalysis/e530/pipo.txt
+++ b/NPAnalysis/e530/pipo.txt
@@ -1,2 +1,2 @@
CalibrationFilePath
- dummy.txt
+ /home/Adrien/Desktop/NPTool/NPTool.dev.CalibMust2/NPAnalysis/e530/calibrrrPipo.txt
diff --git a/NPAnalysis/e530/src/Analysis.cc b/NPAnalysis/e530/src/Analysis.cc
index 555c8163e55ee1d72133708a212f64e76133592b..2561d14e48382f007406e0f1693d80867e5e0aaa 100644
--- a/NPAnalysis/e530/src/Analysis.cc
+++ b/NPAnalysis/e530/src/Analysis.cc
@@ -14,91 +14,62 @@ int main(int argc,char** argv)
return 0;
}
- string reactionfileName = argv[1] ;
+ string reactionfileName = argv[1] ;
string detectorfileName = argv[2] ;
- string calibrationfileName = argv[3] ;
- string runToReadfileName = argv[4] ;
+ string calibrationfileName = argv[3] ;
+ string runToTreatFileName = argv[4] ;
- // First of All instantiate RootInput and Output
- // Detector will be attached later
- RootInput:: getInstance(runToReadfileName) ;
- //RootOutput::getInstance("Analysis/10HeRiken_AnalyzedData", "AnalyzedTree") ;
- RootOutput::getInstance("Analysis/Fe60", "Analysed_Data") ;
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
+ // First of All instantiate RootInput and Output
+ // Detector will be attached later
+ RootInput:: getInstance(runToTreatFileName) ;
+ RootOutput::getInstance("Analysis/60Fe_AnalyzedData", "AnalyzedTree") ;
- // Instantiate some Reaction
+ // Instantiate some Reaction
+ NPL::Reaction* e530Reaction = new Reaction ;
+ e530Reaction -> ReadConfigurationFile(reactionfileName) ;
- NPL::Reaction* Fe60Reaction = new Reaction;
- Fe60Reaction -> ReadConfigurationFile(reactionfileName);
+ // Instantiate the Calibration Manger using a file (WARNING:prior to the detector instantiation)
+ CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ;
- // Instantiate the detector using a file
- NPA::DetectorManager* myDetector = new DetectorManager ;
- myDetector -> ReadConfigurationFile(detectorfileName) ;
-
- // Instantiate the Calibration Manger using a file
- // CalibrationManager* myCalibration = new CalibrationManager(calibrationfileName) ;
- CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName);
-
- // Attach more branch to the output
+ // Instantiate the detector using a file
+ NPA::DetectorManager* myDetector = new DetectorManager ;
+ myDetector -> ReadConfigurationFile(detectorfileName) ;
- double ELab[2], ExcitationEnergy[2] ;
- double ThetaLab[2] , ThetaCM[2] ;
- double X[2] , Y[2] ;
+ // Ask the detector manager to load the parameter added by the detector in the calibrationfileName
+ myCalibration->LoadParameterFromFile() ;
+ /////////////////////////////////////////////////////////////////////////////////////////////////////
- // Exitation Energy
- RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",ExcitationEnergy,"ExcitationEnergy[2]/D") ;
-
- //E lab et Theta lab
- RootOutput::getInstance()->GetTree()->Branch("ThetaCM",ThetaCM,"ThetaCM[2]/D") ;
- RootOutput::getInstance()->GetTree()->Branch("ELab",ELab,"ELab[2]/D") ;
- RootOutput::getInstance()->GetTree()->Branch("ThetaLab",ThetaLab,"ThetaLab[2]/D") ;
- RootOutput::getInstance()->GetTree()->Branch("X",X,"X/D[2]") ;
- RootOutput::getInstance()->GetTree()->Branch("Y",Y,"Y/D[2]") ;
-
-
+ // Attach more branch to the output
+ double ELab=0,ThetaLab=0,ExcitationEnergy=0;
+ RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D") ;
+ RootOutput::getInstance()->GetTree()->Branch("ExcitationEnergy",&ExcitationEnergy,"ExcitationEnergy/D") ;
// Get the formed Chained Tree and Treat it
TChain* Chain = RootInput:: getInstance() -> GetChain() ;
- // Open the ThinSi Branch
- // Chain->SetBranchStatus("InitialConditions",true) ;
- // Chain->SetBranchStatus("fIC_*",true) ;
-
- // TInitialConditions* Init = new TInitialConditions();
- // Chain->SetBranchAddress("InitialConditions" ,&Init );
-
- double XTarget=0 ; double YTarget=0; double BeamTheta = 0 ; double BeamPhi = 0 ; double E=-1000;
-
- // Get Detector Pointer:
TMust2Physics* M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ;
- // TPlasticPhysics* Pl = (TPlasticPhysics*) myDetector -> m_Detector["Plastic"] ;
- // TSSSDPhysics* ThinSi = (TSSSDPhysics*) myDetector -> m_Detector["SSSD"] ;
cout << " ///////// Starting Analysis ///////// "<< endl << endl ;
int i ,N=Chain -> GetEntries();
- N = 1000;
+ //N = 1000;
cout << " Number of Event to be treated : " << N << endl ;
clock_t begin=clock();
clock_t end=begin;
for ( i = 0 ; i < N ; i ++ )
{
- // Clear local branch
- for(int hh = 0 ; hh <2 ; hh++)
- {
- ELab[hh] = -1 ; ExcitationEnergy[hh] = -1 ; ThetaLab[hh] = -1 ;
- X[hh] = -1000 ; Y[hh] = -1000 ; ThetaCM[hh] = -1 ;
- }
-
// Minimum code
- if( i%10000 == 0 && i!=0) {
- cout.precision(5);
- end=clock();
- double TimeElapsed = (end-begin)/CLOCKS_PER_SEC;
- double percent = (double)i/N ;
- double TimeToWait = (TimeElapsed/percent) - TimeElapsed ;
- cout << "\r Progression:" << percent*100
- << " % \t | \t Remaining time : ~"
- << TimeToWait <<"s"<< flush;
+ if( i%10000 == 0 && i!=0) {
+ cout.precision(5);
+ end=clock();
+ double TimeElapsed = (end-begin)/CLOCKS_PER_SEC;
+ double percent = (double)i/N ;
+ double TimeToWait = (TimeElapsed/percent) - TimeElapsed ;
+ cout << "\r Progression:" << percent*100
+ << " % \t | \t Remaining time : ~"
+ << TimeToWait <<"s"<< flush;
}
else if (i==N-1) cout << "\r Progression:"
@@ -110,153 +81,18 @@ int main(int argc,char** argv)
// Build the new one
myDetector -> BuildPhysicalEvent() ;
////
+
-
- // Target (from initial condition)
- //XTarget = Init->GetICPositionX(0);
- // YTarget = Init->GetICPositionY(0);
- // XTarget = RandomEngine.Gaus(Init->GetICPositionX(0),1);
- // YTarget = RandomEngine.Gaus(Init->GetICPositionY(0),1);
- // BeamTheta = Init->GetICIncidentAngleTheta(0)*deg ;
- // BeamPhi = Init->GetICIncidentAnglePhi(0)*deg ;
- TVector3 BeamDirection = TVector3(cos(BeamPhi)*sin(BeamTheta) , sin(BeamPhi)*sin(BeamTheta) , cos(BeamTheta)) ;
- ////
-
- // Must 2 And ThinSi //
- //for(int hit = 0; hit < M2 -> GetEventMultiplicity() ; hit ++)
- for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
- {
- ELab[hit] = -1 ; ThetaLab[hit] = -1;
- // Get Hit Direction
- TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(XTarget,YTarget,0);
-
- // Angle between beam and particle
- ThetaLab[hit] = ThetaCalculation ( HitDirection , BeamDirection ) ;
- // Angle between particule and z axis (target Normal)
- double ThetaN = ThetaCalculation ( HitDirection , TVector3(0,0,1) ) ;
- // ANgle between particule and Must2 Si surface
- double ThetaMM2Surface = ThetaCalculation ( HitDirection , M2 -> GetTelescopeNormal(hit) );
-
- if(M2 -> GetPositionOfInteraction(hit).Z() > 0)
- {
- if( M2 -> CsI_E[hit] == 0 && M2 -> Si_E[hit] > 0)
- {
- ELab[hit] = M2 -> Si_E[hit] ;
-
- /*
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 2*0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- */
- // ELab[hit]= He3StripSi.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- // 20*micrometer , // Target Thickness at 0 degree
- // ThetaMM2Surface );
-
- //if(ThinSi -> Energy.size() > 0)ELab[hit] += ThinSi-> Energy[hit];
-
- /*
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- ThetaN );
-
- ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 1.5*mm , // Target Thickness at 0 degree
- ThetaN );
- */
-
- ThetaCM[hit] = Fe60Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ;
- ExcitationEnergy[hit] = Fe60Reaction -> ReconstructRelativistic( ELab[hit] , ThetaLab[hit] ) ;
- X[hit] = HitDirection . X();
- Y[hit] = HitDirection . Y();
- ThetaLab[hit] = ThetaLab[hit] / deg ;
- }
-
- /*
- else if(M2 ->CsI_E[hit] > 0 && M2 -> Si_E[hit] > 0)
- {
-
- ELab[hit]= M2 ->CsI_E[hit] ;
-
- ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 3*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
- ELab[hit]+= M2 ->Si_E[hit];
-
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- //if(ThinSi -> Energy.size() > 0)ELab[hit] += ThinSi-> Energy[hit];
-
-
- ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- ThetaN );
-
- ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 1.5*mm , // Target Thickness at 0 degree
- ThetaN );
-
- ThetaCM[hit]= Fe60Reaction -> EnergyLabToThetaCM( ELab[hit] ) /deg ;
- ExcitationEnergy[hit] = Fe60Reaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit] ) ;
- X[hit] = HitDirection . X();
- Y[hit] = HitDirection . Y();
- ThetaLab[hit] = ThetaLab[hit] / deg ;
- }
- */
- else {ExcitationEnergy[hit]=-100 ; X[hit] = -100 ; Y[hit]= -100 ; ThetaLab[hit]=-100;ThetaCM[hit]=-100;}
-
- }
-
- /*else if(M2 -> GetPositionOfInteraction(hit).Z()<0)
- {
-
- if(ELab[hit]>-1000 )
- {
- if(ELab[hit]>18)
- {
- ELab[hit]= protonStripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
- }
-
- ELab[hit]= protonStripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 0.4*micrometer , // Target Thickness at 0 degree
- ThetaMM2Surface );
-
- ELab[hit]= protonTargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 15*micrometer , // Target Thickness at 0 degree
- ThetaN );
-
- ELab[hit]= protonTargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle
- 1.5*mm , // Target Thickness at 0 degree
- ThetaN );
- ThetaCM[hit] = myReaction -> EnergyLabToThetaCM( ELab[hit] , 1 ) /deg ;
- ExcitationEnergy[hit] = myReaction -> ReconstructRelativistic( ELab[hit], ThetaLab[hit] ) ;
- X[hit] = HitDirection . X();
- Y[hit] = HitDirection . Y();
- ThetaLab[hit] = ThetaLab[hit] / deg ;
- }
-
- else {ExcitationEnergy[hit]=-100 ; X[hit] = -100 ; Y[hit] = -100 ;ThetaLab[hit]=-100; ThetaCM[hit]=-100 ;}
-
- } */
-
-
- }
+ // Must 2
+ for(int hit = 0; hit < M2 -> Si_E.size() ; hit ++)
+ {
+ ELab = -1 ; ThetaLab = -1;
+ // Get Hit Direction
+ TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(0,0,-40);
+ // Angle between beam and particle
+ ThetaLab = ThetaCalculation ( HitDirection , TVector3(0,0,1) ) ;
+ ELab = M2 -> Si_E[hit] + M2 -> SiLi_E[hit] ;
+ }
RootOutput::getInstance()->GetTree()->Fill() ;
}
diff --git a/NPLib/CalibrationManager/CalibrationManager.cxx b/NPLib/CalibrationManager/CalibrationManager.cxx
index 085c08285c60962b323ad2a89135c9683dd309bb..737c31d5997e7db56b15b7efa1eb08490066d62a 100644
--- a/NPLib/CalibrationManager/CalibrationManager.cxx
+++ b/NPLib/CalibrationManager/CalibrationManager.cxx
@@ -25,6 +25,7 @@
#include <cstdlib>
#include <limits>
#include <cmath>
+#include <sstream>
//////////////////////////////////////////////////////////////////
CalibrationManager* CalibrationManager::instance = 0;
@@ -58,15 +59,17 @@ CalibrationManager::CalibrationManager(string configFileName)
return;
}
- else {
+ else
+ {
+ cout << "Reading list of file from :" << configFileName << endl;
while (!inputConfigFile.eof()) {
getline(inputConfigFile, lineBuffer);
// search for token giving the list of Root files to treat
- if ( lineBuffer.compare(0, 12, "CalibrationFilePath") == 0 ) {
+ if ( lineBuffer.compare(0, 19, "CalibrationFilePath") == 0 ) {
while (!inputConfigFile.eof()) {
inputConfigFile >> dataBuffer;
-
+
// ignore comment Line
if (dataBuffer.compare(0, 1, "%") == 0) {
inputConfigFile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
@@ -76,11 +79,12 @@ CalibrationManager::CalibrationManager(string configFileName)
AddFile(dataBuffer);
cout << "Adding file " << dataBuffer << " to Calibration" << endl;
}
- }
- }
- }
+ }
+ }
+ }
+ }
+ cout << "/////////////////////////////////" << endl;
}
- cout << "/////////////////////////////////" << endl;}
//////////////////////////////////////////////////////////////////
CalibrationManager::~CalibrationManager()
@@ -99,49 +103,56 @@ void CalibrationManager::LoadParameterFromFile()
{
ifstream CalibFile ;
string DataBuffer ;
+ string LineBuffer ;
for(unsigned int i = 0 ; i < fFileList.size() ; i++)
{
CalibFile.open( fFileList[i].c_str() );
- vector<double> Coeff ;
map<string,string>::iterator it ;
if(!CalibFile)
{
cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
- cout << " WARNING: FILE " << fFileList[i] << " IS MISSING " << endl ;
+ cout << " WARNING: FILE " << fFileList[i] << " IS MISSING " << endl ;
cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
-
}
else while( !CalibFile.eof() )
+ {
+ // Read the file Line by line
+ getline(CalibFile, LineBuffer);
+
+ // Create a istringstream to manipulate the line easely
+ istringstream theLine (LineBuffer,istringstream::in);
+ theLine >> DataBuffer ;
+
+ // Comment support, comment symbole is %
+ if(DataBuffer.compare(0, 1, "%") == 0) {
+ CalibFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+
+ // Search word in the token list
+ it=fToken.find(DataBuffer);
+
+ // if the word is find, values are read
+ if( it!=fToken.end() )
{
- CalibFile >> DataBuffer ;
-
- // Search word in the token list
- it=fToken.find(DataBuffer);
-
- // if the word is find, values are read
- if( it!=fToken.end() )
+ vector<double> Coeff ;
+ while( !theLine.eof() )
{
-
- Coeff.clear();
- while(DataBuffer!="\n")
- {
- CalibFile >> DataBuffer ; Coeff.push_back( atof(DataBuffer.c_str()) ) ;
- }
+ theLine >> DataBuffer ; Coeff.push_back( atof(DataBuffer.c_str()) ) ;
+ }
- // Check this parameter is not already define
- if( fCalibrationCoeff.find(it->second) != fCalibrationCoeff.end() ) cout << "WARNING: Parameter " << it->second << " Already found. It will be rewritted " << endl;
+ // Check this parameter is not already define
+ if( fCalibrationCoeff.find(it->second) != fCalibrationCoeff.end() )
+ cout << "WARNING: Parameter " << it->second << " Already found. It will be rewritted " << endl;
- // Add the list of Coeff to the Coeff map using Parameter Path as index
- fCalibrationCoeff[ it->second ] = Coeff ;
- }
-
+ // Add the list of Coeff to the Coeff map using Parameter Path as index
+ fCalibrationCoeff[ it->second ] = Coeff ;
}
- CalibFile.close() ;
+
+ }
+ CalibFile.close() ;
}
-
}
//////////////////////////////////////////////////////////////////
@@ -149,8 +160,12 @@ double CalibrationManager::ApplyCalibration(string ParameterPath , double RawVal
{
double CalibratedValue = 0 ;
map< string , vector<double> >::iterator it ;
+
+ // Find the good parameter in the Map
+ // Using Find method of stl is the fastest way
it = fCalibrationCoeff.find(ParameterPath) ;
+ // If the find methods return the end iterator it's mean the parameter was not found
if(it == fCalibrationCoeff.end() )
{
/* cout << "XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX " << endl ;
@@ -160,8 +175,13 @@ double CalibrationManager::ApplyCalibration(string ParameterPath , double RawVal
return RawValue ;
}
+ // Else we take the second part of the element (first is index, ie: parameter path)
+ // Second is the vector of Coeff
vector<double> Coeff = it->second ;
+ // The vector size give the degree of calibration
+ // We just apply the coeff and returned the calibrated value
+
for(unsigned int i = 0 ; i < Coeff.size() ; i++)
{
CalibratedValue += Coeff[i]*pow(RawValue, (double)i);
diff --git a/NPLib/CalibrationManager/CalibrationManager.h b/NPLib/CalibrationManager/CalibrationManager.h
index 69df3a9f7f5e1ed8ee4519c1dbeaaa488cdd128d..a5fc1950b00a8bcdd5c7c382356ec3826d705f2e 100644
--- a/NPLib/CalibrationManager/CalibrationManager.h
+++ b/NPLib/CalibrationManager/CalibrationManager.h
@@ -33,7 +33,7 @@ using namespace std ;
class CalibrationManager
{
- public: // Constructor and Destructor
+ protected: // Constructor and Destructor are protected because the class is a singleton
CalibrationManager(string configFileName);
~CalibrationManager();
@@ -50,20 +50,20 @@ class CalibrationManager
inline void AddFile(string Path) { fFileList.push_back(Path) ;} ;
- public: // Declaration of Calibration
-
+ public: // Calibration Parameter Related
+
// call like : myCalibrationManager->AddParameter( "MUST2" ,"Telescope5_Si_X38_E", "T5_Si_X38_E" )
// return false if the token is not found in the file list
- bool AddParameter(string DetectorName , string ParameterName , string Token ) ;
+ bool AddParameter(string DetectorName , string ParameterName , string Token) ;
// call like : myCalibrationManager->ApplyCalibration( "MUST2/Telescope5_Si_X38_E" , RawEnergy )
// return the Calibrated value
double ApplyCalibration(string ParameterPath , double RawValue);
- public: // To be called after initialisation
- // Loop over the file list and catch the file used for calibration
- void LoadParameterFromFile();
+ public: // To be called after initialisation
+ // Loop over the file list and catch the file used for calibration
+ void LoadParameterFromFile();
private:
diff --git a/NPLib/MUST2/TMust2Data.h b/NPLib/MUST2/TMust2Data.h
index 654e999cdb441c7e62ae992f3b80cef83acd05fd..511cfc0ca6cc07dad1a38b485611bcacd9b1fd98 100644
--- a/NPLib/MUST2/TMust2Data.h
+++ b/NPLib/MUST2/TMust2Data.h
@@ -151,7 +151,7 @@ class TMust2Data : public TObject {
// CsI
//(E)
- UShort_t GetMMCsIEMult() {return fMM_CsIE_DetectorNbr.size();}
+ UShort_t GetMMCsIEMult() {return fMM_CsIE_DetectorNbr.size();}
UShort_t GetMMCsIEDetectorNbr(Int_t i) {return fMM_CsIE_DetectorNbr.at(i);}
UShort_t GetMMCsIECristalNbr(Int_t i) {return fMM_CsIE_CristalNbr.at(i);}
Double_t GetMMCsIEEnergy(Int_t i) {return fMM_CsIE_Energy.at(i);}
diff --git a/NPLib/MUST2/TMust2Physics.cxx b/NPLib/MUST2/TMust2Physics.cxx
index 13c2d3a9708380cfc761004c39b2753125cce8ad..1cf3262a3f8e76f7f9677310d120a72e59ee606e 100644
--- a/NPLib/MUST2/TMust2Physics.cxx
+++ b/NPLib/MUST2/TMust2Physics.cxx
@@ -16,8 +16,6 @@
* *
*---------------------------------------------------------------------------*
* Comment: *
- * Only multiplicity one and multiplicity 2 are down. *
- * Improvment needed *
* *
*****************************************************************************/
#include "TMust2Physics.h"
@@ -40,17 +38,17 @@ using namespace LOCAL;
ClassImp(TMust2Physics)
///////////////////////////////////////////////////////////////////////////
TMust2Physics::TMust2Physics()
- {
- EventMultiplicity = 0 ;
- EventData = new TMust2Data ;
- EventPhysics = this ;
+ {
+ EventMultiplicity = 0 ;
+ EventData = new TMust2Data ;
+ EventPhysics = this ;
+ NumberOfTelescope = 0 ;
}
///////////////////////////////////////////////////////////////////////////
void TMust2Physics::BuildSimplePhysicalEvent()
{
BuildPhysicalEvent();
-
}
///////////////////////////////////////////////////////////////////////////
@@ -95,14 +93,24 @@ void TMust2Physics::BuildPhysicalEvent()
{
if(EventData->GetMMSiLiEDetectorNbr(j)==N)
{
- // if SiLi energy is above threshold check the compatibility
+ // SiLi energy is above threshold check the compatibility
if( fSiLi_E(EventData , j)>SiLi_E_Threshold )
{
+ // pad vs strip number match
if( Match_Si_SiLi( X, Y , EventData->GetMMSiLiEPadNbr(j) ) )
{
SiLi_N.push_back(EventData->GetMMSiLiEPadNbr(j)) ;
- SiLi_E.push_back(fSiLi_E(EventData , j)) ;
- SiLi_T.push_back(fSiLi_T(EventData , j)) ;
+ SiLi_E.push_back(fSiLi_E(EventData , j)) ;
+
+ // Look for associate energy
+ // Note: in case of use of SiLi "Orsay" time is not coded.
+ for(int k =0 ; k < EventData->GetMMSiLiTMult() ; k ++)
+ {
+ // Same Pad, same Detector
+ if( EventData->GetMMSiLiEPadNbr(j)==EventData->GetMMSiLiEPadNbr(k) && EventData->GetMMSiLiEDetectorNbr(j)==EventData->GetMMSiLiTDetectorNbr(k) )
+ {SiLi_T.push_back(fSiLi_T(EventData , k)) ; break ;}
+ }
+
check_SILI = true ;
}
@@ -114,14 +122,20 @@ void TMust2Physics::BuildPhysicalEvent()
{
if(EventData->GetMMCsIEDetectorNbr(j)==N)
{
- // ifCsI energy is above threshold check the compatibility
+ // CsI energy is above threshold check the compatibility
if( fCsI_T(EventData , j)>CsI_E_Threshold )
{
- if( Match_Si_CsI( X, Y , EventData->GetMMCsIECristalNbr(j) ) )
+ if(Match_Si_CsI( X, Y , EventData->GetMMCsIECristalNbr(j) ) )
{
CsI_N.push_back(EventData->GetMMCsIECristalNbr(j)) ;
- CsI_E.push_back(fCsI_E(EventData , j)) ;
- CsI_T.push_back(fCsI_T(EventData , j)) ;
+ CsI_E.push_back(fCsI_E(EventData , j)) ;
+
+ for(int k =0 ; k < EventData->GetMMCsITMult() ; k ++)
+ {
+ if( EventData->GetMMCsIECristalNbr(j)==EventData->GetMMCsITCristalNbr(k) && EventData->GetMMCsIEDetectorNbr(j)==EventData->GetMMCsITDetectorNbr(k) )
+ {CsI_T.push_back(fCsI_T(EventData , k)) ; break ;}
+ }
+
check_CSI = true ;
}
}
@@ -132,15 +146,15 @@ void TMust2Physics::BuildPhysicalEvent()
if(!check_SILI)
{
SiLi_N.push_back(0) ;
- SiLi_E.push_back(0) ;
- SiLi_T.push_back(0) ;
+ SiLi_E.push_back(-10000) ;
+ SiLi_T.push_back(-10000) ;
}
if(!check_CSI)
{
CsI_N.push_back(0) ;
- CsI_E.push_back(0) ;
- CsI_T.push_back(0) ;
+ CsI_E.push_back(-10000) ;
+ CsI_T.push_back(-10000) ;
}
}
@@ -177,6 +191,11 @@ vector < TVector2 > TMust2Physics :: Match_X_Y()
{
vector < TVector2 > ArrayOfGoodCouple ;
+ // Prevent code from treating very high multiplicity Event
+ // Those event are not physical anyway and that improve speed.
+ if( EventData->GetMMStripXEMult()>6 || EventData->GetMMStripYEMult()>6 )
+ return ArrayOfGoodCouple;
+
for(int i = 0 ; i < EventData->GetMMStripXEMult(); i++)
{
// if X value is above threshold, look at Y value
@@ -191,7 +210,7 @@ vector < TVector2 > TMust2Physics :: Match_X_Y()
// if same detector check energy
if ( EventData->GetMMStripXEDetectorNbr(i) == EventData->GetMMStripYEDetectorNbr(j) )
{
- // Look if energy match
+ // Look if energy match (within 10%)
if( ( fSi_X_E(EventData , i) - fSi_Y_E(EventData , j) ) / fSi_X_E(EventData , i) < 0.1 )
ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ;
}
diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx
index 398198c956ba99c9791339f38dd3d3291b79f9b0..92991cdee1f40dc73a2d35a2e8955d8452f7aac5 100644
--- a/NPLib/VDetector/DetectorManager.cxx
+++ b/NPLib/VDetector/DetectorManager.cxx
@@ -302,6 +302,7 @@ void DetectorManager::InitializeRootOutput()
void DetectorManager::AddDetector(string DetectorName , VDetector* newDetector)
{
m_Detector[DetectorName] = newDetector;
+ newDetector->AddParameterToCalibrationManager();
}