From a87d2e461324e83af412b5b8b8cb1581eec1c936 Mon Sep 17 00:00:00 2001 From: Morfouace <morfouac@ipno.in2p3.fr> Date: Sun, 6 Dec 2015 18:05:54 -0500 Subject: [PATCH] Modifying CsI class to test optical photon --- Inputs/DetectorConfiguration/csi.detector | 109 ++++++------- Inputs/DetectorConfiguration/e552.detector | 2 +- Inputs/EventGenerator/3He.source | 6 +- Inputs/EventGenerator/alpha.source | 6 +- Inputs/EventGenerator/deuton.source | 4 +- Inputs/EventGenerator/proton.source | 4 +- Inputs/EventGenerator/triton.source | 4 +- NPLib/CsI/TCsIData.h | 34 ++-- NPLib/Physics/nubtab03.asc | 18 +-- NPLib/scripts/RootLogon.sh | 2 +- NPSimulation/Core/PhysicsList.cc | 36 +++-- NPSimulation/CsI/CsI.cc | 174 +++++++++++++++------ NPSimulation/CsI/CsI.hh | 9 +- NPSimulation/PhysicsListOption.txt | 4 +- NPSimulation/SSSD/SSSD.hh | 2 +- 15 files changed, 255 insertions(+), 159 deletions(-) diff --git a/Inputs/DetectorConfiguration/csi.detector b/Inputs/DetectorConfiguration/csi.detector index 6abde491b..13e38411f 100644 --- a/Inputs/DetectorConfiguration/csi.detector +++ b/Inputs/DetectorConfiguration/csi.detector @@ -5,8 +5,8 @@ GeneralTarget % Radius in mm % Temperature in K, Pressure in bar Target - THICKNESS= 0.1 - ANGLE= 0 + THICKNESS= 0.001 + ANGLE= 0 RADIUS= 2 MATERIAL= CH2 X= 0 @@ -19,67 +19,68 @@ CsI THETA= 0 PHI= 0 R= 0 - X= 0 - Y= 0 - Z= 450 + X= 0 + Y= 0 + Z= 58 Thickness= 100 - Shape= Trapezoidal - FaceFront= 22 - FaceBack= 27 - Radius= 0 - Scintillator= CsI + Shape= Cylinder + FaceFront= 0 + FaceBack= 0 + Radius= 17 + Scintillator= CsI_Scintillator LeadThickness= 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ScintillatorPlastic %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Plastic - THETA= 0 - PHI= 0 - R= 0 - X= 0 - Y= 0 - Z= 389.25 - Thickness= 1.5 - Shape= Square - Height= 22 - Width= 22 - Radius= 0 - Scintillator= Si - LeadThickness= 0 +%Plastic +% THETA= 0 +% PHI= 0 +% R= 0 +% X= 0 +% Y= 0 +% Z= 389.25 +% Thickness= 1.5 +% Shape= Square +% Height= 22 +% Width= 22 +% Radius= 0 +% Scintillator= Si +% LeadThickness= 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Plastic - THETA= 0 - PHI= 0 - R= 0 - X= 0 - Y= 0 - Z= 505 - Thickness= 10 - Shape= Square - Height= 27 - Width= 27 - Radius= 0 - Scintillator= Epoxy - LeadThickness= 0 +%Plastic +% THETA= 0 +% PHI= 0 +% R= 0 +% X= 0 +% Y= 0 +% +% Z= 505 +% Thickness= 10 +% Shape= Square +% Height= 27 +% Width= 27 +% Radius= 0 +% Scintillator= Epoxy +% LeadThickness= 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -Plastic - THETA= 0 - PHI= 0 - R= 0 - X= 0 - Y= 0 - Z= 510.15 - Thickness= 0.3 - Shape= Square - Height= 27 - Width= 27 - Radius= 0 - Scintillator= Si - LeadThickness= 0 - - +%Plastic +% THETA= 0 +% PHI= 0 +% R= 0 +% X= 0 +% Y= 0 +% Z= 510.15 +% Thickness= 0.3 +% Shape= Square +% Height= 27 +% Width= 27 +% Radius= 0 +% Scintillator= Si +% LeadThickness= 0 +% +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Avaible material: % BC400 diff --git a/Inputs/DetectorConfiguration/e552.detector b/Inputs/DetectorConfiguration/e552.detector index 107753b28..40d822ea5 100755 --- a/Inputs/DetectorConfiguration/e552.detector +++ b/Inputs/DetectorConfiguration/e552.detector @@ -8,7 +8,7 @@ GeneralTarget % Temperature in K, Pressure in bar Target THICKNESS= 3.2 - ANGLE= 28 + ANGLE= 20 RADIUS= 24 MATERIAL= CD2 X= 0 diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source index 22957a57d..66b22a966 100644 --- a/Inputs/EventGenerator/3He.source +++ b/Inputs/EventGenerator/3He.source @@ -5,12 +5,12 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic EnergyLow= 0 - EnergyHigh= 300 + EnergyHigh= 26 HalfOpenAngleMin= 0 - HalfOpenAngleMax= 90 + HalfOpenAngleMax= 45 x0= 0 y0= 0 z0= 0 - particle= He3 + particle= 3He %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Supported particle type: proton, neutron, deuton, triton, He3 , alpha diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source index 8d3a34d18..1fed2da97 100644 --- a/Inputs/EventGenerator/alpha.source +++ b/Inputs/EventGenerator/alpha.source @@ -4,10 +4,10 @@ % Energy are given in MeV , Position in mm % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic - EnergyLow= 5.84 - EnergyHigh= 5.84 + EnergyLow= 5.2 + EnergyHigh= 5.2 HalfOpenAngleMin= 0 - HalfOpenAngleMax= 180 + HalfOpenAngleMax= 45 x0= 0 y0= 0 z0= 0 diff --git a/Inputs/EventGenerator/deuton.source b/Inputs/EventGenerator/deuton.source index 5b70824ee..fd3680952 100644 --- a/Inputs/EventGenerator/deuton.source +++ b/Inputs/EventGenerator/deuton.source @@ -5,9 +5,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic EnergyLow= 0 - EnergyHigh= 100 + EnergyHigh= 9 HalfOpenAngleMin= 0 - HalfOpenAngleMax= 90 + HalfOpenAngleMax= 45 x0= 0 y0= 0 z0= 0 diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source index 052193a50..cf3a61752 100644 --- a/Inputs/EventGenerator/proton.source +++ b/Inputs/EventGenerator/proton.source @@ -5,9 +5,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic EnergyLow= 0 - EnergyHigh= 190 + EnergyHigh= 17 HalfOpenAngleMin= 0 - HalfOpenAngleMax= 100 + HalfOpenAngleMax= 45 x0= 0 y0= 0 z0= 0 diff --git a/Inputs/EventGenerator/triton.source b/Inputs/EventGenerator/triton.source index 191488a05..5a9aa2a5a 100644 --- a/Inputs/EventGenerator/triton.source +++ b/Inputs/EventGenerator/triton.source @@ -5,9 +5,9 @@ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Isotropic EnergyLow= 0 - EnergyHigh= 100 + EnergyHigh= 12 HalfOpenAngleMin= 0 - HalfOpenAngleMax= 90 + HalfOpenAngleMax= 45 x0= 0 y0= 0 z0= 0 diff --git a/NPLib/CsI/TCsIData.h b/NPLib/CsI/TCsIData.h index 45ead2b91..caf462896 100644 --- a/NPLib/CsI/TCsIData.h +++ b/NPLib/CsI/TCsIData.h @@ -28,12 +28,16 @@ using namespace std ; class TCsIData : public TObject { private: - // Energy - vector<short> fCsI_E_Number; - vector<double> fCsI_E_Energy; - // Time - vector<short> fCsI_T_Number; - vector<double> fCsI_T_Time; + // Energy + vector<short> fCsI_E_Number; + vector<double> fCsI_E_Energy; + vector<double> fPhotoDiode_E_Energy; + vector<short> fPhotoDiode_E_Number; + // Time + vector<short> fCsI_T_Number; + vector<double> fCsI_T_Time; + vector<double> fPhotoDiode_T_Time; + vector<short> fPhotoDiode_T_Number; public: TCsIData(); @@ -53,13 +57,17 @@ class TCsIData : public TObject { int GetTNumber(int i) {return fCsI_T_Number[i];} double GetTTime(int i) {return fCsI_T_Time[i];} - ///////////////////// SETTERS //////////////////////// - // Energy - void SetENumber(int N) {fCsI_E_Number.push_back(N);} - void SetEEnergy(double E) {fCsI_E_Energy.push_back(E);} - // time - void SetTNumber(int N) {fCsI_T_Number.push_back(N);} - void SetTTime(double T) {fCsI_T_Time.push_back(T);} + ///////////////////// SETTERS //////////////////////// + // Energy + void SetENumber(int N) {fCsI_E_Number.push_back(N);} + void SetCsIEEnergy(double E) {fCsI_E_Energy.push_back(E);} + void SetPhotoDiodeEnergy(double E) {fPhotoDiode_E_Energy.push_back(E);} + void SetPhotoDiodeEDetectorNbr(int N) {fPhotoDiode_E_Number.push_back(N);} + // time + void SetTNumber(int N) {fCsI_T_Number.push_back(N);} + void SetTTime(double T) {fCsI_T_Time.push_back(T);} + void SetPhotoDiodeTime(double T) {fPhotoDiode_T_Time.push_back(T);} + void SetPhotoDiodeTDetectorNbr(int N) {fPhotoDiode_T_Number.push_back(N);} ClassDef(TCsIData,1) // CsIData structure }; diff --git a/NPLib/Physics/nubtab03.asc b/NPLib/Physics/nubtab03.asc index ff36e90df..b4860b280 100644 --- a/NPLib/Physics/nubtab03.asc +++ b/NPLib/Physics/nubtab03.asc @@ -1338,7 +1338,7 @@ 100 0480 100Cd -74250 100 49.1 s 0.5 0+ 97 B+=100 100 0481 100Cdm -71700 100 2548.6 0.5 60 ns 3 (8)+ 97 IT=100 100 0490 100In -64170 250 5.9 s 0.2 (6,7)+ 97 02Pl03tj B+=100;B+p>3.9 -100 0500 100Sn -56780 710 1.1 s 0.4 0+ 97 B+=100;B+p<17 +100 0500 100Sn -57280 300 1.1 s 0.4 0+ 97 B+=100;B+p<17 101 0370 101Rb -43600 170 32 ms 4 3/2+# 98 B-=100;B-n=28 4 101 0380 101Sr -55410 120 118 ms 3 (5/2-) 98 B-=100;B-n=2.37 14 101 0390 101Y -64910 100 426 ms 20 (5/2+) 98 96Me09t B-=100;B-n=1.94 18 @@ -1357,7 +1357,7 @@ 101 0480 101Cd -75750 150 1.36 m 0.05 (5/2+) 98 B+=100 101 0490 101In -68610# 300# 15.1 s 1.1 9/2+# 98 B+=100;B+p=? 101 0491 101Inm -68060# 320# 550# 100# 10# s 1/2-# B+=95#;IT=5# -101 0500 101Sn -59560# 300# 3 s 1 5/2+# 98 B+=100;B+p=? +101 0500 101Sn -60310 300 3 s 1 5/2+# 98 B+=100;B+p=? 102 0370 102Rb -38310# 500# 37 ms 5 98 B-=100;B-n=18 8 102 0380 102Sr -53080 110 69 ms 6 0+ 98 93Ru01d B-=100;B-n=5.5 15 102 0390 102Y -61890 90 *& 300 ms 10 low 98 B-=100;B-n=4.9 12 @@ -1395,7 +1395,7 @@ 103 0480 103Cd -80649 15 7.3 m 0.1 5/2+ 01 B+=100 103 0490 103In -74599 25 60 s 1 9/2+# 01 97Sz04t B+=100 103 0491 103Inm -73967 25 631.7 0.1 34 s 2 1/2-# 01 97Sz04etd B+=67;IT=33 -103 0500 103Sn -66970# 300# 7 s 3 5/2+# 01 B+=100;B+p=? +103 0500 103Sn -66970 300# 7 s 3 5/2+# 01 B+=100;B+p=? 103 0510 103Sb -56180# 300# 100# ms >1.5us 5/2+# 01 95Ry03i B+ ? 104 0380 104Sr -44400# 700# 30# ms >300ns 0+ 00 97Be70i B- ? 104 0390 104Y -54910# 400# 180 ms 60 00 99Wa09d B-=100;B-n=? @@ -1414,7 +1414,7 @@ 104 0480 104Cd -83975 9 57.7 m 1.0 0+ 00 B+=100 104 0490 104In -76110 80 1.80 m 0.03 5,6(+) 00 B+=100 104 0491 104Inm -76020 80 93.48 0.10 15.7 s 0.5 (3+) 00 IT=80;B+=20 -104 0500 104Sn -71590 100 20.8 s 0.5 0+ 00 B+=100 +104 0500 104Sn -71627 100 20.8 s 0.5 0+ 00 B+=100 104 0510 104Sb -59180# 360# 470 ms 130 00 95Fa.Ad B+=?;B+p<7;p<7;A ? 105 0380 105Sr -38580# 700# 20# ms >300ns 97 97Be70i B- ? 105 0390 105Y -51350# 500# 60# ms >300ns 5/2+# 97 94Be24i B- ? @@ -1431,7 +1431,7 @@ 105 0480 105Cd -84330 12 55.5 m 0.4 5/2+ 93 B+=100 105 0490 105In -79481 17 5.07 m 0.07 9/2+ 93 87Eb02j B+=100 105 0491 105Inm -78807 17 674.1 0.3 48 s 6 (1/2)- 93 IT=?;B+=25# -105 0500 105Sn -73260 80 34 s 1 (5/2+) 93 95Pf01t B+=100;B+p=? +105 0500 105Sn -73338 80 34 s 1 (5/2+) 93 95Pf01t B+=100;B+p=? 105 0510 105Sb -63820 100 1.12 s 0.16 (5/2+) 02 B+ ?;p~1;B+p ? 105 0520 105Te -52500# 500# 1# us 5/2+# A ?;B+ ? 106 0390 106Y -46770# 700# 50# ms >300ns 97 97Be70i B- ? @@ -1448,7 +1448,7 @@ 106 0480 106Cd -87132 6 stbl >410Ey 0+ 94 02Tr04t IS=1.25 6;2B+ ? 106 0490 106In -80606 12 6.2 m 0.1 7+ 94 B+=100 106 0491 106Inm -80577 12 28.6 0.3 5.2 m 0.1 (3+) 94 B+=100 -106 0500 106Sn -77430 50 1.92 m 0.08 0+ 94 B+=100 +106 0500 106Sn -77354 50 1.92 m 0.08 0+ 94 B+=100 106 0510 106Sb -66330# 310# 600 ms 200 (4+) 97 94Se01j B+=100 106 0511 106Sbm -65330# 590# 1000# 500# 220 ns 20 98Li50t IT=100 106 0520 106Te -58210 130 70 us 20 0+ 94 94Pa11t A=100 @@ -1469,7 +1469,7 @@ 107 0480 107Cd -86985 6 6.50 h 0.02 5/2+ 00 B+=100 107 0490 107In -83560 11 32.4 m 0.3 9/2+ 00 B+=100 107 0491 107Inm -82882 11 678.5 0.3 50.4 s 0.6 1/2- 00 IT=100 -107 0500 107Sn -78580 80 2.90 m 0.05 (5/2+) 00 B+=100 +107 0500 107Sn -78512 80 2.90 m 0.05 (5/2+) 00 B+=100 107 0510 107Sb -70650# 300# 4.6 s 0.8 5/2+# 00 B+=100 107 0520 107Te -60540# 300# 3.1 ms 0.1 5/2+# 00 A=70 30;B+=30 30 108 0390 108Y -37740# 800# 20# ms >300ns 00 95Cz.Ai B- ?;B-n ? @@ -1486,7 +1486,7 @@ 108 0480 108Cd -89252 6 stbl >410Py 0+ 02 95Ge14t IS=0.89 3;2B+ ? 108 0490 108In -84116 10 58.0 m 1.2 7+ 00 B+=100 108 0491 108Inm -84086 10 29.75 0.05 39.6 m 0.7 2+ 00 B+=100 -108 0500 108Sn -82041 20 10.30 m 0.08 0+ 00 B+=100 +108 0500 108Sn -82070 20 10.30 m 0.08 0+ 00 B+=100 108 0510 108Sb -72510# 210# 7.4 s 0.3 (4+) 00 B+=100;B+p ? 108 0520 108Te -65720 100 2.1 s 0.1 0+ 00 85Ti02d B+=51 4;A=49 4;B+p=2.4 10;B+A<0.065 108 0530 108I -52650# 360# 36 ms 6 1+# 00 94Pa12d A=?;B+=9#;p<1 @@ -1506,7 +1506,7 @@ 109 0490 109In -86489 6 4.2 h 0.1 9/2+ 99 B+=100 109 0491 109Inm -85839 6 650.1 0.3 1.34 m 0.07 1/2- 99 IT=100 109 0492 109Inn -84387 6 2101.8 0.2 209 ms 6 (19/2+) 99 IT=100 -109 0500 109Sn -82639 10 18.0 m 0.2 5/2(+) 99 B+=100 +109 0500 109Sn -82631 10 18.0 m 0.2 5/2(+) 99 B+=100 109 0510 109Sb -76259 19 17.0 s 0.7 5/2+# 99 B+=100 109 0520 109Te -67610 60 4.6 s 0.3 (5/2+) 99 B+=?;A=3.9 13;B+p=9.4 31;B+A<0.005 109 0530 109I -57610 100 103 us 5 (5/2+) 02 87Gi02j p=100 diff --git a/NPLib/scripts/RootLogon.sh b/NPLib/scripts/RootLogon.sh index cf24507af..7604d53d9 100755 --- a/NPLib/scripts/RootLogon.sh +++ b/NPLib/scripts/RootLogon.sh @@ -42,7 +42,7 @@ then else echo 'File .rootlogon.C was created' - cp scripts/rootlogon_basic.C ~/.rootlogon.C + cp rootlogon_basic.C ~/.rootlogon.C fi diff --git a/NPSimulation/Core/PhysicsList.cc b/NPSimulation/Core/PhysicsList.cc index 2286e7828..37d45597d 100644 --- a/NPSimulation/Core/PhysicsList.cc +++ b/NPSimulation/Core/PhysicsList.cc @@ -84,14 +84,16 @@ PhysicsList::PhysicsList() : G4VModularPhysicsList(){ m_PhysList["HadronPhysicsQGSP_BIC_HP"] = new G4HadronPhysicsQGSP_BIC_HP(); // Optical Photon for scintillator simulation - if(m_OpticalPhysics) - m_PhysList["OpticalPhysics"] = new G4OpticalPhysics(); + if(m_OpticalPhysics) + m_PhysList["OpticalPhysics"] = new G4OpticalPhysics(); + // Decay physics // Add Radioactive decay - // Gamma decay of know states + // Gamma decay of known states if(m_Decay){ + std::cout << "Decay is activated: m_Decay=" << m_Decay << std::endl; decay_List = new G4DecayPhysics(); radioactiveDecay_List = new G4RadioactiveDecayPhysics() ; m_PhysList["decay_list"]= decay_List; @@ -198,21 +200,21 @@ void PhysicsList::ConstructParticle(){ ///////////////////////////////////////////////////////////////////////////// void PhysicsList::ConstructProcess(){ - // Transportation - AddTransportation(); - - // Electromagnetic physics - emPhysicsList -> ConstructProcess(); - em_config.AddModels(); - - // Hadronic physics - std::map<std::string,G4VPhysicsConstructor*>::iterator it; - for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){ - it->second -> ConstructProcess(); - } + // Transportation + AddTransportation(); + + // Electromagnetic physics + emPhysicsList -> ConstructProcess(); + em_config.AddModels(); + + // Hadronic physics + std::map<std::string,G4VPhysicsConstructor*>::iterator it; + for(it = m_PhysList.begin(); it!= m_PhysList.end(); it++){ + it->second -> ConstructProcess(); + } BiasCrossSectionByFactor(m_IonBinaryCascadePhysics); - - return; + SetCuts(); + return; } ///////////////////////////////////////////////////////////////////////////// void PhysicsList::AddStepMax(){ diff --git a/NPSimulation/CsI/CsI.cc b/NPSimulation/CsI/CsI.cc index 7603419ab..7017ec07c 100644 --- a/NPSimulation/CsI/CsI.cc +++ b/NPSimulation/CsI/CsI.cc @@ -46,8 +46,10 @@ #include "ObsoleteGeneralScorers.hh" #include "RootOutput.h" #include "MaterialManager.hh" +#include "SiliconScorers.hh" +#include "CalorimeterScorers.hh" #include "NPSDetectorFactory.hh" -using namespace OBSOLETEGENERALSCORERS ; +//using namespace OBSOLETEGENERALSCORERS ; // CLHEP header #include "CLHEP/Random/RandGauss.h" @@ -57,10 +59,9 @@ using namespace CLHEP; //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... namespace CSI{ - // Energy and time Resolution - const G4double ResoTime = 4.2 ;// = 10ns of Resolution // Unit is MeV/2.35 - //const G4double ResoEnergy = 0.180/2.35 ;// Resolution in keV (sigma) - const G4double ResoEnergy = 0.50 ;// Resolution in % + // Energy and time Resolution + //const G4double ResoTime = 4.2 ;// = 10ns of Resolution // Unit is MeV/2.35 + const G4double ResoCsI = 2.4;//% } using namespace CSI ; @@ -396,32 +397,52 @@ void CsI::VolumeMaker(G4ThreeVector Det_pos, int DetNumber, G4LogicalVolume* wor // Cylindrical Case if(m_CsIRadius[i]!=-1){ - if(m_CsIThickness[i]>0 && m_CsIRadius[i]>0){ - G4Tubs* solidCsI = new G4Tubs( Name , - 0 , - m_CsIRadius[i] , - m_CsIThickness[i]/2 , - 0*deg , - 360*deg); - - G4LogicalVolume* logicCsI = new G4LogicalVolume(solidCsI, CsIMaterial, Name+ "_Scintillator", 0, 0, 0); - logicCsI->SetSensitiveDetector(m_CsIScorer); - - G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1.0, 0.5, 0.0)) ; - logicCsI->SetVisAttributes(CsIVisAtt) ; - - - - new G4PVPlacement(0 , - Det_pos , - logicCsI , - Name + "_Scintillator" , - world , - false , - 0 ); + if(m_CsIThickness[i]>0 && m_CsIRadius[i]>0){ + + // CsI crystal + G4Tubs* solidCsI = new G4Tubs( Name , + 0 , + m_CsIRadius[i] , + m_CsIThickness[i]/2 , + 0*deg , + 360*deg); + + G4LogicalVolume* logicCsI = new G4LogicalVolume(solidCsI, CsIMaterial, Name+ "_Scintillator", 0, 0, 0); + logicCsI->SetSensitiveDetector(m_CsIScorer); + + G4VisAttributes* CsIVisAtt = new G4VisAttributes(G4Colour(1.0, 0.5, 0.0)) ; + logicCsI->SetVisAttributes(CsIVisAtt) ; + + new G4PVPlacement(0 , + Det_pos , + logicCsI , + Name + "_Scintillator" , + world , + false , + 0 ); + + // Photodiode + G4String NamePD = Name+"PhotoDiode"; + + G4Material* PDMaterial = MaterialManager::getInstance()->GetMaterialFromLibrary("Si"); + + G4Box* solidPhotoDiode = new G4Box(NamePD,0.5*PhotoDiodeFace,0.5*PhotoDiodeFace,0.5*PhotoDiodeThickness); + + G4LogicalVolume* logicPD = new G4LogicalVolume(solidPhotoDiode, PDMaterial, NamePD,0,0,0); + logicPD->SetSensitiveDetector(m_PDScorer); + + G4VisAttributes* PDVisAtt = new G4VisAttributes(G4Colour(0.1, 0.2, 0.3)) ; + logicPD->SetVisAttributes(PDVisAtt); + + new G4PVPlacement(0 , + Det_pos+(m_CsIThickness[i]/2+PhotoDiodeThickness/2)*Det_pos.unit() , + logicPD , + NamePD , + world , + false , + 0 ); } - if(m_LeadThickness[i]>0&& m_CsIRadius[i]>0){ G4Tubs* solidLead = new G4Tubs(Name+"_Lead", 0, @@ -499,9 +520,48 @@ void CsI::InitializeRootOutput(){ // Read sensitive part and fill the Root tree. // Called at in the EventAction::EndOfEventAvtion void CsI::ReadSensitive(const G4Event* event){ - G4String DetectorNumber; - m_Event->Clear(); + //G4String DetectorNumber; + m_Event->Clear(); + + // CsI // + G4THitsMap<G4double*>* CsIHitMap; + std::map<G4int, G4double**>::iterator CsI_itr; + G4int CsICollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("CsIScorer/CsI"); + CsIHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(CsICollectionID)); + + // Loop on the CsI map + for (CsI_itr = CsIHitMap->GetMap()->begin() ; CsI_itr != CsIHitMap->GetMap()->end() ; CsI_itr++){ + G4double* Info = *(CsI_itr->second); + double E_CsI = RandGauss::shoot(Info[0],Info[0]*ResoCsI/100); + //cout << "Energy CsI " << endl; + //cout << E_CsI << endl; + m_Event->SetCsIEEnergy(E_CsI); + m_Event->SetENumber(Info[2]); + } + // Clear Map for next event + CsIHitMap->clear(); + + // PhotoDiode // + G4THitsMap<G4double*>* PhotoDiodeHitMap; + std::map<G4int, G4double**>::iterator PhotoDiode_itr; + + G4int PhotoDiodeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("PDScorer/PhotoDiode"); + PhotoDiodeHitMap = (G4THitsMap<G4double*>*)(event->GetHCofThisEvent()->GetHC(PhotoDiodeCollectionID)); + + // Loop on the PhotoDiode map + for (PhotoDiode_itr = PhotoDiodeHitMap->GetMap()->begin() ; PhotoDiode_itr != PhotoDiodeHitMap->GetMap()->end() ; PhotoDiode_itr++){ + G4double* Info = *(PhotoDiode_itr->second); + double E_PhotoDiode = RandGauss::shoot(Info[0],Info[0]*ResoCsI/100); + + m_Event->SetPhotoDiodeEnergy(E_PhotoDiode); + m_Event->SetPhotoDiodeEDetectorNbr(Info[7]); + + m_Event->SetPhotoDiodeTime(Info[1]); + m_Event->SetPhotoDiodeTDetectorNbr(Info[7]); + } + PhotoDiodeHitMap->clear(); + /* ////////////////////////////////////////////////////////////////////////////////////// //////////////////////// Used to Read Event Map of detector ////////////////////////// ////////////////////////////////////////////////////////////////////////////////////// @@ -559,7 +619,7 @@ void CsI::ReadSensitive(const G4Event* event){ G4int ETrackID = Energy_itr->first - N ; G4double E = *(Energy_itr->second) ; if (ETrackID == NTrackID) { - m_Event->SetEEnergy(RandGauss::shoot(E, E*ResoEnergy/100./2.35)) ; + m_Event->SetEEnergy(RandGauss::shoot(E, E*ResoEnergy/100./2.35)) ; //m_Event->SetEEnergy(RandGauss::shoot(E, ResoEnergy)) ; } Energy_itr++; @@ -596,28 +656,48 @@ void CsI::ReadSensitive(const G4Event* event){ TimeHitMap ->clear() ; DetectorNumberHitMap ->clear() ; EnergyHitMap ->clear() ; - PosZHitMap ->clear() ; + PosZHitMap ->clear() ; + */ } //////////////////////////////////////////////////////////////// void CsI::InitializeScorers() { - bool already_exist = false; - m_CsIScorer = CheckScorer("CsIScorer",already_exist) ; + bool already_exist = false; + vector<G4int> NestingLevel; + NestingLevel.push_back(0); + NestingLevel.push_back(1); + m_CsIScorer = CheckScorer("CsIScorer",already_exist) ; + m_PDScorer = CheckScorer("PDScorer",already_exist) ; - if(already_exist) return ; - - G4VPrimitiveScorer* DetNbr = new PSDetectorNumber("CsINumber","CsI", 0) ; - G4VPrimitiveScorer* Energy = new PSEnergy("Energy","CsI", 0) ; - G4VPrimitiveScorer* Time = new PSTOF("Time","CsI", 0) ; - G4VPrimitiveScorer* InteractionCoordinatesZ = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","CsI", 0); - //and register it to the multifunctionnal detector - m_CsIScorer->RegisterPrimitive(DetNbr) ; - m_CsIScorer->RegisterPrimitive(Energy) ; - m_CsIScorer->RegisterPrimitive(Time) ; - m_CsIScorer->RegisterPrimitive(InteractionCoordinatesZ); - G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ; + if(already_exist) return ; + + G4VPrimitiveScorer* CsIScorer= new CALORIMETERSCORERS::PS_Calorimeter("CsI",NestingLevel); + m_CsIScorer->RegisterPrimitive(CsIScorer); + + G4VPrimitiveScorer* PDScorer = new SILICONSCORERS::PS_Silicon_Rectangle("PhotoDiode",0, + PhotoDiodeFace, + PhotoDiodeFace, + 1, + 1); + m_PDScorer->RegisterPrimitive(PDScorer); + + G4SDManager::GetSDMpointer()->AddNewDetector(m_PDScorer) ; + G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ; + + + /*G4VPrimitiveScorer* DetNbr = new PSDetectorNumber("CsINumber","CsI", 0) ; + G4VPrimitiveScorer* Energy = new PSEnergy("Energy","CsI", 0) ; + G4VPrimitiveScorer* Time = new PSTOF("Time","CsI", 0) ; + G4VPrimitiveScorer* InteractionCoordinatesZ = new OBSOLETEGENERALSCORERS::PSInteractionCoordinatesZ("InterCoordZ","CsI", 0); + //and register it to the multifunctionnal detector + m_CsIScorer->RegisterPrimitive(DetNbr) ; + m_CsIScorer->RegisterPrimitive(Energy) ; + m_CsIScorer->RegisterPrimitive(Time) ; + m_CsIScorer->RegisterPrimitive(InteractionCoordinatesZ); + G4SDManager::GetSDMpointer()->AddNewDetector(m_CsIScorer) ; + G4SDManager::GetSDMpointer()->AddNewDetector(m_PDScorer) ;*/ } //////////////////////////////////////////////////////////////// diff --git a/NPSimulation/CsI/CsI.hh b/NPSimulation/CsI/CsI.hh index dbeeb9c13..0072df73c 100644 --- a/NPSimulation/CsI/CsI.hh +++ b/NPSimulation/CsI/CsI.hh @@ -57,6 +57,10 @@ public: //////// Specific Function of this Class /////////// //////////////////////////////////////////////////// public: + + G4double PhotoDiodeFace = 18.;//mm + G4double PhotoDiodeThickness = 3.;//mm + // Cylindric CsI void AddCsI( G4double R , G4double Theta , @@ -101,8 +105,9 @@ public: // Scorer // Initialize all Scorer used by the MUST2Array void InitializeScorers() ; - // Silicon Associate Scorer - G4MultiFunctionalDetector* m_CsIScorer ; + // Silicon Associate Scorer + G4MultiFunctionalDetector* m_CsIScorer ; + G4MultiFunctionalDetector* m_PDScorer; //////////////////////////////////////////////////// ///////////Event class to store Data//////////////// //////////////////////////////////////////////////// diff --git a/NPSimulation/PhysicsListOption.txt b/NPSimulation/PhysicsListOption.txt index 45f0b7032..768664c87 100644 --- a/NPSimulation/PhysicsListOption.txt +++ b/NPSimulation/PhysicsListOption.txt @@ -2,6 +2,6 @@ IonBinaryCascadePhysics 0 EmExtraPhysics 0 HadronElasticPhysics 0 StoppingPhysics 0 -OpticalPhysics 0 +OpticalPhysics 1 HadronPhysicsQGSP_BIC_HP 0 -Decay 0 +Decay 1 diff --git a/NPSimulation/SSSD/SSSD.hh b/NPSimulation/SSSD/SSSD.hh index 3e3021ddf..97bd7d838 100644 --- a/NPSimulation/SSSD/SSSD.hh +++ b/NPSimulation/SSSD/SSSD.hh @@ -44,7 +44,7 @@ using namespace CLHEP; namespace SSSD_LOCAL{ // Energy and time Resolution const G4double ResoTime = 0 ; - const G4double ResoEnergy = 0.064 ;// = 150keV of Resolution // Unit is MeV/2.35 + const G4double ResoEnergy = 0.050 ;// = 150keV of Resolution // Unit is MeV/2.35 const G4double EnergyThreshold = 100*keV; // Geometry const G4double DetectorSize = 68*mm ; -- GitLab