diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector index 2cace566d03d4546598209b313a2100903f7e81a..b9cb10e37b93666dc589fa0891edb6d25c6bf812 100644 --- a/Inputs/DetectorConfiguration/Riken_65mm.detector +++ b/Inputs/DetectorConfiguration/Riken_65mm.detector @@ -7,8 +7,8 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% GeneralTarget %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -%Target - THICKNESS= 1 +Target + THICKNESS= 0.001 RADIUS= 45 MATERIAL= CD2 ANGLE= 45 @@ -16,7 +16,7 @@ GeneralTarget Y= 0 Z= 0 -CryoTarget +%CryoTarget THICKNESS= 3 RADIUS= 45 TEMPERATURE= 26 @@ -86,15 +86,15 @@ CSI= 1 VIS= all %%%%%%% Telescope 6 %%%%%%% %M2Telescope -%X1_Y1= 155.77 50.23 8.18 -%X1_Y128= 155.77 -50.23 8.18 -%X128_Y1= 133.17 50.23 -89.7 -%X128_Y128= 133.17 -50.23 -89.7 -%SI= 1 -%SILI= 1 -%CSI= 0 -%VIS= all -%% Alternate Telescope 6 %% +X1_Y1= 155.77 50.23 8.18 +X1_Y128= 155.77 -50.23 8.18 +X128_Y1= 133.17 50.23 -89.7 +X128_Y128= 133.17 -50.23 -89.7 +SI= 1 +SILI= 1 +CSI= 0 +VIS= all +% Alternate Telescope 6 %% M2Telescope THETA= -130 PHI= 0 @@ -115,77 +115,89 @@ SILI= 1 CSI= 0 VIS= all -%%%%%%%%%%%%%%%%%%%% +%%%%%%% Telescope 8 %%%%%%% +M2Telescope +X1_Y1= 27.07 50.23 -154.49 +X1_Y128= 116.58 50.23 -108.88 +X128_Y1= 27.07 -50.23 -154.49 +X128_Y128= 116.58 -50.23 -108.88 +SI= 1 +SILI= 1 +CSI= 0 +VIS= all + + +%%%%%%%%%%%%%%%%%%%%% AddThinSi -%%%%%%%%% Det 1 %%%%%%%% +%%%%%%%%%% Det 1 %%%%%%%% ThinSi A= 17.61 9.85 104.11 B= 64.48 9.85 85.31 C= 58.66 56.29 70.79 D= 11.79 56.29 89.69 - -%%%%%%%%% Det 2 %%%%%%%% +% +%%%%%%%%%% Det 2 %%%%%%%% ThinSi A= -11.79 56.29 89.59 B= -58.66 56.29 70.79 C= -64.48 9.85 85.31 D= -17.61 9.85 104.11 - -%%%%%%%%% Det 3 %%%%%%%% +% +%%%%%%%%%% Det 3 %%%%%%%% ThinSi A= -17.61 -9.85 104.11 B= -64.48 -9.85 85.31 C= -58.66 -56.29 70.79 D= -11.79 -56.29 89.59 - -%%%%%%%%% Det 4 %%%%%%%% +% +%%%%%%%%%% Det 4 %%%%%%%% ThinSi A= 11.79 -56.29 89.59 B= 58.66 -56.29 70.79 C= 64.48 -9.85 85.31 D= 17.61 -9.85 104.11 - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -ScintillatorPlastic -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Plastic - X= 0 - Y= -15 - Z= 318 - Shape= Square - Height= 70 - Width= 40 - Thickness= 20 - Scintillator= BC400 - LeadThickness= 0 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Plastic - X= 0 - Y= -15 - Z= 343 - Thickness= 30 - Shape= Square - Height= 70 - Width= 40 - Scintillator= BC400 - LeadThickness= 2 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%ScintillatorPlastic +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Plastic - THETA= 0 - PHI= 0 - R= 318 - Thickness= 20 - Radius= 20 - Scintillator= BC400 - LeadThickness= 0 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% X= 0 +% Y= -15 +% Z= 318 +% Shape= Square +% Height= 70 +% Width= 40 +% Thickness= 20 +% Scintillator= BC400 +% LeadThickness= 0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Plastic - THETA= 0 - PHI= 0 - R= 343 - Thickness= 30 - Radius= 20 - Scintillator= BC400 - LeadThickness= 2 -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% X= 0 +% Y= -15 +% Z= 343 +% Thickness= 30 +% Shape= Square +% Height= 70 +% Width= 40 +% Scintillator= BC400 +% LeadThickness= 2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%Plastic +% THETA= 0 +% PHI= 0 +% R= 318 +% Thickness= 20 +% Radius= 20 +% Scintillator= BC400 +% LeadThickness= 0 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%Plastic +% THETA= 0 +% PHI= 0 +% R= 343 +% Thickness= 30 +% Radius= 20 +% Scintillator= BC400 +% LeadThickness= 2 +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% diff --git a/Inputs/EventGenerator/10He.reaction b/Inputs/EventGenerator/10He.reaction index 9b71622344f5ef1e607ba911c66cf308cd444332..ee2cb1ea1089651bdddd99f13061d8086992d5d6 100644 --- a/Inputs/EventGenerator/10He.reaction +++ b/Inputs/EventGenerator/10He.reaction @@ -10,10 +10,10 @@ TransfertToResonance ExcitationEnergy= 1.0 BeamEnergy= 550 BeamEnergySpread= 0 - SigmaThetaX= 0.6921330164 - SigmaPhiY= 0.963142053 - SigmaX= 6.232 - SigmaY= 9.069 + SigmaThetaX= 0.01 + SigmaPhiY= 0.01 + SigmaX= 0.01 + SigmaY= 0.01 ResonanceWidth= 0 ResonanceDecayZ= 2 ResonanceDecayA= 8 diff --git a/NPAnalysis/10He_Riken/configs/ConfigMust2.dat b/NPAnalysis/10He_Riken/configs/ConfigMust2.dat new file mode 100644 index 0000000000000000000000000000000000000000..de05eaf790282a3502cb23d9630fecbc78d05d64 --- /dev/null +++ b/NPAnalysis/10He_Riken/configs/ConfigMust2.dat @@ -0,0 +1,4 @@ +ConfigMust2 + MAX_STRIP_MULTIPLICITY 1 + STRIP_ENERGY_MATCHING_TOLERANCE 10 + MUST2 DISABLE_ALL MM4 diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh index 47e09a3cb2df16b98579a2ef3f7a8f651ef9c268..b161bf02f99c8fba7cedea66dccd00b446e70ebb 100644 --- a/NPAnalysis/10He_Riken/include/ObjectManager.hh +++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh @@ -35,6 +35,7 @@ #include "TMust2Physics.h" #include "TSSSDPhysics.h" #include "TPlasticPhysics.h" +#include "NPOptionManager.h" // Use CLHEP System of unit and Physical Constant #include "CLHEP/Units/GlobalSystemOfUnits.h" diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc index bebe502b09a365729d464274121a6f8c7965f70a..2792d41bf03fdd1e5f386eb50af1f78138773545 100644 --- a/NPAnalysis/10He_Riken/src/Analysis.cc +++ b/NPAnalysis/10He_Riken/src/Analysis.cc @@ -1,26 +1,18 @@ #include "ObjectManager.hh" - using namespace std; int main(int argc,char** argv) { - - if(argc!=4) - { - cout << - "you need to specify both a Reaction file and a Detector file such as : Analysis myReaction.reaction myDetector.detector runToRead.run" - << endl; - return 0; - } - - string detectorfileName = argv[1] ; - string calibrationfileName = argv[2] ; - string runToReadfileName = argv[3] ; + + NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv) ; + string detectorfileName = myOptionManager->GetDetectorFilePath() ; + string calibrationfileName = myOptionManager->GetCalibrationFilePath() ; + string runToReadfileName = myOptionManager->GetRunToReadFilePath() ; // First of All instantiate RootInput and Output // Detector will be attached later RootInput:: getInstance(runToReadfileName) ; - RootOutput::getInstance("Analysis/10HeRiken_AnalyzedData", "AnalyzedTree") ; + RootOutput::getInstance("Analysis/10HeRiken_AnalyzedData", "AnalyzedTree") ; // Instantiate some Reaction NPL::Reaction* He10Reaction = new Reaction ; @@ -78,37 +70,8 @@ int main(int argc,char** argv) // Get Detector Pointer: TMust2Physics* M2 = (TMust2Physics*) myDetector -> m_Detector["MUST2"] ; - TPlasticPhysics* Pl = (TPlasticPhysics*) myDetector -> m_Detector["Plastic"] ; - TSSSDPhysics* ThinSi = (TSSSDPhysics*) myDetector -> m_Detector["SSSD"] ; - - -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); +// TPlasticPhysics* Pl = (TPlasticPhysics*) myDetector -> m_Detector["Plastic"] ; +// TSSSDPhysics* ThinSi = (TSSSDPhysics*) myDetector -> m_Detector["SSSD"] ; cout << " ///////// Starting Analysis ///////// "<< endl << endl ; int i ,N=Chain -> GetEntries(); @@ -183,35 +146,35 @@ RootOutput::getInstance()->GetList()->Add(myHist1D); { 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 ); +// 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]= 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 - 0*micrometer , // Target Thickness at 0 degree - ThetaN ); - - ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle - 3*micrometer , // Target Thickness at 0 degree - ThetaN ); - +// 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 +//// 0*micrometer , // Target Thickness at 0 degree +//// ThetaN ); +//// +//// ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle +//// 3*micrometer , // Target Thickness at 0 degree +//// ThetaN ); +// 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); +// 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(); @@ -223,43 +186,43 @@ RootOutput::getInstance()->GetList()->Add(myHist1D); 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]= 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 ); +// 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]= 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 - 0*micrometer , // Target Thickness at 0 degree - ThetaN ); - - ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle - 3*micrometer , // Target Thickness at 0 degree - ThetaN ); - +// 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 +//// 0*micrometer , // Target Thickness at 0 degree +//// ThetaN ); +//// +//// ELab[hit]= He3TargetGaz.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle +//// 3*micrometer , // 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); +// +// 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(); diff --git a/NPAnalysis/Gaspard/include/ObjectManager.hh b/NPAnalysis/Gaspard/include/ObjectManager.hh index f9e4038e881c11c18621009d880300e536d3c079..6b3d302ca514b2675afc3a2bf34da6d0f1ff2f4f 100644 --- a/NPAnalysis/Gaspard/include/ObjectManager.hh +++ b/NPAnalysis/Gaspard/include/ObjectManager.hh @@ -7,6 +7,7 @@ // NPA #include "DetectorManager.h" +#include "NPOptionManager.h" #include "GaspardTracker.h" // STL C++ diff --git a/NPAnalysis/Paris/include/ObjectManager.hh b/NPAnalysis/Paris/include/ObjectManager.hh index bb19ac44e83382ed4d622ba53fc06fe00793c4c2..6eff0e0b91b53d28d9bed5694082f063dba47209 100644 --- a/NPAnalysis/Paris/include/ObjectManager.hh +++ b/NPAnalysis/Paris/include/ObjectManager.hh @@ -7,6 +7,7 @@ // NPA #include "DetectorManager.h" +#include "NPOptionManager.h" #include "Paris.h" #include "Shield.h" diff --git a/NPAnalysis/Paris/src/Analysis.cc b/NPAnalysis/Paris/src/Analysis.cc index a3347d282d5cb155e0137794951211328c06774f..9114a85f91c48beef3ad6f17de987bfec8bd27b7 100644 --- a/NPAnalysis/Paris/src/Analysis.cc +++ b/NPAnalysis/Paris/src/Analysis.cc @@ -6,19 +6,12 @@ using namespace std; int main(int argc,char** argv) { - // test if number of arguments is correct - if (argc != 4) { - cout << - "you need to specify both a Reaction file and a Detector file such as : Analysis myReaction.reaction myDetector.detector runToRead.run" - << endl; - return 0; - } - - // get arguments - string reactionfileName = argv[1]; - string detectorfileName = argv[2]; - string runToReadfileName = argv[3]; + NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv) ; + string detectorfileName = myOptionManager->GetDetectorFilePath() ; + string reactionfileName = myOptionManager->GetCalibrationFilePath() ; + string runToReadfileName = myOptionManager->GetRunToReadFilePath() ; + // Instantiate RootInput and RootOutput singleton classes RootInput:: getInstance(runToReadfileName); RootOutput::getInstance("Analysis/Paris_AnalyzedData", "AnalyzedTree"); diff --git a/NPAnalysis/Template/include/ObjectManager.hh b/NPAnalysis/Template/include/ObjectManager.hh index be80b3102aab19cf42f067d5cebe23d2dee120dd..bee3f4b78da044dfac6f8bdc1beaa36097976ae2 100644 --- a/NPAnalysis/Template/include/ObjectManager.hh +++ b/NPAnalysis/Template/include/ObjectManager.hh @@ -7,6 +7,7 @@ // NPA #include "DetectorManager.h" +#include "NPOptionManager.h" // STL C++ #include <iostream> diff --git a/NPAnalysis/Template/src/Analysis.cc b/NPAnalysis/Template/src/Analysis.cc index 16f32300f37ea75d7e45962aa128f2c13448657c..e5b365dbe84772d80c01262856538f675e096aa1 100644 --- a/NPAnalysis/Template/src/Analysis.cc +++ b/NPAnalysis/Template/src/Analysis.cc @@ -4,18 +4,12 @@ using namespace std; int main(int argc,char** argv) { - - if(argc!=4) - { - cout << - "you need to specify both a Reaction file and a Detector file such as : Analysis myReaction.reaction myDetector.detector runToRead.run" - << endl; - return 0; - } - - string reactionfileName = argv[1] ; - string detectorfileName = argv[2] ; - string runToReadfileName = argv[3] ; + + NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv) ; + string detectorfileName = myOptionManager->GetDetectorFilePath() ; + string reactionfileName = myOptionManager->GetCalibrationFilePath() ; + string calibrationfileName = myOptionManager->GetCalibrationFilePath() ; + string runToReadfileName = myOptionManager->GetRunToReadFilePath() ; // First of All instantiate RootInput and Output // Detector will be attached later @@ -30,6 +24,9 @@ int main(int argc,char** argv) NPA::DetectorManager* myDetector = new DetectorManager ; myDetector -> ReadConfigurationFile(detectorfileName) ; + // Instantiate the Calibration Manger using a file + CalibrationManager* myCalibration = CalibrationManager::getInstance(calibrationfileName) ; + // Get the formed Chained Tree and Treat it TChain* Chain = RootInput:: getInstance() -> GetChain() ; int i; diff --git a/NPAnalysis/e530/include/ObjectManager.hh b/NPAnalysis/e530/include/ObjectManager.hh index 036e1a34f65df03c02e04424554c6185972b913c..ca9db6818acfce9c028cfbf6628a31cb9474ee64 100644 --- a/NPAnalysis/e530/include/ObjectManager.hh +++ b/NPAnalysis/e530/include/ObjectManager.hh @@ -7,6 +7,7 @@ // NPA #include "DetectorManager.h" +#include "NPOptionManager.h" // STL C++ #include <iostream> diff --git a/NPAnalysis/e530/src/Analysis.cc b/NPAnalysis/e530/src/Analysis.cc index 2561d14e48382f007406e0f1693d80867e5e0aaa..c6558b4d14a62fec56b77ee541a1472fbb37e269 100644 --- a/NPAnalysis/e530/src/Analysis.cc +++ b/NPAnalysis/e530/src/Analysis.cc @@ -4,21 +4,12 @@ using namespace std; int main(int argc,char** argv) { - - if(argc!=5) - { - cout << - "you need to specify both a Reaction file and a Detector file such as : "<< endl; - cout << - "Analysis myReaction.reaction myDetector.detector runToRead.run" << endl; - return 0; - } - - string reactionfileName = argv[1] ; - string detectorfileName = argv[2] ; - string calibrationfileName = argv[3] ; - string runToTreatFileName = argv[4] ; - + NPOptionManager* myOptionManager = NPOptionManager::getInstance(argc,argv) ; + string detectorfileName = myOptionManager->GetDetectorFilePath() ; + string reactionfileName = myOptionManager->GetCalibrationFilePath() ; + string calibrationfileName = myOptionManager->GetCalibrationFilePath() ; + string runToTreatFileName = myOptionManager->GetRunToReadFilePath() ; + ///////////////////////////////////////////////////////////////////////////////////////////////////// // First of All instantiate RootInput and Output // Detector will be attached later diff --git a/NPLib/MUST2/Makefile b/NPLib/MUST2/Makefile index ae103d89957afd3ebd94b76017c51febd7d92c2a..9eae34621d33a888ec754e020fe2b8f0c83787e5 100644 --- a/NPLib/MUST2/Makefile +++ b/NPLib/MUST2/Makefile @@ -288,6 +288,7 @@ TMust2PhysicsDict.cxx: TMust2Physics.h # dependances TMust2Data.o: TMust2Data.cxx TMust2Data.h +TMust2DataDict.o: TMust2DataDict.cxx TMust2DataDict.h TMust2Physics.o: TMust2Physics.cxx TMust2Physics.h ####################################### diff --git a/NPLib/MUST2/TMust2Physics.cxx b/NPLib/MUST2/TMust2Physics.cxx index 6ea861cef697ed75a8c7fa92b31f2d3dd4c77d21..8d48a0e14bdb99bf5411f195abcd11cf6d17a9f4 100644 --- a/NPLib/MUST2/TMust2Physics.cxx +++ b/NPLib/MUST2/TMust2Physics.cxx @@ -9,7 +9,7 @@ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr * * * * Creation Date : febuary 2009 * - * Last update : * + * Last update : october 2010 * *---------------------------------------------------------------------------* * Decription: * * This class hold must2 treated data * @@ -39,141 +39,327 @@ using namespace LOCAL; ClassImp(TMust2Physics) /////////////////////////////////////////////////////////////////////////// TMust2Physics::TMust2Physics() - { + { EventMultiplicity = 0 ; EventData = new TMust2Data ; + PreTreatedData = new TMust2Data ; EventPhysics = this ; NumberOfTelescope = 0 ; + m_MaximumStripMultiplicityAllowed = 0 ; + m_StripEnergyMatchingTolerance = 0 ; +// ReadAnalysisConfig(); } /////////////////////////////////////////////////////////////////////////// void TMust2Physics::BuildSimplePhysicalEvent() { - BuildPhysicalEvent(); + BuildPhysicalEvent(); } /////////////////////////////////////////////////////////////////////////// void TMust2Physics::BuildPhysicalEvent() { + PreTreat(); + bool check_SILI = false ; bool check_CSI = false ; - if( CheckEvent() == 1 ) { vector< TVector2 > couple = Match_X_Y() ; - + EventMultiplicity = couple.size(); for(unsigned int i = 0 ; i < couple.size() ; i++) { check_SILI = false ; check_CSI = false ; - int N = EventData->GetMMStripXEDetectorNbr(couple[i].X()) ; + int N = PreTreatedData->GetMMStripXEDetectorNbr(couple[i].X()) ; - int X = EventData->GetMMStripXEStripNbr(couple[i].X()) ; - int Y = EventData->GetMMStripYEStripNbr(couple[i].Y()) ; + int X = PreTreatedData->GetMMStripXEStripNbr(couple[i].X()) ; + int Y = PreTreatedData->GetMMStripYEStripNbr(couple[i].Y()) ; - double Si_X_E = fSi_X_E(EventData , couple[i].X()) ; - double Si_Y_E = fSi_Y_E(EventData , couple[i].Y()) ; - - double Si_X_T = fSi_X_T(EventData , couple[i].X()) ; - double Si_Y_T = fSi_Y_T(EventData , couple[i].Y()) ; - - Si_X.push_back(X) ; Si_Y.push_back(Y) ; TelescopeNumber.push_back(N) ; + double Si_X_E = PreTreatedData->GetMMStripXEEnergy( couple[i].X() ) ; + double Si_Y_E = PreTreatedData->GetMMStripYEEnergy( couple[i].Y() ) ; - // Take maximum Energy - if(Si_X_E >= Si_Y_E) Si_E.push_back(Si_X_E) ; - else Si_E.push_back(Si_Y_E) ; - // Take minimum Time - if(Si_X_T >= Si_Y_T) Si_T.push_back(Si_Y_T) ; - else Si_T.push_back(Si_X_T) ; + // Search for associate Time + double Si_X_T = -1000 ; + for(unsigned int t = 0 ; t < PreTreatedData->GetMMStripXTMult() ; t++ ) + { + if( PreTreatedData->GetMMStripXTStripNbr( couple[i].X() ) == PreTreatedData->GetMMStripXTStripNbr(t) + ||PreTreatedData->GetMMStripXTDetectorNbr( couple[i].X() ) == PreTreatedData->GetMMStripXTDetectorNbr(t)) + Si_X_T = PreTreatedData->GetMMStripXTTime(t); + } + + double Si_Y_T = -1000 ; + for(unsigned int t = 0 ; t < PreTreatedData->GetMMStripYTMult() ; t++ ) + { + if( PreTreatedData->GetMMStripYTStripNbr( couple[i].Y() ) == PreTreatedData->GetMMStripYTStripNbr(t) + ||PreTreatedData->GetMMStripYTDetectorNbr( couple[i].Y() ) == PreTreatedData->GetMMStripYTDetectorNbr(t)) + Si_Y_T = PreTreatedData->GetMMStripYTTime(t); + } - for(unsigned int j = 0 ; j < EventData->GetMMSiLiEMult() ; j++) + Si_X.push_back(X) ; Si_Y.push_back(Y) ; TelescopeNumber.push_back(N) ; + + // Take X Energy (better resolution, better zero extrapolation) + Si_E.push_back(Si_X_E); + // Take Y Time, better resolution than X. + Si_T.push_back(Si_Y_T) ; + + for(unsigned int j = 0 ; j < PreTreatedData->GetMMSiLiEMult() ; j++) { - if(EventData->GetMMSiLiEDetectorNbr(j)==N) + if(PreTreatedData->GetMMSiLiEDetectorNbr(j)==N) { - // 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) ) ) + if( Match_Si_SiLi( X, Y , PreTreatedData->GetMMSiLiEPadNbr(j) ) ) { - SiLi_N.push_back(EventData->GetMMSiLiEPadNbr(j)) ; - SiLi_E.push_back(fSiLi_E(EventData , j)) ; + SiLi_N.push_back(PreTreatedData->GetMMSiLiEPadNbr(j)) ; + SiLi_E.push_back(PreTreatedData->GetMMSiLiEEnergy(j)) ; - // Look for associate energy + // Look for associate time // Note: in case of use of SiLi "Orsay" time is not coded. - for(int k =0 ; k < EventData->GetMMSiLiTMult() ; k ++) + for(int k =0 ; k < PreTreatedData->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 ;} + if( PreTreatedData->GetMMSiLiEPadNbr(j)==PreTreatedData->GetMMSiLiEPadNbr(k) && PreTreatedData->GetMMSiLiEDetectorNbr(j)==PreTreatedData->GetMMSiLiTDetectorNbr(k) ) + { SiLi_T.push_back( PreTreatedData->GetMMSiLiTTime(k) ) ; break ;} } check_SILI = true ; } - } + } } - for( int j = 0 ; j < EventData->GetMMCsIEMult() ; j++) + for( int j = 0 ; j < PreTreatedData->GetMMCsIEMult() ; j++) { - if(EventData->GetMMCsIEDetectorNbr(j)==N) + + if(PreTreatedData->GetMMCsIEDetectorNbr(j)==N) { - // 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 , PreTreatedData->GetMMCsIECristalNbr(j) ) ) { - CsI_N.push_back(EventData->GetMMCsIECristalNbr(j)) ; - CsI_E.push_back(fCsI_E(EventData , j)) ; + CsI_N.push_back( PreTreatedData->GetMMCsIECristalNbr(j) ) ; + CsI_E.push_back( PreTreatedData->GetMMCsIEEnergy(j) ) ; - for(int k =0 ; k < EventData->GetMMCsITMult() ; k ++) + // Look for associate Time + for(int k =0 ; k < PreTreatedData->GetMMCsITMult() ; k ++) { - if( EventData->GetMMCsIECristalNbr(j)==EventData->GetMMCsITCristalNbr(k) && EventData->GetMMCsIEDetectorNbr(j)==EventData->GetMMCsITDetectorNbr(k) ) - {CsI_T.push_back(fCsI_T(EventData , k)) ; break ;} + // Same Cristal; Same Detector + if( PreTreatedData->GetMMCsIECristalNbr(j)==PreTreatedData->GetMMCsITCristalNbr(k) && PreTreatedData->GetMMCsIEDetectorNbr(j)==PreTreatedData->GetMMCsITDetectorNbr(k) ) + { CsI_T.push_back( PreTreatedData->GetMMCsITTime(j) ) ; break ;} } check_CSI = true ; } - } - } - } - - - if(!check_SILI) + + } + } + + + if(!check_SILI) { - SiLi_N.push_back(0) ; - SiLi_E.push_back(-10000) ; - SiLi_T.push_back(-10000) ; + SiLi_N.push_back(0) ; + SiLi_E.push_back(-1000) ; + SiLi_T.push_back(-1000) ; } if(!check_CSI) { - CsI_N.push_back(0) ; - CsI_E.push_back(-10000) ; - CsI_T.push_back(-10000) ; + CsI_N.push_back(0) ; + CsI_E.push_back(-1000) ; + CsI_T.push_back(-1000) ; } } } - + return; } +/////////////////////////////////////////////////////////////////////////// +void TMust2Physics::PreTreat() + { + ClearPreTreatedData(); + + // X + // E + for(int i = 0 ; i < EventData->GetMMStripXEMult() ; i++) + { + if(IsValidChannel("X", EventData->GetMMStripXEDetectorNbr(i), EventData->GetMMStripXEStripNbr(i))) + { + double EX = fSi_X_E(EventData , i); + if( EX > Si_X_E_Threshold ) + { + PreTreatedData->SetMMStripXEDetectorNbr( EventData->GetMMStripXEDetectorNbr(i) ) ; + PreTreatedData->SetMMStripXEStripNbr( EventData->GetMMStripXEStripNbr(i) ) ; + PreTreatedData->SetMMStripXEEnergy( EX ) ; + } + + } + + + else + { + + } + + } + + // T + for(int i = 0 ; i < EventData->GetMMStripXTMult() ; i++) + { + if(IsValidChannel("X", EventData->GetMMStripXTDetectorNbr(i), EventData->GetMMStripXTStripNbr(i))) + { + PreTreatedData->SetMMStripXTDetectorNbr( EventData->GetMMStripXTDetectorNbr(i) ) ; + PreTreatedData->SetMMStripXTStripNbr( EventData->GetMMStripXTStripNbr(i) ) ; + PreTreatedData->SetMMStripXTTime( fSi_X_T(EventData , i) ) ; + } + + else + { + + } + + } + + + // Y + // E + for(int i = 0 ; i < EventData->GetMMStripYEMult() ; i++) + { + if(IsValidChannel("Y", EventData->GetMMStripYEDetectorNbr(i), EventData->GetMMStripYEStripNbr(i))) + { + double EY = fSi_Y_E(EventData , i); + if( EY > Si_Y_E_Threshold ) + { + PreTreatedData->SetMMStripYEDetectorNbr( EventData->GetMMStripYEDetectorNbr(i) ) ; + PreTreatedData->SetMMStripYEStripNbr( EventData->GetMMStripYEStripNbr(i) ) ; + PreTreatedData->SetMMStripYEEnergy( EY ) ; + } + } + + else + { + + } + + } + + // T + for(int i = 0 ; i < EventData->GetMMStripYTMult() ; i++) + { + if(IsValidChannel("Y", EventData->GetMMStripYTDetectorNbr(i), EventData->GetMMStripYTStripNbr(i))) + { + PreTreatedData->SetMMStripYTDetectorNbr( EventData->GetMMStripYTDetectorNbr(i) ) ; + PreTreatedData->SetMMStripYTStripNbr( EventData->GetMMStripYTStripNbr(i) ) ; + PreTreatedData->SetMMStripYTTime( fSi_Y_T(EventData , i) ) ; + } + + else + { + + } + + } + + + // CsI + // E + for(int i = 0 ; i < EventData->GetMMCsIEMult() ; i++) + { + + if(IsValidChannel("CsI", EventData->GetMMCsIEDetectorNbr(i), EventData->GetMMCsIECristalNbr(i))) + { + double ECsI = fCsI_E(EventData , i); + if( ECsI > CsI_E_Threshold ) + { + PreTreatedData->SetMMCsIEDetectorNbr( EventData->GetMMCsIEDetectorNbr(i) ) ; + PreTreatedData->SetMMCsIECristalNbr( EventData->GetMMCsIECristalNbr(i) ) ; + PreTreatedData->SetMMCsIEEnergy( ECsI ) ; + } + } + + else + { + + } + + } + + // T + for(int i = 0 ; i < EventData->GetMMCsITMult() ; i++) + { + if(IsValidChannel("CsI", EventData->GetMMCsITDetectorNbr(i), EventData->GetMMCsITCristalNbr(i))) + { + PreTreatedData->SetMMCsITDetectorNbr( EventData->GetMMCsITDetectorNbr(i) ) ; + PreTreatedData->SetMMCsITCristalNbr( EventData->GetMMCsITCristalNbr(i) ) ; + PreTreatedData->SetMMCsITTime( fCsI_T(EventData , i) ) ; + } + + else + { + + } + + } + + + // SiLi + // E + for(int i = 0 ; i < EventData->GetMMSiLiEMult() ; i++) + { + if(IsValidChannel("SiLi", EventData->GetMMSiLiEDetectorNbr(i), EventData->GetMMSiLiEPadNbr(i))) + { + double ESiLi = fSiLi_E(EventData , i); + if( ESiLi > SiLi_E_Threshold ) + { + PreTreatedData->SetMMSiLiEDetectorNbr( EventData->GetMMSiLiEDetectorNbr(i) ) ; + PreTreatedData->SetMMSiLiEPadNbr( EventData->GetMMSiLiEPadNbr(i) ) ; + PreTreatedData->SetMMSiLiEEnergy( ESiLi ) ; + } + } + + else + { + + } + + } + + // T + for(int i = 0 ; i < EventData->GetMMSiLiTMult() ; i++) + { + if(IsValidChannel("SiLi", EventData->GetMMSiLiTDetectorNbr(i), EventData->GetMMSiLiTPadNbr(i))) + { + PreTreatedData->SetMMSiLiTDetectorNbr( EventData->GetMMSiLiTDetectorNbr(i) ) ; + PreTreatedData->SetMMSiLiTPadNbr( EventData->GetMMSiLiTPadNbr(i) ) ; + PreTreatedData->SetMMSiLiTTime( fSiLi_T(EventData , i) ) ; + } + + else + { + + } + + } + + + return; + } + + /////////////////////////////////////////////////////////////////////////// int TMust2Physics :: CheckEvent() { // Check the size of the different elements - if( EventData->GetMMStripXEMult() == EventData->GetMMStripYEMult() && EventData->GetMMStripYEMult() == EventData->GetMMStripXTMult() && EventData->GetMMStripXTMult() == EventData->GetMMStripYTMult() ) + if( PreTreatedData->GetMMStripXEMult() == PreTreatedData->GetMMStripYEMult() /*&& PreTreatedData->GetMMStripYEMult() == PreTreatedData->GetMMStripXTMult() && PreTreatedData->GetMMStripXTMult() == PreTreatedData->GetMMStripYTMult()*/ ) return 1 ; // Regular Event - else if( EventData->GetMMStripXEMult() == EventData->GetMMStripYEMult()+1 || EventData->GetMMStripXEMult() == EventData->GetMMStripYEMult()-1 ) + else if( PreTreatedData->GetMMStripXEMult() == PreTreatedData->GetMMStripYEMult()+1 || PreTreatedData->GetMMStripXEMult() == PreTreatedData->GetMMStripYEMult()-1 ) return 2 ; // Pseudo Event, potentially interstrip @@ -194,38 +380,150 @@ vector < TVector2 > TMust2Physics :: Match_X_Y() // 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 ) + if( PreTreatedData->GetMMStripXEMult() > m_MaximumStripMultiplicityAllowed || PreTreatedData->GetMMStripYEMult() > m_MaximumStripMultiplicityAllowed ) return ArrayOfGoodCouple; - for(int i = 0 ; i < EventData->GetMMStripXEMult(); i++) + for(int i = 0 ; i < PreTreatedData->GetMMStripXEMult(); i++) { - // if X value is above threshold, look at Y value - if( fSi_X_E(EventData , i) > Si_X_E_Threshold ) - { - - for(int j = 0 ; j < EventData->GetMMStripYEMult(); j++) + for(int j = 0 ; j < PreTreatedData->GetMMStripYEMult(); j++) { - // if Y value is above threshold look if detector match - if( fSi_Y_E(EventData , j) > Si_Y_E_Threshold ) - { // if same detector check energy - if ( EventData->GetMMStripXEDetectorNbr(i) == EventData->GetMMStripYEDetectorNbr(j) ) + if ( PreTreatedData->GetMMStripXEDetectorNbr(i) == PreTreatedData->GetMMStripYEDetectorNbr(j) ) { // Look if energy match (within 10%) - if( ( fSi_X_E(EventData , i) - fSi_Y_E(EventData , j) ) / fSi_X_E(EventData , i) < 0.1 ) + if( abs(( PreTreatedData->GetMMStripXEEnergy(i) - PreTreatedData->GetMMStripYEEnergy(j) ) / PreTreatedData->GetMMStripXEEnergy(i)) < m_StripEnergyMatchingTolerance/100. ) ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ; } - } } - } } - if( ArrayOfGoodCouple.size() > EventData->GetMMStripXEMult() ) ArrayOfGoodCouple.clear() ; + // Prevent to treat event with ambiguous matchin beetween X and Y + if( ArrayOfGoodCouple.size() > PreTreatedData->GetMMStripXEMult() ) ArrayOfGoodCouple.clear() ; return ArrayOfGoodCouple; } +//////////////////////////////////////////////////////////////////////////// +bool TMust2Physics :: IsValidChannel(string DetectorType, int telescope , int channel) + { + vector<bool>::iterator it ; + if(DetectorType == "X") + return *(XChannelStatus[telescope].begin()+channel); + + else if(DetectorType == "Y") + return *(YChannelStatus[telescope].begin()+channel); + + else if(DetectorType == "SiLi") + return *(SiLiChannelStatus[telescope].begin()+channel); + + else if(DetectorType == "CsI") + return *(CsIChannelStatus[telescope].begin()+channel); + + else return false; + } + +/////////////////////////////////////////////////////////////////////////// +void TMust2Physics::ReadAnalysisConfig() +{ + bool ReadingStatus = false; + bool check_mult = false; + bool check_match = false; + + // path to file + string FileName = "./configs/ConfigMust2.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << " No ConfigMust2.dat found: Default parameter loaded for Analayis " << FileName << endl; + return; + } + cout << " Loading user parameter for Analysis from ConfigMust2.dat " << endl; + + // read analysis config file + string LineBuffer,DataBuffer; + while (!AnalysisConfigFile.eof()) { + // Pick-up next line + getline(AnalysisConfigFile, LineBuffer); + + // search for "header" + if (LineBuffer.compare(0, 11, "ConfigMust2") == 0) ReadingStatus = true; + + // loop on tokens and data + while (ReadingStatus) { + AnalysisConfigFile >> DataBuffer; + + // Search for comment symbol (%) + if (DataBuffer.compare(0, 1, "%") == 0) { + AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' ); + } + else if (DataBuffer.compare(0, 22, "MAX_STRIP_MULTIPLICITY") == 0) { + check_mult = true; + AnalysisConfigFile >> DataBuffer; + m_MaximumStripMultiplicityAllowed = atoi(DataBuffer.c_str() ); + cout << "Maximun strip multiplicity= " << m_MaximumStripMultiplicityAllowed << endl; + } + else if (DataBuffer.compare(0, 31, "STRIP_ENERGY_MATCHING_TOLERANCE") == 0) { + check_match = true; + AnalysisConfigFile >> DataBuffer; + m_StripEnergyMatchingTolerance = atoi(DataBuffer.c_str() ); + cout << "Strip energy matching tolerance= " << m_StripEnergyMatchingTolerance << endl; + } + else if (DataBuffer.compare(0, 5, "MUST2") == 0) { + AnalysisConfigFile >> DataBuffer; + string whatToDo = DataBuffer; + if (whatToDo.compare(0, 11, "DISABLE_ALL") == 0) { + AnalysisConfigFile >> DataBuffer; + cout << whatToDo << " " << DataBuffer << endl; + int telescope = atoi(DataBuffer.substr(2,1).c_str()); + vector< bool > ChannelStatus; + ChannelStatus.resize(128,false); + XChannelStatus[telescope] = ChannelStatus; + YChannelStatus[telescope] = ChannelStatus; + ChannelStatus.resize(16,false); + SiLiChannelStatus[telescope] = ChannelStatus; + CsIChannelStatus[telescope] = ChannelStatus; + } + else if (whatToDo.compare(0, 15, "DISABLE_CHANNEL") == 0) { + AnalysisConfigFile >> DataBuffer; + cout << whatToDo << " " << DataBuffer << endl; + int telescope = atoi(DataBuffer.substr(2,1).c_str()); + int channel = -1; + if (DataBuffer.compare(3,4,"STRX") == 0) { + channel = atoi(DataBuffer.substr(7).c_str()); + *(XChannelStatus[telescope].begin()+channel) = false; + } + else if (DataBuffer.compare(3,4,"STRY") == 0) { + channel = atoi(DataBuffer.substr(7).c_str()); + *(YChannelStatus[telescope].begin()+channel) = false; + } + else if (DataBuffer.compare(3,4,"SILI") == 0) { + channel = atoi(DataBuffer.substr(7).c_str()); + *(SiLiChannelStatus[telescope].begin()+channel) = false; + } + else if (DataBuffer.compare(3,3,"CSI") == 0) { + channel = atoi(DataBuffer.substr(6).c_str()); + *(CsIChannelStatus[telescope].begin()+channel) = false; + } + else { + cout << "Warning: detector type for Must2 unknown!" << endl; + } + } + else { + cout << "Warning: don't know what to do (lost in translation)" << endl; + } + } + else { + ReadingStatus = false; +// cout << "WARNING: Wrong Token Sequence" << endl; + } + } + } +} + /////////////////////////////////////////////////////////////////////////// bool TMust2Physics :: Match_Si_SiLi(int X, int Y , int PadNbr) { @@ -334,107 +632,108 @@ bool TMust2Physics :: Match_Si_SiLi(int X, int Y , int PadNbr) /////////////////////////////////////////////////////////////////////////// bool TMust2Physics :: Match_Si_CsI(int X, int Y , int CristalNbr) - { - if( CristalNbr == 1 - && X<=71 && X>=27 - && Y<=101 && Y>=59 ) + { + if( CristalNbr == 1 + && X<=128 && X>=97 + && Y<=128 && Y>=97 ) return true ; else if( CristalNbr == 2 - && X<=71 && X>=27 - && Y<=128 && Y>=90 ) + && X<=128 && X>=97 + && Y<=96 && Y>=65 ) return true ; else if( CristalNbr == 3 - && X<=35 && X>=1 - && Y<=101 && Y>=59 ) + && X<=128 && X>=97 + && Y<=64 && Y>=33 ) return true ; else if( CristalNbr == 4 - && X<=35 && X>=1 - && Y<=128 && Y>=90 ) + && X<=128 && X>=7 + && Y<=32 && Y>=1 ) return true ; else if( CristalNbr == 5 - && X<=104 && X>=60 - && Y<=71 && Y>=30 ) + && X<=96 && X>=65 + && Y<=32 && Y>=1 ) return true ; else if( CristalNbr == 6 - && X<=104 && X>=60 - && Y<=41 && Y>=1 ) + && X<=96 && X>=65 + && Y<=64 && Y>=33 ) return true ; else if( CristalNbr == 7 - && X<=128 && X>=90 - && Y<=71 && Y>=30 ) + && X<=96 && X>=65 + && Y<=96 && Y>=65 ) return true ; else if( CristalNbr == 8 - && X<=128 && X>=90 - && Y<=41 && Y>=1 ) + && X<=96 && X>=65 + && Y<=128 && Y>=97 ) return true ; else if( CristalNbr == 9 - && X<=71 && X>=27 - && Y<=71 && Y>=40 ) + && X<=64 && X>=33 + && Y<=32 && Y>=1 ) return true ; else if( CristalNbr == 10 - && X<=71 && X>=27 - && Y<=41 && Y>=1 ) + && X<=64 && X>=33 + && Y<=64 && Y>=33 ) return true ; else if( CristalNbr == 11 - && X<=35 && X>=1 - && Y<=71 && Y>=30 ) + && X<=64 && X>=33 + && Y<=96 && Y>=65 ) return true ; else if( CristalNbr == 12 - && X<=35 && X>=1 - && Y<=41 && Y>=1 ) + && X<=64 && X>=33 + && Y<=128 && Y>=97 ) return true ; else if( CristalNbr == 13 - && X<=104 && X>=60 - && Y<=101 && Y>=59 ) + && X<=32 && X>=1 + && Y<=32 && Y>=1 ) return true ; else if( CristalNbr == 14 - && X<=104 && X>=60 - && Y<=128 && Y>=90 ) + && X<=32 && X>=1 + && Y<=64 && Y>=33 ) return true ; else if( CristalNbr == 15 - && X<=128 && X>=90 - && Y<=101 && Y>=59 ) + && X<=32 && X>=1 + && Y<=96 && Y>=65 ) return true ; else if( CristalNbr == 16 - && X<=128 && X>=90 - && Y<=128 && Y>=90 ) + && X<=32 && X>=1 + && Y<=128 && Y>=97 ) return true ; else return false; + } /////////////////////////////////////////////////////////////////////////// @@ -461,14 +760,79 @@ void TMust2Physics::Clear() CsI_E.clear() ; CsI_T.clear() ; CsI_N.clear() ; + + Si_EX.clear() ; + Si_TX.clear() ; + Si_EY.clear() ; + Si_TY.clear() ; + TelescopeNumber_X.clear() ; + TelescopeNumber_Y.clear() ; } /////////////////////////////////////////////////////////////////////////// +void TMust2Physics::ReadCalibrationRun() + { + // X + // E + for(int i = 0 ; i < EventData->GetMMStripXEMult();i++) + { + TelescopeNumber_X.push_back(EventData->GetMMStripXEDetectorNbr(i)); + Si_EX.push_back( fSi_X_E( EventData , i) ) ; + Si_X.push_back(EventData->GetMMStripXEStripNbr(i)); + + } + // T + for(int i = 0 ; i < EventData->GetMMStripXTMult();i++) + { + TelescopeNumber_X.push_back(EventData->GetMMStripXTDetectorNbr(i)); + Si_TX.push_back( fSi_X_T( EventData , i) ) ; + Si_X.push_back(EventData->GetMMStripXTStripNbr(i)); + + } + + // Y + // E + for(int i = 0 ; i < EventData->GetMMStripYEMult();i++) + { + TelescopeNumber_Y.push_back(EventData->GetMMStripYEDetectorNbr(i)); + Si_EY.push_back( fSi_Y_E( EventData , i) ) ; + Si_Y.push_back(EventData->GetMMStripYEStripNbr(i)); + + + } + + // T + for(int i = 0 ; i < EventData->GetMMStripYTMult();i++) + { + TelescopeNumber_Y.push_back(EventData->GetMMStripYTDetectorNbr(i)); + Si_TY.push_back( fSi_Y_T( EventData , i) ) ; + Si_Y.push_back(EventData->GetMMStripYTStripNbr(i)); + + } + + //CsI Energy + for( int j = 0 ; j < EventData->GetMMCsIEMult() ; j++) + { + CsI_N.push_back(EventData->GetMMCsIECristalNbr(j)) ; + CsI_E.push_back(fCsI_E(EventData , j)) ; + + } + + //CsI Time + for( int j = 0 ; j < EventData->GetMMCsITMult() ; j++) + { + //CsI_N.push_back(EventData->GetMMCsITCristalNbr(j)) ; + CsI_T.push_back(fCsI_T(EventData , j)) ; + + } + + } + //// Innherited from VDetector Class //// // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token void TMust2Physics::ReadConfiguration(string Path) -{ +{ ifstream ConfigFile ; ConfigFile.open(Path.c_str()) ; string LineBuffer ; @@ -682,7 +1046,7 @@ void TMust2Physics::ReadConfiguration(string Path) check_C = false ; check_D = false ; - check_Theta = false ; + check_Theta = false ; check_Phi = false ; check_R = false ; check_beta = false ; @@ -690,13 +1054,11 @@ void TMust2Physics::ReadConfiguration(string Path) } } - - - - } + InitializeStandardParameter(); + ReadAnalysisConfig(); - cout << endl << "/////////////////////////////" << endl<<endl; + cout << endl << "/////////////////////////////" << endl << endl; } @@ -728,7 +1090,8 @@ void TMust2Physics::AddParameterToCalibrationManager() Cal->AddParameter("MUST2", "T"+itoa(i+1)+"_CsI"+itoa(j+1)+"_T","MUST2_T"+itoa(i+1)+"_CsI"+itoa(j+1)+"_T") ; } } - + + return; } @@ -737,9 +1100,9 @@ void TMust2Physics::AddParameterToCalibrationManager() void TMust2Physics::InitializeRootInput() { TChain* inputChain = RootInput::getInstance()->GetChain() ; - inputChain->SetBranchStatus( "MUST2" , true ) ; - inputChain->SetBranchStatus( "fMM_*" , true ) ; - inputChain->SetBranchAddress( "MUST2" , &EventData ) ; + inputChain->SetBranchStatus( "MUST2" , true ) ; + inputChain->SetBranchStatus( "fMM_*" , true ) ; + inputChain->SetBranchAddress( "MUST2" , &EventData ) ; } @@ -817,7 +1180,37 @@ void TMust2Physics::AddTelescope( TVector3 C_X1_Y1 , StripPositionZ.push_back( OneTelescopeStripPositionZ ) ; } - + + +void TMust2Physics::InitializeStandardParameter() + { + // Enable all channel + vector< bool > ChannelStatus; + XChannelStatus.clear() ; + YChannelStatus.clear() ; + SiLiChannelStatus.clear() ; + CsIChannelStatus.clear() ; + + ChannelStatus.resize(128,true); + for(int i = 0 ; i < NumberOfTelescope ; i ++) + { + XChannelStatus[i+1] = ChannelStatus; + YChannelStatus[i+1] = ChannelStatus; + } + + ChannelStatus.resize(16,true); + for(int i = 0 ; i < NumberOfTelescope ; i ++) + { + SiLiChannelStatus[i+1] = ChannelStatus; + CsIChannelStatus[i+1] = ChannelStatus; + } + + + m_MaximumStripMultiplicityAllowed = NumberOfTelescope ; + m_StripEnergyMatchingTolerance = 10 ; // % + + return; + } void TMust2Physics::AddTelescope( double theta , double phi , @@ -968,27 +1361,27 @@ namespace LOCAL // X double fSi_X_E(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/" + itoa( EventData->GetMMStripXEDetectorNbr(i) ) + "Si_X" + itoa( EventData->GetMMStripXEStripNbr(i) ) + "_E", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripXEDetectorNbr(i) ) + "_Si_X" + itoa( EventData->GetMMStripXEStripNbr(i) ) + "_E", EventData->GetMMStripXEEnergy(i) ); } double fSi_X_T(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripXTDetectorNbr(i) ) + "Si_X" + itoa( EventData->GetMMStripXTStripNbr(i) ) +"_T", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripXTDetectorNbr(i) ) + "_Si_X" + itoa( EventData->GetMMStripXTStripNbr(i) ) +"_T", EventData->GetMMStripXTTime(i) ); } // Y double fSi_Y_E(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripYEDetectorNbr(i) ) + "Si_Y" + itoa( EventData->GetMMStripYEStripNbr(i) ) +"_E", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripYEDetectorNbr(i) ) + "_Si_Y" + itoa( EventData->GetMMStripYEStripNbr(i) ) +"_E", EventData->GetMMStripYEEnergy(i) ); } double fSi_Y_T(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripYTDetectorNbr(i) ) + "Si_Y" + itoa( EventData->GetMMStripYTStripNbr(i) ) +"_T", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMStripYTDetectorNbr(i) ) + "_Si_Y" + itoa( EventData->GetMMStripYTStripNbr(i) ) +"_T", EventData->GetMMStripYTTime(i) ); } @@ -996,26 +1389,26 @@ namespace LOCAL // SiLi double fSiLi_E(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMSiLiEDetectorNbr(i) ) + "SiLi" + itoa( EventData->GetMMSiLiEPadNbr(i) ) +"_E", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMSiLiEDetectorNbr(i) ) + "_SiLi" + itoa( EventData->GetMMSiLiEPadNbr(i) ) +"_E", EventData->GetMMSiLiEEnergy(i) ); } double fSiLi_T(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMSiLiTDetectorNbr(i) ) + "SiLi" + itoa( EventData->GetMMSiLiTPadNbr(i) )+"_T", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMSiLiTDetectorNbr(i) ) + "_SiLi" + itoa( EventData->GetMMSiLiTPadNbr(i) )+"_T", EventData->GetMMSiLiTTime(i) ); } // CsI double fCsI_E(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMCsIEDetectorNbr(i) ) + "SiLi" + itoa( EventData->GetMMCsIECristalNbr(i) ) +"_E", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMCsIEDetectorNbr(i) ) + "_CsI" + itoa( EventData->GetMMCsIECristalNbr(i) ) +"_E", EventData->GetMMCsIEEnergy(i) ); } double fCsI_T(TMust2Data* EventData , int i) { - return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMCsITDetectorNbr(i) ) + "SiLi" + itoa( EventData->GetMMCsITCristalNbr(i) ) +"_T", + return CalibrationManager::getInstance()->ApplyCalibration( "MUST2/T" + itoa( EventData->GetMMCsITDetectorNbr(i) ) + "_CsI" + itoa( EventData->GetMMCsITCristalNbr(i) ) +"_T", EventData->GetMMCsITTime(i) ); } diff --git a/NPLib/MUST2/TMust2Physics.h b/NPLib/MUST2/TMust2Physics.h index 67cb9c31633b521d3ce59a5f0852baf9a329c7ec..4b497678265dced51b816ce876b53b3035d97eaf 100644 --- a/NPLib/MUST2/TMust2Physics.h +++ b/NPLib/MUST2/TMust2Physics.h @@ -18,8 +18,8 @@ * * *---------------------------------------------------------------------------* * Comment: * - * Only multiplicity one and multiplicity 2 are down. * - * Improvment needed * + * * + * * * * *****************************************************************************/ // STL @@ -64,12 +64,19 @@ class TMust2Physics : public TObject, public NPA::VDetector // Telescope vector<int> TelescopeNumber ; - // Si X - vector<double> Si_E ; - vector<double> Si_T ; - vector<int> Si_X ; - vector<int> Si_Y ; - + // Si + vector<double> Si_E ;//max of Si_EX and Si_EY + vector<double> Si_T ;//min of Si_TX and Si_TY + vector<int> Si_X ; + vector<int> Si_Y ; + + // Use for checking purpose + vector<double> Si_EX ; + vector<double> Si_TX ; + vector<double> Si_EY ; + vector<double> Si_TY ; + vector<int> TelescopeNumber_X ; + vector<int> TelescopeNumber_Y ; // Si(Li) vector<double> SiLi_E ; vector<double> SiLi_T ; @@ -106,7 +113,6 @@ class TMust2Physics : public TObject, public NPA::VDetector // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter. void BuildPhysicalEvent() ; - // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...). // This method aimed to be used for analysis performed during experiment, when speed is requiered. // NB: This method can eventually be the same as BuildPhysicalEvent. @@ -114,23 +120,46 @@ class TMust2Physics : public TObject, public NPA::VDetector // Those two method all to clear the Event Physics or Data void ClearEventPhysics() {Clear();} - void ClearEventData() {EventData->Clear();} + void ClearEventData() {EventData->Clear();} public: // Specific to MUST2 Array + + // Clear The PreTeated object + void ClearPreTreatedData() {PreTreatedData->Clear();} + + // Remove bad channel, calibrate the data and apply threshold + void PreTreat(); + + // Return false if the channel is disabled by user + // Frist argument is either "X","Y","SiLi","CsI" + bool IsValidChannel(string DetectorType, int telescope , int channel); + + + // Initialize the standard parameter for analysis + // ie: all channel enable, maximum multiplicity for strip = number of telescope + void InitializeStandardParameter(); + + // Read the user configuration file; if no file found, load standard one + void ReadAnalysisConfig(); + // Add a Telescope using Corner Coordinate information - void AddTelescope( TVector3 C_X1_Y1 , - TVector3 C_X128_Y1 , - TVector3 C_X1_Y128 , - TVector3 C_X128_Y128 ); + void AddTelescope( TVector3 C_X1_Y1 , + TVector3 C_X128_Y1 , + TVector3 C_X1_Y128 , + TVector3 C_X128_Y128 ); // Add a Telescope using R Theta Phi of Si center information - void AddTelescope( double theta , - double phi , - double distance , - double beta_u , - double beta_v , - double beta_w ); - + void AddTelescope( double theta , + double phi , + double distance , + double beta_u , + double beta_v , + double beta_w ); + + // Use for reading Calibration Run, very simple methods; only apply calibration, no condition + void ReadCalibrationRun(); + + // Use to access the strip position double GetStripPositionX( int N , int X , int Y ) { return StripPositionX[N-1][X-1][Y-1] ; }; double GetStripPositionY( int N , int X , int Y ) { return StripPositionY[N-1][X-1][Y-1] ; }; double GetStripPositionZ( int N , int X , int Y ) { return StripPositionZ[N-1][X-1][Y-1] ; }; @@ -143,14 +172,28 @@ class TMust2Physics : public TObject, public NPA::VDetector double GetEnergyDeposit(int i) { return TotalEnergy[i] ;}; TVector3 GetPositionOfInteraction(int i) ; - TVector3 GetTelescopeNormal(int i) ; + TVector3 GetTelescopeNormal(int i) ; + private: // Parameter used in the analysis + + // Event over this value after pre-treatment are not treated / avoid long treatment time on spurious event + int m_MaximumStripMultiplicityAllowed ;//! + // Give the allowance in percent of the difference in energy between X and Y + double m_StripEnergyMatchingTolerance ; //! + private: // Root Input and Output tree classes TMust2Data* EventData ;//! + TMust2Data* PreTreatedData ;//! TMust2Physics* EventPhysics ;//! + private: // Map of activated channel + map< int, vector<bool> > XChannelStatus ;//! + map< int, vector<bool> > YChannelStatus ;//! + map< int, vector<bool> > SiLiChannelStatus ;//! + map< int, vector<bool> > CsIChannelStatus ;//! + private: // Spatial Position of Strip Calculated on bases of detector position int NumberOfTelescope ;//! @@ -164,12 +207,12 @@ class TMust2Physics : public TObject, public NPA::VDetector namespace LOCAL { - // Threshold - const double Si_X_E_Threshold = 0 ; const double Si_X_T_Threshold = 0 ; - const double Si_Y_E_Threshold = 0 ; const double Si_Y_T_Threshold = 0 ; - const double SiLi_E_Threshold = 0 ; const double SiLi_T_Threshold = 0 ; - const double CsI_E_Threshold = 0 ; const double CsI_T_Threshold = 0 ; - + // Threshold + const double Si_X_E_Threshold = 0 ; + const double Si_Y_E_Threshold = 0 ; + const double SiLi_E_Threshold = 0 ; + const double CsI_E_Threshold = 0 ; + // tranform an integer to a string string itoa(int value); // DSSD diff --git a/NPLib/Makefile b/NPLib/Makefile index 41e9bebb433cd9760fa8882fa8e404fb6e435983..7542417d5c051c374aecbb20d9b81c8a46ea5418 100644 --- a/NPLib/Makefile +++ b/NPLib/Makefile @@ -41,6 +41,7 @@ ROOTLIBS := $(shell $(ROOTCONFIG) --libs) ROOTGLIBS := $(shell $(ROOTCONFIG) --glibs) HASTHREAD := $(shell $(ROOTCONFIG) --has-thread) ROOTDICTTYPE := $(shell $(ROOTCONFIG) --dicttype) +NOSTUBS := $(shell $(ROOTCONFIG) --nostubs) ROOTCINT := rootcint ifeq ($(ARCH),linux) @@ -267,41 +268,189 @@ GLIBS = $(ROOTGLIBS) $(SYSLIBS) INCLUDE = -I$(CLHEP_BASE_DIR)/include #------------------------------------------------------------------------------ -SHARELIB = SharedLib -FILLINC = FillIncludeDir -FILLLIB = FillLibraryDir -LIBLIST = liblistfile +SHARELIB = CalibManager Vdetec InputOutputRoot InitCond InterCoord \ + Must2All GaspardData AnnularS1Data PlasticData DummyDetectorData SSSDData\ + Reaction EnergyLoss ParisData ShieldData -all: $(FILLINC) $(SHARELIB) $(FILLLIB) $(LIBLIST) +all: $(SHARELIB) + rm -f ./include/*Dict.h #------------------------------------------------------------------------------ -############### fillinclib ############## -FillIncludeDir: - ./scripts/fillincdir.sh +############### Calibration ############## -FillLibraryDir: - ./scripts/filllibdir.sh +## CalibrationManager ## +CalibManager: + make -C ./CalibrationManager + cp ./CalibrationManager/*.so ./lib ; cp ./CalibrationManager/*.h ./include ifeq ($(ARCH),macosx) - ./scripts/filllibmacdir.sh + cd lib; ln -sf libCalibrationManager.so libCalibrationManager.dylib endif -############### liblist ############## -liblistfile: - ./scripts/buildliblist.sh +############### Detector ############## -############### sharedlib ############## -SharedLib: - ./scripts/makefile.sh +## VDetector ## +Vdetec: + cp ./VDetector/*.h ./include + make -C ./VDetector + cp ./VDetector/*.so ./lib +ifeq ($(ARCH),macosx) + cd lib; ln -sf libVDetector.so libVDetector.dylib +endif + +## MUST2 ## +Must2All: + make -C ./MUST2 + cp ./MUST2/*.so ./lib ; cp ./MUST2/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libMust2Data.so libMust2Data.dylib + cd lib; ln -sf libMust2Physics.so libMust2Physics.dylib +endif + +## SSSD ## +SSSDData: + make -C ./SSSD + cp ./SSSD/*.so ./lib ; cp ./SSSD/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libSSSDData.so libSSSDData.dylib + cd lib; ln -sf libSSSDPhysics.so libSSSDPhysics.dylib +endif + +## AnnularS1 ## +AnnularS1Data: + make -C ./AnnularS1 + cp ./AnnularS1/*.so ./lib ; cp ./AnnularS1/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libAnnularS1Data.so libAnnularS1Data.dylib +endif + +## Gaspard ## +GaspardData: + make -C ./GASPARD + cp ./GASPARD/*.so ./lib ; cp ./GASPARD/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libGaspardData.so libGaspardData.dylib + cd lib; ln -sf libGaspardPhysics.so libGaspardPhysics.dylib +endif + +## Plastic ## +PlasticData: + make -C ./Plastic + cp ./Plastic/*.so ./lib ; cp ./Plastic/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libPlasticData.so libPlasticData.dylib + cd lib; ln -sf libPlasticPhysics.so libPlasticPhysics.dylib +endif + +## DUMMY Detector ## +DummyDetectorData: + make -C ./DummyDetector + cp ./DummyDetector/*.so ./lib ; cp ./DummyDetector/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libDummyDetectorData.so libDummyDetectorData.dylib +endif + +## Paris Detector ## +ParisData: + make -C ./Paris + cp ./Paris/*.so ./lib ; cp ./Paris/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libParisData.so libParisData.dylib + cd lib; ln -sf libParisPhysics.so libParisPhysics.dylib +endif + + +## Paris Shield Detector ## +ShieldData: + make -C ./Shield + cp ./Shield/*.so ./lib ; cp ./Shield/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libShieldData.so libShieldData.dylib + cd lib; ln -sf libShieldPhysics.so libShieldPhysics.dylib +endif + +############# Simulation ############## + +## InitialConditions ## +InitCond: + make -C ./InitialConditions + cp ./InitialConditions/*.so ./lib ; cp ./InitialConditions/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libInitialConditions.so libInititalConditions.dylib +endif + +## InteractionCoordinates ## +InterCoord: + make -C ./InteractionCoordinates + cp ./InteractionCoordinates/*.so ./lib ; cp ./InteractionCoordinates/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libInteractionCoordinates.so libInteractionCoordinates.dylib +endif + +############# I/O Root File ############ +InputOutputRoot: + make -C ./IORoot + cp ./IORoot/*.so ./lib ; cp ./IORoot/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libIORoot.so libIORoot.dylib +endif + + + +############# Various Tools ############ + +## Reaction ## +Reaction: + make libReaction.so -C ./Tools + cp ./Tools/*.so ./lib ; cp ./Tools/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libReaction.so libReaction.dylib +endif + +EnergyLoss: + make libEnergyLoss.so -C ./Tools + cp ./Tools/*.so ./lib ; cp ./Tools/*.h ./include +ifeq ($(ARCH),macosx) + cd lib; ln -sf libEnergyLoss.so libEnergyLoss.dylib +endif + +####################################### ############# Clean and More ########## clean: - ./scripts/makefile.sh clean + make clean -C ./Tools + make clean -C ./IORoot + make clean -C ./VDetector + make clean -C ./CalibrationManager + make clean -C ./MUST2 + make clean -C ./SSSD + make clean -C ./AnnularS1 + make clean -C ./GASPARD + make clean -C ./InteractionCoordinates + make clean -C ./InitialConditions + make clean -C ./DummyDetector + make clean -C ./Plastic + make clean -C ./Paris + make clean -C ./Shield distclean: - rm -f ./lib/* + rm -f ./lib/*.so +ifeq ($(ARCH),macosx) + rm -f ./lib/*.dylib +endif rm -f ./include/*.h - rm -f liblist - ./scripts/makefile.sh distclean - + make distclean -C ./Tools + make distclean -C ./IORoot + make distclean -C ./VDetector + make distclean -C ./CalibrationManager + make distclean -C ./MUST2 + make distclean -C ./SSSD + make distclean -C ./AnnularS1 + make distclean -C ./GASPARD + make distclean -C ./InteractionCoordinates + make distclean -C ./InitialConditions + make distclean -C ./DummyDetector + make distclean -C ./Plastic + make distclean -C ./Paris + make distclean -C ./Shield .SUFFIXES: .$(SrcSuf) ### diff --git a/NPLib/SSSD/TSSSDPhysics.cxx b/NPLib/SSSD/TSSSDPhysics.cxx index 9ba970a25192fdd0ecb925002a526f0fee21e616..ef7b3ffb88fdf1523dd0e99f7ea9953bdc338a77 100644 --- a/NPLib/SSSD/TSSSDPhysics.cxx +++ b/NPLib/SSSD/TSSSDPhysics.cxx @@ -19,55 +19,58 @@ * * * * *****************************************************************************/ -// NPL +// NPL #include "TSSSDPhysics.h" #include "RootOutput.h" #include "RootInput.h" -// STL +// STL #include <iostream> #include <sstream> #include <fstream> #include <limits> #include <stdlib.h> using namespace std; - -// ROOT +using namespace LOCAL; +// ROOT #include "TChain.h" -// tranform an integer to a string - string itoa(int value) - { - std::ostringstream o; - - if (!(o << value)) - return "" ; - - return o.str(); - } +// tranform an integer to a string + string itoa(int value) + { + std::ostringstream o; + + if (!(o << value)) + return "" ; + + return o.str(); + } ClassImp(TSSSDPhysics) /////////////////////////////////////////////////////////////////////////// TSSSDPhysics::TSSSDPhysics() - { - NumberOfDetector = 0 ; - EventData = new TSSSDData ; - EventPhysics = this ; - } + { + NumberOfDetector = 0 ; + EventData = new TSSSDData ; + PreTreatedData = new TSSSDData ; + EventPhysics = this ; + m_E_Threshold = 0 ; + m_Pedestal_Threshold = 0 ; + } /////////////////////////////////////////////////////////////////////////// TSSSDPhysics::~TSSSDPhysics() - {} + {} /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::Clear() - { - DetectorNumber .clear() ; - StripNumber .clear() ; - Energy .clear() ; - Time .clear() ; - } + { + DetectorNumber .clear() ; + StripNumber .clear() ; + Energy .clear() ; + Time .clear() ; + } /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::ReadConfiguration(string Path) - { + { ifstream ConfigFile ; ConfigFile.open(Path.c_str()) ; string LineBuffer ; @@ -87,218 +90,354 @@ void TSSSDPhysics::ReadConfiguration(string Path) bool ReadingStatus = false ; while (!ConfigFile.eof()) - { + { - getline(ConfigFile, LineBuffer); - - // If line is a Start Up ThinSi bloc, Reading toggle to true - if (LineBuffer.compare(0, 6, "ThinSi") == 0) - { - cout << "Detector found: " << endl ; - ReadingStatus = true ; - } - - // Else don't toggle to Reading Block Status - else ReadingStatus = false ; - - // Reading Block - while(ReadingStatus) - { - // Pickup Next Word - ConfigFile >> DataBuffer ; - - // Comment Line - if (DataBuffer.compare(0, 1, "%") == 0) { ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );} - - // Finding another telescope (safety), toggle out - else if (DataBuffer.compare(0, 6, "ThinSi") == 0) { - cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ; - ReadingStatus = false ; - } - - //Position method - else if (DataBuffer.compare(0, 3, "A=") == 0) { - check_A = true; - ConfigFile >> DataBuffer ; - TLX = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - TLY = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - TLZ = atof(DataBuffer.c_str()) ; - - } - - else if (DataBuffer.compare(0, 3, "B=") == 0) { - check_B = true; - ConfigFile >> DataBuffer ; - BLX = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - BLY = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - BLZ = atof(DataBuffer.c_str()) ; - - } - - else if (DataBuffer.compare(0, 3, "C=") == 0) { - check_C = true; - ConfigFile >> DataBuffer ; - BRX = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - BRY = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - BRZ = atof(DataBuffer.c_str()) ; - - } - - else if (DataBuffer.compare(0, 3, "D=") == 0) { - check_D = true; - ConfigFile >> DataBuffer ; - TRX = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - TRY = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - TRZ = atof(DataBuffer.c_str()) ; - - } - - - //Angle method - else if (DataBuffer.compare(0, 6, "THETA=") == 0) { - check_Theta = true; - ConfigFile >> DataBuffer ; - Theta = atof(DataBuffer.c_str()) ; - } - - else if (DataBuffer.compare(0, 4, "PHI=") == 0) { - check_Phi = true; - ConfigFile >> DataBuffer ; - Phi = atof(DataBuffer.c_str()) ; - } - - else if (DataBuffer.compare(0, 2, "R=") == 0) { - check_R = true; - ConfigFile >> DataBuffer ; - R = atof(DataBuffer.c_str()) ; - } - - - else if (DataBuffer.compare(0, 5, "BETA=") == 0) { - check_beta = true; - ConfigFile >> DataBuffer ; - beta_u = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - beta_v = atof(DataBuffer.c_str()) ; - ConfigFile >> DataBuffer ; - beta_w = atof(DataBuffer.c_str()) ; - } - - /////////////////////////////////////////////////// - // If no Detector Token and no comment, toggle out - else - {ReadingStatus = false; cout << "Wrong Token Sequence: Getting out " << DataBuffer << endl ;} - - ///////////////////////////////////////////////// - // If All necessary information there, toggle out - - if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) ) - { - ReadingStatus = false; - - ///Add The previously define telescope - //With position method - if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R)) { - NumberOfDetector++; - } - - //with angle method - else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D)) { - NumberOfDetector++; - } - - // Reinitialisation of Check Boolean - - check_A = false ; - check_B = false ; - check_C = false ; - check_D = false ; - - check_Theta = false ; - check_Phi = false ; - check_R = false ; - check_beta = false ; - ReadingStatus = false ; - - } - - } - } - + getline(ConfigFile, LineBuffer); + + // If line is a Start Up ThinSi bloc, Reading toggle to true + if (LineBuffer.compare(0, 6, "ThinSi") == 0) + { + cout << "Detector found: " << endl ; + ReadingStatus = true ; + } + + // Else don't toggle to Reading Block Status + else ReadingStatus = false ; + + // Reading Block + while(ReadingStatus) + { + // Pickup Next Word + ConfigFile >> DataBuffer ; + + // Comment Line + if (DataBuffer.compare(0, 1, "%") == 0) { ConfigFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );} + + // Finding another telescope (safety), toggle out + else if (DataBuffer.compare(0, 6, "ThinSi") == 0) { + cout << "WARNING: Another Telescope is find before standard sequence of Token, Error may occured in Telecope definition" << endl ; + ReadingStatus = false ; + } + + //Position method + else if (DataBuffer.compare(0, 3, "A=") == 0) { + check_A = true; + ConfigFile >> DataBuffer ; + TLX = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + TLY = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + TLZ = atof(DataBuffer.c_str()) ; + + } + + else if (DataBuffer.compare(0, 3, "B=") == 0) { + check_B = true; + ConfigFile >> DataBuffer ; + BLX = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + BLY = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + BLZ = atof(DataBuffer.c_str()) ; + + } + + else if (DataBuffer.compare(0, 3, "C=") == 0) { + check_C = true; + ConfigFile >> DataBuffer ; + BRX = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + BRY = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + BRZ = atof(DataBuffer.c_str()) ; + + } + + else if (DataBuffer.compare(0, 3, "D=") == 0) { + check_D = true; + ConfigFile >> DataBuffer ; + TRX = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + TRY = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + TRZ = atof(DataBuffer.c_str()) ; + + } + + + //Angle method + else if (DataBuffer.compare(0, 6, "THETA=") == 0) { + check_Theta = true; + ConfigFile >> DataBuffer ; + Theta = atof(DataBuffer.c_str()) ; + } + + else if (DataBuffer.compare(0, 4, "PHI=") == 0) { + check_Phi = true; + ConfigFile >> DataBuffer ; + Phi = atof(DataBuffer.c_str()) ; + } + + else if (DataBuffer.compare(0, 2, "R=") == 0) { + check_R = true; + ConfigFile >> DataBuffer ; + R = atof(DataBuffer.c_str()) ; + } + + + else if (DataBuffer.compare(0, 5, "BETA=") == 0) { + check_beta = true; + ConfigFile >> DataBuffer ; + beta_u = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + beta_v = atof(DataBuffer.c_str()) ; + ConfigFile >> DataBuffer ; + beta_w = atof(DataBuffer.c_str()) ; + } + + /////////////////////////////////////////////////// + // If no Detector Token and no comment, toggle out + else + {ReadingStatus = false; cout << "Wrong Token Sequence: Getting out " << DataBuffer << endl ;} + + ///////////////////////////////////////////////// + // If All necessary information there, toggle out + + if ( (check_A && check_B && check_C && check_D) || (check_Theta && check_Phi && check_R && check_beta) ) + { + ReadingStatus = false; + + ///Add The previously define telescope + //With position method + if ((check_A && check_B && check_C && check_D) || !(check_Theta && check_Phi && check_R)) { + NumberOfDetector++; + } + + //with angle method + else if ((check_Theta && check_Phi && check_R) || !(check_A && check_B && check_C && check_D)) { + NumberOfDetector++; + } + + // Reinitialisation of Check Boolean + + check_A = false ; + check_B = false ; + check_C = false ; + check_D = false ; + + check_Theta = false ; + check_Phi = false ; + check_R = false ; + check_beta = false ; + ReadingStatus = false ; + + } + + } + } + + InitializeStandardParameter() ; + ReadAnalysisConfig() ; } /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::AddParameterToCalibrationManager() - { - CalibrationManager* Cal = CalibrationManager::getInstance(); - - for(int i = 0 ; i < NumberOfDetector ; i++) - { - - for( int j = 0 ; j < 16 ; j++) - { - Cal->AddParameter("SSSD", "Detector"+itoa(i+1)+"_Strip"+itoa(j+1)+"_E","SSSD_Detector"+itoa(i+1)+"_Strip"+itoa(j+1)+"_E") ; - Cal->AddParameter("SSSD", "Detector"+itoa(i+1)+"_Strip"+itoa(j+1)+"_T","SSSD_Detector"+itoa(i+1)+"_Strip"+itoa(j+1)+"_T") ; - } - - } - } - + { + CalibrationManager* Cal = CalibrationManager::getInstance(); + + for(int i = 0 ; i < NumberOfDetector ; i++) + { + + for( int j = 0 ; j < 16 ; j++) + { + Cal->AddParameter("SSSD", "Detector"+itoa(i+1)+"_Strip"+itoa(j+1)+"_E","SSSD_DETECTOR_"+itoa(i+1)+"_STRIP_"+itoa(j+1)+"_E") ; + Cal->AddParameter("SSSD", "Detector"+itoa(i+1)+"_Strip"+itoa(j+1)+"_T","SSSD_DETECTOR_"+itoa(i+1)+"_STRIP_"+itoa(j+1)+"_T") ; + } + + } + } + /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::InitializeRootInput() - { - TChain* inputChain = RootInput::getInstance()->GetChain() ; - inputChain->SetBranchStatus ( "ThinSi" , true ) ; - inputChain->SetBranchStatus ( "fSSSD_*" , true ) ; - inputChain->SetBranchAddress( "ThinSi" , &EventData ) ; - } + { + TChain* inputChain = RootInput::getInstance()->GetChain() ; + inputChain->SetBranchStatus ( "SSSD" , true ) ; + inputChain->SetBranchStatus ( "fSSSD_*" , true ) ; + inputChain->SetBranchAddress( "SSSD" , &EventData ) ; + } /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::InitializeRootOutput() - { - TTree* outputTree = RootOutput::getInstance()->GetTree() ; - outputTree->Branch( "SSSD" , "TSSSDPhysics" , &EventPhysics ) ; - } + { + TTree* outputTree = RootOutput::getInstance()->GetTree() ; + outputTree->Branch( "SSSD" , "TSSSDPhysics" , &EventPhysics ) ; + } /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::BuildPhysicalEvent() - { - BuildSimplePhysicalEvent() ; - } + { + BuildSimplePhysicalEvent() ; + } /////////////////////////////////////////////////////////////////////////// void TSSSDPhysics::BuildSimplePhysicalEvent() - { - if( EventData->GetEnergyMult() == EventData->GetTimeMult() ) - { - for(unsigned int i = 0 ; i < EventData->GetEnergyMult() ; i++) - { - - DetectorNumber .push_back( EventData->GetEnergyDetectorNbr(i) ) ; - StripNumber .push_back( EventData->GetEnergyStripNbr(i) ) ; - - Energy .push_back( - CalibrationManager::getInstance()->ApplyCalibration( "SSSD/Detector" + itoa( EventData->GetEnergyDetectorNbr(i) ) + "_Strip" + itoa( EventData->GetEnergyStripNbr(i) ) +"_E", - EventData->GetEnergy(i) ) - ) ; - - Time .push_back( - CalibrationManager::getInstance()->ApplyCalibration( "SSSD/Detector" + itoa( EventData->GetEnergyDetectorNbr(i) ) + "_Strip" + itoa( EventData->GetEnergyStripNbr(i) ) +"_T", - EventData->GetTime(i) ) - ) ; - } - } - - else return ; + { + PreTreat(); + + for(unsigned int i = 0 ; i < PreTreatedData->GetEnergyMult() ; i++) + { + + DetectorNumber .push_back( PreTreatedData->GetEnergyDetectorNbr(i) ) ; + StripNumber .push_back( PreTreatedData->GetEnergyStripNbr(i) ) ; + Energy .push_back( PreTreatedData->GetEnergy(i) ) ; + // Look For associate Time + for(unsigned int j = 0 ; j < PreTreatedData->GetTimeMult() ; j ++ ) + { + if(PreTreatedData->GetEnergyDetectorNbr(i) == PreTreatedData->GetTimeDetectorNbr(j) && PreTreatedData->GetEnergyStripNbr(i)==PreTreatedData->GetTimeStripNbr(j)) + Time.push_back(PreTreatedData->GetTime(j)); + } + } + return; + } + +/////////////////////////////////////////////////////////////////////////// +void TSSSDPhysics::PreTreat() + { + ClearPreTreatedData(); + + // E + for(int i = 0 ; i < EventData->GetEnergyMult() ; i++) + { + if(ChannelStatus[EventData->GetEnergyDetectorNbr(i)][EventData->GetEnergyStripNbr(i)]) + { + double E = fSi_E(EventData , i); + if( E > m_E_Threshold && EventData->GetEnergy(i) > m_Pedestal_Threshold) + { + PreTreatedData->SetEnergyDetectorNbr( EventData->GetEnergyDetectorNbr(i) ) ; + PreTreatedData->SetEnergyStripNbr( EventData->GetEnergyStripNbr(i) ) ; + PreTreatedData->SetEnergy( E ) ; + } + } + } - } + // T + for(int i = 0 ; i < EventData->GetTimeMult() ; i++) + { + if(ChannelStatus[EventData->GetEnergyDetectorNbr(i)][EventData->GetEnergyStripNbr(i)]) + { + PreTreatedData->SetTimeDetectorNbr( EventData->GetTimeDetectorNbr(i) ) ; + PreTreatedData->SetTimeStripNbr( EventData->GetTimeStripNbr(i) ) ; + PreTreatedData->SetTime( fSi_T(EventData , i) ) ; + } + } + } +/////////////////////////////////////////////////////////////////////////// +void TSSSDPhysics::InitializeStandardParameter() + { + // Enable all channel + vector<bool> TempChannelStatus; + ChannelStatus.clear(); + TempChannelStatus.resize(16,true); + for(int i = 0 ; i < NumberOfDetector ; i ++) + ChannelStatus[i+1] = TempChannelStatus; + } +/////////////////////////////////////////////////////////////////////////// +void TSSSDPhysics::ReadAnalysisConfig() +{ + bool ReadingStatus = false; + + // path to file + string FileName = "./configs/ConfigSSSD.dat"; + + // open analysis config file + ifstream AnalysisConfigFile; + AnalysisConfigFile.open(FileName.c_str()); + + if (!AnalysisConfigFile.is_open()) { + cout << " No ConfigMust2.dat found: Default parameter loaded for Analayis " << FileName << endl; + return; + } + cout << " Loading user parameter for Analysis from ConfigSSSD.dat " << endl; + + // read analysis config file + string LineBuffer,DataBuffer; + while (!AnalysisConfigFile.eof()) { + // Pick-up next line + getline(AnalysisConfigFile, LineBuffer); + + // search for "header" + if (LineBuffer.compare(0, 11, "ConfigSSSD") == 0) ReadingStatus = true; + + // loop on tokens and data + while (ReadingStatus) { + AnalysisConfigFile >> DataBuffer; + + // Search for comment symbol (%) + if (DataBuffer.compare(0, 1, "%") == 0) { + AnalysisConfigFile.ignore(numeric_limits<streamsize>::max(), '\n' ); + } + + else if (DataBuffer.compare(0, 18, "PEDESTAL_THRESHOLD") == 0) { + AnalysisConfigFile >> DataBuffer; + m_Pedestal_Threshold = atoi(DataBuffer.c_str() ); + cout << "Pedestal threshold = " << m_Pedestal_Threshold << endl; + } + + else if (DataBuffer.compare(0, 4, "SSSD") == 0) { + AnalysisConfigFile >> DataBuffer; + string whatToDo = DataBuffer; + if (whatToDo.compare(0, 11, "DISABLE_ALL") == 0) { + AnalysisConfigFile >> DataBuffer; + cout << whatToDo << " " << DataBuffer << endl; + int Detector = atoi(DataBuffer.substr(2,1).c_str()); + vector< bool > ChannelStatusBuffer; + ChannelStatusBuffer.resize(16,false); + ChannelStatus[Detector] = ChannelStatusBuffer; + } + + else if (whatToDo.compare(0, 15, "DISABLE_CHANNEL") == 0) { + AnalysisConfigFile >> DataBuffer; + cout << whatToDo << " " << DataBuffer << endl; + int telescope = atoi(DataBuffer.substr(2,1).c_str()); + int channel = -1; + if (DataBuffer.compare(3,3,"STR") == 0) { + channel = atoi(DataBuffer.substr(6).c_str()); + *(ChannelStatus[telescope].begin()+channel) = false; + } + + else { + cout << "Warning: detector type for Must2 unknown!" << endl; + } + } + else { + cout << "Warning: don't know what to do (lost in translation)" << endl; + } + } + else { + ReadingStatus = false; +// cout << "WARNING: Wrong Token Sequence" << endl; + } + } + } +} + + + + /////////////////////////////////////////////////////////////////////////// + double LOCAL::fSi_E( TSSSDData* EventData , int i ) + { + return CalibrationManager::getInstance()->ApplyCalibration( "SSSD/Detector" + itoa( EventData->GetEnergyDetectorNbr(i) ) + "_Strip" + itoa( EventData->GetEnergyStripNbr(i) ) +"_E", + EventData->GetEnergy(i) ); + } + + + double LOCAL::fSi_T( TSSSDData* EventData , int i ) + { + return CalibrationManager::getInstance()->ApplyCalibration( "SSSD/Detector" + itoa( EventData->GetEnergyDetectorNbr(i) ) + "_Strip" + itoa( EventData->GetEnergyStripNbr(i) ) +"_T", + EventData->GetTime(i) ); + } + + + diff --git a/NPLib/SSSD/TSSSDPhysics.h b/NPLib/SSSD/TSSSDPhysics.h index fa341dc97e52f53da1a3e3d424a379b76a656cc0..ff10938ec17892b03782e4dce082befe59a14be1 100644 --- a/NPLib/SSSD/TSSSDPhysics.h +++ b/NPLib/SSSD/TSSSDPhysics.h @@ -33,60 +33,90 @@ using namespace std ; #include "../include/VDetector.h" #include "../include/CalibrationManager.h" + class TSSSDPhysics : public TObject, public NPA::VDetector { - public: // Constructor and Destructor - TSSSDPhysics(); - ~TSSSDPhysics(); - - - public: // Calibrated Data + public: // Constructor and Destructor + TSSSDPhysics(); + ~TSSSDPhysics(); - vector<UShort_t> DetectorNumber ; - vector<UShort_t> StripNumber ; - vector<Double_t> Energy ; - vector<Double_t> Time ; + public: // Calibrated Data + vector<UShort_t> DetectorNumber; + vector<UShort_t> StripNumber; + vector<Double_t> Energy; + vector<Double_t> Time; - public: // inherrited from VDetector - // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token - void ReadConfiguration(string) ; + public: // inherrited from VDetector + // Read stream at ConfigFile to pick-up parameters of detector (Position,...) using Token + void ReadConfiguration(string); - - // Add Parameter to the CalibrationManger - void AddParameterToCalibrationManager() ; + // Add Parameter to the CalibrationManger + void AddParameterToCalibrationManager(); - // Activated associated Branches and link it to the private member DetectorData address - // In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated - void InitializeRootInput() ; + // Activated associated Branches and link it to the private member DetectorData address + // In this method mother Branches (Detector) AND daughter leaf (fDetector_parameter) have to be activated + void InitializeRootInput(); + // Create associated branches and associated private member DetectorPhysics address + void InitializeRootOutput(); - // Create associated branches and associated private member DetectorPhysics address - void InitializeRootOutput() ; - + // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter. + void BuildPhysicalEvent(); - // This method is called at each event read from the Input Tree. Aime is to build treat Raw dat in order to extract physical parameter. - void BuildPhysicalEvent() ; + // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...). + // This method aimed to be used for analysis performed during experiment, when speed is requiered. + // NB: This method can eventually be the same as BuildPhysicalEvent. + void BuildSimplePhysicalEvent(); + + // Those two method all to clear the Event Physics or Data + void ClearEventPhysics() {Clear();} + void ClearEventData() {EventData->Clear();} + + + public: // Specific to SSSD + // Clear The PreTeated object + void ClearPreTreatedData() {PreTreatedData->Clear();} - // Same as above, but only the simplest event and/or simple method are used (low multiplicity, faster algorythm but less efficient ...). - // This method aimed to be used for analysis performed during experiment, when speed is requiered. - // NB: This method can eventually be the same as BuildPhysicalEvent. - void BuildSimplePhysicalEvent() ; - - // Those two method all to clear the Event Physics or Data - void ClearEventPhysics() {Clear();} - void ClearEventData() {EventData->Clear();} - - private: // Data not writted in the tree - int NumberOfDetector ;//! - TSSSDData* EventData ;//! - TSSSDPhysics* EventPhysics ;//! - - void Clear(); - void Clear(const Option_t*) {}; + // Remove bad channel, calibrate the data and apply threshold + void PreTreat(); + + // Initialize the standard parameter for analysis + // ie: all channel enable, maximum multiplicity for strip = number of telescope + void InitializeStandardParameter(); + + // Read the user configuration file; if no file found, load standard one + void ReadAnalysisConfig(); + + + private: // Data not written in the tree + int NumberOfDetector; //! + TSSSDData* EventData; //! + TSSSDData* PreTreatedData; //! + TSSSDPhysics* EventPhysics; //! + + double m_E_Threshold; //! + double m_Pedestal_Threshold; //! + + + private: // Map of activated Channel + map< int, vector<bool> > ChannelStatus;//! + + public: // Return True if the channel is activated + // bool IsValidChannel(int DetectorNbr, int StripNbr) ; + + void Clear(); + void Clear(const Option_t*) {}; - ClassDef(TSSSDPhysics,1) // SSSDPhysics structure + ClassDef(TSSSDPhysics,1) // SSSDPhysics structure }; + +namespace LOCAL +{ + double fSi_E( TSSSDData* EventData , int i ); + double fSi_T( TSSSDData* EventData , int i ); +} + #endif diff --git a/NPLib/Tools/Makefile b/NPLib/Tools/Makefile index f3e167bf51ad70fc998a71f5f7c8b8406d890400..54693b981ed50e041752d0a348297251011941cf 100644 --- a/NPLib/Tools/Makefile +++ b/NPLib/Tools/Makefile @@ -267,7 +267,7 @@ GLIBS = $(ROOTGLIBS) $(SYSLIBS) INCLUDE = -I$(CLHEP_BASE_DIR)/include #------------------------------------------------------------------------------ -SHARELIB = libReaction.so libEnergyLoss.so +SHARELIB = libReaction.so libEnergyLoss.so libTagManager.so libOptionManager.so all: $(SHARELIB) #------------------------------------------------------------------------------ @@ -281,10 +281,23 @@ libReaction.so: NPReaction.o NPNucleus.o libEnergyLoss.so: NPEnergyLoss.o $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ +## TAGManager ## +libTagManager.so: NPTagManager.o NPTagManagerDict.o + $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ + +NPTagManagerDict.cxx: NPTagManager.h + rootcint -f $@ -c $^ + +## OptionManager ## +libOptionManager.so: NPOptionManager.o + $(LD) $(SOFLAGS) $^ $(OutPutOpt) $@ + # dependances NPReaction.o:NPReaction.cxx NPReaction.h NPNucleus.o: NPNucleus.cxx NPNucleus.h NPEnergyLoss.o:NPEnergyLoss.cxx NPEnergyLoss.h +NPTagManager.o:NPTagManager.cxx NPTagManager.h +NPOptionManager.o:NPOptionManager.cxx NPOptionManager.h ####################################### ############# Clean and More ########## diff --git a/NPLib/Tools/NPEnergyLoss.cxx b/NPLib/Tools/NPEnergyLoss.cxx index afb7567cd56efcc8c5b0f1551488fbe2b8841a95..85582a9a0266769bffa99da15267e4c9292840c0 100644 --- a/NPLib/Tools/NPEnergyLoss.cxx +++ b/NPLib/Tools/NPEnergyLoss.cxx @@ -64,20 +64,28 @@ EnergyLoss::EnergyLoss(string Path , string Source, int NumberOfSlice=100 , int fNumberOfMass = NumberOfMass ; string globalPath = getenv("NPTOOL"); - Path = globalPath + "/Inputs/EnergyLoss/" + Path; + string StandardPath = globalPath + "/Inputs/EnergyLoss/" + Path; cout << "///////////////////////////////// " << endl ; cout << "Initialising an EnergyLoss object " << endl ; ifstream TableFile ; - TableFile.open(Path.c_str()) ; + TableFile.open(StandardPath.c_str()) ; // Opening dE/dX file - if(!TableFile) cout << "ERROR: TABLE FILE NOT FOUND" << endl; + if(TableFile.is_open()) cout << "Reading Energy Loss File: " << Path << endl ; + // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path. + else + { + TableFile.open( Path.c_str() ); + if(TableFile.is_open()) { cout << "Reading Energy Loss File: " << Path << endl ;} + + else { cout << "ERROR: TABLE FILE NOT FOUND" << endl; return; } + } + if (Source == "G4Table") { - cout << "Reading Energy Loss File: " << Path << endl ; // Reading Data double energy, total; string dummy; @@ -230,25 +238,53 @@ double EnergyLoss::Slow( double Energy , // Energy of the detected particle double Angle ) // Particle Angle const { - TargetThickness = TargetThickness / cos(Angle) ; - double SliceThickness = TargetThickness / (double)fNumberOfSlice ; - - Interpolator* s = new Interpolator( fEnergy , fdEdX_Total ) ; +// // Lise file are given in MeV/u +// // For SRIM and geant4 file fNumberOfMass = 1 whatever is the nucleus, file are given in MeV +// Energy = Energy / (double) fNumberOfMass ; +// +// if (Angle > halfpi) Angle = pi-Angle ; +// +// TargetThickness = TargetThickness / cos(Angle) ; +// double SliceThickness = TargetThickness / (double)fNumberOfSlice ; +// +// //Interpolator* s = new Interpolator( fEnergy , fdEdX_Total ) ; +// +//// double InitialEnergy = Energy ; +//// //double slow = 0. ; +// +// for (int i = 0; i < fNumberOfSlice ; i++) +// { +// // double de = s->Eval(Energy) * SliceThickness; +// double de = fInter->Eval(Energy) * SliceThickness ; +// // slow += de ; +// Energy -= de/fNumberOfMass ; +// // If ion do not cross the target +//// if (Energy < 0) {slow = InitialEnergy; break;} +// if (Energy < 0) {Energy=0; break;} +// } +// +// // delete s ; +// return slow ; + + // Lise file are given in MeV/u + // For SRIM and geant4 file fNumberOfMass = 1 whatever is the nucleus, file are given in MeV + Energy = Energy / (double) fNumberOfMass ; + + if (Angle > halfpi) Angle = pi-Angle ; + TargetThickness = TargetThickness / ( cos(Angle) ) ; - double InitialEnergy = Energy ; - double slow = 0. ; - + double SliceThickness = TargetThickness / (double)fNumberOfSlice ; + for (int i = 0; i < fNumberOfSlice ; i++) { - double de = s->Eval(Energy) * SliceThickness; - slow += de ; - Energy -= de ; - // If ion do not cross the target - if (Energy < 0) {slow = InitialEnergy; break;} + double de = fInter->Eval(Energy) * SliceThickness ; + Energy -= de/fNumberOfMass ; + + if(Energy<0) {Energy=0;break;} } - - delete s ; - return slow ; + + return (Energy*fNumberOfMass) ; + } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... double EnergyLoss::EvaluateInitialEnergy( double Energy , // Energy of the detected particle diff --git a/NPLib/Tools/NPOptionManager.cxx b/NPLib/Tools/NPOptionManager.cxx new file mode 100644 index 0000000000000000000000000000000000000000..18cf752ae7786dd1b899c2c03c314dee549c1fcb --- /dev/null +++ b/NPLib/Tools/NPOptionManager.cxx @@ -0,0 +1,79 @@ +/***************************************************************************** + * Copyright (C) 2009-2010 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: A. MATTA contact address: matta@ipno.in2p3.fr * + * * + * Creation Date : 21/07/09 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class is a singleton class which deals with input * + * arguments of the different NPTool programm (NPS and NPA) * + *---------------------------------------------------------------------------* + * Comment: The singleton form allow users to call the object from anywhere * + * in the code * + * * + * * + *****************************************************************************/ + +#include "NPOptionManager.h" + + +NPOptionManager* NPOptionManager::instance = 0 ; + +NPOptionManager* NPOptionManager::getInstance(int argc,char** argv) + { + if( instance == 0 ) instance = new NPOptionManager(argc,argv); + + return instance ; + } + + +////////////////////////////////////////////////////// +NPOptionManager::NPOptionManager(int argc,char** argv) + { + // Default Setting + fReactionFileName = "myReaction.reaction" ; + fDetectorFileName = "myDetector.detector" ; + fOutputFileName = "myResult.root" ; + fRunToReadFileName = "RunToRead.txt" ; + fCalibrationFileName = "Calibration.txt" ; + + for (int i = 0 ; i < argc ; i++) + { + string argument = argv[i]; + + if(argument == "-H" || argument == "-h" || argument == "-help") + DisplayHelp(); + + else if(argument == "-R" && argc>=i+1 ) fReactionFileName = argv[i+1] ; + + else if(argument == "-C" && argc>=i+1 ) fDetectorFileName = argv[i+1] ; + + else if(argument == "-O" && argc>=i+1 ) fOutputFileName = argv[i+1] ; + + else if(argument == "-run" && argc>=i+1 ) fRunToReadFileName = argv[i+1] ; + + else if(argument == "-cal" && argc>=i+1 ) fCalibrationFileName = argv[i+1] ; + + else ; + } + + } + +/////////////////////////////////////////////////// +void NPOptionManager::DisplayHelp() + { + cout << "----NPOptionManager Help----" << endl ; + cout << "List of Option " << endl ; + cout << "\t -C <arg>\t \t \t \t \t Set arg as the detector configuration file" << endl ; + cout << "\t -cal <arg>\t \t \t \t \t Set arg as the calibration file list" << endl ; + cout << "\t -H -h -help\t \t \t \t \t Display this help message" << endl ; + cout << "\t -O <arg>\t \t \t \t \t Set arg as the Output File Name (output tree)" << endl ; + cout << "\t -run <arg>\t \t \t \t \t Set arg as the run to read file list" << endl ; + cout << endl << endl ; + } diff --git a/NPLib/Tools/NPOptionManager.h b/NPLib/Tools/NPOptionManager.h new file mode 100644 index 0000000000000000000000000000000000000000..1c638a3e4b2842aac107470fc0164a8a8bd11c0e --- /dev/null +++ b/NPLib/Tools/NPOptionManager.h @@ -0,0 +1,80 @@ +#ifndef NPOPTIONMANAGER_HH +#define NPOPTIONMANAGER_HH + +/***************************************************************************** + * Copyright (C) 2009-2010 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * Original Author: A. MATTA contact address: matta@ipno.in2p3.fr * + * * + * Creation Date : 21/07/09 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: This class is a singleton class which deals with input * + * arguments of the different NPTool programm (NPS and NPA) * + *---------------------------------------------------------------------------* + * Comment: The singleton form allow users to call the object from anywhere * + * in the code * + * * + * * + *****************************************************************************/ + +// STL headers +#include <iostream> +#include <string> +using namespace std; + +class NPOptionManager + { + public: + // The analysis class is designed to be a singleton (i.e. only one instance + // can exist). A member function called Instance is defined, which allows + // the user to get a pointer to the existing instance or to create it if + // it does not yet exist: + // (see the constructor for an explanation of the arguments) + static NPOptionManager* getInstance(int argc=0,char** argv=NULL); + + // The analysis class instance can be deleted by calling the Destroy + // method (NOTE: The class destructor is protected, and can thus not be + // called directly): + static void Destroy(); + + protected: + // Constructor (protected) + NPOptionManager(int argc,char** argv); + + // Destructor (protected) + ~NPOptionManager() {}; + + // Prevent copying + NPOptionManager(const NPOptionManager& only); + const NPOptionManager& operator=(const NPOptionManager& only); + + private: + // The static instance of the NPOptionManager class: + static NPOptionManager* instance; + + private: + void DisplayHelp(); + + public: + string GetReactionFilePath() { return fReactionFileName ; } ; + string GetDetectorFilePath() { return fDetectorFileName ; } ; + string GetRunToReadFilePath() { return fRunToReadFileName ; } ; + string GetCalibrationFilePath() { return fCalibrationFileName ; } ; + string GetOutputFilePath() { return fOutputFileName ; } ; + + private: + string fReactionFileName ; + string fDetectorFileName ; + string fRunToReadFileName ; + string fCalibrationFileName ; + string fOutputFileName ; + + }; + +#endif diff --git a/NPLib/Tools/NPReaction.cxx b/NPLib/Tools/NPReaction.cxx index 140776985584bc51f8edf0efd4e200a838f9108f..1697f58285b51eae16cda769179f34954f3a5b4a 100644 --- a/NPLib/Tools/NPReaction.cxx +++ b/NPLib/Tools/NPReaction.cxx @@ -92,12 +92,20 @@ void Reaction::SetEveryThing(string name1, string name2, string name3, string na ///Read the differential cross section string GlobalPath = getenv("NPTOOL"); - Path = GlobalPath + "/Inputs/CrossSection/" + Path; + string StandardPath = GlobalPath + "/Inputs/CrossSection/" + Path; ifstream CSFile; - CSFile.open( Path.c_str() ); + CSFile.open( StandardPath.c_str() ); - if(CSFile.is_open()) { cout << "Reading Cross Section File " << Path << endl;} - else {cout << "Cross Section File " << Path << " not found" << endl;return;} + if(CSFile.is_open()) cout << "Reading Cross Section File " << Path << endl; + + // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path. + else + { + CSFile.open( Path.c_str() ); + if(CSFile.is_open()) { cout << "Reading Cross Section File " << Path << endl;} + + else {cout << "Cross Section File " << Path << " not found" << endl;return;} + } double CSBuffer,AngleBuffer; vector<double> CrossSectionBuffer ; @@ -262,15 +270,19 @@ void Reaction::ReadConfigurationFile(string Path) ////////////////////////////////////////////////////////////////////////////////////////// ifstream ReactionFile; string GlobalPath = getenv("NPTOOL"); - Path = GlobalPath + "/Inputs/EventGenerator/" + Path; - ReactionFile.open(Path.c_str()); + string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path; + ReactionFile.open(StandardPath.c_str()); if (ReactionFile.is_open()) {cout << "Reading Reaction File " << Path << endl ;} - else - { - cout << "Reaction File " << Path << " not Found! " << endl ; - return; - } + + // In case the file is not found in the standard path, the programm try to interpret the file name as an absolute or relative file path. + else + { + ReactionFile.open( Path.c_str() ); + if(ReactionFile.is_open()) { cout << "Reading Reaction File " << Path << endl;} + + else {cout << "Reaction File " << Path << " not found" << endl;return;} + } while (!ReactionFile.eof()) { //Pick-up next line diff --git a/NPLib/Tools/NPTagManager.cxx b/NPLib/Tools/NPTagManager.cxx new file mode 100644 index 0000000000000000000000000000000000000000..51be7299bfb19d5ddd119e73dcb1233dce37884a --- /dev/null +++ b/NPLib/Tools/NPTagManager.cxx @@ -0,0 +1,90 @@ +/***************************************************************************** + * Copyright (C) 2009 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * * + * Original Author : Adrien MATTA contact address: matta@ipno.in2p3.fr * + * * + * Creation Date : November 2010 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class will held a set of string, that can be used as a TAG manager * + * Users can write macro and add different TAG to that object based on users * + * condition. Then the TAG branch can be open and close alone to select event* + * without loading the whole tree. * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ + +#include "NPTagManager.h" + +ClassImp(NPTagManager) + +//////////////////////////////////////////// +NPTagManager::NPTagManager() {} +NPTagManager::~NPTagManager(){} + + +//////////////////////////////////////////// +bool NPTagManager::Is(string condition) + { + // return True is the element is find, false other wise + return !( fTAG.find(condition)==fTAG.end() ); + } +//////////////////////////////////////////// +void NPTagManager::AddCondition(string condition) + { + fTAG.insert(condition); + } + +//////////////////////////////////////////// +void NPTagManager::PrintCondition() + { + set<string>::iterator it ; + + cout << "------------------ Event Condition ------------------" << endl ; + + for ( it=fTAG.begin() ; it!=fTAG.end() ; it++ ) + { + cout << " " << *it << endl ; + } + + cout << "-------------------------------------------------------" << endl ; + + } + +//////////////////////////////////////////// +void NPTagManager::PrintConditionToFile(string filename) + { + + ofstream file; + file.open(filename.c_str()); + + if(!file) {cout << "Warning: file " << filename << " not found " << endl ; return ;} + + else + { + set<string>::iterator it ; + + file << "------------------ Event Condition ------------------" << endl ; + + for ( it=fTAG.begin() ; it!=fTAG.end() ; it++ ) + { + file << " " << *it << endl ; + } + + file << "-------------------------------------------------------" << endl ; + + } + + } + + diff --git a/NPLib/Tools/NPTagManager.h b/NPLib/Tools/NPTagManager.h new file mode 100644 index 0000000000000000000000000000000000000000..8e12aef99ebc02d9e720eff5a029c0bdc324a99a --- /dev/null +++ b/NPLib/Tools/NPTagManager.h @@ -0,0 +1,67 @@ +#ifndef _TAG_ +#define _TAG_ +/***************************************************************************** + * Copyright (C) 2009 this file is part of the NPTool Project * + * * + * For the licensing terms see $NPTOOL/Licence/NPTool_Licence * + * For the list of contributors see $NPTOOL/Licence/Contributors * + *****************************************************************************/ + +/***************************************************************************** + * * + * Original Author : Adrien MATTA contact address: matta@ipno.in2p3.fr * + * * + * Creation Date : November 2010 * + * Last update : * + *---------------------------------------------------------------------------* + * Decription: * + * This class will held a set of string, that can be used as a TAG manager * + * Users can write macro and add different TAG to that object based on users * + * condition. Then the TAG branch can be open and close alone to select event* + * without loading the whole tree. * + * * + *---------------------------------------------------------------------------* + * Comment: * + * * + * * + *****************************************************************************/ +// C++ header +#include <iostream> +#include <fstream> +#include <string> +#include <set> +using namespace std; + +// ROOT header +#include "TObject.h" + +class NPTagManager : public TObject + { + public: + NPTagManager(); + ~NPTagManager(); + + private: + set<string> fTAG; + + public: + // Return True if condition is in the set fTAG + bool Is(string condition); + + // Add condition to the set fTAG + void AddCondition(string condition); + + // Print all the condition that exist in fTAG + void PrintCondition(); + inline void Print() {PrintCondition();}; + + // Print to File filename all the condition that exist in fTAG + void PrintConditionToFile(string filename); + + // Clear all the fTAG + inline void Clear() {fTAG.clear();} ; + + ClassDef(NPTagManager,1) + }; + +#endif diff --git a/NPLib/VDetector/DetectorManager.cxx b/NPLib/VDetector/DetectorManager.cxx index a1cd19be3844e516e5a678780bf621630e46bb56..3041dcf8398ba37b8692e50ff2acc43441c9096c 100644 --- a/NPLib/VDetector/DetectorManager.cxx +++ b/NPLib/VDetector/DetectorManager.cxx @@ -49,18 +49,27 @@ void DetectorManager::ReadConfigurationFile(string Path) ////////////////////////////////////////////////////////////////////////////////////////// string GlobalPath = getenv("NPTOOL"); - Path = GlobalPath + "/Inputs/DetectorConfiguration/" + Path; + string StandardPath = GlobalPath + "/Inputs/DetectorConfiguration/" + Path; ifstream ConfigFile; - ConfigFile.open(Path.c_str()); + ConfigFile.open(StandardPath.c_str()); - if (ConfigFile.is_open()) { + if (ConfigFile.is_open()) + { cout << "/////////////////////////////" << endl; cout << " Configuration file " << Path << " loading " << endl; - } - else { - cout << " Error, no configuration file" << Path << " found" << endl; - return; - } + Path = StandardPath; + } + + else + { + ConfigFile.open( Path.c_str() ); + if(ConfigFile.is_open()) { + cout << "/////////////////////////////" << endl; + cout << " Configuration file " << Path << " loading " << endl; + } + + else {cout << "Configuration File " << Path << " not found" << endl;return;} + } while (!ConfigFile.eof()) { diff --git a/NPSimulation/Simulation.cc b/NPSimulation/Simulation.cc index 1c4d652cb23cfb52e65b9a18a897a16699146db0..0a49c6b354ba7bb22383cb046cae93fe782afdfb 100644 --- a/NPSimulation/Simulation.cc +++ b/NPSimulation/Simulation.cc @@ -13,103 +13,100 @@ #ifdef G4VIS_USE #include "G4VisExecutive.hh" #endif -// NPS Source +// NPS headers #include "EventAction.hh" #include "VDetector.hh" +//NPL headers +#include "NPOptionManager.h" #include "RootOutput.h" -// ROOT Source +// ROOT headers #include "TTree.h" #include "TFile.h" // STL #include <vector> int main(int argc, char** argv) -{ - - if (argc != 3) { - cout << "you need to specify both a Reaction file and a Detector file such as : Simulation myReaction.reaction myDetector.detector" << endl ; - return 0; - } - - // Getting arguments - G4String EventGeneratorFileName = argv[1]; - G4String DetectorFileName = argv[2]; - - //my Verbose output class - G4VSteppingVerbose::SetInstance(new SteppingVerbose); - - //Construct the default run manager - G4RunManager* runManager = new G4RunManager; - - //set mandatory initialization classes - DetectorConstruction* detector = new DetectorConstruction(); - runManager->SetUserInitialization(detector); - - PhysicsList* physics = new PhysicsList(); - runManager->SetUserInitialization(physics); - PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector); - - //Initialize Geant4 kernel - runManager->Initialize(); - physics->MyOwnConstruction(); - - /////////////////////////////////////////////////////////////// - ///////////////// Initializing the Root Output //////////////// - /////////////////////////////////////////////////////////////// - RootOutput::getInstance("Simulation/mySimul"); - - /////////////////////////////////////////////////////////////// - ////////////// Reading Detector Configuration ///////////////// - /////////////////////////////////////////////////////////////// - detector->ReadConfigurationFile(DetectorFileName); - - /////////////////////////////////////////////////////////////// - ////////////////////// Reading Reaction /////////////////////// - /////////////////////////////////////////////////////////////// - primary->ReadEventGeneratorFile(EventGeneratorFileName); - runManager->SetUserAction(primary); - - /////////////////////////////////////////////////////////////// - ////////////////// Starting the Event Action ////////////////// - /////////////////////////////////////////////////////////////// - EventAction* event_action = new EventAction() ; - event_action->SetDetector(detector) ; - runManager->SetUserAction(event_action) ; - - /////////////////////////////////////////////////////////////// - /////// Get the pointer to the User Interface manager //////// - /////////////////////////////////////////////////////////////// - G4UImanager* UI = G4UImanager::GetUIpointer(); - - /////////////////////////////////////////////////////////////// - /////////// Define UI terminal for interactive mode /////////// - /////////////////////////////////////////////////////////////// -#ifdef G4VIS_USE - G4VisManager* visManager = new G4VisExecutive; - visManager->Initialize(); -#endif - - G4UIsession* session = 0; - -#ifdef G4UI_USE_TCSH - session = new G4UIterminal(new G4UItcsh); -#else - session = new G4UIterminal(); -#endif - - UI->ApplyCommand("/control/execute vis.mac"); - session->SessionStart(); - delete session; - -#ifdef G4VIS_USE - delete visManager; -#endif - - /////////////////////////////////////////////////////////////// - ////////////////////// Job termination //////////////////////// - /////////////////////////////////////////////////////////////// - RootOutput::getInstance()->Destroy(); - - delete runManager; - return 0; -} + { + NPOptionManager* OptionManager = NPOptionManager::getInstance(argc,argv); + G4String EventGeneratorFileName = OptionManager->GetReactionFilePath(); + G4String DetectorFileName = OptionManager->GetDetectorFilePath(); + + + //my Verbose output class + G4VSteppingVerbose::SetInstance(new SteppingVerbose); + + //Construct the default run manager + G4RunManager* runManager = new G4RunManager; + + //set mandatory initialization classes + DetectorConstruction* detector = new DetectorConstruction(); + runManager->SetUserInitialization(detector); + + PhysicsList* physics = new PhysicsList(); + runManager->SetUserInitialization(physics); + PrimaryGeneratorAction* primary = new PrimaryGeneratorAction(detector); + + //Initialize Geant4 kernel + runManager->Initialize(); + physics->MyOwnConstruction(); + + /////////////////////////////////////////////////////////////// + ///////////////// Initializing the Root Output //////////////// + /////////////////////////////////////////////////////////////// + RootOutput::getInstance("Simulation/"+OptionManager->GetOutputFilePath()); + + /////////////////////////////////////////////////////////////// + ////////////// Reading Detector Configuration ///////////////// + /////////////////////////////////////////////////////////////// + detector->ReadConfigurationFile(DetectorFileName); + + /////////////////////////////////////////////////////////////// + ////////////////////// Reading Reaction /////////////////////// + /////////////////////////////////////////////////////////////// + primary->ReadEventGeneratorFile(EventGeneratorFileName); + runManager->SetUserAction(primary); + + /////////////////////////////////////////////////////////////// + ////////////////// Starting the Event Action ////////////////// + /////////////////////////////////////////////////////////////// + EventAction* event_action = new EventAction() ; + event_action->SetDetector(detector) ; + runManager->SetUserAction(event_action) ; + + /////////////////////////////////////////////////////////////// + /////// Get the pointer to the User Interface manager //////// + /////////////////////////////////////////////////////////////// + G4UImanager* UI = G4UImanager::GetUIpointer(); + + /////////////////////////////////////////////////////////////// + /////////// Define UI terminal for interactive mode /////////// + /////////////////////////////////////////////////////////////// + #ifdef G4VIS_USE + G4VisManager* visManager = new G4VisExecutive; + visManager->Initialize(); + #endif + + G4UIsession* session = 0; + + #ifdef G4UI_USE_TCSH + session = new G4UIterminal(new G4UItcsh); + #else + session = new G4UIterminal(); + #endif + + UI->ApplyCommand("/control/execute vis.mac"); + session->SessionStart(); + delete session; + + #ifdef G4VIS_USE + delete visManager; + #endif + + /////////////////////////////////////////////////////////////// + ////////////////////// Job termination //////////////////////// + /////////////////////////////////////////////////////////////// + RootOutput::getInstance()->Destroy(); + + delete runManager; + return 0; + } diff --git a/NPSimulation/src/DetectorConstruction.cc b/NPSimulation/src/DetectorConstruction.cc index cbe0026f6e6fc155eea2638102c82c18b4ca7883..8ca7279ed95935fcf913c79e8160ee61b014da84 100644 --- a/NPSimulation/src/DetectorConstruction.cc +++ b/NPSimulation/src/DetectorConstruction.cc @@ -170,17 +170,25 @@ void DetectorConstruction::ReadConfigurationFile(string Path) ////////////////////////////////////////////////////////////////////////////////////////// // added by Nicolas [07/05/09] string GlobalPath = getenv("NPTOOL"); - Path = GlobalPath + "/Inputs/DetectorConfiguration/" + Path; + string StandardPath = GlobalPath + "/Inputs/DetectorConfiguration/" + Path; ifstream ConfigFile; - ConfigFile.open(Path.c_str()); + ConfigFile.open(StandardPath.c_str()); if (ConfigFile.is_open()) - cout << " Configuration file " << Path << " loading " << endl; - else { - cout << " Error, no configuration file" << Path << " found" << endl; - return; - } + { + cout << " Configuration file " << Path << " loading " << endl; + Path=StandardPath; + } + else + { + ConfigFile.open( Path.c_str() ); + if(ConfigFile.is_open()) { + cout << " Configuration file " << Path << " loading " << endl; + } + + else { cout << " Error, no configuration file" << Path << " found" << endl;return;} + } while (!ConfigFile.eof()) { //Pick-up next line diff --git a/NPSimulation/src/GaspardTrackerDummyShape.cc b/NPSimulation/src/GaspardTrackerDummyShape.cc index 4c4f1c79195cc4fc5767704dd21ef77dbf26bbd9..e4dcf07c007ae19acd419c3172edefb1ba5fe622 100644 --- a/NPSimulation/src/GaspardTrackerDummyShape.cc +++ b/NPSimulation/src/GaspardTrackerDummyShape.cc @@ -555,7 +555,6 @@ void GaspardTrackerDummyShape::ConstructDetector(G4LogicalVolume* world) MMv = m_X1_Y128[i] - m_X1_Y1[i]; MMv = MMv.unit(); - MMw = MMu.cross(MMv); MMw = MMw.unit(); diff --git a/NPSimulation/src/GaspardTrackerSquare.cc b/NPSimulation/src/GaspardTrackerSquare.cc index c8ba9e5a9983ff4f46c6aaa2a455dc8973b2f508..9d820d70293acc8f5e05ac9c6570a3b62185eecd 100644 --- a/NPSimulation/src/GaspardTrackerSquare.cc +++ b/NPSimulation/src/GaspardTrackerSquare.cc @@ -111,14 +111,14 @@ void GaspardTrackerSquare::AddModule(G4ThreeVector X1_Y1 , //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void GaspardTrackerSquare::AddModule(G4double R , - G4double Theta , - G4double Phi , - G4double beta_u , - G4double beta_v , - G4double beta_w , - bool wFirstStage , - bool wSecondStage , - bool wThirdStage) + G4double Theta , + G4double Phi , + G4double beta_u , + G4double beta_v , + G4double beta_w , + bool wFirstStage , + bool wSecondStage , + bool wThirdStage) { G4ThreeVector empty = G4ThreeVector(0, 0, 0); @@ -188,7 +188,7 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber, // Al density = 2.702 * g / cm3; a = 26.98 * g / mole; -// G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density); + G4Material* Aluminium = new G4Material("Aluminium", z = 13., a, density); // Iron // density = 7.874 * g / cm3; @@ -228,15 +228,37 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber, // Little trick to avoid warning in compilation: Use a PVPlacement "buffer". // If don't you will have a Warning unused variable 'myPVP' G4PVPlacement* PVPBuffer ; - G4String Name = "GPDSquare" + DetectorNumber; - G4Box* solidGPDSquare = new G4Box(Name, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length); - G4LogicalVolume* logicGPDSquare = new G4LogicalVolume(solidGPDSquare, Vacuum, Name, 0, 0, 0); + G4Trd* solidMM = new G4Trd("GPDSquare" + DetectorNumber, 0.5*FaceFront, 0.5*FaceFront, 0.5*FaceFront, 0.5*FaceFront, 0.5*Length); +// G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, Iron, "GPDSquare" + DetectorNumber, 0, 0, 0) ; + G4LogicalVolume* logicMM = new G4LogicalVolume(solidMM, Vacuum, "GPDSquare" + DetectorNumber, 0, 0, 0) ; + + G4String Name = "GPDSquare" + DetectorNumber ; + PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot, MMpos) , + logicMM , + Name , + world , + false , + 0); + + logicMM->SetVisAttributes(G4VisAttributes::Invisible); + if (m_non_sensitive_part_visiualisation) logicMM->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90))); + + G4ThreeVector positionVacBox = G4ThreeVector(0, 0, VacBox_PosZ); - PVPBuffer = new G4PVPlacement(G4Transform3D(*MMrot, MMpos), logicGPDSquare, Name, world, false, 0); + G4Trd* solidVacBox = new G4Trd("solidVacBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconFace, 0.5*VacBoxThickness); + G4LogicalVolume* logicVacBox = new G4LogicalVolume(solidVacBox, Vacuum, "logicVacBox", 0, 0, 0); - logicGPDSquare->SetVisAttributes(G4VisAttributes::Invisible); - if (m_non_sensitive_part_visiualisation) logicGPDSquare->SetVisAttributes(G4VisAttributes(G4Colour(0.90, 0.90, 0.90))); + PVPBuffer = new G4PVPlacement(0, positionVacBox, logicVacBox, "G" + DetectorNumber + "VacBox", logicMM, false, 0); + + logicVacBox->SetVisAttributes(G4VisAttributes::Invisible); + + // Add a degrader plate between Si and CsI: + /* + G4Box* Degrader = new G4Box("Degrader" , 50*mm , 50*mm , 0.1*mm ); + G4LogicalVolume* logicDegrader = new G4LogicalVolume( Degrader , Harvar, "logicDegrader",0,0,0); + PVPBuffer = new G4PVPlacement(0,G4ThreeVector(0,0,0),logicDegrader,"Degrader",logicVacBox,false,0) ; + */ //Place two marker to identify the u and v axis on silicon face: //marker are placed a bit before the silicon itself so they don't perturbate simulation @@ -262,55 +284,166 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber, */ //////////////////////////////////////////////////////////////// - //////////////// First Stage Construction ////////////////////// + /////////////////Si Strip Construction////////////////////////// //////////////////////////////////////////////////////////////// if (wFirstStage) { + // Aluminium dead layers + G4ThreeVector positionAluStripFront = G4ThreeVector(0, 0, AluStripFront_PosZ); + G4ThreeVector positionAluStripBack = G4ThreeVector(0, 0, AluStripBack_PosZ); + + G4Box* solidAluStrip = new G4Box("AluBox", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*AluStripThickness); +// G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, Aluminium, "logicAluStrip", 0, 0, 0); + G4LogicalVolume* logicAluStrip = new G4LogicalVolume(solidAluStrip, Vacuum, "logicAluStrip", 0, 0, 0); + + PVPBuffer = new G4PVPlacement(0, positionAluStripFront, logicAluStrip, "G" + DetectorNumber + "AluStripFront", logicMM, false, 0); + PVPBuffer = new G4PVPlacement(0, positionAluStripBack, logicAluStrip, "G" + DetectorNumber + "AluStripBack", logicMM, false, 0); + + logicAluStrip->SetVisAttributes(G4VisAttributes::Invisible); + // Silicon detector itself - G4ThreeVector positionFirstStage = G4ThreeVector(0, 0, FirstStage_PosZ); + G4ThreeVector positionSilicon = G4ThreeVector(0, 0, Silicon_PosZ); + + G4Box* solidSilicon = new G4Box("solidSilicon", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*SiliconThickness); + G4LogicalVolume* logicSilicon = new G4LogicalVolume(solidSilicon, Silicon, "logicSilicon", 0, 0, 0); - G4Box* solidFirstStage = new G4Box("solidFirstStage", 0.5*FirstStageFace, 0.5*FirstStageFace, 0.5*FirstStageThickness); - G4LogicalVolume* logicFirstStage = new G4LogicalVolume(solidFirstStage, Silicon, "logicFirstStage", 0, 0, 0); + PVPBuffer = new G4PVPlacement(0, positionSilicon, logicSilicon, Name + "_Silicon", logicMM, false, 0); - PVPBuffer = new G4PVPlacement(0, - positionFirstStage, - logicFirstStage, - Name + "_FirstStage", - logicGPDSquare, - false, - 0); // Set First Stage sensible - logicFirstStage->SetSensitiveDetector(m_FirstStageScorer); + logicSilicon->SetSensitiveDetector(m_FirstStageScorer); - ///Visualisation of FirstStage Strip - G4VisAttributes* FirstStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.0, 0.9)); // blue - logicFirstStage->SetVisAttributes(FirstStageVisAtt); + ///Visualisation of Silicon Strip + G4VisAttributes* SiliconVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ; + logicSilicon->SetVisAttributes(SiliconVisAtt) ; } //////////////////////////////////////////////////////////////// - //////////////////// Second Stage Construction //////////////// + //////////////////// SiLi Construction //////////////////////// //////////////////////////////////////////////////////////////// if (wSecondStage) { - // Second stage silicon detector - G4ThreeVector positionSecondStage = G4ThreeVector(0, 0, SecondStage_PosZ); - - G4Box* solidSecondStage = new G4Box("solidSecondStage", 0.5*SecondStageFace, 0.5*SecondStageFace, 0.5*SecondStageThickness); - G4LogicalVolume* logicSecondStage = new G4LogicalVolume(solidSecondStage, Silicon, "logicSecondStage", 0, 0, 0); - - PVPBuffer = new G4PVPlacement(0, - positionSecondStage, - logicSecondStage, - Name + "_SecondStage", - logicGPDSquare, - false, - 0); + G4double SiLiSpace = 8 * mm; + + G4RotationMatrix* rotSiLi = Rotation(0., 0., 0.); + + G4ThreeVector positionSiLi = G4ThreeVector(-0.25 * SiliconFace - 0.5 * SiLiSpace, 0, 0); + G4ThreeVector positionSiLi2 = G4ThreeVector(0.25 * SiliconFace + 0.5 * SiLiSpace, 0, 0); + + G4Box* solidSiLi = new G4Box("SiLi", 0.5*SiLiFaceX, 0.5*SiLiFaceY, 0.5*SiLiThickness); + + G4LogicalVolume* logicSiLi = new G4LogicalVolume(solidSiLi, Aluminium, "SiLi" + DetectorNumber, 0, 0, 0); + + // First Si(Li) 2 time 4 detectore + PVPBuffer = new G4PVPlacement(G4Transform3D(*rotSiLi, positionSiLi) , + logicSiLi , + "SiLi" + DetectorNumber , + logicVacBox , + false , + 0); + + // Second Si(Li) 2 time 4 detectore + PVPBuffer = new G4PVPlacement(G4Transform3D(*rotSiLi, positionSiLi2) , + logicSiLi , + "SiLi" + DetectorNumber , + logicVacBox , + false , + 1); + + // SiLi are placed inside of the VacBox... + // Left/Right define when looking to detector from Si to CsI + G4double SiLi_HighY_Upper = 19.86 * mm; + G4double SiLi_HighY_Center = 25.39 * mm ; + G4double SiLi_WidthX_Left = 22.85 * mm; + G4double SiLi_WidthX_Right = 24.9 * mm ; + G4double SiLi_ShiftX = 0.775 * mm ; + + // SiLi : left side of SiLi detector + G4Box* solidSiLi_LT = new G4Box("SiLi_LT" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); + G4Box* solidSiLi_RT = new G4Box("SiLi_RT" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); + G4Box* solidSiLi_LC1 = new G4Box("SiLi_LC1" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); + G4Box* solidSiLi_RC1 = new G4Box("SiLi_RC1" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); + G4Box* solidSiLi_LB = new G4Box("SiLi_LB" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); + G4Box* solidSiLi_RB = new G4Box("SiLi_RB" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Upper , 0.5*SiLiThickness); + G4Box* solidSiLi_LC2 = new G4Box("SiLi_LC2" , 0.5*SiLi_WidthX_Left , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); + G4Box* solidSiLi_RC2 = new G4Box("SiLi_RC2" , 0.5*SiLi_WidthX_Right , 0.5*SiLi_HighY_Center , 0.5*SiLiThickness); + + G4LogicalVolume* logicSiLi_LT = new G4LogicalVolume(solidSiLi_LT , Silicon , "SiLi_LT" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_RT = new G4LogicalVolume(solidSiLi_RT , Silicon , "SiLi_RT" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_LC1 = new G4LogicalVolume(solidSiLi_LC1 , Silicon , "SiLi_LC1" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_RC1 = new G4LogicalVolume(solidSiLi_RC1 , Silicon , "SiLi_RC1" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_LB = new G4LogicalVolume(solidSiLi_LB , Silicon , "SiLi_LB" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_RB = new G4LogicalVolume(solidSiLi_RB , Silicon , "SiLi_RB" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_LC2 = new G4LogicalVolume(solidSiLi_LC2 , Silicon , "SiLi_LC2" , 0 , 0 , 0); + G4LogicalVolume* logicSiLi_RC2 = new G4LogicalVolume(solidSiLi_RC2 , Silicon , "SiLi_RC2" , 0 , 0 , 0); + + G4double interSiLi = 0.5 * mm; + + // Top + G4ThreeVector positionSiLi_LT = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi , + 0); + + G4ThreeVector positionSiLi_RT = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX , + 0.5 * SiLi_HighY_Upper + SiLi_HighY_Center + 1.5 * interSiLi , + 0); + + G4ThreeVector positionSiLi_LC1 = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi , + 0); + + G4ThreeVector positionSiLi_RC1 = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX , + 0.5 * SiLi_HighY_Center + 0.5 * interSiLi , + 0); + + // Bottom + G4ThreeVector positionSiLi_LB = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , + 0); + + G4ThreeVector positionSiLi_RB = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX , + -0.5 * SiLi_HighY_Upper - SiLi_HighY_Center - 1.5 * interSiLi , + 0); + + G4ThreeVector positionSiLi_LC2 = G4ThreeVector(-0.5 * SiLi_WidthX_Left - interSiLi - SiLi_ShiftX , + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi , + 0); + + G4ThreeVector positionSiLi_RC2 = G4ThreeVector(0.5 * SiLi_WidthX_Right - SiLi_ShiftX , + -0.5 * SiLi_HighY_Center - 0.5 * interSiLi , + 0); + + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LT , logicSiLi_LT , Name + "_SiLi_LT" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RT , logicSiLi_RT , Name + "_SiLi_RT" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LC1 , logicSiLi_LC1 , Name + "_SiLi_LC1" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RC1 , logicSiLi_RC1 , Name + "_SiLi_RC1" , logicSiLi , false , 0) ; + + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LB , logicSiLi_LB , Name + "_SiLi_LB" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RB , logicSiLi_RB , Name + "_SiLi_RB" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LC2 , logicSiLi_LC2 , Name + "_SiLi_LC2" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RC2 , logicSiLi_RC2 , Name + "_SiLi_RC2" , logicSiLi , false , 0) ; + + logicSiLi->SetVisAttributes(G4VisAttributes(G4Colour(1, 1., 1.))); // Set Second Stage sensible - logicSecondStage->SetSensitiveDetector(m_SecondStageScorer); - - ///Visualisation of SecondStage Strip - G4VisAttributes* SecondStageVisAtt = new G4VisAttributes(G4Colour(0.5, 0.5, 0.5)) ; - logicSecondStage->SetVisAttributes(SecondStageVisAtt) ; + logicSiLi_LT->SetSensitiveDetector(m_SecondStageScorer); + logicSiLi_RT->SetSensitiveDetector(m_SecondStageScorer); + logicSiLi_LC1->SetSensitiveDetector(m_SecondStageScorer); + logicSiLi_RC1->SetSensitiveDetector(m_SecondStageScorer); + + logicSiLi_LB->SetSensitiveDetector(m_SecondStageScorer); + logicSiLi_RB->SetSensitiveDetector(m_SecondStageScorer); + logicSiLi_LC2->SetSensitiveDetector(m_SecondStageScorer); + logicSiLi_RC2->SetSensitiveDetector(m_SecondStageScorer); + + // Mark blue a SiLi to see telescope orientation + logicSiLi_LT->SetVisAttributes(G4VisAttributes(G4Colour(0, 0., 1.))); + logicSiLi_RT->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); + logicSiLi_LC1->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); + logicSiLi_RC1->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); + + logicSiLi_LB->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); + logicSiLi_RB->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); + logicSiLi_LC2->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); + logicSiLi_RC2->SetVisAttributes(G4VisAttributes(G4Colour(0, 1., 0))); } //////////////////////////////////////////////////////////////// @@ -320,23 +453,19 @@ void GaspardTrackerSquare::VolumeMaker(G4int TelescopeNumber, // Third stage silicon detector G4ThreeVector positionThirdStage = G4ThreeVector(0, 0, ThirdStage_PosZ); - G4Box* solidThirdStage = new G4Box("solidThirdStage", 0.5*ThirdStageFace, 0.5*ThirdStageFace, 0.5*ThirdStageThickness); +// G4Box* solidThirdStage = new G4Box("solidThirdStage", 0.5*SiliconFace, 0.5*SiliconFace, 0.5*ThirdStageThickness); + G4Box* solidThirdStage = new G4Box("solidThirdStage", 0.5*FaceFront, 0.5*FaceFront, 0.5*ThirdStageThickness); G4LogicalVolume* logicThirdStage = new G4LogicalVolume(solidThirdStage, Silicon, "logicThirdStage", 0, 0, 0); - PVPBuffer = new G4PVPlacement(0, - positionThirdStage, - logicThirdStage, - Name + "_ThirdStage", - logicGPDSquare, - false, - 0); + PVPBuffer = new G4PVPlacement(0, positionThirdStage, logicThirdStage, Name + "_ThirdStage", logicMM, false, 0); + + ///Visualisation of Third Stage + G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.7, 0.7, 0.7)) ; + logicThirdStage->SetVisAttributes(ThirdStageVisAtt) ; +// logicThirdStage->SetVisAttributes(G4VisAttributes::Invisible); // Set Third Stage sensible logicThirdStage->SetSensitiveDetector(m_ThirdStageScorer); - - ///Visualisation of Third Stage - G4VisAttributes* ThirdStageVisAtt = new G4VisAttributes(G4Colour(0.0, 0.9, 0.0)); // green - logicThirdStage->SetVisAttributes(ThirdStageVisAtt); } } @@ -582,14 +711,18 @@ void GaspardTrackerSquare::ReadConfiguration(string Path) void GaspardTrackerSquare::ConstructDetector(G4LogicalVolume* world) { G4RotationMatrix* MMrot = NULL ; - G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ; +/* G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ; G4ThreeVector MMu = G4ThreeVector(0, 0, 0) ; G4ThreeVector MMv = G4ThreeVector(0, 0, 0) ; - G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ; + G4ThreeVector MMw = G4ThreeVector(0, 0, 0) ;*/ + MMpos = G4ThreeVector(0, 0, 0) ; + MMu = G4ThreeVector(0, 0, 0) ; + MMv = G4ThreeVector(0, 0, 0) ; + MMw = G4ThreeVector(0, 0, 0) ; G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ; - bool FirstStage = true; - bool SecondStage = true; - bool ThirdStage = true; + bool FirstStage = true ; + bool SecondStage = true ; + bool ThirdStage = true ; G4int NumberOfTelescope = m_DefinitionType.size() ; @@ -599,21 +732,31 @@ void GaspardTrackerSquare::ConstructDetector(G4LogicalVolume* world) // (u,v,w) unitary vector associated to telescope referencial // (u,v) // to silicon plan // w perpendicular to (u,v) plan and pointing ThirdStage - MMu = m_X128_Y1[i] - m_X1_Y1[i]; - MMu = MMu.unit(); + G4cout << "############ Gaspard " << i << " #############" << G4endl; + MMu = m_X128_Y1[i] - m_X1_Y1[i] ; + G4cout << "MMu: X = " << MMu(0) << " , Y = " << MMu(1) << " , Z = " << MMu(2) << G4endl; + MMu = MMu.unit() ; + G4cout << "Norm MMu: X = " << MMu(0) << " , Y = " << MMu(1) << " , Z = " << MMu(2) << G4endl; - MMv = m_X1_Y128[i] - m_X1_Y1[i]; - MMv = MMv.unit(); + MMv = m_X1_Y128[i] - m_X1_Y1[i] ; + G4cout << "MMv X = " << MMv(0) << " , Y = " << MMv(1) << " , Z = " << MMv(2) << G4endl; + MMv = MMv.unit() ; + G4cout << "Norm MMv X = " << MMv(0) << " , Y = " << MMv(1) << " , Z = " << MMv(2) << G4endl; - MMw = MMu.cross(MMv); - MMw = MMw.unit(); + G4ThreeVector MMscal = MMu.dot(MMv); + G4cout << "Norm MMu.MMv X = " << MMv(0) << " , Y = " << MMv(1) << " , Z = " << MMv(2) << G4endl; - MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4; + MMw = MMu.cross(MMv) ; +// if (MMw.z() > 0) MMw = MMv.cross(MMu) ; + MMw = MMw.unit() ; + + MMCenter = (m_X1_Y1[i] + m_X1_Y128[i] + m_X128_Y1[i] + m_X128_Y128[i]) / 4 ; // Passage Matrix from Lab Referential to Telescope Referential - MMrot = new G4RotationMatrix(MMu, MMv, MMw); + // MUST2 + MMrot = new G4RotationMatrix(MMu, MMv, MMw) ; // translation to place Telescope - MMpos = MMw * Length * 0.5 + MMCenter; + MMpos = MMw * Length * 0.5 + MMCenter ; } // By Angle @@ -626,13 +769,13 @@ void GaspardTrackerSquare::ConstructDetector(G4LogicalVolume* world) // w perpendicular to (u,v) plan and pointing ThirdStage // Phi is angle between X axis and projection in (X,Y) plan // Theta is angle between position vector and z axis - G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad); - G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad); - G4double wZ = m_R[i] * cos(Theta / rad); - MMw = G4ThreeVector(wX, wY, wZ); + G4double wX = m_R[i] * sin(Theta / rad) * cos(Phi / rad) ; + G4double wY = m_R[i] * sin(Theta / rad) * sin(Phi / rad) ; + G4double wZ = m_R[i] * cos(Theta / rad) ; + MMw = G4ThreeVector(wX, wY, wZ) ; // vector corresponding to the center of the module - MMCenter = MMw; + CT = MMw; // vector parallel to one axis of silicon plane G4double ii = cos(Theta / rad) * cos(Phi / rad); @@ -654,7 +797,7 @@ void GaspardTrackerSquare::ConstructDetector(G4LogicalVolume* world) MMrot->rotate(m_beta_v[i], MMv); MMrot->rotate(m_beta_w[i], MMw); // translation to place Telescope - MMpos = MMw * Length * 0.5 + MMCenter; + MMpos = MMw * Length * 0.5 + CT ; } FirstStage = m_wFirstStage[i] ; diff --git a/NPSimulation/src/MUST2Array.cc b/NPSimulation/src/MUST2Array.cc index b8c5f43281a739a3ab1ec0d4bc16fcc661fbe14a..e644c5950b1859f5bd080eb362d7aff080be8042 100644 --- a/NPSimulation/src/MUST2Array.cc +++ b/NPSimulation/src/MUST2Array.cc @@ -164,11 +164,11 @@ void MUST2Array::VolumeMaker(G4int TelescopeNumber , PVPBuffer = new G4PVPlacement( G4Transform3D(*MMrot, MMpos) , - logicMM , - Name , - world , - false , - 0); + logicMM , + Name , + world , + false , + 0 ); diff --git a/NPSimulation/src/Must2Scorers.cc b/NPSimulation/src/Must2Scorers.cc index 1d9d6bc813ab1b3a4b362744fefb6c2b6f525d20..0493519a3b5f7a65833b5b64b81c89a82857ca7c 100644 --- a/NPSimulation/src/Must2Scorers.cc +++ b/NPSimulation/src/Must2Scorers.cc @@ -9,7 +9,7 @@ * Original Author: Adrien MATTA contact address: matta@ipno.in2p3.fr * * * * Creation Date : January 2009 * - * Last update : * + * Last update : october 2010 * *---------------------------------------------------------------------------* * Decription: * * File old the scorer specific to the MUST2 Detector * @@ -52,7 +52,7 @@ G4bool PSStripNumberX::ProcessHits(G4Step* aStep, G4TouchableHistory*) G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); - G4double temp = (POS(0) + m_StripPlaneSize / 2.) / m_StripPitch ; + G4double temp = (POS(1) + m_StripPlaneSize / 2.) / m_StripPitch ; G4int X = int(temp) + 1 ; //Rare case where particle is close to edge of silicon plan if (X == m_NumberOfStrip+1) X = m_NumberOfStrip; @@ -118,7 +118,7 @@ G4bool PSStripNumberY::ProcessHits(G4Step* aStep, G4TouchableHistory*) G4ThreeVector POS = aStep->GetPreStepPoint()->GetPosition(); POS = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().TransformPoint(POS); - G4double temp = (POS(1) + m_StripPlaneSize / 2.) / m_StripPitch ; + G4double temp = (POS(0) + m_StripPlaneSize / 2.) / m_StripPitch ; G4int temp2 = temp ; G4int Y = temp2 + 1 ; //Rare case where particle is close to edge of silicon plan diff --git a/NPSimulation/src/ParisCluster.cc b/NPSimulation/src/ParisCluster.cc index 9f35358a8ae2650c5e9baf3c79f25051c4dacffb..3ff1f1c10c6c8c3816e8b96f0033b6d584c9d305 100644 --- a/NPSimulation/src/ParisCluster.cc +++ b/NPSimulation/src/ParisCluster.cc @@ -564,8 +564,6 @@ void ParisCluster::ConstructDetector(G4LogicalVolume* world) MMv = m_X1_Y128[i] - m_X1_Y1[i]; MMv = MMv.unit(); - G4ThreeVector MMscal = MMu.dot(MMv); - MMw = MMu.cross(MMv); // if (MMw.z() > 0) MMw = MMv.cross(MMu) ; MMw = MMw.unit(); diff --git a/NPSimulation/src/ParisPhoswich.cc b/NPSimulation/src/ParisPhoswich.cc index b32609a0ad86d1fcef427efab9edd55fcaecb48f..a0a956a8114f1989945a5223c3f295ad594f4c87 100644 --- a/NPSimulation/src/ParisPhoswich.cc +++ b/NPSimulation/src/ParisPhoswich.cc @@ -475,8 +475,6 @@ void ParisPhoswich::ConstructDetector(G4LogicalVolume* world) MMv = m_X1_Y128[i] - m_X1_Y1[i]; MMv = MMv.unit(); - G4ThreeVector MMscal = MMu.dot(MMv); - MMw = MMu.cross(MMv); // if (MMw.z() > 0) MMw = MMv.cross(MMu) ; MMw = MMw.unit(); diff --git a/NPSimulation/src/PrimaryGeneratorAction.cc b/NPSimulation/src/PrimaryGeneratorAction.cc index c6a852cdb4fe334a5a68e3425e3bd3908571b835..4e61398936dfc381717afac5f09f29980325c3bd 100644 --- a/NPSimulation/src/PrimaryGeneratorAction.cc +++ b/NPSimulation/src/PrimaryGeneratorAction.cc @@ -79,10 +79,10 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path) { // added by Nicolas [07/05/09] string GlobalPath = getenv("NPTOOL"); - Path = GlobalPath + "/Inputs/EventGenerator/" + Path; + string StandardPath = GlobalPath + "/Inputs/EventGenerator/" + Path; string LineBuffer; ifstream EventGeneratorFile; - EventGeneratorFile.open(Path.c_str()); + EventGeneratorFile.open(StandardPath.c_str()); bool check_TransfertToResonance = false ; bool check_PhaseSpace = false ; @@ -90,11 +90,22 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path) bool check_Transfert = false ; bool check_Beam = false ; - if (EventGeneratorFile.is_open()) cout << " Event Generator file " << Path << " loading " << endl ; - else { - cout << " Error, Event Generator file " << Path << " found" << endl ; - return; - } + if (EventGeneratorFile.is_open()) + { + cout << " Event Generator file " << Path << " loading " << endl ; + Path = StandardPath; + } + + else + { + EventGeneratorFile.open( Path.c_str() ); + if(EventGeneratorFile.is_open()) { + cout << " Event Generator file " << Path << " loading " << endl ; + } + + else { cout << " Error, Event Generator file " << Path << " found" << endl ; return;} + } + while (!EventGeneratorFile.eof()) { //Pick-up next line diff --git a/NPSimulation/src/ShieldClParis.cc b/NPSimulation/src/ShieldClParis.cc index 2b77daa1872f600c4d441b9d0551885fcc843bc7..f3fe6a464dddf0450840d05d4fa8e398faea9632 100644 --- a/NPSimulation/src/ShieldClParis.cc +++ b/NPSimulation/src/ShieldClParis.cc @@ -463,8 +463,6 @@ void ShieldClParis::ConstructDetector(G4LogicalVolume* world) //MMv = -0.5 * (m_X1_Y1[i] + m_X128_Y128[i] - m_X1_Y128[i] - m_X128_Y1[i]); MMv = MMv.unit(); - G4ThreeVector MMscal = MMu.dot(MMv); - MMw = MMu.cross(MMv); // if (MMw.z() > 0) MMw = MMv.cross(MMu) ; MMw = MMw.unit(); diff --git a/NPSimulation/src/ShieldPhParis.cc b/NPSimulation/src/ShieldPhParis.cc index 0f6b140dcf54dbd26f18d4ed95396aea0ea61cc1..8d8edc55db700584e415c66ac2189649a034cdba 100644 --- a/NPSimulation/src/ShieldPhParis.cc +++ b/NPSimulation/src/ShieldPhParis.cc @@ -473,8 +473,6 @@ void ShieldPhParis::ConstructDetector(G4LogicalVolume* world) //MMv = -0.5 * (m_X1_Y1[i] + m_X128_Y128[i] - m_X1_Y128[i] - m_X128_Y1[i]); MMv = MMv.unit(); - G4ThreeVector MMscal = MMu.dot(MMv); - MMw = MMu.cross(MMv); // if (MMw.z() > 0) MMw = MMv.cross(MMu) ; MMw = MMw.unit(); diff --git a/NPSimulation/src/Target.cc b/NPSimulation/src/Target.cc index 6e8d1d6b4b46bb7ab898c5ddf8309e0c8ab16f80..3b3b32aa591a7e37608bbe219e51bbdce7c15dd2 100644 --- a/NPSimulation/src/Target.cc +++ b/NPSimulation/src/Target.cc @@ -216,6 +216,16 @@ G4Material* Target::GetMaterialFromLibrary(G4String MaterialName, G4double Tempe return myMaterial; } + else if (MaterialName == "Si") { + G4Material* myMaterial = new G4Material("Si", 14., 28.0855 * g / mole, 2.321 * g / cm3); + return myMaterial; + } + + else if (MaterialName == "Al") { + G4Material* myMaterial = new G4Material("Aluminium", 13., 26.98 * g / mole, 2.702 * g / cm3); + return myMaterial; + } + else { G4cout << "No Matching Material in the Target Library Default is Vacuum" << G4endl; G4Element* N = new G4Element("Nitrogen", "N", 7., 14.01*g / mole);