diff --git a/Inputs/DetectorConfiguration/Riken_65mm.detector b/Inputs/DetectorConfiguration/Riken_65mm.detector index 31d9cbda56cb74ee9702fc462440060d3930e6b1..4e86c36cce8c2061491072f5dfcdecce16767416 100644 --- a/Inputs/DetectorConfiguration/Riken_65mm.detector +++ b/Inputs/DetectorConfiguration/Riken_65mm.detector @@ -115,7 +115,7 @@ CSI= 0 VIS= all %%%%%%%%%%%%%%%%%%%% -AddThinSi +%AddThinSi %%%%%%%%% Det 1 %%%%%%%% ThinSi A= 17.61 9.85 104.11 diff --git a/NPAnalysis/10He_Riken/Analysis b/NPAnalysis/10He_Riken/Analysis index 59f527f9d00c3ebb0de4ab83fe2eb345fcc3f2f7..f0612a01b99c8a28060b8aa7d3b0266620d96a3a 100755 Binary files a/NPAnalysis/10He_Riken/Analysis and b/NPAnalysis/10He_Riken/Analysis differ diff --git a/NPAnalysis/10He_Riken/include/ObjectManager.hh b/NPAnalysis/10He_Riken/include/ObjectManager.hh index 9f84a86bd3dbeee5fc2b356fa15e4b145bb623b1..7898e3e60608ce3982bfd338cf645181faa36077 100644 --- a/NPAnalysis/10He_Riken/include/ObjectManager.hh +++ b/NPAnalysis/10He_Riken/include/ObjectManager.hh @@ -109,7 +109,7 @@ namespace ENERGYLOSS // 3He Energy Loss EnergyLoss He3TargetWind = EnergyLoss ( "3He_Mylar.txt" , 1000 , - 3 , + 3 , 3 ); EnergyLoss He3TargetGaz = EnergyLoss ( "3He_D2solid_1b_26K.txt" , diff --git a/NPAnalysis/10He_Riken/src/Analysis.cc b/NPAnalysis/10He_Riken/src/Analysis.cc index e6b193ec156b4b8a399d53eff8d9ffae66b23916..4009505cb2da39ef76bd6f432c6c1ad4833ee3b0 100644 --- a/NPAnalysis/10He_Riken/src/Analysis.cc +++ b/NPAnalysis/10He_Riken/src/Analysis.cc @@ -73,7 +73,8 @@ int main(int argc,char** argv) // Get Must2 Pointer: MUST2Array* M2 = (MUST2Array*) myDetector -> m_Detector["MUST2"] ; - + // Allow directe acces to MUST2 Physics event + TMust2Physics* M2Physics = M2 -> GetPhysics(); cout << " ///////// Starting Analysis ///////// "<< endl << endl ; int i ,N=Chain -> GetEntries(); @@ -115,6 +116,8 @@ int main(int argc,char** argv) myDetector -> BuildPhysicalEvent() ; //// + + // Target (from initial condition) XTarget = Init->GetICPositionX(0); YTarget = Init->GetICPositionY(0); @@ -127,7 +130,7 @@ int main(int argc,char** argv) // Must 2 And ThinSi // for(int hit = 0; hit < M2 -> GetEventMultiplicity() ; hit ++) { - ELab[hit] = M2 -> GetEnergyDeposit(hit); + // ELab[hit] = M2 -> GetEnergyDeposit(hit); TVector3 HitDirection = M2 -> GetPositionOfInteraction(hit) - TVector3(XTarget,YTarget,0); // Angle between beam and particle @@ -140,10 +143,12 @@ int main(int argc,char** argv) if(M2 -> GetPositionOfInteraction(hit).Z()>0) { if(ELab[hit]>-1000 && ThinSi_E>0 ) + if(M2Physics.CsI.size) { ELab[hit]= He3StripAl.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle 2*0.4*micrometer , // Target Thickness at 0 degree ThetaMM2Surface ); + ELab[hit]= He3StripSi.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle 20*micrometer , // Target Thickness at 0 degree @@ -172,7 +177,8 @@ int main(int argc,char** argv) else if(ELab[hit]>-1000 ) { - if(ELab[hit]>21.66)//CsI are inside a Mylar foil, plus rear alu strip + + if(ELab[hit]>21.66)//CsI are inside a Mylar foil, plus rear alu strip { ELab[hit]= He3TargetWind.EvaluateInitialEnergy( ELab[hit] , // Energy of the detected particle 3*micrometer , // Target Thickness at 0 degree @@ -180,8 +186,7 @@ int main(int argc,char** argv) 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 ); @@ -259,7 +264,6 @@ int main(int argc,char** argv) return 0 ; } - double ThetaCalculation (TVector3 A , TVector3 B) { double Theta = acos( (A.Dot(B)) / (A.Mag()*B.Mag()) ) ; diff --git a/NPLib/DummyDetector/TDUMMYDetectorDataDict.cxx b/NPLib/DummyDetector/TDUMMYDetectorDataDict.cxx index f0f60697fba77a65fb00eca509cec583b65b0d74..00fdadd5b61ea753aab2d008b99d740eea487897 100644 --- a/NPLib/DummyDetector/TDUMMYDetectorDataDict.cxx +++ b/NPLib/DummyDetector/TDUMMYDetectorDataDict.cxx @@ -1,5 +1,5 @@ // -// File generated by rootcint at Fri Oct 16 12:56:33 2009 +// File generated by rootcint at Tue Oct 20 12:17:59 2009 // Do NOT change. Changes will be lost next time file is generated // diff --git a/NPLib/InitialConditions/TInitialConditions.h b/NPLib/InitialConditions/TInitialConditions.h index 96d8e9cec26be5fa2b1a2b918fd80e2ac7d31c01..468e941d2ea92152d55524ceccf0cda2e1089a47 100644 --- a/NPLib/InitialConditions/TInitialConditions.h +++ b/NPLib/InitialConditions/TInitialConditions.h @@ -28,7 +28,6 @@ #include <vector> #include "TObject.h" - using namespace std ; diff --git a/NPLib/MUST2/Must2Array.h b/NPLib/MUST2/Must2Array.h index 6e91fcab2b40baec023f74e9321dd3c0eeef18ab..d8babeea6bfb4596faa3f2566ab857ef198bcd26 100644 --- a/NPLib/MUST2/Must2Array.h +++ b/NPLib/MUST2/Must2Array.h @@ -98,11 +98,14 @@ class MUST2Array : public NPA::VDetector TVector3 GetPositionOfInteraction(int i) ; TVector3 GetTelescopeNormal(int i) ; + + TMust2Physics* GetPhysics() {return EventPhysics ;}; + void Print() ; private: // Root Input and Output tree classes - TMust2Data* EventData ; + TMust2Data* EventData ; TMust2Physics* EventPhysics ; diff --git a/NPLib/MUST2/TMust2Data.h b/NPLib/MUST2/TMust2Data.h index 04393789731e33d124dab9a8bdb698438b1f35cf..654e999cdb441c7e62ae992f3b80cef83acd05fd 100644 --- a/NPLib/MUST2/TMust2Data.h +++ b/NPLib/MUST2/TMust2Data.h @@ -139,12 +139,12 @@ class TMust2Data : public TObject { // SiLi //(E) - UShort_t GetMMSiLiEMult() {return fMM_SiLiE_DetectorNbr.size();} + UShort_t GetMMSiLiEMult() {return fMM_SiLiE_DetectorNbr.size();} UShort_t GetMMSiLiEDetectorNbr(Int_t i) {return fMM_SiLiE_DetectorNbr.at(i);} UShort_t GetMMSiLiEPadNbr(Int_t i) {return fMM_SiLiE_PadNbr.at(i);} Double_t GetMMSiLiEEnergy(Int_t i) {return fMM_SiLiE_Energy.at(i);} //(T) - UShort_t GetMMSiLiTMult() {return fMM_SiLiT_DetectorNbr.size();} + UShort_t GetMMSiLiTMult() {return fMM_SiLiT_DetectorNbr.size();} UShort_t GetMMSiLiTDetectorNbr(Int_t i) {return fMM_SiLiT_DetectorNbr.at(i);} UShort_t GetMMSiLiTPadNbr(Int_t i) {return fMM_SiLiT_PadNbr.at(i);} Double_t GetMMSiLiTTime(Int_t i) {return fMM_SiLiT_Time.at(i);} diff --git a/NPLib/MUST2/TMust2Physics.cxx b/NPLib/MUST2/TMust2Physics.cxx index 137e8a86a7293f5063e8c45048785a6744ad6244..20d7240bede8fe74a335c3b65c2f55b5d340c93e 100644 --- a/NPLib/MUST2/TMust2Physics.cxx +++ b/NPLib/MUST2/TMust2Physics.cxx @@ -24,6 +24,9 @@ #include <iostream> #include <cmath> +// ROOT +#include "TVector2.h" + ClassImp(TMust2Physics) TMust2Physics::TMust2Physics() @@ -71,13 +74,15 @@ void TMust2Physics::BuildPhysicalEvent(TMust2Data* Data) TelescopeNumber.push_back(Data->GetMMStripXEDetectorNbr(0)) ; // Data->Get Max Energy - if(Data->GetMMStripXEEnergy(0) > Data->GetMMStripYEEnergy(0)) Si_E.push_back( Data->GetMMStripXEEnergy(0) ) ; - else Si_E.push_back( Data->GetMMStripYEEnergy(0) ) ; - + // if(Data->GetMMStripXEEnergy(0) > Data->GetMMStripYEEnergy(0)) Si_E.push_back( Data->GetMMStripXEEnergy(0) ) ; + // else Si_E.push_back( Data->GetMMStripYEEnergy(0) ) ; + Si_E.push_back( Data->GetMMStripXEEnergy(0) ) ; + // Data->Get Min Time - if(Data->GetMMStripXTTime(0) < Data->GetMMStripYTTime(0)) Si_T.push_back( Data->GetMMStripXTTime(0) ) ; - else Si_T.push_back( Data->GetMMStripYTTime(0) ) ; - + // if(Data->GetMMStripXTTime(0) < Data->GetMMStripYTTime(0)) Si_T.push_back( Data->GetMMStripXTTime(0) ) ; + // else Si_T.push_back( Data->GetMMStripYTTime(0) ) ; + Si_T.push_back( Data->GetMMStripXTTime(0) ) ; + Si_X.push_back( Data->GetMMStripXEStripNbr(0) ) ; Si_Y.push_back( Data->GetMMStripYEStripNbr(0) ) ; @@ -125,11 +130,17 @@ void TMust2Physics::BuildPhysicalEvent(TMust2Data* Data) } - if (!Check_SiLi && !Check_CsI ) TotalEnergy.push_back( Si_E.at(0) ); - else if (Check_SiLi && !Check_CsI ) TotalEnergy.push_back( + Si_E.at(0) + SiLi_E.at(0) ); - else if (Check_CsI && !Check_SiLi) TotalEnergy.push_back( CsI_E .at(0) + Si_E.at(0) ); - else if (Check_CsI && Check_SiLi) TotalEnergy.push_back( CsI_E .at(0) + Si_E.at(0) + SiLi_E.at(0) ); + if(!Check_SiLi && !Check_CsI ) { TotalEnergy.push_back( Si_E.at(0) ) ; + CsI_E.push_back(0) ; CsI_T.push_back(0) ; + SiLi_E.push_back(0) ; SiLi_T.push_back(0) ; } + + else if (Check_SiLi && !Check_CsI ) { TotalEnergy.push_back( Si_E.at(0) + SiLi_E.at(0) ) ; + CsI_E.push_back(0) ; CsI_T.push_back(0) ; } + + else if (Check_CsI && !Check_SiLi) { TotalEnergy.push_back( CsI_E .at(0) + Si_E.at(0) ) ; + SiLi_E.push_back(0) ; SiLi_T.push_back(0) ; } + TotalEnergy.push_back( CsI_E[0] + Si_E[0] + SiLi_E[0] ); return; } @@ -217,10 +228,9 @@ void TMust2Physics::BuildPhysicalEvent(TMust2Data* Data) && Data->GetMMCsIEEnergy(i) > CsI_E_Threshold ) { Check_CsI = true ; - CsI_E.push_back(Data->GetMMCsIEEnergy(i)) ; -// cout << Data->GetMMCsIEEnergy(i) << endl; + CsI_E.push_back(Data->GetMMCsIEEnergy(i)) ; CsI_N.push_back(Data->GetMMCsIECristalNbr(i)) ; - CsI_T.push_back(Data->GetMMCsITTime(i)) ; + CsI_T.push_back(Data->GetMMCsITTime(i)) ; } } @@ -232,9 +242,9 @@ void TMust2Physics::BuildPhysicalEvent(TMust2Data* Data) CsI_N.push_back(-2) ; } - TotalEnergy.push_back(Si_E.at(jj)) ; - if (Check_SiLi) TotalEnergy.at(jj) += SiLi_E.at(jj) ; - if (Check_CsI) TotalEnergy.at(jj) += CsI_E.at(jj) ; + TotalEnergy.push_back(Si_E[jj])) ; + if (Check_SiLi) TotalEnergy[jj] += SiLi_E[jj] ; + if (Check_CsI) TotalEnergy[jj] += CsI_E[jj] ; } } return; @@ -284,13 +294,279 @@ void TMust2Physics::BuildPhysicalEvent(TMust2Data* Data) } +///////////////////////////////////////////// +int TMust2Physics :: CheckEvent() + { + // Check the size of the different elements + if( Data->GetMMStripXEMult() == Data->GetMMStripYEMult() && Data->GetMMStripYEMult() == Data->GetMMStripXTMult() && Data->GetMMStripXTMult() == Data->GetMMStripYTMult() ) + + return 1 ; // Regular Event + + else if( Data->GetMMStripXEMult() == Data->GetMMStripYEMult()+1 || Data->GetMMStripXEMult() == Data->GetMMStripYEMult()-1 ) + + return 2 ; // Pseudo Event, potentially interstrip + + else return -1 ; // Rejected Event + + + return -1 ; // Rejected Event + } + +///////////////////////////////////////////// +bool TMust2Physics :: ResolvePseudoEvent() + { + + return false; + } + +///////////////////////////////////////////// +vector < TVector2 > TMust2Physics :: Match_X_Y() + { + vector < TVector2 > ArrayOfGoodCouple ; + + for(int i = 0 ; i < Data->GetMMStripXEMult(); i++) + { + for(int j = 0 ; j < Data->GetMMStripYEMult(); j++) + { + if ( Data->GetMMStripXEDetectorNbr(i) == Data->GetMMStripYEDetectorNbr(j) ) + { + if( fabs(Data->GetMMStripXEEnergy(i) - Data->GetMMStripYEEnergy(j))/Data->GetMMStripXEEnergy(i) < 0.1 ) + ArrayOfGoodCouple . push_back ( TVector2(i,j) ) ; + } + } + } + + if( ArrayOfGoodCouple.size() > Data->GetMMStripXEMult() ) ArrayOfGoodCouple.clear() ; + + return ArrayOfGoodCouple; + } + + +///////////////////////////////////////////// +bool TMust2Physics :: Match_Si_SiLi(TVector2 SiCouple) + { + + for(int i = 0 ; i < Data->GetMMSiLiEMult() ; i++ ) + { + if( Data->GetMMSiLiEPadNbr(i) == 1 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + + else if( Data->GetMMSiLiEPadNbr(i) == 2 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + + else if( Data->GetMMSiLiEPadNbr(i) == 3 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 4 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 5 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 6 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 7 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 8 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 9 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 10 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 11 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 12 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 13 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 14 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 15 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMSiLiEPadNbr(i) == 16 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + } + + return false; + } + + +///////////////////////////////////////////// +vector < TVector3 > TMust2Physics :: Match_Si_CsI() + { + for(int i = 0 ; i < Data->GetMMCsIEMult() ; i++ ) + { + if( Data->GetMMCsIECristalNbr(i) == 1 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + + else if( Data->GetMMCsIECristalNbr(i) == 2 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + + else if( Data->GetMMCsIECristalNbr(i) == 3 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 4 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 5 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 6 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 7 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 8 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 9 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 10 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 11 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 12 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 13 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 14 + && Data->GetMMStripXEStripNbr(Int_t i)<1 && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)<1 && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 15 + && Data->GetMMStripXEStripNbr(Int_t i)< && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)< && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + + else if( Data->GetMMCsIECristalNbr(i) == 16 + && Data->GetMMStripXEStripNbr(Int_t i)< && Data->GetMMStripXEStripNbr(Int_t i)>1 + && Data->GetMMStripYEStripNbr(Int_t i)< && Data->GetMMStripYEStripNbr(Int_t i)>1 ) + + return true ; + } + + + return false; + } + +///////////////////////////////////////////// void TMust2Physics::Clear() { EventMultiplicity= 0 ; + TelescopeNumber .clear() ; - EventType .clear() ; - TotalEnergy .clear() ; + EventType .clear() ; + TotalEnergy .clear() ; // Si X Si_E.clear() ; @@ -309,3 +585,5 @@ void TMust2Physics::Clear() CsI_N.clear() ; } + + diff --git a/NPLib/MUST2/TMust2Physics.h b/NPLib/MUST2/TMust2Physics.h index 22bd7b3b518042138485d27d55fb8e14b2925728..d2b52c6c149013b435169c1eaa54d7d54d43b29b 100644 --- a/NPLib/MUST2/TMust2Physics.h +++ b/NPLib/MUST2/TMust2Physics.h @@ -36,11 +36,12 @@ class TMust2Physics : public TObject public: void Clear() ; - void Clear(const Option_t*) {}; + void Clear(const Option_t*) {}; void BuildPhysicalEvent(TMust2Data* Data) ; void BuildSimplePhysicalEvent(TMust2Data* Data) ; public: + // Provide Physical Multiplicity Int_t EventMultiplicity ; diff --git a/NPSimulation/include/Must2Scorers.hh b/NPSimulation/include/Must2Scorers.hh index 098cdbd38ba39127e3066f06921cb74c52bba236..b9aac0754c98eac26bbf3427ff18c855cac55e04 100644 --- a/NPSimulation/include/Must2Scorers.hh +++ b/NPSimulation/include/Must2Scorers.hh @@ -78,6 +78,30 @@ namespace MUST2 { G4THitsMap<G4double>* EvtMap; }; + + class PSPadOrCristalNumber : public G4VPrimitiveScorer + { + + public: // with description + PSPadOrCristalNumber(G4String name, G4int depth = 0, G4String Type = "XXX"); + virtual ~PSPadOrCristalNumber(); + + protected: // with description + virtual G4bool ProcessHits(G4Step*, G4TouchableHistory*); + + public: + virtual void Initialize(G4HCofThisEvent*); + virtual void EndOfEvent(G4HCofThisEvent*); + virtual void clear(); + virtual void DrawAll(); + virtual void PrintAll(); + + private: + bool m_type; + G4int HCID; + G4THitsMap<G4double>* EvtMap; + }; + } #endif diff --git a/NPSimulation/src/EventGeneratorTransfert.cc b/NPSimulation/src/EventGeneratorTransfert.cc index d0b3b32b6098f17e72b93da3461e5abc02937c5d..ef6a7a4e2104ac3a331b6b551b26438ba5e0ef11 100644 --- a/NPSimulation/src/EventGeneratorTransfert.cc +++ b/NPSimulation/src/EventGeneratorTransfert.cc @@ -350,9 +350,9 @@ void EventGeneratorTransfert::GenerateEvent(G4Event* anEvent , G4ParticleGun* pa G4double x0 = InterCoord.x(); G4double y0 = InterCoord.y(); G4double z0 = InterCoord.z(); - m_InitConditions->SetICPositionX(x0); - m_InitConditions->SetICPositionY(y0); - m_InitConditions->SetICPositionZ(z0); + m_InitConditions->SetICPositionX(x0);// + m_InitConditions->SetICPositionY(y0);// + m_InitConditions->SetICPositionZ(z0);// // write emittance angles to ROOT file m_InitConditions->SetICIncidentEmittanceTheta(Beam_thetaX / deg); diff --git a/NPSimulation/src/EventGeneratorTransfertToResonance.cc b/NPSimulation/src/EventGeneratorTransfertToResonance.cc index e1229e51fa4cab33b29da769b4b6f00a3bd82724..8454e90d98970da5029ec465609762ae1ea9823a 100644 --- a/NPSimulation/src/EventGeneratorTransfertToResonance.cc +++ b/NPSimulation/src/EventGeneratorTransfertToResonance.cc @@ -410,6 +410,7 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part else InterCoord = G4ThreeVector(0,0,0); + // write vertex position to ROOT file G4double x0 = InterCoord.x(); G4double y0 = InterCoord.y(); @@ -472,17 +473,6 @@ void EventGeneratorTransfertToResonance::GenerateEvent(G4Event* anEvent , G4Part sin(ThetaHeavy) * sin(phi), cos(ThetaHeavy)); - - // Move to the target -if( m_Target !=0) { - x0 += m_Target->GetTargetX() ; - y0 += m_Target->GetTargetY() ; - z0 += m_Target->GetTargetZ() ; - } - // write vertex position to ROOT file - m_InitConditions->SetICPositionX(x0); - m_InitConditions->SetICPositionY(y0); - m_InitConditions->SetICPositionZ(z0); ////////////////////////////////////////////////// ///////// Set up everything for shooting ///////// ////////////////////////////////////////////////// diff --git a/NPSimulation/src/GeneralScorers.cc b/NPSimulation/src/GeneralScorers.cc index 7c260c6df2921628d7b9305a522b7847bb38aa40..4ab30937e05ca8c07907f75ff5507b90f0a40d12 100644 --- a/NPSimulation/src/GeneralScorers.cc +++ b/NPSimulation/src/GeneralScorers.cc @@ -36,7 +36,6 @@ int GENERALSCORERS::PickUpDetectorNumber(G4Step* aStep, std::string DetName) { std::string name = aStep->GetTrack()->GetVolume()->GetName(); std::string nbr ; - size_t start, end ; start = name.find(DetName) + DetName.length(); end = name.find("_"); diff --git a/NPSimulation/src/MUST2Array.cc b/NPSimulation/src/MUST2Array.cc index 83c65039c50c7cee9ee430608ef34dee3d4b57fc..d8420df91dd0ad962e8938b4ad661c75d5231004 100644 --- a/NPSimulation/src/MUST2Array.cc +++ b/NPSimulation/src/MUST2Array.cc @@ -342,15 +342,15 @@ void MUST2Array::VolumeMaker(G4int TelescopeNumber , -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_LT , logicSiLi_LT , Name + "_SiLi_Pad1" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RT , logicSiLi_RT , Name + "_SiLi_Pad2" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LC1 , logicSiLi_LC1 , Name + "_SiLi_Pad3" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RC1 , logicSiLi_RC1 , Name + "_SiLi_Pad4" , 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) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LB , logicSiLi_LB , Name + "_SiLi_Pad5" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RB , logicSiLi_RB , Name + "_SiLi_Pad6" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_LC2 , logicSiLi_LC2 , Name + "_SiLi_Pad7" , logicSiLi , false , 0) ; + PVPBuffer = new G4PVPlacement(0 , positionSiLi_RC2 , logicSiLi_RC2 , Name + "_SiLi_Pad8" , logicSiLi , false , 0) ; logicSiLi->SetVisAttributes(G4VisAttributes(G4Colour(1, 1., 1.))); @@ -812,15 +812,15 @@ void MUST2Array::ReadConfiguration(string Path) // Called After DetecorConstruction::AddDetector Method void MUST2Array::ConstructDetector(G4LogicalVolume* world) { - G4RotationMatrix* MMrot = NULL ; - G4ThreeVector MMpos = G4ThreeVector(0, 0, 0) ; + G4RotationMatrix* MMrot = NULL ; + 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 MMCenter = G4ThreeVector(0, 0, 0) ; - bool Si = true ; - bool SiLi = true ; - bool CsI = true ; + G4ThreeVector MMCenter = G4ThreeVector(0, 0, 0) ; + bool Si = true ; + bool SiLi = true ; + bool CsI = true ; G4int NumberOfTelescope = m_DefinitionType.size() ; @@ -923,102 +923,115 @@ void MUST2Array::ReadSensitive(const G4Event* event) // Si std::map<G4int, G4int*>::iterator DetectorNumber_itr ; - std::map<G4int, G4double*>::iterator Energy_itr ; - std::map<G4int, G4double*>::iterator Time_itr ; - std::map<G4int, G4double*>::iterator X_itr ; - std::map<G4int, G4double*>::iterator Y_itr ; - std::map<G4int, G4double*>::iterator Pos_X_itr ; - std::map<G4int, G4double*>::iterator Pos_Y_itr ; - std::map<G4int, G4double*>::iterator Pos_Z_itr ; + std::map<G4int, G4double*>::iterator Energy_itr ; + std::map<G4int, G4double*>::iterator Time_itr ; + std::map<G4int, G4double*>::iterator X_itr ; + std::map<G4int, G4double*>::iterator Y_itr ; + std::map<G4int, G4double*>::iterator Pos_X_itr ; + std::map<G4int, G4double*>::iterator Pos_Y_itr ; + std::map<G4int, G4double*>::iterator Pos_Z_itr ; std::map<G4int, G4double*>::iterator Ang_Theta_itr ; std::map<G4int, G4double*>::iterator Ang_Phi_itr ; - G4THitsMap<G4int>* DetectorNumberHitMap ; - G4THitsMap<G4double>* EnergyHitMap ; - G4THitsMap<G4double>* TimeHitMap ; - G4THitsMap<G4double>* XHitMap ; - G4THitsMap<G4double>* YHitMap ; - G4THitsMap<G4double>* PosXHitMap ; - G4THitsMap<G4double>* PosYHitMap ; - G4THitsMap<G4double>* PosZHitMap ; - G4THitsMap<G4double>* AngThetaHitMap ; - G4THitsMap<G4double>* AngPhiHitMap ; - -// NULL pointer are given to avoid warning at compilation + G4THitsMap<G4int>* DetectorNumberHitMap ; + G4THitsMap<G4double>* EnergyHitMap ; + G4THitsMap<G4double>* TimeHitMap ; + G4THitsMap<G4double>* XHitMap ; + G4THitsMap<G4double>* YHitMap ; + G4THitsMap<G4double>* PosXHitMap ; + G4THitsMap<G4double>* PosYHitMap ; + G4THitsMap<G4double>* PosZHitMap ; + G4THitsMap<G4double>* AngThetaHitMap ; + G4THitsMap<G4double>* AngPhiHitMap ; // Si(Li) std::map<G4int, G4double*>::iterator SiLiEnergy_itr ; - G4THitsMap<G4double>* SiLiEnergyHitMap ; + std::map<G4int, G4double*>::iterator SiLiPadNbr_itr ; + G4THitsMap<G4double>* SiLiEnergyHitMap ; + G4THitsMap<G4double>* SiLiPadNbrHitMap ; + // CsI std::map<G4int, G4double*>::iterator CsIEnergy_itr ; - G4THitsMap<G4double>* CsIEnergyHitMap ; + std::map<G4int, G4double*>::iterator CsICristalNbr_itr ; + G4THitsMap<G4double>* CsIEnergyHitMap ; + G4THitsMap<G4double>* CsICristalNbrHitMap ; ////////////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////// // Read the Scorer associate to the Silicon Strip //Detector Number - G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/DetectorNumber") ; - DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ; - DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ; + G4int StripDetCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/DetectorNumber") ; + DetectorNumberHitMap = (G4THitsMap<G4int>*)(event->GetHCofThisEvent()->GetHC(StripDetCollectionID)) ; + DetectorNumber_itr = DetectorNumberHitMap->GetMap()->begin() ; //Energy - G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripEnergy") ; - EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ; - Energy_itr = EnergyHitMap->GetMap()->begin() ; + G4int StripEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripEnergy") ; + EnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripEnergyCollectionID)) ; + Energy_itr = EnergyHitMap->GetMap()->begin() ; //Time of Flight - G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripTime") ; - TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ; - Time_itr = TimeHitMap->GetMap()->begin() ; + G4int StripTimeCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripTime") ; + TimeHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripTimeCollectionID)) ; + Time_itr = TimeHitMap->GetMap()->begin() ; //Strip Number X - G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberX") ; - XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ; - X_itr = XHitMap->GetMap()->begin() ; + G4int StripXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberX") ; + XHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripXCollectionID)) ; + X_itr = XHitMap->GetMap()->begin() ; //Strip Number Y - G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberY") ; - YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ; - Y_itr = YHitMap->GetMap()->begin() ; + G4int StripYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/StripNumberY") ; + YHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(StripYCollectionID)) ; + Y_itr = YHitMap->GetMap()->begin() ; //Interaction Coordinate X G4int InterCoordXCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordX") ; - PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ; - Pos_X_itr = PosXHitMap->GetMap()->begin() ; + PosXHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordXCollectionID)) ; + Pos_X_itr = PosXHitMap->GetMap()->begin() ; //Interaction Coordinate Y G4int InterCoordYCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordY") ; - PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ; - Pos_Y_itr = PosYHitMap->GetMap()->begin() ; + PosYHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordYCollectionID)) ; + Pos_Y_itr = PosYHitMap->GetMap()->begin() ; //Interaction Coordinate Z G4int InterCoordZCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordZ") ; - PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ; - Pos_Z_itr = PosXHitMap->GetMap()->begin() ; + PosZHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordZCollectionID)) ; + Pos_Z_itr = PosXHitMap->GetMap()->begin() ; //Interaction Coordinate Angle Theta - G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngTheta") ; + G4int InterCoordAngThetaCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngTheta"); AngThetaHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngThetaCollectionID)) ; - Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ; + Ang_Theta_itr = AngThetaHitMap->GetMap()->begin() ; //Interaction Coordinate Angle Phi - G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngPhi") ; - AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)) ; - Ang_Phi_itr = AngPhiHitMap->GetMap()->begin() ; + G4int InterCoordAngPhiCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_StripScorer/InterCoordAngPhi") ; + AngPhiHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(InterCoordAngPhiCollectionID)) ; + Ang_Phi_itr = AngPhiHitMap->GetMap()->begin() ; // Read the Scorer associate to the SiLi //Energy - G4int SiLiEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiEnergy") ; - SiLiEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiEnergyCollectionID)) ; - SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ; - + G4int SiLiEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiEnergy") ; + SiLiEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiEnergyCollectionID)) ; + SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ; + G4int SiLiPadNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_SiLiScorer/SiLiPadNbr") ; + SiLiPadNbrHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(SiLiPadNbrCollectionID)) ; + SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin() ; + // Read the Scorer associate to the CsI crystal //Energy - G4int CsIEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIEnergy") ; - CsIEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIEnergyCollectionID)) ; - CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ; + G4int CsIEnergyCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsIEnergy") ; + CsIEnergyHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsIEnergyCollectionID)) ; + CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ; + + G4int CsICristalNbrCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID("MUST2_CsIScorer/CsICristalNbr") ; + CsICristalNbrHitMap = (G4THitsMap<G4double>*)(event->GetHCofThisEvent()->GetHC(CsICristalNbrCollectionID)) ; + CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin() ; + + +///////////////////// G4int sizeN = DetectorNumberHitMap->entries() ; G4int sizeE = EnergyHitMap->entries() ; @@ -1153,39 +1166,48 @@ void MUST2Array::ReadSensitive(const G4Event* event) // Si(Li) SiLiEnergy_itr = SiLiEnergyHitMap->GetMap()->begin() ; + SiLiPadNbr_itr = SiLiPadNbrHitMap->GetMap()->begin() ; + for (G4int h = 0 ; h < SiLiEnergyHitMap->entries() ; h++) { G4int SiLiEnergyTrackID = SiLiEnergy_itr->first -N ; G4double SiLiEnergy = *(SiLiEnergy_itr->second) ; + G4int PadNbr = *(SiLiPadNbr_itr->second) ; + if (SiLiEnergyTrackID == NTrackID) { m_Event->SetMMSiLiEEnergy(RandGauss::shoot(SiLiEnergy, ResoSiLi)) ; - m_Event->SetMMSiLiEPadNbr(1); - m_Event->SetMMSiLiTPadNbr(1); + m_Event->SetMMSiLiEPadNbr(PadNbr); + m_Event->SetMMSiLiTPadNbr(PadNbr); m_Event->SetMMSiLiTTime(1); m_Event->SetMMSiLiTDetectorNbr(N); m_Event->SetMMSiLiEDetectorNbr(N); } - SiLiEnergy_itr++; + SiLiEnergy_itr++ ; + SiLiPadNbr_itr++ ; } // CsI - CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ; + CsIEnergy_itr = CsIEnergyHitMap->GetMap()->begin() ; + CsICristalNbr_itr = CsICristalNbrHitMap->GetMap()->begin() ; + for (G4int h = 0 ; h < CsIEnergyHitMap->entries() ; h++) { - G4int CsIEnergyTrackID = CsIEnergy_itr->first -N ; - G4double CsIEnergy = *(CsIEnergy_itr->second) ; - + G4int CsIEnergyTrackID = CsIEnergy_itr->first-N ; + G4double CsIEnergy = *(CsIEnergy_itr->second) ; + + G4int CristalNumber = *(CsICristalNbr_itr->second) ; if (CsIEnergyTrackID == NTrackID) { m_Event->SetMMCsIEEnergy(RandGauss::shoot(CsIEnergy, ResoCsI*sqrt(CsIEnergy))); - m_Event->SetMMCsIECristalNbr(1); - m_Event->SetMMCsITCristalNbr(1); + m_Event->SetMMCsIECristalNbr(CristalNumber); + m_Event->SetMMCsITCristalNbr(CristalNumber); m_Event->SetMMCsITTime(1); m_Event->SetMMCsITDetectorNbr(N); m_Event->SetMMCsIEDetectorNbr(N); } - CsIEnergy_itr++; + CsIEnergy_itr++ ; + CsICristalNbr_itr++ ; } @@ -1201,8 +1223,10 @@ void MUST2Array::ReadSensitive(const G4Event* event) TimeHitMap ->clear() ; XHitMap ->clear() ; YHitMap ->clear() ; - CsIEnergyHitMap ->clear() ; SiLiEnergyHitMap ->clear() ; + SiLiPadNbrHitMap ->clear() ; + CsIEnergyHitMap ->clear() ; + CsICristalNbrHitMap ->clear() ; PosXHitMap ->clear() ; PosYHitMap ->clear() ; PosZHitMap ->clear() ; @@ -1247,12 +1271,16 @@ void MUST2Array::InitializeScorers() // SiLi Associate Scorer m_SiLiScorer = new G4MultiFunctionalDetector("MUST2_SiLiScorer") ; G4VPrimitiveScorer* SiLiEnergy = new GENERALSCORERS::PSEnergy("SiLiEnergy","MUST2Telescope", 0) ; + G4VPrimitiveScorer* SiLiPadNbr = new PSPadOrCristalNumber("SiLiPadNbr",0,"SiLi") ; m_SiLiScorer->RegisterPrimitive(SiLiEnergy) ; + m_SiLiScorer->RegisterPrimitive(SiLiPadNbr) ; // CsI Associate Scorer m_CsIScorer = new G4MultiFunctionalDetector("MUST2_CsIScorer") ; G4VPrimitiveScorer* CsIEnergy = new GENERALSCORERS::PSEnergy("CsIEnergy","MUST2Telescope", 0) ; + G4VPrimitiveScorer* CsICristalNbr = new PSPadOrCristalNumber("CsICristalNbr",0,"CsI") ; m_CsIScorer->RegisterPrimitive(CsIEnergy) ; + m_CsIScorer->RegisterPrimitive(CsICristalNbr) ; // Add All Scorer to the Global Scorer Manager G4SDManager::GetSDMpointer()->AddNewDetector(m_StripScorer) ; diff --git a/NPSimulation/src/Must2Scorers.cc b/NPSimulation/src/Must2Scorers.cc index ea0c60b30498b34ae1588fa2788af9bdc76ab46a..39c58d7590e4e2e831bba4aac13d711300f9d54b 100644 --- a/NPSimulation/src/Must2Scorers.cc +++ b/NPSimulation/src/Must2Scorers.cc @@ -171,3 +171,103 @@ void PSStripNumberY::PrintAll() G4cout << " PrimitiveScorer " << GetName() << G4endl ; G4cout << " Number of entries " << EvtMap->entries() << G4endl ; } + + + + + + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//CsI Cristal / SiLi Pad Number Scorer +// +PSPadOrCristalNumber::PSPadOrCristalNumber(G4String name, G4int depth,G4String type) + : G4VPrimitiveScorer(name, depth), HCID(-1) +{ + if (type=="SiLi") m_type = true ; + else if (type=="CsI" ) m_type = false ; + else G4cout << "Problem in MUST2 Scorer definition: Type should be SiLi or CsI" << G4endl ; +} + +PSPadOrCristalNumber::~PSPadOrCristalNumber() +{ + ; +} + +G4bool PSPadOrCristalNumber::ProcessHits(G4Step* aStep, G4TouchableHistory*) +{ + std::string name = aStep->GetTrack()->GetVolume()->GetName(); + std::string nbr ; + unsigned int numberOfCharacterInDetectorNumber ; + + if(m_type)// 24 character before pad number MUST2Telescope4_SiLi_PadXX + { + numberOfCharacterInDetectorNumber = name.length() - 24 ; + for(unsigned int i = 24 ; i < 24 + numberOfCharacterInDetectorNumber ; i++ ) + nbr += name[i] ; + } + + else // 27 character before cristal number : MUST2Telescope4_CsI_CristalXX + { + numberOfCharacterInDetectorNumber = name.length() - 27 ; + for(unsigned int i = 27 ; i < 27 + numberOfCharacterInDetectorNumber ; i++ ) + nbr += name[i] ; + } + + + double VolumeNumber = atoi( nbr.c_str() ); + + int DetNbr = GENERALSCORERS::PickUpDetectorNumber(aStep, "MUST2Telescope"); + + G4double edep = aStep->GetTotalEnergyDeposit(); + if (edep < 100*keV) return FALSE; + G4int index = aStep->GetTrack()->GetTrackID(); + EvtMap->set(index+DetNbr, VolumeNumber); + return TRUE; +} + +void PSPadOrCristalNumber::Initialize(G4HCofThisEvent* HCE) +{ + EvtMap = new G4THitsMap<G4double>(GetMultiFunctionalDetector()->GetName(), GetName()); + if (HCID < 0) { + HCID = GetCollectionID(0); + } + HCE->AddHitsCollection(HCID, (G4VHitsCollection*)EvtMap); +} + +void PSPadOrCristalNumber::EndOfEvent(G4HCofThisEvent*) +{ + ; +} + +void PSPadOrCristalNumber::clear() +{ + EvtMap->clear(); +} + +void PSPadOrCristalNumber::DrawAll() +{ + ; +} + +void PSPadOrCristalNumber::PrintAll() +{ + G4cout << " MultiFunctionalDet " << detector->GetName() << G4endl ; + G4cout << " PrimitiveScorer " << GetName() << G4endl ; + G4cout << " Number of entries " << EvtMap->entries() << G4endl ; +} + + + + + + + + + + + + + + +