From 1a22548383c5945c9d4185eb113e1064e815fd82 Mon Sep 17 00:00:00 2001 From: sparzych Date: Thu, 12 Jan 2023 11:20:34 +0100 Subject: [PATCH 01/15] Update for saving information about prompt scatters in phantom and crystals --- source/digits_hits/include/GateToRoot.hh | 6 +- source/digits_hits/src/GateAnalysis.cc | 47 +++++ source/digits_hits/src/GateToRoot.cc | 185 ++++++++++-------- .../src/GateTrajectoryNavigator.cc | 65 +++++- 4 files changed, 221 insertions(+), 82 deletions(-) diff --git a/source/digits_hits/include/GateToRoot.hh b/source/digits_hits/include/GateToRoot.hh index 530ca17d5..ab9fa8a52 100644 --- a/source/digits_hits/include/GateToRoot.hh +++ b/source/digits_hits/include/GateToRoot.hh @@ -62,9 +62,9 @@ class GateTrajectoryNavigator; class ComptonRayleighData { public: - G4int photon1_phantom_Rayleigh, photon2_phantom_Rayleigh; - G4int photon1_phantom_compton, photon2_phantom_compton; - Char_t theComptonVolumeName1[60], theComptonVolumeName2[60], theRayleighVolumeName1[60], theRayleighVolumeName2[60]; + G4int photon1_phantom_Rayleigh,photon2_phantom_Rayleigh,photon3_phantom_Rayleigh; + G4int photon1_phantom_compton , photon2_phantom_compton,photon3_phantom_compton; + Char_t theComptonVolumeName1[60] , theComptonVolumeName2[60], theComptonVolumeName3[60], theRayleighVolumeName1[60], theRayleighVolumeName2[60], theRayleighVolumeName3[60]; ComptonRayleighData(ComptonRayleighData &); diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 7e0a53715..c04e8e00a 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -149,20 +149,26 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) //G4int positronID = 0; // no more needed G4int photon1ID = 0; G4int photon2ID = 0; + G4int photon3ID = 0; G4int rootID = 0; G4int primaryID = 0; G4int photon1_phantom_compton = 0; G4int photon2_phantom_compton = 0; + G4int photon3_phantom_compton = 0; + G4int photon1_crystal_compton = 0; G4int photon2_crystal_compton = 0; + G4int photon3_crystal_compton = 0; G4int photon1_phantom_Rayleigh = 0; G4int photon2_phantom_Rayleigh = 0; + G4int photon3_phantom_Rayleigh = 0; G4int photon1_crystal_Rayleigh = 0; G4int photon2_crystal_Rayleigh = 0; + G4int photon3_crystal_Rayleigh = 0; G4int septalNb = 0; // HDS : septal penetration @@ -197,6 +203,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) photon1ID = photonIDVec[0]; photon2ID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; + photon3ID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; } if (photon1ID == 0) @@ -208,6 +215,11 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (nVerboseLevel > 1) G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0\n"; } + if (photon3ID == 0) { + if (nVerboseLevel > 1) G4cout + << "GateAnalysis::RecordEndOfEvent : WARNING : photon3ID == 0\n"; + } + if (nVerboseLevel > 1) G4cout << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon1ID << " photon2ID : " << photon2ID << Gateendl; @@ -221,10 +233,12 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4String theComptonVolumeName("NULL"); G4String theComptonVolumeName1("NULL"); G4String theComptonVolumeName2("NULL"); + G4String theComptonVolumeName3("NULL"); G4String theRayleighVolumeName("NULL"); G4String theRayleighVolumeName1("NULL"); G4String theRayleighVolumeName2("NULL"); + G4String theRayleighVolumeName3("NULL"); for (G4int iPHit=0;iPHit 0) G4cout << "GateAnalysis::RecordEndOfEvent : photon2_phantom_compton : " << photon2_phantom_compton << Gateendl; } + if (phantomTrackID == photon3ID) + { + photon3_phantom_compton++; + theComptonVolumeName3 = theComptonVolumeName; + if (nVerboseLevel > 0) G4cout + << "GateAnalysis::RecordEndOfEvent : photon3_phantom_compton : " << photon3_phantom_compton << Gateendl; + } } // Counting Rayleigh scatter in phantom @@ -313,6 +334,13 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (nVerboseLevel > 0) G4cout << "GateAnalysis::RecordEndOfEvent : photon2_phantom_Rayleigh : " << photon2_phantom_Rayleigh << Gateendl; } + if (phantomTrackID == photon3ID) + { + photon3_phantom_Rayleigh++; + theRayleighVolumeName3 = theRayleighVolumeName; + if (nVerboseLevel > 0) G4cout + << "GateAnalysis::RecordEndOfEvent : photon3_phantom_Rayleigh : " << photon3_phantom_Rayleigh << Gateendl; + } } } // end loop NpHits @@ -325,12 +353,16 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) ComptonRayleighData aCRData; aCRData.photon1_phantom_Rayleigh = photon1_phantom_Rayleigh; aCRData.photon2_phantom_Rayleigh = photon2_phantom_Rayleigh; + aCRData.photon3_phantom_Rayleigh = photon3_phantom_Rayleigh; aCRData.photon1_phantom_compton = photon1_phantom_compton; aCRData.photon2_phantom_compton = photon2_phantom_compton; + aCRData.photon3_phantom_compton = photon3_phantom_compton; strcpy(aCRData.theComptonVolumeName1 , theComptonVolumeName1.c_str() ); strcpy(aCRData.theComptonVolumeName2 , theComptonVolumeName2.c_str() ); + strcpy(aCRData.theComptonVolumeName3 , theComptonVolumeName3.c_str() ); strcpy(aCRData.theRayleighVolumeName1 , theRayleighVolumeName1.c_str() ); strcpy(aCRData.theRayleighVolumeName2 , theRayleighVolumeName2.c_str() ); + strcpy(aCRData.theRayleighVolumeName3 , theRayleighVolumeName3.c_str() ); gateToRoot->RecordPHData( aCRData ); // return; } @@ -344,8 +376,10 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) photon1_phantom_Rayleigh += aCRData.photon1_phantom_Rayleigh; photon2_phantom_Rayleigh += aCRData.photon2_phantom_Rayleigh; + photon3_phantom_Rayleigh += aCRData.photon3_phantom_Rayleigh; photon1_phantom_compton += aCRData.photon1_phantom_compton; photon2_phantom_compton += aCRData.photon2_phantom_compton; + photon3_phantom_compton += aCRData.photon3_phantom_compton; /* if( theComptonVolumeName1 == G4String("NULL") ) {theComptonVolumeName1 = aCRData.theComptonVolumeName1;} if( theComptonVolumeName2 == G4String("NULL") ) {theComptonVolumeName2 = aCRData.theComptonVolumeName2;} @@ -354,8 +388,10 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) */ theComptonVolumeName1 = aCRData.theComptonVolumeName1; theComptonVolumeName2 = aCRData.theComptonVolumeName2; + theComptonVolumeName3 = aCRData.theComptonVolumeName3; theRayleighVolumeName1 = aCRData.theRayleighVolumeName1; theRayleighVolumeName2 = aCRData.theRayleighVolumeName2; + theRayleighVolumeName3 = aCRData.theRayleighVolumeName3; } @@ -381,6 +417,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (crystalTrackID == photon1ID) photon1_crystal_compton++; if (crystalTrackID == photon2ID) photon2_crystal_compton++; + if (crystalTrackID == photon3ID) photon3_crystal_compton++; } // Counting Rayleigh scatter in crystal @@ -389,6 +426,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (crystalTrackID == photon1ID) photon1_crystal_Rayleigh++; if (crystalTrackID == photon2ID) photon2_crystal_Rayleigh++; + if (crystalTrackID == photon3ID) photon3_crystal_Rayleigh++; } G4int PDGEncoding = (*CHC)[iHit]->GetPDGEncoding(); @@ -439,6 +477,15 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) nCrystalCompton = photon2_crystal_compton; nCrystalRayleigh = photon2_crystal_Rayleigh; } + else if (photonID == 3) + { + nPhantomCompton = photon3_phantom_compton; + nPhantomRayleigh = photon3_phantom_Rayleigh; + theComptonVolumeName = theComptonVolumeName3; + theRayleighVolumeName = theRayleighVolumeName3; + nCrystalCompton = photon3_crystal_compton; + nCrystalRayleigh = photon3_crystal_Rayleigh; + } // search the primary that originated the track primaryID = m_trajectoryNavigator->FindPrimaryID(trackID); diff --git a/source/digits_hits/src/GateToRoot.cc b/source/digits_hits/src/GateToRoot.cc index 15a7976aa..c13ea3813 100644 --- a/source/digits_hits/src/GateToRoot.cc +++ b/source/digits_hits/src/GateToRoot.cc @@ -70,24 +70,32 @@ ComptonRayleighData::ComptonRayleighData() { ; } ComptonRayleighData::ComptonRayleighData(ComptonRayleighData &aCRData) { photon1_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; - photon2_phantom_Rayleigh = aCRData.photon2_phantom_Rayleigh; - photon1_phantom_compton = aCRData.photon1_phantom_compton; - photon2_phantom_compton = aCRData.photon2_phantom_compton; - strcpy(theComptonVolumeName1, aCRData.theComptonVolumeName1); - strcpy(theComptonVolumeName2, aCRData.theComptonVolumeName2); - strcpy(theRayleighVolumeName1, aCRData.theRayleighVolumeName1); - strcpy(theRayleighVolumeName2, aCRData.theRayleighVolumeName2); + photon2_phantom_Rayleigh = aCRData.photon2_phantom_Rayleigh; + photon3_phantom_Rayleigh = aCRData.photon3_phantom_Rayleigh; + photon1_phantom_compton = aCRData.photon1_phantom_compton; + photon2_phantom_compton = aCRData.photon2_phantom_compton; + photon3_phantom_compton = aCRData.photon3_phantom_compton; + strcpy(theComptonVolumeName1,aCRData.theComptonVolumeName1); + strcpy(theComptonVolumeName2,aCRData.theComptonVolumeName2); + strcpy(theComptonVolumeName3,aCRData.theComptonVolumeName3); + strcpy(theRayleighVolumeName1,aCRData.theRayleighVolumeName1); + strcpy(theRayleighVolumeName2,aCRData.theRayleighVolumeName2); + strcpy(theRayleighVolumeName3,aCRData.theRayleighVolumeName3); } ComptonRayleighData &ComptonRayleighData::operator=(const ComptonRayleighData &aCR) { photon1_phantom_Rayleigh = aCR.photon1_phantom_Rayleigh; - photon2_phantom_Rayleigh = aCR.photon2_phantom_Rayleigh; - photon1_phantom_compton = aCR.photon1_phantom_compton; - photon2_phantom_compton = aCR.photon2_phantom_compton; - strcpy(theComptonVolumeName1, aCR.theComptonVolumeName1); - strcpy(theComptonVolumeName2, aCR.theComptonVolumeName2); - strcpy(theRayleighVolumeName1, aCR.theRayleighVolumeName1); - strcpy(theRayleighVolumeName2, aCR.theRayleighVolumeName2); + photon2_phantom_Rayleigh = aCR.photon2_phantom_Rayleigh; + photon3_phantom_Rayleigh = aCR.photon3_phantom_Rayleigh; + photon1_phantom_compton = aCR.photon1_phantom_compton; + photon2_phantom_compton = aCR.photon2_phantom_compton; + photon3_phantom_compton = aCR.photon3_phantom_compton; + strcpy(theComptonVolumeName1 , aCR.theComptonVolumeName1); + strcpy(theComptonVolumeName2 , aCR.theComptonVolumeName2); + strcpy(theComptonVolumeName3 , aCR.theComptonVolumeName3); + strcpy(theRayleighVolumeName1 , aCR.theRayleighVolumeName1); + strcpy(theRayleighVolumeName2 , aCR.theRayleighVolumeName2); + strcpy(theRayleighVolumeName3 , aCR.theRayleighVolumeName3); return *this; } @@ -442,14 +450,18 @@ void GateToRoot::RecordBeginOfAcquisition() { m_RecStepTree->Branch(G4String("dxg2").c_str(), &dxg2, "dxg2/D"); m_RecStepTree->Branch(G4String("dyg2").c_str(), &dyg2, "dyg2/D"); m_RecStepTree->Branch(G4String("dzg2").c_str(), &dzg2, "dzg2/D"); - m_RecStepTree->Branch(G4String("photon1PhR").c_str(), &theCRData.photon1_phantom_Rayleigh, "photon1PhR/I"); - m_RecStepTree->Branch(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh, "photon2PhR/I"); - m_RecStepTree->Branch(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton, "photon1PhC/I"); - m_RecStepTree->Branch(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton, "photon2PhC/I"); - m_RecStepTree->Branch(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1, "ComptVol1/C"); - m_RecStepTree->Branch(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2, "ComptVol2/C"); - m_RecStepTree->Branch(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1, "RaylVol1/C"); - m_RecStepTree->Branch(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2, "RaylVol2/C"); + m_RecStepTree->Branch(G4String("photon1PhR").c_str(), &theCRData.photon1_phantom_Rayleigh,"photon1PhR/I"); + m_RecStepTree->Branch(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh,"photon2PhR/I"); + m_RecStepTree->Branch(G4String("photon3PhR").c_str(), &theCRData.photon3_phantom_Rayleigh,"photon3PhR/I"); + m_RecStepTree->Branch(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton,"photon1PhC/I"); + m_RecStepTree->Branch(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton,"photon2PhC/I"); + m_RecStepTree->Branch(G4String("photon3PhC").c_str(), &theCRData.photon3_phantom_compton,"photon3PhC/I"); + m_RecStepTree->Branch(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1,"ComptVol1/C"); + m_RecStepTree->Branch(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2,"ComptVol2/C"); + m_RecStepTree->Branch(G4String("ComptVol3").c_str(), &theCRData.theComptonVolumeName3,"ComptVol3/C"); + m_RecStepTree->Branch(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1,"RaylVol1/C"); + m_RecStepTree->Branch(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2,"RaylVol2/C"); + m_RecStepTree->Branch(G4String("RaylVol3").c_str(), &theCRData.theRayleighVolumeName3,"RaylVol3/C"); m_RecStepTree->Branch(G4String("EventID").c_str(), &m_RSEventID, "EventID/I"); m_RecStepTree->Branch(G4String("RunID").c_str(), &m_RSRunID, "RunID/I"); //////////////////////// @@ -620,16 +632,17 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) { /*PY Descourt 08/09/2009 */ theCRData.photon1_phantom_Rayleigh = 0; - theCRData.photon2_phantom_Rayleigh = 0; - theCRData.photon1_phantom_compton = 0; - theCRData.photon2_phantom_compton = 0; - strcpy(theCRData.theComptonVolumeName1, G4String("NULL").c_str()); - strcpy(theCRData.theComptonVolumeName2, G4String("NULL").c_str()); - strcpy(theCRData.theRayleighVolumeName1, G4String("NULL").c_str()); - strcpy(theCRData.theRayleighVolumeName2, G4String("NULL").c_str()); - m_ionDecayPos = G4ThreeVector(0., 0., 0.); - m_positronGenerationPos = G4ThreeVector(0., 0., 0.); - m_positronAnnihilPos = G4ThreeVector(0., 0., 0.); + theCRData.photon2_phantom_Rayleigh = 0; + theCRData.photon3_phantom_Rayleigh = 0; + theCRData.photon1_phantom_compton = 0; + theCRData.photon2_phantom_compton = 0; + theCRData.photon3_phantom_compton = 0; + strcpy( theCRData.theComptonVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName3, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName3, G4String("NULL").c_str() ); dxg1 = 0.; dyg1 = 0.; dzg1 = 0.; @@ -638,13 +651,17 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) { dzg2 = 0.; theCRData.photon1_phantom_Rayleigh = 0; - theCRData.photon2_phantom_Rayleigh = 0; - theCRData.photon1_phantom_compton = 0; - theCRData.photon2_phantom_compton = 0; - strcpy(theCRData.theComptonVolumeName1, G4String("NULL").c_str()); - strcpy(theCRData.theComptonVolumeName2, G4String("NULL").c_str()); - strcpy(theCRData.theRayleighVolumeName1, G4String("NULL").c_str()); - strcpy(theCRData.theRayleighVolumeName2, G4String("NULL").c_str()); + theCRData.photon2_phantom_Rayleigh = 0; + theCRData.photon3_phantom_Rayleigh = 0; + theCRData.photon1_phantom_compton = 0; + theCRData.photon2_phantom_compton = 0; + theCRData.photon3_phantom_compton = 0; + strcpy( theCRData.theComptonVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName3, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName3, G4String("NULL").c_str() ); TrackingMode theMode = ((GateSteppingAction *) (GateRunManager::GetRunManager()->GetUserSteppingAction()))->GetMode(); @@ -1266,14 +1283,18 @@ void GateToRoot::PrintRecStep() { G4cout << "dxg2 = " << dxg2_copy << Gateendl; G4cout << "dyg2 = " << dyg2_copy << Gateendl; G4cout << "dzg2 = " << dzg2_copy << Gateendl; - G4cout << "photon1_phantom_Rayleigh = " << theCRData_copy.photon2_phantom_Rayleigh << Gateendl; - G4cout << "photon2_phantom_Rayleigh = " << theCRData_copy.photon2_phantom_Rayleigh << Gateendl; - G4cout << "photon1_phantom_compton = " << theCRData_copy.photon1_phantom_compton << Gateendl; - G4cout << "photon2_phantom_compton = " << theCRData_copy.photon2_phantom_compton << Gateendl; - G4cout << "theComptonVolumeName1 = " << theCRData_copy.theComptonVolumeName1 << Gateendl; - G4cout << "theComptonVolumeName2 = " << theCRData_copy.theComptonVolumeName2 << Gateendl; - G4cout << "theRayleighVolumeName1 = " << theCRData_copy.theRayleighVolumeName1 << Gateendl; - G4cout << "theRayleighVolumeName2 = " << theCRData_copy.theRayleighVolumeName2 << Gateendl; + G4cout << "photon1_phantom_Rayleigh = "<< theCRData_copy.photon1_phantom_Rayleigh<< Gateendl; + G4cout << "photon2_phantom_Rayleigh = "<< theCRData_copy.photon2_phantom_Rayleigh << Gateendl; + G4cout << "photon3_phantom_Rayleigh = "<< theCRData_copy.photon3_phantom_Rayleigh << Gateendl; + G4cout << "photon1_phantom_compton = "<< theCRData_copy.photon1_phantom_compton << Gateendl; + G4cout <<"photon2_phantom_compton = "<< theCRData_copy.photon2_phantom_compton << Gateendl; + G4cout <<"photon3_phantom_compton = "<< theCRData_copy.photon3_phantom_compton << Gateendl; + G4cout <<"theComptonVolumeName1 = "<< theCRData_copy.theComptonVolumeName1<< Gateendl; + G4cout <<"theComptonVolumeName2 = "<< theCRData_copy.theComptonVolumeName2<< Gateendl; + G4cout <<"theComptonVolumeName3 = "<< theCRData_copy.theComptonVolumeName3<< Gateendl; + G4cout <<"theRayleighVolumeName1 = "<< theCRData_copy.theRayleighVolumeName1<< Gateendl; + G4cout <<"theRayleighVolumeName2 = "<< theCRData_copy.theRayleighVolumeName2<< Gateendl; + G4cout <<"theRayleighVolumeName3 = "<< theCRData_copy.theRayleighVolumeName3<< Gateendl; } /// OPEN ROOT TRACKS DATA FILE IN READ MODE @@ -1378,21 +1399,19 @@ void GateToRoot::OpenTracksFile() { tracksTuple->SetBranchAddress(G4String("PPName").c_str(), &m_parentparticleName); tracksTuple->SetBranchAddress(G4String("LogAtVertex").c_str(), &m_volumeName); - if (currentN > 0) { + if (currentN > 0) + { // store the content of the RecStep data /// check that the first event ID entry is the same as the one of the preceding file // if not we should close this file and open it next time we get here tracksTuple->GetEntry(0); - if (EventID != lastEventID) { // if ( m_verboseLevel > 3 ) { - G4cout << " GateToRoot::OpenTracksFile() :: last Event ID was " << lastEventID - << " --- current one read from last Tracks Root File " << fTracksFN << " is " << EventID - << Gateendl; + if (EventID != lastEventID) + { // if ( m_verboseLevel > 3 ) { + G4cout << " GateToRoot::OpenTracksFile() :: last Event ID was " << lastEventID<< " --- current one read from last Tracks Root File " << fTracksFN << " is " << EventID<< Gateendl; G4cout << " SAVING RecStep Data to be used in GateToRoot::RecordBeginOfEvent \n"; - G4cout << " GateToRoot::OpenTracksFile() :: LAST event ID read from RecStep Data in file " << previousFN - << " is " << m_RSEventID << Gateendl; - G4cout << " GateToRoot::OpenTracksFile() :: LAST Event ID read from Tracks Data in previous file " - << previousFN << " is " << lastEventID << Gateendl; + G4cout << " GateToRoot::OpenTracksFile() :: LAST event ID read from RecStep Data in file " << previousFN<< " is " << m_RSEventID << Gateendl; + G4cout << " GateToRoot::OpenTracksFile() :: LAST Event ID read from Tracks Data in previous file "<SetBranchAddress(G4String("dyg2").c_str(), &dyg2); m_RecStepTree->SetBranchAddress(G4String("dzg2").c_str(), &dzg2); m_RecStepTree->SetBranchAddress(G4String("photon1PhR").c_str(), &theCRData.photon1_phantom_Rayleigh); - m_RecStepTree->SetBranchAddress(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh); - m_RecStepTree->SetBranchAddress(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton); - m_RecStepTree->SetBranchAddress(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton); - m_RecStepTree->SetBranchAddress(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1); - m_RecStepTree->SetBranchAddress(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2); - m_RecStepTree->SetBranchAddress(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1); - m_RecStepTree->SetBranchAddress(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2); + m_RecStepTree->SetBranchAddress(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh); + m_RecStepTree->SetBranchAddress(G4String("photon3PhR").c_str(), &theCRData.photon3_phantom_Rayleigh); + m_RecStepTree->SetBranchAddress(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton); + m_RecStepTree->SetBranchAddress(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton); + m_RecStepTree->SetBranchAddress(G4String("photon3PhC").c_str(), &theCRData.photon3_phantom_compton); + m_RecStepTree->SetBranchAddress(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1); + m_RecStepTree->SetBranchAddress(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2); + m_RecStepTree->SetBranchAddress(G4String("ComptVol3").c_str(), &theCRData.theComptonVolumeName3); + m_RecStepTree->SetBranchAddress(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1); + m_RecStepTree->SetBranchAddress(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2); + m_RecStepTree->SetBranchAddress(G4String("RaylVol3").c_str(), &theCRData.theRayleighVolumeName3); m_RecStepTree->SetBranchAddress(G4String("RunID").c_str(), &m_RSRunID); // rewind counters @@ -1431,24 +1454,32 @@ void GateToRoot::OpenTracksFile() { void GateToRoot::RecordPHData(ComptonRayleighData aCRData) { theCRData.photon1_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; - theCRData.photon2_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; - theCRData.photon1_phantom_compton = aCRData.photon1_phantom_compton; - theCRData.photon2_phantom_compton = aCRData.photon2_phantom_compton; - strcpy(theCRData.theComptonVolumeName1, aCRData.theComptonVolumeName1); - strcpy(theCRData.theComptonVolumeName2, aCRData.theComptonVolumeName2); - strcpy(theCRData.theRayleighVolumeName1, aCRData.theRayleighVolumeName1); - strcpy(theCRData.theRayleighVolumeName2, aCRData.theRayleighVolumeName2); + theCRData.photon2_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; + theCRData.photon3_phantom_Rayleigh = aCRData.photon3_phantom_Rayleigh; + theCRData.photon1_phantom_compton = aCRData.photon1_phantom_compton; + theCRData.photon2_phantom_compton = aCRData.photon2_phantom_compton; + theCRData.photon3_phantom_compton = aCRData.photon3_phantom_compton; + strcpy( theCRData.theComptonVolumeName1 , aCRData.theComptonVolumeName1 ); + strcpy( theCRData.theComptonVolumeName2 , aCRData.theComptonVolumeName2 ); + strcpy( theCRData.theComptonVolumeName3 , aCRData.theComptonVolumeName3 ); + strcpy( theCRData.theRayleighVolumeName1 , aCRData.theRayleighVolumeName1 ); + strcpy( theCRData.theRayleighVolumeName2 , aCRData.theRayleighVolumeName2 ); + strcpy( theCRData.theRayleighVolumeName3 , aCRData.theRayleighVolumeName3 ); } void GateToRoot::GetPHData(ComptonRayleighData &aCRData) { aCRData.photon1_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; - aCRData.photon2_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; - aCRData.photon1_phantom_compton = theCRData.photon1_phantom_compton; - aCRData.photon2_phantom_compton = theCRData.photon2_phantom_compton; - strcpy(aCRData.theComptonVolumeName1, theCRData.theComptonVolumeName1); - strcpy(aCRData.theComptonVolumeName2, theCRData.theComptonVolumeName2); - strcpy(aCRData.theRayleighVolumeName1, theCRData.theRayleighVolumeName1); - strcpy(aCRData.theRayleighVolumeName2, theCRData.theRayleighVolumeName2); + aCRData.photon2_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; + aCRData.photon3_phantom_Rayleigh = theCRData.photon3_phantom_Rayleigh; + aCRData.photon1_phantom_compton = theCRData.photon1_phantom_compton; + aCRData.photon2_phantom_compton = theCRData.photon2_phantom_compton; + aCRData.photon3_phantom_compton = theCRData.photon3_phantom_compton; + strcpy( aCRData.theComptonVolumeName1 , theCRData.theComptonVolumeName1 ); + strcpy( aCRData.theComptonVolumeName2 , theCRData.theComptonVolumeName2 ); + strcpy( aCRData.theComptonVolumeName3 , theCRData.theComptonVolumeName3 ); + strcpy( aCRData.theRayleighVolumeName1 , theCRData.theRayleighVolumeName1 ); + strcpy( aCRData.theRayleighVolumeName2 , theCRData.theRayleighVolumeName2 ); + strcpy( aCRData.theRayleighVolumeName3 , theCRData.theRayleighVolumeName3 ); } diff --git a/source/digits_hits/src/GateTrajectoryNavigator.cc b/source/digits_hits/src/GateTrajectoryNavigator.cc index 115899d58..a4eb35a29 100644 --- a/source/digits_hits/src/GateTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateTrajectoryNavigator.cc @@ -183,17 +183,24 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() // considered by GateAnalysis G4int i1; G4int i2; + G4int i3; + bool only2gamma = true; + for (G4int j1=0; j1= 0) && (i2 >= 0)) { + i3 = photonIndices[j3]; + if ((i1 >= 0) && (i2 >= 0) && (i3 >= 0)) { // both gammas were not already taken G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); + G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); + G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); // in detector mode the vertex position is stored at the last trajectory point // not the first one if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh @@ -202,6 +209,8 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); n_points =trj2->GetPointEntries(); vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj3->GetPointEntries(); + vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint( n_points - 1 )))->GetPosition(); } G4double dist = (vert1-vert2).mag(); @@ -215,14 +224,63 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() // we add both photons to the vertex m_photonIDVec.push_back(trj1->GetTrackID()); m_photonIDVec.push_back(trj2->GetTrackID()); + if ((vert2-vert3).mag()/mm < 1E-7) + { + m_photonIDVec.push_back(trj3->GetTrackID()); + photonIndices[j3] = -1; + only2gamma = false; + // std::cout << "\n\n\n\n\n\n\n\n\n" << photonIndices[j3] << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n" << std::endl; + } // we cancel the 2nd photon from the list to avoid double counting later photonIndices[j2] = -1; } } } + + } + } + + if(only2gamma){ + for (G4int j1=0; j1= 0) && (i2 >= 0)) { + // both gammas were not already taken + G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); + G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); + + G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); + G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); + // in detector mode the vertex position is stored at the last trajectory point + // not the first one + if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh + { + G4int n_points =trj1->GetPointEntries(); + vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj2->GetPointEntries(); + vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); + } + + G4double dist = (vert1-vert2).mag(); + if (nVerboseLevel > 2) + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << Gateendl; + if (dist/mm < 1E-7) { + if (nVerboseLevel > 1) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :" + << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << Gateendl; + } + // we add both photons to the vertex + m_photonIDVec.push_back(trj1->GetTrackID()); + m_photonIDVec.push_back(trj2->GetTrackID()); + // we cancel the 2nd photon from the list to avoid double counting later + photonIndices[j2] = -1; } } } +} + } } + } /*if ((m_positronTrackID != 0) && ((m_photonIDVec.size() > 2) || (m_photonIDVec.size() == 1))) { // strange case: we would expect to see a true annihilation, thus only 2 gammas backto back from the positron @@ -273,11 +331,12 @@ G4int GateTrajectoryNavigator::FindPhotonID(G4int trackID) G4int photon1ID = m_photonIDVec[0]; // in case there are 2 gammas OK, if not we put 0 (compatible with subsequent analysis) G4int photon2ID = (m_photonIDVec.size() >= 2) ? m_photonIDVec[1] : 0; + G4int photon3ID = (m_photonIDVec.size() >= 3) ? m_photonIDVec[2] : 0; G4int rootID = 0; // we go up and up, starting from the present trackID, to the parentID, the parentID, ecc until // we find that the ID of the track is equal to the ID of: one of the photons, or rootID(==0) - while (!((photonID==photon1ID)||(photonID==photon2ID)||(photonID==rootID))) { + while (!((photonID==photon1ID)||(photonID==photon2ID)||(photonID==photon3ID)||(photonID==rootID))) { found = false; G4int n_trajectories = m_trajectoryContainer->entries(); for (G4int iTrj=0; (iTrj Date: Thu, 12 Jan 2023 16:33:23 +0100 Subject: [PATCH 02/15] Small visual fixes like indentations --- source/digits_hits/include/GateToRoot.hh | 4 +- source/digits_hits/src/GateToRoot.cc | 161 +++++++++--------- .../src/GateTrajectoryNavigator.cc | 8 +- 3 files changed, 86 insertions(+), 87 deletions(-) diff --git a/source/digits_hits/include/GateToRoot.hh b/source/digits_hits/include/GateToRoot.hh index ab9fa8a52..ba237d43c 100644 --- a/source/digits_hits/include/GateToRoot.hh +++ b/source/digits_hits/include/GateToRoot.hh @@ -63,8 +63,8 @@ class GateTrajectoryNavigator; class ComptonRayleighData { public: G4int photon1_phantom_Rayleigh,photon2_phantom_Rayleigh,photon3_phantom_Rayleigh; - G4int photon1_phantom_compton , photon2_phantom_compton,photon3_phantom_compton; - Char_t theComptonVolumeName1[60] , theComptonVolumeName2[60], theComptonVolumeName3[60], theRayleighVolumeName1[60], theRayleighVolumeName2[60], theRayleighVolumeName3[60]; + G4int photon1_phantom_compton , photon2_phantom_compton,photon3_phantom_compton; + Char_t theComptonVolumeName1[60] , theComptonVolumeName2[60], theComptonVolumeName3[60], theRayleighVolumeName1[60], theRayleighVolumeName2[60], theRayleighVolumeName3[60]; ComptonRayleighData(ComptonRayleighData &); diff --git a/source/digits_hits/src/GateToRoot.cc b/source/digits_hits/src/GateToRoot.cc index c13ea3813..c67d8cff6 100644 --- a/source/digits_hits/src/GateToRoot.cc +++ b/source/digits_hits/src/GateToRoot.cc @@ -69,7 +69,7 @@ ComptonRayleighData::ComptonRayleighData() { ; } ComptonRayleighData::ComptonRayleighData(ComptonRayleighData &aCRData) { - photon1_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; + photon1_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; photon2_phantom_Rayleigh = aCRData.photon2_phantom_Rayleigh; photon3_phantom_Rayleigh = aCRData.photon3_phantom_Rayleigh; photon1_phantom_compton = aCRData.photon1_phantom_compton; @@ -84,7 +84,7 @@ ComptonRayleighData::ComptonRayleighData(ComptonRayleighData &aCRData) { } ComptonRayleighData &ComptonRayleighData::operator=(const ComptonRayleighData &aCR) { - photon1_phantom_Rayleigh = aCR.photon1_phantom_Rayleigh; + photon1_phantom_Rayleigh = aCR.photon1_phantom_Rayleigh; photon2_phantom_Rayleigh = aCR.photon2_phantom_Rayleigh; photon3_phantom_Rayleigh = aCR.photon3_phantom_Rayleigh; photon1_phantom_compton = aCR.photon1_phantom_compton; @@ -96,7 +96,7 @@ ComptonRayleighData &ComptonRayleighData::operator=(const ComptonRayleighData &a strcpy(theRayleighVolumeName1 , aCR.theRayleighVolumeName1); strcpy(theRayleighVolumeName2 , aCR.theRayleighVolumeName2); strcpy(theRayleighVolumeName3 , aCR.theRayleighVolumeName3); - return *this; + return *this; } //-------------------------------------------------------------------------- @@ -451,17 +451,17 @@ void GateToRoot::RecordBeginOfAcquisition() { m_RecStepTree->Branch(G4String("dyg2").c_str(), &dyg2, "dyg2/D"); m_RecStepTree->Branch(G4String("dzg2").c_str(), &dzg2, "dzg2/D"); m_RecStepTree->Branch(G4String("photon1PhR").c_str(), &theCRData.photon1_phantom_Rayleigh,"photon1PhR/I"); - m_RecStepTree->Branch(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh,"photon2PhR/I"); - m_RecStepTree->Branch(G4String("photon3PhR").c_str(), &theCRData.photon3_phantom_Rayleigh,"photon3PhR/I"); - m_RecStepTree->Branch(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton,"photon1PhC/I"); - m_RecStepTree->Branch(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton,"photon2PhC/I"); - m_RecStepTree->Branch(G4String("photon3PhC").c_str(), &theCRData.photon3_phantom_compton,"photon3PhC/I"); - m_RecStepTree->Branch(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1,"ComptVol1/C"); - m_RecStepTree->Branch(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2,"ComptVol2/C"); - m_RecStepTree->Branch(G4String("ComptVol3").c_str(), &theCRData.theComptonVolumeName3,"ComptVol3/C"); - m_RecStepTree->Branch(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1,"RaylVol1/C"); - m_RecStepTree->Branch(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2,"RaylVol2/C"); - m_RecStepTree->Branch(G4String("RaylVol3").c_str(), &theCRData.theRayleighVolumeName3,"RaylVol3/C"); + m_RecStepTree->Branch(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh,"photon2PhR/I"); + m_RecStepTree->Branch(G4String("photon3PhR").c_str(), &theCRData.photon3_phantom_Rayleigh,"photon3PhR/I"); + m_RecStepTree->Branch(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton,"photon1PhC/I"); + m_RecStepTree->Branch(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton,"photon2PhC/I"); + m_RecStepTree->Branch(G4String("photon3PhC").c_str(), &theCRData.photon3_phantom_compton,"photon3PhC/I"); + m_RecStepTree->Branch(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1,"ComptVol1/C"); + m_RecStepTree->Branch(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2,"ComptVol2/C"); + m_RecStepTree->Branch(G4String("ComptVol3").c_str(), &theCRData.theComptonVolumeName3,"ComptVol3/C"); + m_RecStepTree->Branch(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1,"RaylVol1/C"); + m_RecStepTree->Branch(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2,"RaylVol2/C"); + m_RecStepTree->Branch(G4String("RaylVol3").c_str(), &theCRData.theRayleighVolumeName3,"RaylVol3/C"); m_RecStepTree->Branch(G4String("EventID").c_str(), &m_RSEventID, "EventID/I"); m_RecStepTree->Branch(G4String("RunID").c_str(), &m_RSRunID, "RunID/I"); //////////////////////// @@ -632,17 +632,17 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) { /*PY Descourt 08/09/2009 */ theCRData.photon1_phantom_Rayleigh = 0; - theCRData.photon2_phantom_Rayleigh = 0; - theCRData.photon3_phantom_Rayleigh = 0; - theCRData.photon1_phantom_compton = 0; - theCRData.photon2_phantom_compton = 0; - theCRData.photon3_phantom_compton = 0; - strcpy( theCRData.theComptonVolumeName1, G4String("NULL").c_str() ); - strcpy( theCRData.theComptonVolumeName2, G4String("NULL").c_str() ); - strcpy( theCRData.theComptonVolumeName3, G4String("NULL").c_str() ); - strcpy( theCRData.theRayleighVolumeName1, G4String("NULL").c_str() ); - strcpy( theCRData.theRayleighVolumeName2, G4String("NULL").c_str() ); - strcpy( theCRData.theRayleighVolumeName3, G4String("NULL").c_str() ); + theCRData.photon2_phantom_Rayleigh = 0; + theCRData.photon3_phantom_Rayleigh = 0; + theCRData.photon1_phantom_compton = 0; + theCRData.photon2_phantom_compton = 0; + theCRData.photon3_phantom_compton = 0; + strcpy( theCRData.theComptonVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName3, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName3, G4String("NULL").c_str() ); dxg1 = 0.; dyg1 = 0.; dzg1 = 0.; @@ -651,17 +651,17 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) { dzg2 = 0.; theCRData.photon1_phantom_Rayleigh = 0; - theCRData.photon2_phantom_Rayleigh = 0; - theCRData.photon3_phantom_Rayleigh = 0; - theCRData.photon1_phantom_compton = 0; - theCRData.photon2_phantom_compton = 0; - theCRData.photon3_phantom_compton = 0; - strcpy( theCRData.theComptonVolumeName1, G4String("NULL").c_str() ); - strcpy( theCRData.theComptonVolumeName2, G4String("NULL").c_str() ); - strcpy( theCRData.theComptonVolumeName3, G4String("NULL").c_str() ); - strcpy( theCRData.theRayleighVolumeName1, G4String("NULL").c_str() ); - strcpy( theCRData.theRayleighVolumeName2, G4String("NULL").c_str() ); - strcpy( theCRData.theRayleighVolumeName3, G4String("NULL").c_str() ); + theCRData.photon2_phantom_Rayleigh = 0; + theCRData.photon3_phantom_Rayleigh = 0; + theCRData.photon1_phantom_compton = 0; + theCRData.photon2_phantom_compton = 0; + theCRData.photon3_phantom_compton = 0; + strcpy( theCRData.theComptonVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theComptonVolumeName3, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName1, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName2, G4String("NULL").c_str() ); + strcpy( theCRData.theRayleighVolumeName3, G4String("NULL").c_str() ); TrackingMode theMode = ((GateSteppingAction *) (GateRunManager::GetRunManager()->GetUserSteppingAction()))->GetMode(); @@ -1284,17 +1284,17 @@ void GateToRoot::PrintRecStep() { G4cout << "dyg2 = " << dyg2_copy << Gateendl; G4cout << "dzg2 = " << dzg2_copy << Gateendl; G4cout << "photon1_phantom_Rayleigh = "<< theCRData_copy.photon1_phantom_Rayleigh<< Gateendl; - G4cout << "photon2_phantom_Rayleigh = "<< theCRData_copy.photon2_phantom_Rayleigh << Gateendl; - G4cout << "photon3_phantom_Rayleigh = "<< theCRData_copy.photon3_phantom_Rayleigh << Gateendl; - G4cout << "photon1_phantom_compton = "<< theCRData_copy.photon1_phantom_compton << Gateendl; - G4cout <<"photon2_phantom_compton = "<< theCRData_copy.photon2_phantom_compton << Gateendl; - G4cout <<"photon3_phantom_compton = "<< theCRData_copy.photon3_phantom_compton << Gateendl; - G4cout <<"theComptonVolumeName1 = "<< theCRData_copy.theComptonVolumeName1<< Gateendl; - G4cout <<"theComptonVolumeName2 = "<< theCRData_copy.theComptonVolumeName2<< Gateendl; - G4cout <<"theComptonVolumeName3 = "<< theCRData_copy.theComptonVolumeName3<< Gateendl; - G4cout <<"theRayleighVolumeName1 = "<< theCRData_copy.theRayleighVolumeName1<< Gateendl; - G4cout <<"theRayleighVolumeName2 = "<< theCRData_copy.theRayleighVolumeName2<< Gateendl; - G4cout <<"theRayleighVolumeName3 = "<< theCRData_copy.theRayleighVolumeName3<< Gateendl; + G4cout << "photon2_phantom_Rayleigh = "<< theCRData_copy.photon2_phantom_Rayleigh << Gateendl; + G4cout << "photon3_phantom_Rayleigh = "<< theCRData_copy.photon3_phantom_Rayleigh << Gateendl; + G4cout << "photon1_phantom_compton = "<< theCRData_copy.photon1_phantom_compton << Gateendl; + G4cout <<"photon2_phantom_compton = "<< theCRData_copy.photon2_phantom_compton << Gateendl; + G4cout <<"photon3_phantom_compton = "<< theCRData_copy.photon3_phantom_compton << Gateendl; + G4cout <<"theComptonVolumeName1 = "<< theCRData_copy.theComptonVolumeName1<< Gateendl; + G4cout <<"theComptonVolumeName2 = "<< theCRData_copy.theComptonVolumeName2<< Gateendl; + G4cout <<"theComptonVolumeName3 = "<< theCRData_copy.theComptonVolumeName3<< Gateendl; + G4cout <<"theRayleighVolumeName1 = "<< theCRData_copy.theRayleighVolumeName1<< Gateendl; + G4cout <<"theRayleighVolumeName2 = "<< theCRData_copy.theRayleighVolumeName2<< Gateendl; + G4cout <<"theRayleighVolumeName3 = "<< theCRData_copy.theRayleighVolumeName3<< Gateendl; } /// OPEN ROOT TRACKS DATA FILE IN READ MODE @@ -1430,17 +1430,17 @@ void GateToRoot::OpenTracksFile() { m_RecStepTree->SetBranchAddress(G4String("dyg2").c_str(), &dyg2); m_RecStepTree->SetBranchAddress(G4String("dzg2").c_str(), &dzg2); m_RecStepTree->SetBranchAddress(G4String("photon1PhR").c_str(), &theCRData.photon1_phantom_Rayleigh); - m_RecStepTree->SetBranchAddress(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh); - m_RecStepTree->SetBranchAddress(G4String("photon3PhR").c_str(), &theCRData.photon3_phantom_Rayleigh); - m_RecStepTree->SetBranchAddress(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton); - m_RecStepTree->SetBranchAddress(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton); - m_RecStepTree->SetBranchAddress(G4String("photon3PhC").c_str(), &theCRData.photon3_phantom_compton); - m_RecStepTree->SetBranchAddress(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1); - m_RecStepTree->SetBranchAddress(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2); - m_RecStepTree->SetBranchAddress(G4String("ComptVol3").c_str(), &theCRData.theComptonVolumeName3); - m_RecStepTree->SetBranchAddress(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1); - m_RecStepTree->SetBranchAddress(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2); - m_RecStepTree->SetBranchAddress(G4String("RaylVol3").c_str(), &theCRData.theRayleighVolumeName3); + m_RecStepTree->SetBranchAddress(G4String("photon2PhR").c_str(), &theCRData.photon2_phantom_Rayleigh); + m_RecStepTree->SetBranchAddress(G4String("photon3PhR").c_str(), &theCRData.photon3_phantom_Rayleigh); + m_RecStepTree->SetBranchAddress(G4String("photon1PhC").c_str(), &theCRData.photon1_phantom_compton); + m_RecStepTree->SetBranchAddress(G4String("photon2PhC").c_str(), &theCRData.photon2_phantom_compton); + m_RecStepTree->SetBranchAddress(G4String("photon3PhC").c_str(), &theCRData.photon3_phantom_compton); + m_RecStepTree->SetBranchAddress(G4String("ComptVol1").c_str(), &theCRData.theComptonVolumeName1); + m_RecStepTree->SetBranchAddress(G4String("ComptVol2").c_str(), &theCRData.theComptonVolumeName2); + m_RecStepTree->SetBranchAddress(G4String("ComptVol3").c_str(), &theCRData.theComptonVolumeName3); + m_RecStepTree->SetBranchAddress(G4String("RaylVol1").c_str(), &theCRData.theRayleighVolumeName1); + m_RecStepTree->SetBranchAddress(G4String("RaylVol2").c_str(), &theCRData.theRayleighVolumeName2); + m_RecStepTree->SetBranchAddress(G4String("RaylVol3").c_str(), &theCRData.theRayleighVolumeName3); m_RecStepTree->SetBranchAddress(G4String("RunID").c_str(), &m_RSRunID); // rewind counters @@ -1454,33 +1454,32 @@ void GateToRoot::OpenTracksFile() { void GateToRoot::RecordPHData(ComptonRayleighData aCRData) { theCRData.photon1_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; - theCRData.photon2_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; - theCRData.photon3_phantom_Rayleigh = aCRData.photon3_phantom_Rayleigh; - theCRData.photon1_phantom_compton = aCRData.photon1_phantom_compton; - theCRData.photon2_phantom_compton = aCRData.photon2_phantom_compton; - theCRData.photon3_phantom_compton = aCRData.photon3_phantom_compton; - strcpy( theCRData.theComptonVolumeName1 , aCRData.theComptonVolumeName1 ); - strcpy( theCRData.theComptonVolumeName2 , aCRData.theComptonVolumeName2 ); - strcpy( theCRData.theComptonVolumeName3 , aCRData.theComptonVolumeName3 ); - strcpy( theCRData.theRayleighVolumeName1 , aCRData.theRayleighVolumeName1 ); - strcpy( theCRData.theRayleighVolumeName2 , aCRData.theRayleighVolumeName2 ); - strcpy( theCRData.theRayleighVolumeName3 , aCRData.theRayleighVolumeName3 ); + theCRData.photon2_phantom_Rayleigh = aCRData.photon1_phantom_Rayleigh; + theCRData.photon3_phantom_Rayleigh = aCRData.photon3_phantom_Rayleigh; + theCRData.photon1_phantom_compton = aCRData.photon1_phantom_compton; + theCRData.photon2_phantom_compton = aCRData.photon2_phantom_compton; + theCRData.photon3_phantom_compton = aCRData.photon3_phantom_compton; + strcpy( theCRData.theComptonVolumeName1 , aCRData.theComptonVolumeName1 ); + strcpy( theCRData.theComptonVolumeName2 , aCRData.theComptonVolumeName2 ); + strcpy( theCRData.theComptonVolumeName3 , aCRData.theComptonVolumeName3 ); + strcpy( theCRData.theRayleighVolumeName1 , aCRData.theRayleighVolumeName1 ); + strcpy( theCRData.theRayleighVolumeName2 , aCRData.theRayleighVolumeName2 ); + strcpy( theCRData.theRayleighVolumeName3 , aCRData.theRayleighVolumeName3 ); } void GateToRoot::GetPHData(ComptonRayleighData &aCRData) { aCRData.photon1_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; - aCRData.photon2_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; - aCRData.photon3_phantom_Rayleigh = theCRData.photon3_phantom_Rayleigh; - aCRData.photon1_phantom_compton = theCRData.photon1_phantom_compton; - aCRData.photon2_phantom_compton = theCRData.photon2_phantom_compton; - aCRData.photon3_phantom_compton = theCRData.photon3_phantom_compton; - strcpy( aCRData.theComptonVolumeName1 , theCRData.theComptonVolumeName1 ); - strcpy( aCRData.theComptonVolumeName2 , theCRData.theComptonVolumeName2 ); - strcpy( aCRData.theComptonVolumeName3 , theCRData.theComptonVolumeName3 ); - strcpy( aCRData.theRayleighVolumeName1 , theCRData.theRayleighVolumeName1 ); - strcpy( aCRData.theRayleighVolumeName2 , theCRData.theRayleighVolumeName2 ); - strcpy( aCRData.theRayleighVolumeName3 , theCRData.theRayleighVolumeName3 ); - + aCRData.photon2_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; + aCRData.photon3_phantom_Rayleigh = theCRData.photon3_phantom_Rayleigh; + aCRData.photon1_phantom_compton = theCRData.photon1_phantom_compton; + aCRData.photon2_phantom_compton = theCRData.photon2_phantom_compton; + aCRData.photon3_phantom_compton = theCRData.photon3_phantom_compton; + strcpy( aCRData.theComptonVolumeName1 , theCRData.theComptonVolumeName1 ); + strcpy( aCRData.theComptonVolumeName2 , theCRData.theComptonVolumeName2 ); + strcpy( aCRData.theComptonVolumeName3 , theCRData.theComptonVolumeName3 ); + strcpy( aCRData.theRayleighVolumeName1 , theCRData.theRayleighVolumeName1 ); + strcpy( aCRData.theRayleighVolumeName2 , theCRData.theRayleighVolumeName2 ); + strcpy( aCRData.theRayleighVolumeName3 , theCRData.theRayleighVolumeName3 ); } void GateToRoot::GetCurrentRecStepData(const G4Event *evt) { diff --git a/source/digits_hits/src/GateTrajectoryNavigator.cc b/source/digits_hits/src/GateTrajectoryNavigator.cc index a4eb35a29..4ecdaf4ad 100644 --- a/source/digits_hits/src/GateTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateTrajectoryNavigator.cc @@ -191,16 +191,16 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() for(G4int j3=j2+1; j3= 0) && (i2 >= 0) && (i3 >= 0)) { // both gammas were not already taken G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); - G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); + G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); - G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); + G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); // in detector mode the vertex position is stored at the last trajectory point // not the first one if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh @@ -267,7 +267,7 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << Gateendl; if (dist/mm < 1E-7) { if (nVerboseLevel > 1) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :" + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :" << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << Gateendl; } // we add both photons to the vertex From 04cdaaeb655eaab36def719a4ad1408968c09ae3 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 19:33:41 +0100 Subject: [PATCH 03/15] style:code is more readable without commented out code --- source/digits_hits/src/GateAnalysis.cc | 32 ++------------------------ 1 file changed, 2 insertions(+), 30 deletions(-) diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index c04e8e00a..911e65b14 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -124,9 +124,6 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4int eventID = event->GetEventID(); G4int runID = GateRunManager::GetRunManager()->GetCurrentRun()->GetRunID(); - //G4cout << "GateAnalysis::EventID et RunID : " <entries(); - //G4cout << " " << NbHits << " hits are stored in essaiCrystalHitsCollection.\n"; - //G4int ionID = 1; // the primary vertex particle - //G4int positronID = 0; // no more needed G4int photon1ID = 0; G4int photon2ID = 0; G4int photon3ID = 0; @@ -172,20 +166,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4int septalNb = 0; // HDS : septal penetration - //////////// - // search the positron - //positronID = // No more needed m_trajectoryNavigator->FindPositronTrackID(); - - /*if (positronID == 0) - { - if (nVerboseLevel > 0) G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : positronID == 0\n"; - } - - if (nVerboseLevel > 1) G4cout << "GateAnalysis::RecordEndOfEvent : positronID : " << positronID << Gateendl; - */ - - //////////// //search the two gammas std::vector photonIDVec = m_trajectoryNavigator->FindAnnihilationGammasTrackID(); @@ -348,7 +329,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if ( theMode == TrackingMode::kTracker ) // in tracker mode we store the infos about the number of compton and rayleigh - { // G4cout << " GateAnalysis eventID "<GetModule("root")); ComptonRayleighData aCRData; aCRData.photon1_phantom_Rayleigh = photon1_phantom_Rayleigh; @@ -364,7 +345,6 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) strcpy(aCRData.theRayleighVolumeName2 , theRayleighVolumeName2.c_str() ); strcpy(aCRData.theRayleighVolumeName3 , theRayleighVolumeName3.c_str() ); gateToRoot->RecordPHData( aCRData ); - // return; } if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh @@ -380,12 +360,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) photon1_phantom_compton += aCRData.photon1_phantom_compton; photon2_phantom_compton += aCRData.photon2_phantom_compton; photon3_phantom_compton += aCRData.photon3_phantom_compton; - /* - if( theComptonVolumeName1 == G4String("NULL") ) {theComptonVolumeName1 = aCRData.theComptonVolumeName1;} - if( theComptonVolumeName2 == G4String("NULL") ) {theComptonVolumeName2 = aCRData.theComptonVolumeName2;} - if( theRayleighVolumeName1 == G4String("NULL") ) {theRayleighVolumeName1 = aCRData.theRayleighVolumeName1;} - if( theRayleighVolumeName2 == G4String("NULL") ){theRayleighVolumeName2 = aCRData.theRayleighVolumeName2;} - */ + theComptonVolumeName1 = aCRData.theComptonVolumeName1; theComptonVolumeName2 = aCRData.theComptonVolumeName2; theComptonVolumeName3 = aCRData.theComptonVolumeName3; @@ -411,7 +386,6 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4int crystalTrackID = (*CHC)[iHit]->GetTrackID(); G4String processName = (*CHC)[iHit]->GetProcess(); // Counting Compton in the Crystal - // if (processName.find("ompt") != G4String::npos || processName.find("Rayleigh") != G4String::npos) { if (processName.find("ompt") != G4String::npos) { @@ -438,7 +412,6 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) // fill in values with the branch with C struct G4int trackID = (*CHC)[iHit]->GetTrackID(); - //G4int parentID = (*CHC)[iHit]->GetParentID(); G4int photonID = 0; G4int nPhantomCompton = 0; @@ -447,7 +420,6 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4int nPhantomRayleigh = 0; G4int nCrystalRayleigh = 0; - // if ((photon1ID != 0) && (photon2ID != 0)) { if (photon1ID != 0) { // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT // search the gamma from which this hit comes --> photonID From 8da55dbafa6296ae5751930eb27be890dd5a87f9 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 20:07:15 +0100 Subject: [PATCH 04/15] style:code is more readable with corrected tabulation and added brackets --- source/digits_hits/src/GateAnalysis.cc | 645 ++++++++++++------------- 1 file changed, 299 insertions(+), 346 deletions(-) diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 911e65b14..09d36f983 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -115,370 +115,323 @@ void GateAnalysis::RecordBeginOfEvent(const G4Event* ) void GateAnalysis::RecordEndOfEvent(const G4Event* event) { if (nVerboseLevel > 2) - G4cout << "GateAnalysis::RecordEndOfEvent "<< Gateendl; + G4cout << "GateAnalysis::RecordEndOfEvent" << G4endl; G4TrajectoryContainer* trajectoryContainer = event->GetTrajectoryContainer(); - if (trajectoryContainer) + if (trajectoryContainer) { m_trajectoryNavigator->SetTrajectoryContainer(trajectoryContainer); + } G4int eventID = event->GetEventID(); G4int runID = GateRunManager::GetRunManager()->GetCurrentRun()->GetRunID(); - if (!trajectoryContainer) - { - if (nVerboseLevel > 0) - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : G4TrajectoryContainer not found\n"; + if (!trajectoryContainer) { + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : G4TrajectoryContainer not found" << G4endl; } - else - { - GateCrystalHitsCollection* CHC = GetOutputMgr()->GetCrystalHitCollection(); - - G4int NbHits = 0; - G4int NpHits = 0; - - if (CHC) - { - NbHits = CHC->entries(); - - G4int photon1ID = 0; - G4int photon2ID = 0; - G4int photon3ID = 0; - G4int rootID = 0; - G4int primaryID = 0; - - G4int photon1_phantom_compton = 0; - G4int photon2_phantom_compton = 0; - G4int photon3_phantom_compton = 0; - - - G4int photon1_crystal_compton = 0; - G4int photon2_crystal_compton = 0; - G4int photon3_crystal_compton = 0; - - G4int photon1_phantom_Rayleigh = 0; - G4int photon2_phantom_Rayleigh = 0; - G4int photon3_phantom_Rayleigh = 0; - - G4int photon1_crystal_Rayleigh = 0; - G4int photon2_crystal_Rayleigh = 0; - G4int photon3_crystal_Rayleigh = 0; - - G4int septalNb = 0; // HDS : septal penetration - - m_trajectoryNavigator->FindPositronTrackID(); - //search the two gammas - - std::vector photonIDVec = m_trajectoryNavigator->FindAnnihilationGammasTrackID(); - if (photonIDVec.size() == 0) - { - // no gamma coming from a positron or an ion, or shooted as primary - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : WARNING : photonIDs not found\n"; + } else { + GateCrystalHitsCollection* CHC = GetOutputMgr()->GetCrystalHitCollection(); + G4int NbHits = 0; + G4int NpHits = 0; + + if (CHC) { + NbHits = CHC->entries(); + G4int photon1ID = 0; + G4int photon2ID = 0; + G4int photon3ID = 0; + G4int rootID = 0; + G4int primaryID = 0; + G4int photon1_phantom_compton = 0; + G4int photon2_phantom_compton = 0; + G4int photon3_phantom_compton = 0; + G4int photon1_crystal_compton = 0; + G4int photon2_crystal_compton = 0; + G4int photon3_crystal_compton = 0; + G4int photon1_phantom_Rayleigh = 0; + G4int photon2_phantom_Rayleigh = 0; + G4int photon3_phantom_Rayleigh = 0; + G4int photon1_crystal_Rayleigh = 0; + G4int photon2_crystal_Rayleigh = 0; + G4int photon3_crystal_Rayleigh = 0; + G4int septalNb = 0; // HDS : septal penetration + + m_trajectoryNavigator->FindPositronTrackID(); + //search the two gammas + std::vector photonIDVec = m_trajectoryNavigator->FindAnnihilationGammasTrackID(); + if (photonIDVec.size() == 0) { + // no gamma coming from a positron or an ion, or shooted as primary + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonIDs not found" << G4endl; + } + } else { + // This warning is somewhat irrelevant with 124I + if (nVerboseLevel > 0 && photonIDVec.size() > 2) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonID vector size > 2" << G4endl; + } + photon1ID = photonIDVec[0]; + photon2ID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; + photon3ID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; + } + + if (nVerboseLevel > 0) { + if (photon1ID == 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon1ID == 0" << G4endl; + } + if (photon2ID == 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0" << G4endl; + } + if (photon3ID == 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon3ID == 0" << G4endl; + } + if (nVerboseLevel > 1) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon1ID; + G4cout << " photon2ID : " << photon2ID << G4endl; + } + } + + // analysis of the phantom hits to count the comptons, etc. + + GatePhantomHitsCollection* PHC = GetOutputMgr()->GetPhantomHitCollection(); + NpHits = PHC->entries(); + G4String theComptonVolumeName("NULL"); + G4String theComptonVolumeName1("NULL"); + G4String theComptonVolumeName2("NULL"); + G4String theComptonVolumeName3("NULL"); + G4String theRayleighVolumeName("NULL"); + G4String theRayleighVolumeName1("NULL"); + G4String theRayleighVolumeName2("NULL"); + G4String theRayleighVolumeName3("NULL"); + + for (G4int iPHit=0;iPHitGetPhysVolName() == m_septalPhysVolumeName) { + ++septalNb; + } + } + G4int phantomTrackID = (*PHC)[iPHit]->GetTrackID(); + G4String processName = (*PHC)[iPHit]->GetProcess(); + G4int PDGcode = (*PHC)[iPHit]->GetPDGEncoding(); + G4ThreeVector hitPos = (*PHC)[iPHit]->GetPos(); + theComptonVolumeName = G4String("NULL"); + theRayleighVolumeName = G4String("NULL"); + + if (nVerboseLevel > 2) { + G4cout << "GateAnalysis::RecordEndOfEvent : GatePhantomHitsCollection : trackID : " << std::setw(5) << phantomTrackID; + G4cout << " PDG code : " << std::setw(5) << PDGcode << " processName : <" << processName << G4endl; + } + + // Modif by DS and LS on Oct 4, 2002: we need to be able to recognise both 'compt' + // and 'LowEnCompt", hence the find on 'ompt' + // if (processName.find("ompt") != G4String::npos || processName.find("Rayleigh") != G4String::npos) { + // modif. by CJG to separate Compton and Rayleigh photons + if (processName.find("ompt") != G4String::npos) { + if ((phantomTrackID == photon1ID)||(phantomTrackID == photon2ID)) { + G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4ThreeVector null(0.,0.,0.); + G4ThreeVector *ptr = &null; + theComptonVolumeName = gNavigator->LocateGlobalPointAndSetup(hitPos,ptr,false)->GetName(); + if (nVerboseLevel > 1) { + G4cout << "GateAnalysis::RecordEndOfEvent : theComptonVolumeName: " << theComptonVolumeName << G4endl; } - else - { - // This warning is somewhat irrelevant with 124I - if (nVerboseLevel > 0 && photonIDVec.size() > 2) - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonID vector size > 2\n"; - - photon1ID = photonIDVec[0]; - photon2ID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; - photon3ID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; + } + if (phantomTrackID == photon1ID) { + photon1_phantom_compton++; + theComptonVolumeName1 = theComptonVolumeName; + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon1_phantom_compton : " << photon1_phantom_compton << G4endl; } - - if (photon1ID == 0) - { - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : WARNING : photon1ID == 0\n"; + } + if (phantomTrackID == photon2ID) { + photon2_phantom_compton++; + theComptonVolumeName2 = theComptonVolumeName; + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon2_phantom_compton : " << photon2_phantom_compton << G4endl; } - if (photon2ID == 0) { - if (nVerboseLevel > 1) G4cout - << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0\n"; } - if (photon3ID == 0) { - if (nVerboseLevel > 1) G4cout - << "GateAnalysis::RecordEndOfEvent : WARNING : photon3ID == 0\n"; + if (phantomTrackID == photon3ID) { + photon3_phantom_compton++; + theComptonVolumeName3 = theComptonVolumeName; + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon3_phantom_compton : " << photon3_phantom_compton << G4endl; + } } - - if (nVerboseLevel > 1) G4cout - << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon1ID - << " photon2ID : " << photon2ID << Gateendl; - - - // analysis of the phantom hits to count the comptons, etc. - - GatePhantomHitsCollection* PHC = GetOutputMgr()->GetPhantomHitCollection(); - NpHits = PHC->entries(); - - G4String theComptonVolumeName("NULL"); - G4String theComptonVolumeName1("NULL"); - G4String theComptonVolumeName2("NULL"); - G4String theComptonVolumeName3("NULL"); - - G4String theRayleighVolumeName("NULL"); - G4String theRayleighVolumeName1("NULL"); - G4String theRayleighVolumeName2("NULL"); - G4String theRayleighVolumeName3("NULL"); - - for (G4int iPHit=0;iPHitGetPhysVolName() == m_septalPhysVolumeName) { - ++septalNb; - } - } - // - - G4int phantomTrackID = (*PHC)[iPHit]->GetTrackID(); - G4String processName = (*PHC)[iPHit]->GetProcess(); - G4int PDGcode = (*PHC)[iPHit]->GetPDGEncoding(); - G4ThreeVector hitPos = (*PHC)[iPHit]->GetPos(); - - if (nVerboseLevel > 2) - G4cout << "GateAnalysis::RecordEndOfEvent : GatePhantomHitsCollection : trackID : " << std::setw(5) << phantomTrackID - << " PDG code : " << std::setw(5) << PDGcode << " processName : <" << processName << ">\n"; - theComptonVolumeName = G4String("NULL"); - - if (nVerboseLevel > 2) G4cout - << "GateAnalysis::RecordEndOfEvent : GatePhantomHitsCollection : trackID : " << std::setw(5) << phantomTrackID - << " PDG code : " << std::setw(5) << PDGcode << " processName : <" << processName << ">\n"; - theRayleighVolumeName = G4String("NULL"); - - // Modif by DS and LS on Oct 4, 2002: we need to be able to recognise both 'compt' - // and 'LowEnCompt", hence the find on 'ompt' - // if (processName.find("ompt") != G4String::npos || processName.find("Rayleigh") != G4String::npos) { - // modif. by CJG to separate Compton and Rayleigh photons - if (processName.find("ompt") != G4String::npos) - { - if ((phantomTrackID == photon1ID)||(phantomTrackID == photon2ID)) - { - G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4ThreeVector null(0.,0.,0.); - G4ThreeVector *ptr; - ptr = &null; - theComptonVolumeName = gNavigator->LocateGlobalPointAndSetup(hitPos,ptr,false)->GetName(); - if (nVerboseLevel > 1) - G4cout << "GateAnalysis::RecordEndOfEvent : theComptonVolumeName: " - << theComptonVolumeName << Gateendl; - } - if (phantomTrackID == photon1ID) - { - photon1_phantom_compton++; - theComptonVolumeName1 = theComptonVolumeName; - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : photon1_phantom_compton : " << photon1_phantom_compton << Gateendl; - } - if (phantomTrackID == photon2ID) - { - photon2_phantom_compton++; - theComptonVolumeName2 = theComptonVolumeName; - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : photon2_phantom_compton : " << photon2_phantom_compton << Gateendl; - } - if (phantomTrackID == photon3ID) - { - photon3_phantom_compton++; - theComptonVolumeName3 = theComptonVolumeName; - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : photon3_phantom_compton : " << photon3_phantom_compton << Gateendl; - } - } - - // Counting Rayleigh scatter in phantom - if (processName.find("Rayl") != G4String::npos) - { - if ((phantomTrackID == photon1ID)||(phantomTrackID == photon2ID)) - { - G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4ThreeVector null(0.,0.,0.); - G4ThreeVector *ptr; - ptr = &null; - theRayleighVolumeName = gNavigator->LocateGlobalPointAndSetup(hitPos,ptr,false)->GetName(); - if (nVerboseLevel > 1) - G4cout << "GateAnalysis::RecordEndOfEvent : theRayleighVolumeName: " - << theRayleighVolumeName << Gateendl; - } - if (phantomTrackID == photon1ID) - { - photon1_phantom_Rayleigh++; - theRayleighVolumeName1 = theRayleighVolumeName; - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : photon1_phantom_Rayleigh : " << photon1_phantom_Rayleigh << Gateendl; - } - if (phantomTrackID == photon2ID) - { - photon2_phantom_Rayleigh++; - theRayleighVolumeName2 = theRayleighVolumeName; - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : photon2_phantom_Rayleigh : " << photon2_phantom_Rayleigh << Gateendl; - } - if (phantomTrackID == photon3ID) - { - photon3_phantom_Rayleigh++; - theRayleighVolumeName3 = theRayleighVolumeName; - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : photon3_phantom_Rayleigh : " << photon3_phantom_Rayleigh << Gateendl; - } - } - } // end loop NpHits - - TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); - - - if ( theMode == TrackingMode::kTracker ) // in tracker mode we store the infos about the number of compton and rayleigh - { - GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); - ComptonRayleighData aCRData; - aCRData.photon1_phantom_Rayleigh = photon1_phantom_Rayleigh; - aCRData.photon2_phantom_Rayleigh = photon2_phantom_Rayleigh; - aCRData.photon3_phantom_Rayleigh = photon3_phantom_Rayleigh; - aCRData.photon1_phantom_compton = photon1_phantom_compton; - aCRData.photon2_phantom_compton = photon2_phantom_compton; - aCRData.photon3_phantom_compton = photon3_phantom_compton; - strcpy(aCRData.theComptonVolumeName1 , theComptonVolumeName1.c_str() ); - strcpy(aCRData.theComptonVolumeName2 , theComptonVolumeName2.c_str() ); - strcpy(aCRData.theComptonVolumeName3 , theComptonVolumeName3.c_str() ); - strcpy(aCRData.theRayleighVolumeName1 , theRayleighVolumeName1.c_str() ); - strcpy(aCRData.theRayleighVolumeName2 , theRayleighVolumeName2.c_str() ); - strcpy(aCRData.theRayleighVolumeName3 , theRayleighVolumeName3.c_str() ); - gateToRoot->RecordPHData( aCRData ); + } + + // Counting Rayleigh scatter in phantom + if (processName.find("Rayl") != G4String::npos) { + if ((phantomTrackID == photon1ID)||(phantomTrackID == photon2ID)) { + G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4ThreeVector null(0.,0.,0.); + G4ThreeVector *ptr = &null; + theRayleighVolumeName = gNavigator->LocateGlobalPointAndSetup(hitPos,ptr,false)->GetName(); + if (nVerboseLevel > 1) { + G4cout << "GateAnalysis::RecordEndOfEvent : theRayleighVolumeName: " << theRayleighVolumeName << G4endl; } - - if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh - { - // we are in detector mode - GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); - ComptonRayleighData aCRData; - gateToRoot->GetPHData( aCRData); - - photon1_phantom_Rayleigh += aCRData.photon1_phantom_Rayleigh; - photon2_phantom_Rayleigh += aCRData.photon2_phantom_Rayleigh; - photon3_phantom_Rayleigh += aCRData.photon3_phantom_Rayleigh; - photon1_phantom_compton += aCRData.photon1_phantom_compton; - photon2_phantom_compton += aCRData.photon2_phantom_compton; - photon3_phantom_compton += aCRData.photon3_phantom_compton; - - theComptonVolumeName1 = aCRData.theComptonVolumeName1; - theComptonVolumeName2 = aCRData.theComptonVolumeName2; - theComptonVolumeName3 = aCRData.theComptonVolumeName3; - theRayleighVolumeName1 = aCRData.theRayleighVolumeName1; - theRayleighVolumeName2 = aCRData.theRayleighVolumeName2; - theRayleighVolumeName3 = aCRData.theRayleighVolumeName3; - + } + if (phantomTrackID == photon1ID) { + photon1_phantom_Rayleigh++; + theRayleighVolumeName1 = theRayleighVolumeName; + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon1_phantom_Rayleigh : " << photon1_phantom_Rayleigh << G4endl; + } + } + if (phantomTrackID == photon2ID) { + photon2_phantom_Rayleigh++; + theRayleighVolumeName2 = theRayleighVolumeName; + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon2_phantom_Rayleigh : " << photon2_phantom_Rayleigh << G4endl; + } + } + if (phantomTrackID == photon3ID) { + photon3_phantom_Rayleigh++; + theRayleighVolumeName3 = theRayleighVolumeName; + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon3_phantom_Rayleigh : " << photon3_phantom_Rayleigh << G4endl; } + } + } + } // end loop NpHits + + TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); + + if (theMode == TrackingMode::kTracker) { + // in tracker mode we store the infos about the number of compton and rayleigh + GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); + ComptonRayleighData aCRData; + aCRData.photon1_phantom_Rayleigh = photon1_phantom_Rayleigh; + aCRData.photon2_phantom_Rayleigh = photon2_phantom_Rayleigh; + aCRData.photon3_phantom_Rayleigh = photon3_phantom_Rayleigh; + aCRData.photon1_phantom_compton = photon1_phantom_compton; + aCRData.photon2_phantom_compton = photon2_phantom_compton; + aCRData.photon3_phantom_compton = photon3_phantom_compton; + strcpy(aCRData.theComptonVolumeName1 , theComptonVolumeName1.c_str() ); + strcpy(aCRData.theComptonVolumeName2 , theComptonVolumeName2.c_str() ); + strcpy(aCRData.theComptonVolumeName3 , theComptonVolumeName3.c_str() ); + strcpy(aCRData.theRayleighVolumeName1 , theRayleighVolumeName1.c_str() ); + strcpy(aCRData.theRayleighVolumeName2 , theRayleighVolumeName2.c_str() ); + strcpy(aCRData.theRayleighVolumeName3 , theRayleighVolumeName3.c_str() ); + gateToRoot->RecordPHData( aCRData ); + } else if ( theMode == TrackingMode::kDetector ) { + // in tracker mode we store the infos about the number of compton and rayleigh - we are in detector mode + GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); + ComptonRayleighData aCRData; + gateToRoot->GetPHData( aCRData); + photon1_phantom_Rayleigh += aCRData.photon1_phantom_Rayleigh; + photon2_phantom_Rayleigh += aCRData.photon2_phantom_Rayleigh; + photon3_phantom_Rayleigh += aCRData.photon3_phantom_Rayleigh; + photon1_phantom_compton += aCRData.photon1_phantom_compton; + photon2_phantom_compton += aCRData.photon2_phantom_compton; + photon3_phantom_compton += aCRData.photon3_phantom_compton; + + theComptonVolumeName1 = aCRData.theComptonVolumeName1; + theComptonVolumeName2 = aCRData.theComptonVolumeName2; + theComptonVolumeName3 = aCRData.theComptonVolumeName3; + theRayleighVolumeName1 = aCRData.theRayleighVolumeName1; + theRayleighVolumeName2 = aCRData.theRayleighVolumeName2; + theRayleighVolumeName3 = aCRData.theRayleighVolumeName3; + } + + // Source info + // DS : if gate source is not used (with /gate/EnableGeneralParticleSource) there are no GateSource, so skip + if ((GateSourceMgr::GetInstance())->GetSourcesForThisEvent().size() == 0) { + return; + } + + G4int sourceID = (((GateSourceMgr::GetInstance())->GetSourcesForThisEvent())[0])->GetSourceID(); + G4ThreeVector sourceVertex = m_trajectoryNavigator->FindSourcePosition(); + + // Hits loop + for (G4int iHit=0;iHitGetTrackID(); + G4String processName = (*CHC)[iHit]->GetProcess(); + // Counting Compton in the Crystal + if (processName.find("ompt") != G4String::npos) { + if (crystalTrackID == photon1ID) { + photon1_crystal_compton++; + } else if (crystalTrackID == photon2ID) { + photon2_crystal_compton++; + } else if (crystalTrackID == photon3ID) { + photon3_crystal_compton++; + } + } + // Counting Rayleigh scatter in crystal + if (processName.find("Rayl") != G4String::npos) { + if (crystalTrackID == photon1ID) { + photon1_crystal_Rayleigh++; + } else if (crystalTrackID == photon2ID) { + photon2_crystal_Rayleigh++; + } else if (crystalTrackID == photon3ID) { + photon3_crystal_Rayleigh++; + } + } + + G4int PDGEncoding = (*CHC)[iHit]->GetPDGEncoding(); + if (nVerboseLevel > 2) { + G4cout << "GateAnalysis::RecordEndOfEvent : CrystalHitsCollection: processName : <" << processName << G4endl; + G4cout << "> Particls PDG code : " << PDGEncoding << G4endl; + } + if ((*CHC)[iHit]->GoodForAnalysis()) { + // fill in values with the branch with C struct + G4int trackID = (*CHC)[iHit]->GetTrackID(); + G4int photonID = 0; + G4int nPhantomCompton = 0; + G4int nCrystalCompton = 0; + G4int nPhantomRayleigh = 0; + G4int nCrystalRayleigh = 0; + + if (photon1ID != 0) { + // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT + // search the gamma from which this hit comes --> photonID + photonID = m_trajectoryNavigator->FindPhotonID(trackID); + if (nVerboseLevel > 2 && photonID == rootID) { + G4cout << "GateAnalysis::RecordEndOfEvent : trackID: " << trackID << " photonID = " << rootID << G4endl; + } + } + if (photonID == 1) { + nPhantomCompton = photon1_phantom_compton; + nPhantomRayleigh = photon1_phantom_Rayleigh; + theComptonVolumeName = theComptonVolumeName1; + theRayleighVolumeName = theRayleighVolumeName1; + nCrystalCompton = photon1_crystal_compton; + nCrystalRayleigh = photon1_crystal_Rayleigh; + } else if (photonID == 2) { + nPhantomCompton = photon2_phantom_compton; + nPhantomRayleigh = photon2_phantom_Rayleigh; + theComptonVolumeName = theComptonVolumeName2; + theRayleighVolumeName = theRayleighVolumeName2; + nCrystalCompton = photon2_crystal_compton; + nCrystalRayleigh = photon2_crystal_Rayleigh; + } else if (photonID == 3) { + nPhantomCompton = photon3_phantom_compton; + nPhantomRayleigh = photon3_phantom_Rayleigh; + theComptonVolumeName = theComptonVolumeName3; + theRayleighVolumeName = theRayleighVolumeName3; + nCrystalCompton = photon3_crystal_compton; + nCrystalRayleigh = photon3_crystal_Rayleigh; + } - ///////// - // Source info - // DS : if gate source is not used (with /gate/EnableGeneralParticleSource) - // there are no GateSource, so skip - if ((GateSourceMgr::GetInstance())->GetSourcesForThisEvent().size() == 0) return; - - G4int sourceID = (((GateSourceMgr::GetInstance())->GetSourcesForThisEvent())[0])->GetSourceID(); - G4ThreeVector sourceVertex = m_trajectoryNavigator->FindSourcePosition(); - - // Hits loop - for (G4int iHit=0;iHitGetTrackID(); - G4String processName = (*CHC)[iHit]->GetProcess(); - // Counting Compton in the Crystal - if (processName.find("ompt") != G4String::npos) - { - - if (crystalTrackID == photon1ID) photon1_crystal_compton++; - if (crystalTrackID == photon2ID) photon2_crystal_compton++; - if (crystalTrackID == photon3ID) photon3_crystal_compton++; - } - - // Counting Rayleigh scatter in crystal - if (processName.find("Rayl") != G4String::npos) - { - - if (crystalTrackID == photon1ID) photon1_crystal_Rayleigh++; - if (crystalTrackID == photon2ID) photon2_crystal_Rayleigh++; - if (crystalTrackID == photon3ID) photon3_crystal_Rayleigh++; - } - - G4int PDGEncoding = (*CHC)[iHit]->GetPDGEncoding(); - if (nVerboseLevel > 2) - G4cout << "GateAnalysis::RecordEndOfEvent : CrystalHitsCollection: processName : <" << processName - << "> Particls PDG code : " << PDGEncoding << Gateendl; - if ((*CHC)[iHit]->GoodForAnalysis()) - { - // fill in values with the branch with C struct - - G4int trackID = (*CHC)[iHit]->GetTrackID(); - - G4int photonID = 0; - G4int nPhantomCompton = 0; - G4int nCrystalCompton = 0; - - G4int nPhantomRayleigh = 0; - G4int nCrystalRayleigh = 0; - - if (photon1ID != 0) - { // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT - // search the gamma from which this hit comes --> photonID - photonID = m_trajectoryNavigator->FindPhotonID(trackID); - if (photonID == rootID) - { - if (nVerboseLevel > 2) G4cout - << "GateAnalysis::RecordEndOfEvent : trackID: " << trackID << " photonID = " << rootID << Gateendl; - } - } - - if (photonID == 1) - { - nPhantomCompton = photon1_phantom_compton; - nPhantomRayleigh = photon1_phantom_Rayleigh; - theComptonVolumeName = theComptonVolumeName1; - theRayleighVolumeName = theRayleighVolumeName1; - nCrystalCompton = photon1_crystal_compton; - nCrystalRayleigh = photon1_crystal_Rayleigh; - } - else if (photonID == 2) - { - nPhantomCompton = photon2_phantom_compton; - nPhantomRayleigh = photon2_phantom_Rayleigh; - theComptonVolumeName = theComptonVolumeName2; - theRayleighVolumeName = theRayleighVolumeName2; - nCrystalCompton = photon2_crystal_compton; - nCrystalRayleigh = photon2_crystal_Rayleigh; - } - else if (photonID == 3) - { - nPhantomCompton = photon3_phantom_compton; - nPhantomRayleigh = photon3_phantom_Rayleigh; - theComptonVolumeName = theComptonVolumeName3; - theRayleighVolumeName = theRayleighVolumeName3; - nCrystalCompton = photon3_crystal_compton; - nCrystalRayleigh = photon3_crystal_Rayleigh; - } - - // search the primary that originated the track - primaryID = m_trajectoryNavigator->FindPrimaryID(trackID); - - (*CHC)[iHit]->SetSourceID (sourceID); - (*CHC)[iHit]->SetSourcePosition (sourceVertex); - (*CHC)[iHit]->SetNPhantomCompton (nPhantomCompton); - (*CHC)[iHit]->SetNPhantomRayleigh (nPhantomRayleigh); - (*CHC)[iHit]->SetComptonVolumeName (theComptonVolumeName); - (*CHC)[iHit]->SetRayleighVolumeName (theRayleighVolumeName); - (*CHC)[iHit]->SetPhotonID (photonID); - (*CHC)[iHit]->SetPrimaryID (primaryID); - (*CHC)[iHit]->SetEventID (eventID); - (*CHC)[iHit]->SetRunID (runID); - (*CHC)[iHit]->SetNCrystalCompton (nCrystalCompton); - (*CHC)[iHit]->SetNCrystalRayleigh (nCrystalRayleigh); - (*CHC)[iHit]->SetNSeptal (septalNb); // HDS : septal penetration - } - } - } // end if (CHC) - } // end if (!trajectoryContainer) + // search the primary that originated the track + primaryID = m_trajectoryNavigator->FindPrimaryID(trackID); + (*CHC)[iHit]->SetSourceID(sourceID); + (*CHC)[iHit]->SetSourcePosition(sourceVertex); + (*CHC)[iHit]->SetNPhantomCompton(nPhantomCompton); + (*CHC)[iHit]->SetNPhantomRayleigh(nPhantomRayleigh); + (*CHC)[iHit]->SetComptonVolumeName(theComptonVolumeName); + (*CHC)[iHit]->SetRayleighVolumeName(theRayleighVolumeName); + (*CHC)[iHit]->SetPhotonID(photonID); + (*CHC)[iHit]->SetPrimaryID(primaryID); + (*CHC)[iHit]->SetEventID(eventID); + (*CHC)[iHit]->SetRunID(runID); + (*CHC)[iHit]->SetNCrystalCompton(nCrystalCompton); + (*CHC)[iHit]->SetNCrystalRayleigh(nCrystalRayleigh); + (*CHC)[iHit]->SetNSeptal(septalNb); // HDS : septal penetration + } + } + } // end if (CHC) + } // end if (!trajectoryContainer) } // end function //-------------------------------------------------------------------------------------------------- From fb3f26e1aea53474ae490d8e56f0cfeca3406cb2 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 20:57:47 +0100 Subject: [PATCH 05/15] feat:new structure simplifies the RecordEndOfEvent method code --- source/digits_hits/include/GateAnalysis.hh | 10 + source/digits_hits/src/GateAnalysis.cc | 234 +++++++-------------- 2 files changed, 81 insertions(+), 163 deletions(-) diff --git a/source/digits_hits/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index b2f8133a3..ed131b10e 100644 --- a/source/digits_hits/include/GateAnalysis.hh +++ b/source/digits_hits/include/GateAnalysis.hh @@ -26,6 +26,16 @@ class GateAnalysis : public GateVOutputModule { public: + struct PhotonScatterings { + G4int nPhantomCompton = 0; + G4int nPhantomRayleigh = 0; + G4int nCrystalCompton = 0; + G4int nCrystalRayleigh = 0; + G4int photonID = 0; + G4String theComptonVolumeName = G4String("NULL"); + G4String theRayleighVolumeName = G4String("NULL"); + }; + GateAnalysis(const G4String& name, GateOutputMgr* outputMgr,DigiMode digiMode); virtual ~GateAnalysis(); const G4String& GiveNameOfFile(); diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 09d36f983..2c0bcaf34 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -137,23 +137,9 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (CHC) { NbHits = CHC->entries(); - G4int photon1ID = 0; - G4int photon2ID = 0; - G4int photon3ID = 0; + std::vector photon_scatterings(3); G4int rootID = 0; G4int primaryID = 0; - G4int photon1_phantom_compton = 0; - G4int photon2_phantom_compton = 0; - G4int photon3_phantom_compton = 0; - G4int photon1_crystal_compton = 0; - G4int photon2_crystal_compton = 0; - G4int photon3_crystal_compton = 0; - G4int photon1_phantom_Rayleigh = 0; - G4int photon2_phantom_Rayleigh = 0; - G4int photon3_phantom_Rayleigh = 0; - G4int photon1_crystal_Rayleigh = 0; - G4int photon2_crystal_Rayleigh = 0; - G4int photon3_crystal_Rayleigh = 0; G4int septalNb = 0; // HDS : septal penetration m_trajectoryNavigator->FindPositronTrackID(); @@ -169,24 +155,24 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (nVerboseLevel > 0 && photonIDVec.size() > 2) { G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonID vector size > 2" << G4endl; } - photon1ID = photonIDVec[0]; - photon2ID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; - photon3ID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; + photon_scatterings[0].photonID = photonIDVec[0]; + photon_scatterings[1].photonID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; + photon_scatterings[2].photonID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; } if (nVerboseLevel > 0) { - if (photon1ID == 0) { + if (photon_scatterings[0].photonID == 0) { G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon1ID == 0" << G4endl; } - if (photon2ID == 0) { + if (photon_scatterings[1].photonID == 0) { G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0" << G4endl; } - if (photon3ID == 0) { + if (photon_scatterings[2].photonID == 0) { G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon3ID == 0" << G4endl; } if (nVerboseLevel > 1) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon1ID; - G4cout << " photon2ID : " << photon2ID << G4endl; + G4cout << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon_scatterings[0].photonID; + G4cout << " photon2ID : " << photon_scatterings[1].photonID << G4endl; } } @@ -195,13 +181,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) GatePhantomHitsCollection* PHC = GetOutputMgr()->GetPhantomHitCollection(); NpHits = PHC->entries(); G4String theComptonVolumeName("NULL"); - G4String theComptonVolumeName1("NULL"); - G4String theComptonVolumeName2("NULL"); - G4String theComptonVolumeName3("NULL"); G4String theRayleighVolumeName("NULL"); - G4String theRayleighVolumeName1("NULL"); - G4String theRayleighVolumeName2("NULL"); - G4String theRayleighVolumeName3("NULL"); for (G4int iPHit=0;iPHitGetProcess(); G4int PDGcode = (*PHC)[iPHit]->GetPDGEncoding(); G4ThreeVector hitPos = (*PHC)[iPHit]->GetPos(); - theComptonVolumeName = G4String("NULL"); - theRayleighVolumeName = G4String("NULL"); if (nVerboseLevel > 2) { G4cout << "GateAnalysis::RecordEndOfEvent : GatePhantomHitsCollection : trackID : " << std::setw(5) << phantomTrackID; G4cout << " PDG code : " << std::setw(5) << PDGcode << " processName : <" << processName << G4endl; } - - // Modif by DS and LS on Oct 4, 2002: we need to be able to recognise both 'compt' - // and 'LowEnCompt", hence the find on 'ompt' - // if (processName.find("ompt") != G4String::npos || processName.find("Rayleigh") != G4String::npos) { - // modif. by CJG to separate Compton and Rayleigh photons - if (processName.find("ompt") != G4String::npos) { - if ((phantomTrackID == photon1ID)||(phantomTrackID == photon2ID)) { - G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4ThreeVector null(0.,0.,0.); - G4ThreeVector *ptr = &null; - theComptonVolumeName = gNavigator->LocateGlobalPointAndSetup(hitPos,ptr,false)->GetName(); - if (nVerboseLevel > 1) { - G4cout << "GateAnalysis::RecordEndOfEvent : theComptonVolumeName: " << theComptonVolumeName << G4endl; - } - } - if (phantomTrackID == photon1ID) { - photon1_phantom_compton++; - theComptonVolumeName1 = theComptonVolumeName; - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon1_phantom_compton : " << photon1_phantom_compton << G4endl; - } - } - if (phantomTrackID == photon2ID) { - photon2_phantom_compton++; - theComptonVolumeName2 = theComptonVolumeName; - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon2_phantom_compton : " << photon2_phantom_compton << G4endl; - } - } - if (phantomTrackID == photon3ID) { - photon3_phantom_compton++; - theComptonVolumeName3 = theComptonVolumeName; - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon3_phantom_compton : " << photon3_phantom_compton << G4endl; - } - } - } - - // Counting Rayleigh scatter in phantom - if (processName.find("Rayl") != G4String::npos) { - if ((phantomTrackID == photon1ID)||(phantomTrackID == photon2ID)) { + + if ((phantomTrackID == photon_scatterings[0].photonID)||(phantomTrackID == photon_scatterings[1].photonID)) { + //Modif by DS and LS on Oct 4, 2002: we need to be able to recognise both 'compt' + //and 'LowEnCompt", hence the find on 'ompt' modif. by CJG to separate Compton and Rayleigh photons + const bool isComptonInPhantom = processName.find("ompt") != G4String::npos; + const bool isRayleighInPhantom = processName.find("Rayl") != G4String::npos; + if (isComptonInPhantom || isRayleighInPhantom) { G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4ThreeVector null(0.,0.,0.); - G4ThreeVector *ptr = &null; - theRayleighVolumeName = gNavigator->LocateGlobalPointAndSetup(hitPos,ptr,false)->GetName(); - if (nVerboseLevel > 1) { - G4cout << "GateAnalysis::RecordEndOfEvent : theRayleighVolumeName: " << theRayleighVolumeName << G4endl; - } - } - if (phantomTrackID == photon1ID) { - photon1_phantom_Rayleigh++; - theRayleighVolumeName1 = theRayleighVolumeName; - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon1_phantom_Rayleigh : " << photon1_phantom_Rayleigh << G4endl; - } + G4ThreeVector zero_vector(0.,0.,0.); + G4String volume_name = gNavigator->LocateGlobalPointAndSetup(hitPos,&zero_vector,false)->GetName(); + theComptonVolumeName = isComptonInPhantom ? volume_name : theComptonVolumeName; + theRayleighVolumeName = isRayleighInPhantom ? volume_name : theRayleighVolumeName; } - if (phantomTrackID == photon2ID) { - photon2_phantom_Rayleigh++; - theRayleighVolumeName2 = theRayleighVolumeName; - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon2_phantom_Rayleigh : " << photon2_phantom_Rayleigh << G4endl; - } - } - if (phantomTrackID == photon3ID) { - photon3_phantom_Rayleigh++; - theRayleighVolumeName3 = theRayleighVolumeName; - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon3_phantom_Rayleigh : " << photon3_phantom_Rayleigh << G4endl; + auto found_photon_scatterings = find_if(photon_scatterings.begin(),photon_scatterings.end(), [&phantomTrackID](PhotonScatterings& ps) { return ps.photonID == phantomTrackID;}); + if (found_photon_scatterings != photon_scatterings.end()) { + if (isComptonInPhantom) { + found_photon_scatterings->nPhantomCompton += 1; + found_photon_scatterings->theComptonVolumeName = theComptonVolumeName; + } else { + found_photon_scatterings->nPhantomRayleigh += 1; + found_photon_scatterings->theRayleighVolumeName = theComptonVolumeName; } } } @@ -300,37 +231,37 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) // in tracker mode we store the infos about the number of compton and rayleigh GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); ComptonRayleighData aCRData; - aCRData.photon1_phantom_Rayleigh = photon1_phantom_Rayleigh; - aCRData.photon2_phantom_Rayleigh = photon2_phantom_Rayleigh; - aCRData.photon3_phantom_Rayleigh = photon3_phantom_Rayleigh; - aCRData.photon1_phantom_compton = photon1_phantom_compton; - aCRData.photon2_phantom_compton = photon2_phantom_compton; - aCRData.photon3_phantom_compton = photon3_phantom_compton; - strcpy(aCRData.theComptonVolumeName1 , theComptonVolumeName1.c_str() ); - strcpy(aCRData.theComptonVolumeName2 , theComptonVolumeName2.c_str() ); - strcpy(aCRData.theComptonVolumeName3 , theComptonVolumeName3.c_str() ); - strcpy(aCRData.theRayleighVolumeName1 , theRayleighVolumeName1.c_str() ); - strcpy(aCRData.theRayleighVolumeName2 , theRayleighVolumeName2.c_str() ); - strcpy(aCRData.theRayleighVolumeName3 , theRayleighVolumeName3.c_str() ); + aCRData.photon1_phantom_Rayleigh = photon_scatterings[0].nPhantomRayleigh; + aCRData.photon2_phantom_Rayleigh = photon_scatterings[1].nPhantomRayleigh; + aCRData.photon3_phantom_Rayleigh = photon_scatterings[2].nPhantomRayleigh; + aCRData.photon1_phantom_compton = photon_scatterings[0].nPhantomCompton; + aCRData.photon2_phantom_compton = photon_scatterings[1].nPhantomCompton; + aCRData.photon3_phantom_compton = photon_scatterings[2].nPhantomCompton; + strcpy(aCRData.theComptonVolumeName1 , photon_scatterings[0].theComptonVolumeName.c_str() ); + strcpy(aCRData.theComptonVolumeName2 , photon_scatterings[1].theComptonVolumeName.c_str() ); + strcpy(aCRData.theComptonVolumeName3 , photon_scatterings[2].theComptonVolumeName.c_str() ); + strcpy(aCRData.theRayleighVolumeName1 , photon_scatterings[0].theRayleighVolumeName.c_str() ); + strcpy(aCRData.theRayleighVolumeName2 , photon_scatterings[1].theRayleighVolumeName.c_str() ); + strcpy(aCRData.theRayleighVolumeName3 , photon_scatterings[2].theRayleighVolumeName.c_str() ); gateToRoot->RecordPHData( aCRData ); } else if ( theMode == TrackingMode::kDetector ) { // in tracker mode we store the infos about the number of compton and rayleigh - we are in detector mode GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); ComptonRayleighData aCRData; gateToRoot->GetPHData( aCRData); - photon1_phantom_Rayleigh += aCRData.photon1_phantom_Rayleigh; - photon2_phantom_Rayleigh += aCRData.photon2_phantom_Rayleigh; - photon3_phantom_Rayleigh += aCRData.photon3_phantom_Rayleigh; - photon1_phantom_compton += aCRData.photon1_phantom_compton; - photon2_phantom_compton += aCRData.photon2_phantom_compton; - photon3_phantom_compton += aCRData.photon3_phantom_compton; - - theComptonVolumeName1 = aCRData.theComptonVolumeName1; - theComptonVolumeName2 = aCRData.theComptonVolumeName2; - theComptonVolumeName3 = aCRData.theComptonVolumeName3; - theRayleighVolumeName1 = aCRData.theRayleighVolumeName1; - theRayleighVolumeName2 = aCRData.theRayleighVolumeName2; - theRayleighVolumeName3 = aCRData.theRayleighVolumeName3; + photon_scatterings[0].nPhantomRayleigh += aCRData.photon1_phantom_Rayleigh; + photon_scatterings[1].nPhantomRayleigh += aCRData.photon2_phantom_Rayleigh; + photon_scatterings[2].nPhantomRayleigh += aCRData.photon3_phantom_Rayleigh; + photon_scatterings[0].nPhantomCompton += aCRData.photon1_phantom_compton; + photon_scatterings[1].nPhantomCompton += aCRData.photon2_phantom_compton; + photon_scatterings[2].nPhantomCompton += aCRData.photon3_phantom_compton; + + photon_scatterings[0].theComptonVolumeName = aCRData.theComptonVolumeName1; + photon_scatterings[1].theComptonVolumeName = aCRData.theComptonVolumeName2; + photon_scatterings[2].theComptonVolumeName = aCRData.theComptonVolumeName3; + photon_scatterings[0].theRayleighVolumeName = aCRData.theRayleighVolumeName1; + photon_scatterings[1].theRayleighVolumeName = aCRData.theRayleighVolumeName2; + photon_scatterings[2].theRayleighVolumeName = aCRData.theRayleighVolumeName3; } // Source info @@ -346,24 +277,14 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) for (G4int iHit=0;iHitGetTrackID(); G4String processName = (*CHC)[iHit]->GetProcess(); - // Counting Compton in the Crystal - if (processName.find("ompt") != G4String::npos) { - if (crystalTrackID == photon1ID) { - photon1_crystal_compton++; - } else if (crystalTrackID == photon2ID) { - photon2_crystal_compton++; - } else if (crystalTrackID == photon3ID) { - photon3_crystal_compton++; - } - } - // Counting Rayleigh scatter in crystal - if (processName.find("Rayl") != G4String::npos) { - if (crystalTrackID == photon1ID) { - photon1_crystal_Rayleigh++; - } else if (crystalTrackID == photon2ID) { - photon2_crystal_Rayleigh++; - } else if (crystalTrackID == photon3ID) { - photon3_crystal_Rayleigh++; + auto found_photon_scatterings_crystal = find_if(photon_scatterings.begin(),photon_scatterings.end(), [&crystalTrackID](PhotonScatterings& ps) { return ps.photonID == crystalTrackID;}); + if (found_photon_scatterings_crystal != photon_scatterings.end()) { + if (processName.find("ompt") != G4String::npos) { + // Counting Compton in the Crystal + found_photon_scatterings_crystal->nCrystalCompton += 1; + } else if (processName.find("Rayl") != G4String::npos) { + // Counting Rayleigh scatter in crystal + found_photon_scatterings_crystal->nCrystalRayleigh += 1; } } @@ -381,7 +302,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4int nPhantomRayleigh = 0; G4int nCrystalRayleigh = 0; - if (photon1ID != 0) { + if (photon_scatterings[0].photonID != 0) { // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT // search the gamma from which this hit comes --> photonID photonID = m_trajectoryNavigator->FindPhotonID(trackID); @@ -389,28 +310,15 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4cout << "GateAnalysis::RecordEndOfEvent : trackID: " << trackID << " photonID = " << rootID << G4endl; } } - - if (photonID == 1) { - nPhantomCompton = photon1_phantom_compton; - nPhantomRayleigh = photon1_phantom_Rayleigh; - theComptonVolumeName = theComptonVolumeName1; - theRayleighVolumeName = theRayleighVolumeName1; - nCrystalCompton = photon1_crystal_compton; - nCrystalRayleigh = photon1_crystal_Rayleigh; - } else if (photonID == 2) { - nPhantomCompton = photon2_phantom_compton; - nPhantomRayleigh = photon2_phantom_Rayleigh; - theComptonVolumeName = theComptonVolumeName2; - theRayleighVolumeName = theRayleighVolumeName2; - nCrystalCompton = photon2_crystal_compton; - nCrystalRayleigh = photon2_crystal_Rayleigh; - } else if (photonID == 3) { - nPhantomCompton = photon3_phantom_compton; - nPhantomRayleigh = photon3_phantom_Rayleigh; - theComptonVolumeName = theComptonVolumeName3; - theRayleighVolumeName = theRayleighVolumeName3; - nCrystalCompton = photon3_crystal_compton; - nCrystalRayleigh = photon3_crystal_Rayleigh; + + if (photonID > 0 && photonID < 4) { + const G4int index = photonID - 1; + nPhantomCompton = photon_scatterings[index].nPhantomCompton; + nPhantomRayleigh = photon_scatterings[index].nPhantomRayleigh; + theComptonVolumeName = photon_scatterings[index].theComptonVolumeName; + theRayleighVolumeName = photon_scatterings[index].theRayleighVolumeName; + nCrystalCompton = photon_scatterings[index].nCrystalCompton; + nCrystalRayleigh = photon_scatterings[index].nCrystalRayleigh; } // search the primary that originated the track From be84582f9df98ae4493ed5090c475987d6c7020d Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 21:08:11 +0100 Subject: [PATCH 06/15] feat:separation of phantom scatterings collection from the RecordEndOfEvent method makes RecordEndOfEvent method easier to understand --- source/digits_hits/include/GateAnalysis.hh | 2 + source/digits_hits/src/GateAnalysis.cc | 100 +++++++++++---------- 2 files changed, 53 insertions(+), 49 deletions(-) diff --git a/source/digits_hits/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index ed131b10e..408a02cb5 100644 --- a/source/digits_hits/include/GateAnalysis.hh +++ b/source/digits_hits/include/GateAnalysis.hh @@ -67,6 +67,8 @@ public: inline G4bool GetRecordSeptalFlag() const { return m_recordSeptalFlag; } inline void SetSeptalPhysVolumeName(const G4String& name) { m_septalPhysVolumeName = name; } inline void SetRecordSeptalFlag(G4bool flag) { m_recordSeptalFlag = flag; } + + void CollectPhantomScatterings(std::vector& photon_scatterings, G4int& septalNb); private: diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 2c0bcaf34..cee07806e 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -110,7 +110,54 @@ void GateAnalysis::RecordBeginOfEvent(const G4Event* ) } //-------------------------------------------------------------------------------------------------- - +void GateAnalysis::CollectPhantomScatterings(std::vector& photon_scatterings, G4int& septalNb) { + GatePhantomHitsCollection* PHC = GetOutputMgr()->GetPhantomHitCollection(); + G4int NpHits = PHC->entries(); + G4String theComptonVolumeName("NULL"); + G4String theRayleighVolumeName("NULL"); + + for (G4int iPHit=0;iPHitGetPhysVolName() == m_septalPhysVolumeName) { + ++septalNb; + } + } + G4int phantomTrackID = (*PHC)[iPHit]->GetTrackID(); + G4String processName = (*PHC)[iPHit]->GetProcess(); + G4int PDGcode = (*PHC)[iPHit]->GetPDGEncoding(); + G4ThreeVector hitPos = (*PHC)[iPHit]->GetPos(); + + if (nVerboseLevel > 2) { + G4cout << "GateAnalysis::RecordEndOfEvent : GatePhantomHitsCollection : trackID : " << std::setw(5) << phantomTrackID; + G4cout << " PDG code : " << std::setw(5) << PDGcode << " processName : <" << processName << G4endl; + } + + if ((phantomTrackID == photon_scatterings[0].photonID)||(phantomTrackID == photon_scatterings[1].photonID)) { + //Modif by DS and LS on Oct 4, 2002: we need to be able to recognise both 'compt' + //and 'LowEnCompt", hence the find on 'ompt' modif. by CJG to separate Compton and Rayleigh photons + const bool isComptonInPhantom = processName.find("ompt") != G4String::npos; + const bool isRayleighInPhantom = processName.find("Rayl") != G4String::npos; + if (isComptonInPhantom || isRayleighInPhantom) { + G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); + G4ThreeVector zero_vector(0.,0.,0.); + G4String volume_name = gNavigator->LocateGlobalPointAndSetup(hitPos,&zero_vector,false)->GetName(); + theComptonVolumeName = isComptonInPhantom ? volume_name : theComptonVolumeName; + theRayleighVolumeName = isRayleighInPhantom ? volume_name : theRayleighVolumeName; + } + auto found_photon_scatterings = find_if(photon_scatterings.begin(),photon_scatterings.end(), [&phantomTrackID](PhotonScatterings& ps) { return ps.photonID == phantomTrackID;}); + if (found_photon_scatterings != photon_scatterings.end()) { + if (isComptonInPhantom) { + found_photon_scatterings->nPhantomCompton += 1; + found_photon_scatterings->theComptonVolumeName = theComptonVolumeName; + } else { + found_photon_scatterings->nPhantomRayleigh += 1; + found_photon_scatterings->theRayleighVolumeName = theComptonVolumeName; + } + } + } + } +} //-------------------------------------------------------------------------------------------------- void GateAnalysis::RecordEndOfEvent(const G4Event* event) { @@ -133,7 +180,6 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) } else { GateCrystalHitsCollection* CHC = GetOutputMgr()->GetCrystalHitCollection(); G4int NbHits = 0; - G4int NpHits = 0; if (CHC) { NbHits = CHC->entries(); @@ -177,53 +223,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) } // analysis of the phantom hits to count the comptons, etc. - - GatePhantomHitsCollection* PHC = GetOutputMgr()->GetPhantomHitCollection(); - NpHits = PHC->entries(); - G4String theComptonVolumeName("NULL"); - G4String theRayleighVolumeName("NULL"); - - for (G4int iPHit=0;iPHitGetPhysVolName() == m_septalPhysVolumeName) { - ++septalNb; - } - } - G4int phantomTrackID = (*PHC)[iPHit]->GetTrackID(); - G4String processName = (*PHC)[iPHit]->GetProcess(); - G4int PDGcode = (*PHC)[iPHit]->GetPDGEncoding(); - G4ThreeVector hitPos = (*PHC)[iPHit]->GetPos(); - - if (nVerboseLevel > 2) { - G4cout << "GateAnalysis::RecordEndOfEvent : GatePhantomHitsCollection : trackID : " << std::setw(5) << phantomTrackID; - G4cout << " PDG code : " << std::setw(5) << PDGcode << " processName : <" << processName << G4endl; - } - - if ((phantomTrackID == photon_scatterings[0].photonID)||(phantomTrackID == photon_scatterings[1].photonID)) { - //Modif by DS and LS on Oct 4, 2002: we need to be able to recognise both 'compt' - //and 'LowEnCompt", hence the find on 'ompt' modif. by CJG to separate Compton and Rayleigh photons - const bool isComptonInPhantom = processName.find("ompt") != G4String::npos; - const bool isRayleighInPhantom = processName.find("Rayl") != G4String::npos; - if (isComptonInPhantom || isRayleighInPhantom) { - G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4ThreeVector zero_vector(0.,0.,0.); - G4String volume_name = gNavigator->LocateGlobalPointAndSetup(hitPos,&zero_vector,false)->GetName(); - theComptonVolumeName = isComptonInPhantom ? volume_name : theComptonVolumeName; - theRayleighVolumeName = isRayleighInPhantom ? volume_name : theRayleighVolumeName; - } - auto found_photon_scatterings = find_if(photon_scatterings.begin(),photon_scatterings.end(), [&phantomTrackID](PhotonScatterings& ps) { return ps.photonID == phantomTrackID;}); - if (found_photon_scatterings != photon_scatterings.end()) { - if (isComptonInPhantom) { - found_photon_scatterings->nPhantomCompton += 1; - found_photon_scatterings->theComptonVolumeName = theComptonVolumeName; - } else { - found_photon_scatterings->nPhantomRayleigh += 1; - found_photon_scatterings->theRayleighVolumeName = theComptonVolumeName; - } - } - } - } // end loop NpHits + CollectPhantomScatterings(photon_scatterings, septalNb); TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); @@ -301,6 +301,8 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) G4int nCrystalCompton = 0; G4int nPhantomRayleigh = 0; G4int nCrystalRayleigh = 0; + G4String theComptonVolumeName("NULL"); + G4String theRayleighVolumeName("NULL"); if (photon_scatterings[0].photonID != 0) { // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT From 8d719fab295eb7a4cd66bfc0c40f82182d8f4f13 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 21:27:55 +0100 Subject: [PATCH 07/15] fix:class GateSteppingAction has const method GetMode --- source/digits_hits/include/GateActions.hh | 2 +- source/digits_hits/src/GateActions.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/digits_hits/include/GateActions.hh b/source/digits_hits/include/GateActions.hh index 564defd9d..54dbf7513 100644 --- a/source/digits_hits/include/GateActions.hh +++ b/source/digits_hits/include/GateActions.hh @@ -142,7 +142,7 @@ public : void StopOnBoundary(G4int aI); void StopAndKill(G4String aString); void SetMode( TrackingMode aMode); - TrackingMode GetMode(); + TrackingMode GetMode() const; void SetTxtOut(G4String aString); G4int GetTxtOn() { return TxtOn;}; void SetEnergyThreshold(G4double); diff --git a/source/digits_hits/src/GateActions.cc b/source/digits_hits/src/GateActions.cc index 4ec4a4dc6..c8ea93556 100644 --- a/source/digits_hits/src/GateActions.cc +++ b/source/digits_hits/src/GateActions.cc @@ -589,7 +589,7 @@ void GateSteppingAction::SetMode( TrackingMode aMode) m_trackingMode = aMode; } -TrackingMode GateSteppingAction::GetMode() +TrackingMode GateSteppingAction::GetMode() const { return m_trackingMode;} //----------------------------------------------------------------------------- From 414ebabb186282e969829cc8ee91350d95dc82aa Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 21:29:06 +0100 Subject: [PATCH 08/15] feat:separation of ComptonRayleighData updates from the RecordEndOfEvent make code more readable --- source/digits_hits/include/GateAnalysis.hh | 3 + source/digits_hits/src/GateAnalysis.cc | 90 ++++++++++++---------- 2 files changed, 54 insertions(+), 39 deletions(-) diff --git a/source/digits_hits/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index 408a02cb5..4d794bf5e 100644 --- a/source/digits_hits/include/GateAnalysis.hh +++ b/source/digits_hits/include/GateAnalysis.hh @@ -69,6 +69,9 @@ public: inline void SetRecordSeptalFlag(G4bool flag) { m_recordSeptalFlag = flag; } void CollectPhantomScatterings(std::vector& photon_scatterings, G4int& septalNb); + void UpdateComptonRayleighDataFromScatterings(const std::vector& photon_scatterings); + void UpdateScatteringsFromComptonRayleighData(std::vector& photon_scatterings); + void MakeComptonRayleighDataUpdates(std::vector& photon_scatterings); private: diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index cee07806e..fe5a4855b 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -33,6 +33,7 @@ #include "GateVVolume.hh" #include "GateActions.hh" #include "GateToRoot.hh" +#include "GateActions.hh" //-------------------------------------------------------------------------------------------------- GateAnalysis::GateAnalysis(const G4String& name, GateOutputMgr* outputMgr,DigiMode digiMode) : GateVOutputModule(name,outputMgr,digiMode) @@ -158,6 +159,55 @@ void GateAnalysis::CollectPhantomScatterings(std::vector& pho } } } + +void GateAnalysis::UpdateComptonRayleighDataFromScatterings(const std::vector& photon_scatterings) { + GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); + ComptonRayleighData aCRData; + aCRData.photon1_phantom_Rayleigh = photon_scatterings.at(0).nPhantomRayleigh; + aCRData.photon2_phantom_Rayleigh = photon_scatterings.at(1).nPhantomRayleigh; + aCRData.photon3_phantom_Rayleigh = photon_scatterings.at(2).nPhantomRayleigh; + aCRData.photon1_phantom_compton = photon_scatterings.at(0).nPhantomCompton; + aCRData.photon2_phantom_compton = photon_scatterings.at(1).nPhantomCompton; + aCRData.photon3_phantom_compton = photon_scatterings.at(2).nPhantomCompton; + strcpy(aCRData.theComptonVolumeName1 , photon_scatterings.at(0).theComptonVolumeName.c_str() ); + strcpy(aCRData.theComptonVolumeName2 , photon_scatterings.at(1).theComptonVolumeName.c_str() ); + strcpy(aCRData.theComptonVolumeName3 , photon_scatterings.at(2).theComptonVolumeName.c_str() ); + strcpy(aCRData.theRayleighVolumeName1 , photon_scatterings.at(0).theRayleighVolumeName.c_str() ); + strcpy(aCRData.theRayleighVolumeName2 , photon_scatterings.at(1).theRayleighVolumeName.c_str() ); + strcpy(aCRData.theRayleighVolumeName3 , photon_scatterings.at(2).theRayleighVolumeName.c_str() ); + gateToRoot->RecordPHData( aCRData ); +} + +void GateAnalysis::UpdateScatteringsFromComptonRayleighData(std::vector& photon_scatterings) { + GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); + ComptonRayleighData aCRData; + gateToRoot->GetPHData( aCRData); + photon_scatterings[0].nPhantomRayleigh += aCRData.photon1_phantom_Rayleigh; + photon_scatterings[1].nPhantomRayleigh += aCRData.photon2_phantom_Rayleigh; + photon_scatterings[2].nPhantomRayleigh += aCRData.photon3_phantom_Rayleigh; + photon_scatterings[0].nPhantomCompton += aCRData.photon1_phantom_compton; + photon_scatterings[1].nPhantomCompton += aCRData.photon2_phantom_compton; + photon_scatterings[2].nPhantomCompton += aCRData.photon3_phantom_compton; + + photon_scatterings[0].theComptonVolumeName = aCRData.theComptonVolumeName1; + photon_scatterings[1].theComptonVolumeName = aCRData.theComptonVolumeName2; + photon_scatterings[2].theComptonVolumeName = aCRData.theComptonVolumeName3; + photon_scatterings[0].theRayleighVolumeName = aCRData.theRayleighVolumeName1; + photon_scatterings[1].theRayleighVolumeName = aCRData.theRayleighVolumeName2; + photon_scatterings[2].theRayleighVolumeName = aCRData.theRayleighVolumeName3; +} + +void GateAnalysis::MakeComptonRayleighDataUpdates(std::vector& photon_scatterings) { + auto action = dynamic_cast(GateRunManager::GetRunManager()->GetUserSteppingAction()); + TrackingMode theMode = action->GetMode(); + if (theMode == TrackingMode::kTracker) { + // in tracker mode we store the infos about the number of compton and rayleigh + UpdateComptonRayleighDataFromScatterings(photon_scatterings); + } else if ( theMode == TrackingMode::kDetector ) { + // in tracker mode we store the infos about the number of compton and rayleigh - we are in detector mode + UpdateScatteringsFromComptonRayleighData(photon_scatterings); + } +} //-------------------------------------------------------------------------------------------------- void GateAnalysis::RecordEndOfEvent(const G4Event* event) { @@ -224,45 +274,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) // analysis of the phantom hits to count the comptons, etc. CollectPhantomScatterings(photon_scatterings, septalNb); - - TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); - - if (theMode == TrackingMode::kTracker) { - // in tracker mode we store the infos about the number of compton and rayleigh - GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); - ComptonRayleighData aCRData; - aCRData.photon1_phantom_Rayleigh = photon_scatterings[0].nPhantomRayleigh; - aCRData.photon2_phantom_Rayleigh = photon_scatterings[1].nPhantomRayleigh; - aCRData.photon3_phantom_Rayleigh = photon_scatterings[2].nPhantomRayleigh; - aCRData.photon1_phantom_compton = photon_scatterings[0].nPhantomCompton; - aCRData.photon2_phantom_compton = photon_scatterings[1].nPhantomCompton; - aCRData.photon3_phantom_compton = photon_scatterings[2].nPhantomCompton; - strcpy(aCRData.theComptonVolumeName1 , photon_scatterings[0].theComptonVolumeName.c_str() ); - strcpy(aCRData.theComptonVolumeName2 , photon_scatterings[1].theComptonVolumeName.c_str() ); - strcpy(aCRData.theComptonVolumeName3 , photon_scatterings[2].theComptonVolumeName.c_str() ); - strcpy(aCRData.theRayleighVolumeName1 , photon_scatterings[0].theRayleighVolumeName.c_str() ); - strcpy(aCRData.theRayleighVolumeName2 , photon_scatterings[1].theRayleighVolumeName.c_str() ); - strcpy(aCRData.theRayleighVolumeName3 , photon_scatterings[2].theRayleighVolumeName.c_str() ); - gateToRoot->RecordPHData( aCRData ); - } else if ( theMode == TrackingMode::kDetector ) { - // in tracker mode we store the infos about the number of compton and rayleigh - we are in detector mode - GateToRoot* gateToRoot = (GateToRoot*) (GateOutputMgr::GetInstance()->GetModule("root")); - ComptonRayleighData aCRData; - gateToRoot->GetPHData( aCRData); - photon_scatterings[0].nPhantomRayleigh += aCRData.photon1_phantom_Rayleigh; - photon_scatterings[1].nPhantomRayleigh += aCRData.photon2_phantom_Rayleigh; - photon_scatterings[2].nPhantomRayleigh += aCRData.photon3_phantom_Rayleigh; - photon_scatterings[0].nPhantomCompton += aCRData.photon1_phantom_compton; - photon_scatterings[1].nPhantomCompton += aCRData.photon2_phantom_compton; - photon_scatterings[2].nPhantomCompton += aCRData.photon3_phantom_compton; - - photon_scatterings[0].theComptonVolumeName = aCRData.theComptonVolumeName1; - photon_scatterings[1].theComptonVolumeName = aCRData.theComptonVolumeName2; - photon_scatterings[2].theComptonVolumeName = aCRData.theComptonVolumeName3; - photon_scatterings[0].theRayleighVolumeName = aCRData.theRayleighVolumeName1; - photon_scatterings[1].theRayleighVolumeName = aCRData.theRayleighVolumeName2; - photon_scatterings[2].theRayleighVolumeName = aCRData.theRayleighVolumeName3; - } + MakeComptonRayleighDataUpdates(photon_scatterings); // Source info // DS : if gate source is not used (with /gate/EnableGeneralParticleSource) there are no GateSource, so skip From 15e30df5ed527e312fcdbbe7ce02d1cbe93b5540 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 22:27:10 +0100 Subject: [PATCH 09/15] feat:separation of crystal scatterings collection from the RecordEndOfEvent method makes this method easier to understand --- source/digits_hits/include/GateAnalysis.hh | 11 ++ source/digits_hits/src/GateAnalysis.cc | 171 +++++++++++---------- 2 files changed, 103 insertions(+), 79 deletions(-) diff --git a/source/digits_hits/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index 4d794bf5e..3064c56e0 100644 --- a/source/digits_hits/include/GateAnalysis.hh +++ b/source/digits_hits/include/GateAnalysis.hh @@ -14,6 +14,7 @@ See LICENSE.md for further details #define GateAnalysis_H #include "GateVOutputModule.hh" +#include "GateCrystalHit.hh" class GateTrajectoryNavigator; @@ -72,6 +73,16 @@ public: void UpdateComptonRayleighDataFromScatterings(const std::vector& photon_scatterings); void UpdateScatteringsFromComptonRayleighData(std::vector& photon_scatterings); void MakeComptonRayleighDataUpdates(std::vector& photon_scatterings); + void SetCrystalScatterings(std::vector& photon_scatterings, GateCrystalHit* hit); + void UpdateHitDataForAnalysis( + std::vector& photon_scatterings, GateCrystalHit* hit, + const G4int septalNb, + const G4int sourceID, + const G4int eventID, + const G4int runID, + const G4ThreeVector& sourceVertex + ); + void CollectCrystalScatterings(std::vector& photon_scatterings, const G4int septalNb, GateCrystalHitsCollection* CHC, const G4int eventID); private: diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index fe5a4855b..56e3f9343 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -208,9 +208,100 @@ void GateAnalysis::MakeComptonRayleighDataUpdates(std::vector UpdateScatteringsFromComptonRayleighData(photon_scatterings); } } + +void GateAnalysis::SetCrystalScatterings(std::vector& photon_scatterings, GateCrystalHit* hit) { + G4int crystalTrackID = hit->GetTrackID(); + G4String processName = hit->GetProcess(); + auto found_photon_scatterings_crystal = find_if(photon_scatterings.begin(),photon_scatterings.end(), [&crystalTrackID](PhotonScatterings& ps) { return ps.photonID == crystalTrackID;}); + if (found_photon_scatterings_crystal != photon_scatterings.end()) { + if (processName.find("ompt") != G4String::npos) { + // Counting Compton in the Crystal + found_photon_scatterings_crystal->nCrystalCompton += 1; + } else if (processName.find("Rayl") != G4String::npos) { + // Counting Rayleigh scatter in crystal + found_photon_scatterings_crystal->nCrystalRayleigh += 1; + } + } +} + +void GateAnalysis::UpdateHitDataForAnalysis( + std::vector& photon_scatterings, GateCrystalHit* hit, + const G4int septalNb, + const G4int sourceID, + const G4int eventID, + const G4int runID, + const G4ThreeVector& sourceVertex +) { + // fill in values with the branch with C struct + G4int trackID = hit->GetTrackID(); + G4int primaryID = m_trajectoryNavigator->FindPrimaryID(trackID); + G4int photonID = 0; + G4int nPhantomCompton = 0; + G4int nCrystalCompton = 0; + G4int nPhantomRayleigh = 0; + G4int nCrystalRayleigh = 0; + G4String theComptonVolumeName("NULL"); + G4String theRayleighVolumeName("NULL"); + const G4int rootID = 0; + + if (photon_scatterings[0].photonID != 0) { + // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT + // search the gamma from which this hit comes --> photonID + photonID = m_trajectoryNavigator->FindPhotonID(trackID); + if (nVerboseLevel > 2 && photonID == rootID) { + G4cout << "GateAnalysis::RecordEndOfEvent : trackID: " << trackID << " photonID = " << rootID << G4endl; + } + } + + if (photonID > 0 && photonID < 4) { + const G4int index = photonID - 1; + nPhantomCompton = photon_scatterings[index].nPhantomCompton; + nPhantomRayleigh = photon_scatterings[index].nPhantomRayleigh; + theComptonVolumeName = photon_scatterings[index].theComptonVolumeName; + theRayleighVolumeName = photon_scatterings[index].theRayleighVolumeName; + nCrystalCompton = photon_scatterings[index].nCrystalCompton; + nCrystalRayleigh = photon_scatterings[index].nCrystalRayleigh; + } + + // search the primary that originated the track + hit->SetSourceID(sourceID); + hit->SetSourcePosition(sourceVertex); + hit->SetNPhantomCompton(nPhantomCompton); + hit->SetNPhantomRayleigh(nPhantomRayleigh); + hit->SetComptonVolumeName(theComptonVolumeName); + hit->SetRayleighVolumeName(theRayleighVolumeName); + hit->SetPhotonID(photonID); + hit->SetPrimaryID(primaryID); + hit->SetEventID(eventID); + hit->SetRunID(runID); + hit->SetNCrystalCompton(nCrystalCompton); + hit->SetNCrystalRayleigh(nCrystalRayleigh); + hit->SetNSeptal(septalNb); // HDS : septal penetration +} + +void GateAnalysis::CollectCrystalScatterings(std::vector& photon_scatterings, const G4int septalNb, GateCrystalHitsCollection* CHC, const G4int eventID) { + G4int NbHits = CHC->entries();; + G4int sourceID = (((GateSourceMgr::GetInstance())->GetSourcesForThisEvent())[0])->GetSourceID(); + G4int runID = GateRunManager::GetRunManager()->GetCurrentRun()->GetRunID(); + G4ThreeVector sourceVertex = m_trajectoryNavigator->FindSourcePosition(); + // Hits loop + for (G4int iHit=0;iHit 2) { + G4cout << "GateAnalysis::RecordEndOfEvent : CrystalHitsCollection: processName : <" << hit->GetProcess() << G4endl; + G4cout << "> Particls PDG code : " << hit->GetPDGEncoding() << G4endl; + } + if (hit->GoodForAnalysis()) { + UpdateHitDataForAnalysis(photon_scatterings,hit, septalNb, sourceID, eventID,runID,sourceVertex); + } + } +} + //-------------------------------------------------------------------------------------------------- void GateAnalysis::RecordEndOfEvent(const G4Event* event) { + const G4int eventID = event->GetEventID(); if (nVerboseLevel > 2) G4cout << "GateAnalysis::RecordEndOfEvent" << G4endl; @@ -220,22 +311,14 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) m_trajectoryNavigator->SetTrajectoryContainer(trajectoryContainer); } - G4int eventID = event->GetEventID(); - G4int runID = GateRunManager::GetRunManager()->GetCurrentRun()->GetRunID(); - if (!trajectoryContainer) { if (nVerboseLevel > 0) { G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : G4TrajectoryContainer not found" << G4endl; } } else { GateCrystalHitsCollection* CHC = GetOutputMgr()->GetCrystalHitCollection(); - G4int NbHits = 0; - if (CHC) { - NbHits = CHC->entries(); std::vector photon_scatterings(3); - G4int rootID = 0; - G4int primaryID = 0; G4int septalNb = 0; // HDS : septal penetration m_trajectoryNavigator->FindPositronTrackID(); @@ -281,77 +364,7 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if ((GateSourceMgr::GetInstance())->GetSourcesForThisEvent().size() == 0) { return; } - - G4int sourceID = (((GateSourceMgr::GetInstance())->GetSourcesForThisEvent())[0])->GetSourceID(); - G4ThreeVector sourceVertex = m_trajectoryNavigator->FindSourcePosition(); - - // Hits loop - for (G4int iHit=0;iHitGetTrackID(); - G4String processName = (*CHC)[iHit]->GetProcess(); - auto found_photon_scatterings_crystal = find_if(photon_scatterings.begin(),photon_scatterings.end(), [&crystalTrackID](PhotonScatterings& ps) { return ps.photonID == crystalTrackID;}); - if (found_photon_scatterings_crystal != photon_scatterings.end()) { - if (processName.find("ompt") != G4String::npos) { - // Counting Compton in the Crystal - found_photon_scatterings_crystal->nCrystalCompton += 1; - } else if (processName.find("Rayl") != G4String::npos) { - // Counting Rayleigh scatter in crystal - found_photon_scatterings_crystal->nCrystalRayleigh += 1; - } - } - - G4int PDGEncoding = (*CHC)[iHit]->GetPDGEncoding(); - if (nVerboseLevel > 2) { - G4cout << "GateAnalysis::RecordEndOfEvent : CrystalHitsCollection: processName : <" << processName << G4endl; - G4cout << "> Particls PDG code : " << PDGEncoding << G4endl; - } - if ((*CHC)[iHit]->GoodForAnalysis()) { - // fill in values with the branch with C struct - G4int trackID = (*CHC)[iHit]->GetTrackID(); - G4int photonID = 0; - G4int nPhantomCompton = 0; - G4int nCrystalCompton = 0; - G4int nPhantomRayleigh = 0; - G4int nCrystalRayleigh = 0; - G4String theComptonVolumeName("NULL"); - G4String theRayleighVolumeName("NULL"); - - if (photon_scatterings[0].photonID != 0) { - // this means that at least 1 photon has been found, requiring 2 is wrong for SPECT - // search the gamma from which this hit comes --> photonID - photonID = m_trajectoryNavigator->FindPhotonID(trackID); - if (nVerboseLevel > 2 && photonID == rootID) { - G4cout << "GateAnalysis::RecordEndOfEvent : trackID: " << trackID << " photonID = " << rootID << G4endl; - } - } - - if (photonID > 0 && photonID < 4) { - const G4int index = photonID - 1; - nPhantomCompton = photon_scatterings[index].nPhantomCompton; - nPhantomRayleigh = photon_scatterings[index].nPhantomRayleigh; - theComptonVolumeName = photon_scatterings[index].theComptonVolumeName; - theRayleighVolumeName = photon_scatterings[index].theRayleighVolumeName; - nCrystalCompton = photon_scatterings[index].nCrystalCompton; - nCrystalRayleigh = photon_scatterings[index].nCrystalRayleigh; - } - - // search the primary that originated the track - primaryID = m_trajectoryNavigator->FindPrimaryID(trackID); - (*CHC)[iHit]->SetSourceID(sourceID); - (*CHC)[iHit]->SetSourcePosition(sourceVertex); - (*CHC)[iHit]->SetNPhantomCompton(nPhantomCompton); - (*CHC)[iHit]->SetNPhantomRayleigh(nPhantomRayleigh); - (*CHC)[iHit]->SetComptonVolumeName(theComptonVolumeName); - (*CHC)[iHit]->SetRayleighVolumeName(theRayleighVolumeName); - (*CHC)[iHit]->SetPhotonID(photonID); - (*CHC)[iHit]->SetPrimaryID(primaryID); - (*CHC)[iHit]->SetEventID(eventID); - (*CHC)[iHit]->SetRunID(runID); - (*CHC)[iHit]->SetNCrystalCompton(nCrystalCompton); - (*CHC)[iHit]->SetNCrystalRayleigh(nCrystalRayleigh); - (*CHC)[iHit]->SetNSeptal(septalNb); // HDS : septal penetration - } - } + CollectCrystalScatterings(photon_scatterings,septalNb,CHC,eventID); } // end if (CHC) } // end if (!trajectoryContainer) } // end function From a1ead1c0b128f3aa9b354caf7b9f21cef1c48555 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 23:49:17 +0100 Subject: [PATCH 10/15] feat:separation of setting photonIDs from the RecordEndOfEvent method makes this method easier to understand --- source/digits_hits/include/GateAnalysis.hh | 1 + source/digits_hits/src/GateAnalysis.cc | 75 +++++++++++----------- 2 files changed, 39 insertions(+), 37 deletions(-) diff --git a/source/digits_hits/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index 3064c56e0..741746fbb 100644 --- a/source/digits_hits/include/GateAnalysis.hh +++ b/source/digits_hits/include/GateAnalysis.hh @@ -83,6 +83,7 @@ public: const G4ThreeVector& sourceVertex ); void CollectCrystalScatterings(std::vector& photon_scatterings, const G4int septalNb, GateCrystalHitsCollection* CHC, const G4int eventID); + void SetPhotonIDs(std::vector& photon_scatterings); private: diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 56e3f9343..991139f4b 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -298,6 +298,42 @@ void GateAnalysis::CollectCrystalScatterings(std::vector& pho } } +void GateAnalysis::SetPhotonIDs(std::vector& photon_scatterings) { + m_trajectoryNavigator->FindPositronTrackID(); + //search the two gammas + std::vector photonIDVec = m_trajectoryNavigator->FindAnnihilationGammasTrackID(); + if (photonIDVec.size() == 0) { + // no gamma coming from a positron or an ion, or shooted as primary + if (nVerboseLevel > 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonIDs not found" << G4endl; + } + } else { + // This warning is somewhat irrelevant with 124I + if (nVerboseLevel > 0 && photonIDVec.size() > 2) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonID vector size > 2" << G4endl; + } + photon_scatterings[0].photonID = photonIDVec[0]; + photon_scatterings[1].photonID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; + photon_scatterings[2].photonID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; + } + + if (nVerboseLevel > 0) { + if (photon_scatterings[0].photonID == 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon1ID == 0" << G4endl; + } + if (photon_scatterings[1].photonID == 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0" << G4endl; + } + if (photon_scatterings[2].photonID == 0) { + G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon3ID == 0" << G4endl; + } + if (nVerboseLevel > 1) { + G4cout << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon_scatterings[0].photonID; + G4cout << " photon2ID : " << photon_scatterings[1].photonID << G4endl; + } + } +} + //-------------------------------------------------------------------------------------------------- void GateAnalysis::RecordEndOfEvent(const G4Event* event) { @@ -319,46 +355,11 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) GateCrystalHitsCollection* CHC = GetOutputMgr()->GetCrystalHitCollection(); if (CHC) { std::vector photon_scatterings(3); - G4int septalNb = 0; // HDS : septal penetration - - m_trajectoryNavigator->FindPositronTrackID(); - //search the two gammas - std::vector photonIDVec = m_trajectoryNavigator->FindAnnihilationGammasTrackID(); - if (photonIDVec.size() == 0) { - // no gamma coming from a positron or an ion, or shooted as primary - if (nVerboseLevel > 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonIDs not found" << G4endl; - } - } else { - // This warning is somewhat irrelevant with 124I - if (nVerboseLevel > 0 && photonIDVec.size() > 2) { - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonID vector size > 2" << G4endl; - } - photon_scatterings[0].photonID = photonIDVec[0]; - photon_scatterings[1].photonID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; - photon_scatterings[2].photonID = (photonIDVec.size() >= 3) ? photonIDVec[2] : 0; - } - - if (nVerboseLevel > 0) { - if (photon_scatterings[0].photonID == 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon1ID == 0" << G4endl; - } - if (photon_scatterings[1].photonID == 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0" << G4endl; - } - if (photon_scatterings[2].photonID == 0) { - G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photon3ID == 0" << G4endl; - } - if (nVerboseLevel > 1) { - G4cout << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon_scatterings[0].photonID; - G4cout << " photon2ID : " << photon_scatterings[1].photonID << G4endl; - } - } - + SetPhotonIDs(photon_scatterings); // analysis of the phantom hits to count the comptons, etc. + G4int septalNb = 0; // HDS : septal penetration CollectPhantomScatterings(photon_scatterings, septalNb); MakeComptonRayleighDataUpdates(photon_scatterings); - // Source info // DS : if gate source is not used (with /gate/EnableGeneralParticleSource) there are no GateSource, so skip if ((GateSourceMgr::GetInstance())->GetSourcesForThisEvent().size() == 0) { From 1840f79ae699adde0b2709166e5738f45bec868d Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Wed, 22 Feb 2023 23:50:44 +0100 Subject: [PATCH 11/15] style:removed commented out code,edited tabulation and new brackets make code more readable --- .../include/GateTrajectoryNavigator.hh | 35 +- .../src/GateTrajectoryNavigator.cc | 354 ++++++++---------- 2 files changed, 165 insertions(+), 224 deletions(-) diff --git a/source/digits_hits/include/GateTrajectoryNavigator.hh b/source/digits_hits/include/GateTrajectoryNavigator.hh index b062fce25..bf9f5915c 100644 --- a/source/digits_hits/include/GateTrajectoryNavigator.hh +++ b/source/digits_hits/include/GateTrajectoryNavigator.hh @@ -27,43 +27,42 @@ public: virtual ~GateTrajectoryNavigator(); - G4int FindSourceIndex(); + G4int FindSourceIndex(); G4ThreeVector FindSourcePosition(); - G4int FindPositronTrackID(); + G4int FindPositronTrackID(); std::vector FindAnnihilationGammasTrackID(); - G4int FindPhotonID(G4int trackID); + G4int FindPhotonID(G4int trackID); - G4int FindPrimaryID(G4int trackID); + G4int FindPrimaryID(G4int trackID); - void Initialize(); + void Initialize(); - void SetTrajectoryContainer(G4TrajectoryContainer* trajectoryContainer); + void SetTrajectoryContainer(G4TrajectoryContainer* trajectoryContainer); - inline G4TrajectoryContainer* GetTrajectoryContainer() { return m_trajectoryContainer; }; + inline G4TrajectoryContainer* GetTrajectoryContainer() { return m_trajectoryContainer; }; - inline std::vector GetPhotonIDVec() { return m_photonIDVec; }; + inline std::vector GetPhotonIDVec() { return m_photonIDVec; }; - inline G4int GetPositronTrackID() { return m_positronTrackID; }; + inline G4int GetPositronTrackID() { return m_positronTrackID; }; - void SetIonID(); + void SetIonID(); - void SetVerboseLevel(G4int val) { nVerboseLevel = val; }; + void SetVerboseLevel(G4int val) { nVerboseLevel = val; }; protected: private: - G4TrajectoryContainer* m_trajectoryContainer; + G4TrajectoryContainer* m_trajectoryContainer = NULL; - std::vector m_photonIDVec; - G4int m_positronTrackID; - G4Trajectory* m_positronTrj; - G4int m_ionID; - - G4int nVerboseLevel; + std::vector m_photonIDVec; + G4int m_positronTrackID = 0; + G4Trajectory* m_positronTrj = NULL; + G4int m_ionID = 0; + G4int nVerboseLevel = 0; }; diff --git a/source/digits_hits/src/GateTrajectoryNavigator.cc b/source/digits_hits/src/GateTrajectoryNavigator.cc index 4ecdaf4ad..c4716c827 100644 --- a/source/digits_hits/src/GateTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateTrajectoryNavigator.cc @@ -19,25 +19,22 @@ #include "GateMessageManager.hh" #include "GateRunManager.hh" -GateTrajectoryNavigator::GateTrajectoryNavigator() : m_trajectoryContainer(NULL), m_positronTrackID(0), m_positronTrj(NULL), m_ionID(0), nVerboseLevel(0) -{ -} +GateTrajectoryNavigator::GateTrajectoryNavigator() {} -GateTrajectoryNavigator::~GateTrajectoryNavigator() -{ -} +GateTrajectoryNavigator::~GateTrajectoryNavigator() {} G4ThreeVector GateTrajectoryNavigator::FindSourcePosition() { - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::FindSourcePosition\n"; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::FindSourcePosition" << G4endl; + } G4ThreeVector sourcePosition = G4ThreeVector(); G4int sourceIndex = FindSourceIndex(); if ((sourceIndex < 0) || ((unsigned int)sourceIndex >= m_trajectoryContainer->entries())) { - G4cout << "GateTrajectoryNavigator::FindSourcePosition : WARNING : sourceIndex out of range: " << sourceIndex << Gateendl; + G4cout << "GateTrajectoryNavigator::FindSourcePosition : WARNING : sourceIndex out of range: " << sourceIndex << G4endl; } else { G4Trajectory* trjSource = (G4Trajectory*)((*m_trajectoryContainer)[sourceIndex]); sourcePosition = ((G4TrajectoryPoint*)(trjSource->GetPoint(0)))->GetPosition(); @@ -48,8 +45,9 @@ G4ThreeVector GateTrajectoryNavigator::FindSourcePosition() G4int GateTrajectoryNavigator::FindSourceIndex() { - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::FindSourceIndex\n"; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::FindSourceIndex"<< G4endl; + } G4int sourceIndex = -1; G4bool found = false; G4int n_trajectories = m_trajectoryContainer->entries(); @@ -74,10 +72,6 @@ void GateTrajectoryNavigator::SetIonID() for (G4int iTrj=0; iTrjGetCharge() > 2) m_ionID = 1; - - /* G4cout << "GateTrajectoryNavigator::SetIonID:GetTrackID " << trj->GetTrackID() - << " Name <" << trj->GetParticleName() << ">\n"; - G4cout << "GateTrajectoryNavigator::SetIonID:GetCharge " << trj->GetCharge() << Gateendl;*/ } } @@ -86,85 +80,64 @@ void GateTrajectoryNavigator::SetIonID() G4int GateTrajectoryNavigator::FindPositronTrackID() { - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::FindPositronTrackID\n"; - + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::FindPositronTrackID" << G4endl; + } m_positronTrackID = 0; - SetIonID(); - - if (!m_trajectoryContainer) { - G4cout << "GateTrajectoryNavigator::FindPositronTrackID : ERROR : NULL trajectoryContainer\n"; + G4cout << "GateTrajectoryNavigator::FindPositronTrackID : ERROR : NULL trajectoryContainer" << G4endl; } else { G4int n_trajectories = m_trajectoryContainer->entries(); G4bool found = false; for (G4int iTrj=0; (iTrjGetParentID() == m_ionID)&&(trj->GetPDGEncoding() == -11)) { // -11 == Positron - - - m_positronTrackID = trj->GetTrackID(); - m_positronTrj = trj; - found = true; + //PDGEncoding = -11 is a positron + if ((trj->GetParentID() == m_ionID)&&(trj->GetPDGEncoding() == -11)) { + m_positronTrackID = trj->GetTrackID(); + m_positronTrj = trj; + found = true; } } } - return m_positronTrackID; } /* modifs pour cas du mode detector : PY Descourt 08/09/2009 */ -std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() -{ - +std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() { TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID\n"; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID" << G4endl; + } std::vector photonIndices; photonIndices.clear(); - SetIonID(); if (!m_trajectoryContainer) { - G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID : ERROR : NULL trajectoryContainer\n"; + G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID : ERROR : NULL trajectoryContainer" << G4endl; } else { - // prepare the list of gammas for later analysis G4int n_trajectories = m_trajectoryContainer->entries(); - for (G4int iTrj=0; iTrjGetParentID() == m_positronTrackID) - && (trj->GetPDGEncoding() == 22) ) { // 22 == Gamma - photonIndices.push_back(iTrj); - } + // in case the positronTrackID has been found, we put in the list only + // the gammas generated by the positron + if ((trj->GetParentID() == m_positronTrackID) && (trj->GetPDGEncoding() == 22) ) { // 22 == Gamma + photonIndices.push_back(iTrj); + } } else { - // in case the positron has not been found, we accept all the photons - // either coming from the ion (assuming m_ionID==1) or shooted as primary - // (both single and back-to-back pair) - //if (((trj->GetParentID() != 0) ||(trj->GetParentID() == 0) || (trj->GetParentID() == m_ionID)) - // && (trj->GetPDGEncoding() == 22)) - if ((trj->GetParentID() >= 0) && (trj->GetPDGEncoding() == 22)) - { // 22 == Gamma - photonIndices.push_back(iTrj); - - } + // in case the positron has not been found, we accept all the photons + // either coming from the ion (assuming m_ionID==1) or shooted as primary + // (both single and back-to-back pair) + if ((trj->GetParentID() >= 0) && (trj->GetPDGEncoding() == 22)) { // 22 == Gamma + photonIndices.push_back(iTrj); + } } } - // we start the analysis of the gammas in the list // in both cases (in case a positron has been found or not) we look for @@ -185,127 +158,98 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() G4int i2; G4int i3; bool only2gamma = true; - + for (G4int j1=0; j1= 0) && (i2 >= 0) && (i3 >= 0)) { - // both gammas were not already taken - G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); - G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); - G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); - - G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); - G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); - G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); - // in detector mode the vertex position is stored at the last trajectory point - // not the first one - if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh - { + for (G4int j2=j1+1; j2= 0) && (i2 >= 0) && (i3 >= 0)) { + // both gammas were not already taken + G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); + G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); + G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); + G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); + G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); + G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); + // in detector mode the vertex position is stored at the last trajectory point + // not the first one + if ( theMode == TrackingMode::kDetector ) { + // in tracker mode we store the infos about the number of compton and rayleigh G4int n_points =trj1->GetPointEntries(); vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); n_points =trj2->GetPointEntries(); vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); - n_points =trj3->GetPointEntries(); - vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj3->GetPointEntries(); + vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint( n_points - 1 )))->GetPosition(); } - - G4double dist = (vert1-vert2).mag(); - if (nVerboseLevel > 2) - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << Gateendl; - if (dist/mm < 1E-7) { - if (nVerboseLevel > 1) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :" - << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << Gateendl; + G4double dist = (vert1-vert2).mag(); + if (nVerboseLevel > 2) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << G4endl; } - // we add both photons to the vertex - m_photonIDVec.push_back(trj1->GetTrackID()); - m_photonIDVec.push_back(trj2->GetTrackID()); - if ((vert2-vert3).mag()/mm < 1E-7) - { - m_photonIDVec.push_back(trj3->GetTrackID()); - photonIndices[j3] = -1; - only2gamma = false; - // std::cout << "\n\n\n\n\n\n\n\n\n" << photonIndices[j3] << "\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n" << std::endl; - } - // we cancel the 2nd photon from the list to avoid double counting later - photonIndices[j2] = -1; - } - } - } - - } - } - - if(only2gamma){ - for (G4int j1=0; j1= 0) && (i2 >= 0)) { - // both gammas were not already taken - G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); - G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); - - G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); - G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); - // in detector mode the vertex position is stored at the last trajectory point - // not the first one - if ( theMode == TrackingMode::kDetector ) // in tracker mode we store the infos about the number of compton and rayleigh - { - G4int n_points =trj1->GetPointEntries(); - vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); - n_points =trj2->GetPointEntries(); - vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); + if (dist/mm < 1E-7) { + if (nVerboseLevel > 1) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; + G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; + } + // we add both photons to the vertex + m_photonIDVec.push_back(trj1->GetTrackID()); + m_photonIDVec.push_back(trj2->GetTrackID()); + if ((vert2-vert3).mag()/mm < 1E-7) { + m_photonIDVec.push_back(trj3->GetTrackID()); + photonIndices[j3] = -1; + only2gamma = false; + } + // we cancel the 2nd photon from the list to avoid double counting later + photonIndices[j2] = -1; + } + } + } } + } + if(only2gamma){ + for (G4int j1=0; j1= 0) && (i2 >= 0)) { + // both gammas were not already taken + G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); + G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); + G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); + G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); + // in detector mode the vertex position is stored at the last trajectory point + // not the first one + if ( theMode == TrackingMode::kDetector ) { + // in tracker mode we store the infos about the number of compton and rayleigh + G4int n_points =trj1->GetPointEntries(); + vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj2->GetPointEntries(); + vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); + } - G4double dist = (vert1-vert2).mag(); - if (nVerboseLevel > 2) - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << Gateendl; - if (dist/mm < 1E-7) { - if (nVerboseLevel > 1) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :" - << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << Gateendl; - } - // we add both photons to the vertex - m_photonIDVec.push_back(trj1->GetTrackID()); - m_photonIDVec.push_back(trj2->GetTrackID()); - // we cancel the 2nd photon from the list to avoid double counting later - photonIndices[j2] = -1; - } + G4double dist = (vert1-vert2).mag(); + if (nVerboseLevel > 2) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << G4endl; + } + if (dist/mm < 1E-7) { + if (nVerboseLevel > 1) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; + G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; + } + // we add both photons to the vertex + m_photonIDVec.push_back(trj1->GetTrackID()); + m_photonIDVec.push_back(trj2->GetTrackID()); + // we cancel the 2nd photon from the list to avoid double counting later + photonIndices[j2] = -1; + } + } + } + } + } } } -} - } } - } - - /*if ((m_positronTrackID != 0) && ((m_photonIDVec.size() > 2) || (m_photonIDVec.size() == 1))) { - // strange case: we would expect to see a true annihilation, thus only 2 gammas backto back from the positron - // 0 can be normal (positron escaped and didn't annihilate), 1 is strange again. - G4cout - << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] WARNING : positron found and photonID vector size: " << m_photonIDVec.size() << Gateendl; - } */ - - /*if ((m_positronTrackID != 0) && (m_photonIDVec.size() == 0)) { - // still, it's not really what we want to see the positron going out of the scene without doing anything... - // if (nVerboseLevel > 0) - if (nVerboseLevel > 0) - G4cout - << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] WARNING : positron found and photonID vector size: " << m_photonIDVec.size() << Gateendl; - G4Navigator *gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); - G4ThreeVector null(0.,0.,0.); - G4ThreeVector *ptr; - ptr = &null; - G4ThreeVector lastPoint = ((G4TrajectoryPoint*)(m_positronTrj->GetPoint(m_positronTrj->GetPointEntries()-1)))->GetPosition(); - G4String theLastVolumeName = gNavigator->LocateGlobalPointAndSetup(lastPoint,ptr,false)->GetName(); - if (nVerboseLevel > 0) - G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID: last position of the positron track has been found in volume: " << theLastVolumeName << Gateendl; - }*/ - - return m_photonIDVec; } @@ -313,18 +257,17 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() // search the gamma from which this track comes --> photonID G4int GateTrajectoryNavigator::FindPhotonID(G4int trackID) { - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::FindPhotonID \n"; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::FindPhotonID " << G4endl; + } G4int photonID = 0; if (!m_trajectoryContainer) { - G4cout << "GateTrajectoryNavigator::FindPhotonID : ERROR : NULL trajectoryContainer\n"; + G4cout << "GateTrajectoryNavigator::FindPhotonID : ERROR : NULL trajectoryContainer" << G4endl; } else { - if (m_photonIDVec.size() == 0) { - G4cout << "GateTrajectoryNavigator::FindPhotonID : m_photonIDVec.size() == 0\n"; + G4cout << "GateTrajectoryNavigator::FindPhotonID : m_photonIDVec.size() == 0" << G4endl; } else { - G4bool found; // search the gamma related to this hit --> photonID photonID = trackID; @@ -333,32 +276,31 @@ G4int GateTrajectoryNavigator::FindPhotonID(G4int trackID) G4int photon2ID = (m_photonIDVec.size() >= 2) ? m_photonIDVec[1] : 0; G4int photon3ID = (m_photonIDVec.size() >= 3) ? m_photonIDVec[2] : 0; G4int rootID = 0; - // we go up and up, starting from the present trackID, to the parentID, the parentID, ecc until // we find that the ID of the track is equal to the ID of: one of the photons, or rootID(==0) while (!((photonID==photon1ID)||(photonID==photon2ID)||(photonID==photon3ID)||(photonID==rootID))) { - found = false; - G4int n_trajectories = m_trajectoryContainer->entries(); - for (G4int iTrj=0; (iTrjGetTrackID()) { - photonID = trj->GetParentID(); - found = true; - } - } + found = false; + G4int n_trajectories = m_trajectoryContainer->entries(); + for (G4int iTrj=0; (iTrjGetTrackID()) { + photonID = trj->GetParentID(); + found = true; + } + } } if (photonID == rootID) { - if (nVerboseLevel > 2) G4cout - << "GateTrajectoryNavigator::FindPhotonID : trackID: " << trackID << " photonID = " << rootID << Gateendl; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::FindPhotonID : trackID: " << trackID << " photonID = " << rootID << G4endl; + } } if (photonID == photon1ID) { - photonID = 1; + photonID = 1; } else if (photonID == photon2ID) { - photonID = 2; + photonID = 2; } else if (photonID == photon3ID) { - photonID = 0; + photonID = 0; } - } } @@ -371,10 +313,9 @@ G4int GateTrajectoryNavigator::FindPrimaryID(G4int trackID) G4int primaryID = 0; if (!m_trajectoryContainer) { - G4cout << "GateTrajectoryNavigator::FindPrimaryID : ERROR : NULL trajectoryContainer\n"; + G4cout << "GateTrajectoryNavigator::FindPrimaryID : ERROR : NULL trajectoryContainer" << G4endl; } else { G4int tempParentID = trackID; - //G4int primaryIndex = -1; G4bool found; G4int n_trajectories = m_trajectoryContainer->entries(); // we go up and up starting from the trackID, via the parentID's, until @@ -383,12 +324,11 @@ G4int GateTrajectoryNavigator::FindPrimaryID(G4int trackID) primaryID = tempParentID; found = false; for (G4int iTrj=0; (iTrjGetTrackID()) { - tempParentID = trj->GetParentID(); - //primaryIndex = iTrj; - found = true; - } + G4Trajectory* trj = (G4Trajectory*)((*m_trajectoryContainer)[iTrj]); + if (primaryID == trj->GetTrackID()) { + tempParentID = trj->GetParentID(); + found = true; + } } } while (tempParentID != 0); } @@ -399,11 +339,12 @@ G4int GateTrajectoryNavigator::FindPrimaryID(G4int trackID) void GateTrajectoryNavigator::SetTrajectoryContainer(G4TrajectoryContainer* trajectoryContainer) { - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::SetTrajectoryContainer\n"; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::SetTrajectoryContainer" << G4endl; + } if (!trajectoryContainer) { - G4cout << "GateTrajectoryNavigator::SetTrajectoryContainer : ERROR : NULL trajectoryContainer\n"; + G4cout << "GateTrajectoryNavigator::SetTrajectoryContainer : ERROR : NULL trajectoryContainer" << G4endl; } else { m_trajectoryContainer = trajectoryContainer; } @@ -413,8 +354,9 @@ void GateTrajectoryNavigator::SetTrajectoryContainer(G4TrajectoryContainer* traj void GateTrajectoryNavigator::Initialize() { - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::Initialize\n"; + if (nVerboseLevel > 2) { + G4cout << "GateTrajectoryNavigator::Initialize" << G4endl; + } m_photonIDVec.clear(); m_positronTrackID = -1; } From 18b1d86af872ceee8febf50936d42809ce7cfb6b Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Thu, 23 Feb 2023 00:54:37 +0100 Subject: [PATCH 12/15] feat:allows to determine whether we are working with the positronium model and which --- .../include/GatePositroniumDecayInfo.hh | 28 +++++++ .../src/GatePositroniumDecayInfo.cc | 83 +++++++++++++++++++ 2 files changed, 111 insertions(+) create mode 100644 source/digits_hits/include/GatePositroniumDecayInfo.hh create mode 100644 source/digits_hits/src/GatePositroniumDecayInfo.cc diff --git a/source/digits_hits/include/GatePositroniumDecayInfo.hh b/source/digits_hits/include/GatePositroniumDecayInfo.hh new file mode 100644 index 000000000..b6026ce26 --- /dev/null +++ b/source/digits_hits/include/GatePositroniumDecayInfo.hh @@ -0,0 +1,28 @@ +#ifndef POSITRONIUM_DECAY_INFO_HH +#define POSITRONIUM_DECAY_INFO_HH +#include "GateCrystalHit.hh" +#include + +namespace GatePositroniumDecay { + enum class GammaType {Unknown = 0,SinglePhoton = 1,AnnihilationPhoton = 2,PromptPhoton = 3}; + enum class SourceType {Unknown = 0, SingleGamaEmitter = 1, ParaPositronium = 2, OrthoPositronium = 3 }; + enum class DecayType {Unknown = 0, StandardChannel = 1, DeexciatationChannel = 2}; + + GammaType IntToGammaType(const G4int value); + SourceType IntToSourceType(const G4int value); + DecayType IntToDecayType(const G4int value); + G4int GetMaxPhotonsNumber(const SourceType source_type,const DecayType decay_type); + + class GatePositroniumDecayInfo { + public: + void Fill(const GateCrystalHitsCollection* collection); + bool HasGamma(const GammaType gamma_type) const; + bool HasSource(const SourceType source_type) const; + bool HasDecay(const DecayType decay_type) const; + private: + std::set gammaTypes; + std::set sourceTypes; + std::set decayTypes; + }; +} +#endif \ No newline at end of file diff --git a/source/digits_hits/src/GatePositroniumDecayInfo.cc b/source/digits_hits/src/GatePositroniumDecayInfo.cc new file mode 100644 index 000000000..eda6e26c2 --- /dev/null +++ b/source/digits_hits/src/GatePositroniumDecayInfo.cc @@ -0,0 +1,83 @@ +#include "GatePositroniumDecayInfo.hh" +namespace GatePositroniumDecay { + GammaType IntToGammaType(const G4int value) { + switch (value) { + case 0: + return GammaType::Unknown; + case 1: + return GammaType::SinglePhoton; + case 2: + return GammaType::AnnihilationPhoton; + case 3: + return GammaType::PromptPhoton; + default: + return GammaType::Unknown; + } + } + + SourceType IntToSourceType(const G4int value) { + switch (value) { + case 0: + return SourceType::Unknown; + case 1: + return SourceType::SingleGamaEmitter; + case 2: + return SourceType::ParaPositronium; + case 3: + return SourceType::OrthoPositronium; + default: + return SourceType::Unknown; + } + } + + DecayType IntToDecayType(const G4int value) { + switch (value) { + case 0: + return DecayType::Unknown; + case 1: + return DecayType::StandardChannel; + case 2: + return DecayType::DeexciatationChannel; + default: + return DecayType::Unknown; + } + } + + G4int GetMaxPhotonsNumber(const SourceType source_type,const DecayType decay_type) { + bool unknown_model = source_type == SourceType::Unknown || decay_type == DecayType::Unknown; + if (unknown_model){ + return 0; + } + if (source_type == SourceType::SingleGamaEmitter) { + return 1; + } + if (source_type == SourceType::ParaPositronium) { + return decay_type == DecayType::DeexciatationChannel ? 3 : 2; + } + //OrthoPositronium source + return decay_type == DecayType::DeexciatationChannel ? 4 : 3; + } + + void GatePositroniumDecayInfo::Fill(const GateCrystalHitsCollection* collection) { + const G4int n_hits = collection->GetSize(); + for (G4int i = 0; i < n_hits; ++i) { + const GateCrystalHit* hit = (*collection)[i]; + gammaTypes.emplace(IntToGammaType(hit->GetGammaType())); + sourceTypes.emplace(IntToSourceType(hit->GetSourceType())); + decayTypes.emplace(IntToDecayType(hit->GetDecayType())); + } + } + + bool GatePositroniumDecayInfo::HasGamma(const GammaType gamma_type) const { + return gammaTypes.find(gamma_type) != gammaTypes.cend(); + } + + bool GatePositroniumDecayInfo::HasSource(const SourceType source_type) const { + return sourceTypes.find(source_type) != sourceTypes.cend(); + } + + bool GatePositroniumDecayInfo::HasDecay(const DecayType decay_type) const { + return decayTypes.find(decay_type) != decayTypes.cend(); + } + +} \ No newline at end of file From 8dbca0fa17ae0c2b9c6f985b81e418a14c17c448 Mon Sep 17 00:00:00 2001 From: BlurredChoise Date: Fri, 3 Mar 2023 05:32:45 +0100 Subject: [PATCH 13/15] feat:methods to find photons IDs are in the separated methods --- .../include/GateTrajectoryNavigator.hh | 4 + .../src/GateTrajectoryNavigator.cc | 245 +++++++++--------- 2 files changed, 133 insertions(+), 116 deletions(-) diff --git a/source/digits_hits/include/GateTrajectoryNavigator.hh b/source/digits_hits/include/GateTrajectoryNavigator.hh index bf9f5915c..4737995fa 100644 --- a/source/digits_hits/include/GateTrajectoryNavigator.hh +++ b/source/digits_hits/include/GateTrajectoryNavigator.hh @@ -54,6 +54,10 @@ public: void SetVerboseLevel(G4int val) { nVerboseLevel = val; }; protected: + + std::vector GetPhotonIndices(); + void FillPhotonIDsForThreePhotons(std::vector& photonIndices, bool& only2gamma); + void FillPhotonIDsForTwoPhotons(std::vector& photonIndices); private: G4TrajectoryContainer* m_trajectoryContainer = NULL; diff --git a/source/digits_hits/src/GateTrajectoryNavigator.cc b/source/digits_hits/src/GateTrajectoryNavigator.cc index c4716c827..faf2f1d9d 100644 --- a/source/digits_hits/src/GateTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateTrajectoryNavigator.cc @@ -103,41 +103,143 @@ G4int GateTrajectoryNavigator::FindPositronTrackID() return m_positronTrackID; } +std::vector GateTrajectoryNavigator::GetPhotonIndices() { + std::vector photonIndices; + G4int n_trajectories = m_trajectoryContainer->entries(); + for (G4int iTrj=0; iTrjGetParentID() == m_positronTrackID) && (trj->GetPDGEncoding() == 22) ) { // 22 == Gamma + photonIndices.push_back(iTrj); + } + } else { + // in case the positron has not been found, we accept all the photons + // either coming from the ion (assuming m_ionID==1) or shooted as primary + // (both single and back-to-back pair) + if ((trj->GetParentID() >= 0) && (trj->GetPDGEncoding() == 22)) { // 22 == Gamma + photonIndices.push_back(iTrj); + } + } + } + return photonIndices; +} -/* modifs pour cas du mode detector : PY Descourt 08/09/2009 */ -std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() { +void GateTrajectoryNavigator::FillPhotonIDsForThreePhotons(std::vector& photonIndices, bool& only2gamma) { TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); + const G4int nPh = photonIndices.size(); + G4int i1; + G4int i2; + G4int i3; + for (G4int j1=0; j1= 0) && (i2 >= 0) && (i3 >= 0)) { + // both gammas were not already taken + G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); + G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); + G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); + G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); + G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); + G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); + // in detector mode the vertex position is stored at the last trajectory point + // not the first one + if ( theMode == TrackingMode::kDetector ) { + // in tracker mode we store the infos about the number of compton and rayleigh + G4int n_points =trj1->GetPointEntries(); + vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj2->GetPointEntries(); + vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj3->GetPointEntries(); + vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint( n_points - 1 )))->GetPosition(); + } + G4double dist = (vert1-vert2).mag(); + if (nVerboseLevel > 2) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << G4endl; + } + if (dist/mm < 1E-7) { + if (nVerboseLevel > 1) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; + G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; + } + // we add both photons to the vertex + m_photonIDVec.push_back(trj1->GetTrackID()); + m_photonIDVec.push_back(trj2->GetTrackID()); + if ((vert2-vert3).mag()/mm < 1E-7) { + m_photonIDVec.push_back(trj3->GetTrackID()); + photonIndices[j3] = -1; + only2gamma = false; + } + // we cancel the 2nd photon from the list to avoid double counting later + photonIndices[j2] = -1; + } + } + } + } + } +} + +void GateTrajectoryNavigator::FillPhotonIDsForTwoPhotons(std::vector& photonIndices) { + TrackingMode theMode =( (GateSteppingAction *)(GateRunManager::GetRunManager()->GetUserSteppingAction() ) )->GetMode(); + const G4int nPh = photonIndices.size(); + G4int i1; + G4int i2; + for (G4int j1=0; j1= 0) && (i2 >= 0)) { + // both gammas were not already taken + G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); + G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); + G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); + G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); + // in detector mode the vertex position is stored at the last trajectory point + // not the first one + if ( theMode == TrackingMode::kDetector ) { + // in tracker mode we store the infos about the number of compton and rayleigh + G4int n_points =trj1->GetPointEntries(); + vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); + n_points =trj2->GetPointEntries(); + vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); + } + G4double dist = (vert1-vert2).mag(); + if (nVerboseLevel > 2) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << G4endl; + } + if (dist/mm < 1E-7) { + if (nVerboseLevel > 1) { + G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; + G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; + } + // we add both photons to the vertex + m_photonIDVec.push_back(trj1->GetTrackID()); + m_photonIDVec.push_back(trj2->GetTrackID()); + // we cancel the 2nd photon from the list to avoid double counting later + photonIndices[j2] = -1; + } + } + } + } +} + +/* modifs pour cas du mode detector : PY Descourt 08/09/2009 */ +std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() { if (nVerboseLevel > 2) { G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID" << G4endl; } - - std::vector photonIndices; - photonIndices.clear(); SetIonID(); if (!m_trajectoryContainer) { G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID : ERROR : NULL trajectoryContainer" << G4endl; } else { // prepare the list of gammas for later analysis - G4int n_trajectories = m_trajectoryContainer->entries(); - for (G4int iTrj=0; iTrjGetParentID() == m_positronTrackID) && (trj->GetPDGEncoding() == 22) ) { // 22 == Gamma - photonIndices.push_back(iTrj); - } - } else { - // in case the positron has not been found, we accept all the photons - // either coming from the ion (assuming m_ionID==1) or shooted as primary - // (both single and back-to-back pair) - if ((trj->GetParentID() >= 0) && (trj->GetPDGEncoding() == 22)) { // 22 == Gamma - photonIndices.push_back(iTrj); - } - } - } + std::vector photonIndices = GetPhotonIndices(); // we start the analysis of the gammas in the list // in both cases (in case a positron has been found or not) we look for @@ -154,100 +256,11 @@ std::vector GateTrajectoryNavigator::FindAnnihilationGammasTrackID() { // if more than 1 gamma, we select those coming from a common vertex // this design should be open to more than 2 gammas, even if for the moment not // considered by GateAnalysis - G4int i1; - G4int i2; - G4int i3; bool only2gamma = true; - - for (G4int j1=0; j1= 0) && (i2 >= 0) && (i3 >= 0)) { - // both gammas were not already taken - G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); - G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); - G4Trajectory* trj3 = (G4Trajectory*)((*m_trajectoryContainer)[i3]); - G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); - G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); - G4ThreeVector vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint(0)))->GetPosition(); - // in detector mode the vertex position is stored at the last trajectory point - // not the first one - if ( theMode == TrackingMode::kDetector ) { - // in tracker mode we store the infos about the number of compton and rayleigh - G4int n_points =trj1->GetPointEntries(); - vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); - n_points =trj2->GetPointEntries(); - vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); - n_points =trj3->GetPointEntries(); - vert3 = ((G4TrajectoryPoint*)(trj3->GetPoint( n_points - 1 )))->GetPosition(); - } - G4double dist = (vert1-vert2).mag(); - if (nVerboseLevel > 2) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << G4endl; - } - if (dist/mm < 1E-7) { - if (nVerboseLevel > 1) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; - G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; - } - // we add both photons to the vertex - m_photonIDVec.push_back(trj1->GetTrackID()); - m_photonIDVec.push_back(trj2->GetTrackID()); - if ((vert2-vert3).mag()/mm < 1E-7) { - m_photonIDVec.push_back(trj3->GetTrackID()); - photonIndices[j3] = -1; - only2gamma = false; - } - // we cancel the 2nd photon from the list to avoid double counting later - photonIndices[j2] = -1; - } - } - } - } - } - if(only2gamma){ - for (G4int j1=0; j1= 0) && (i2 >= 0)) { - // both gammas were not already taken - G4Trajectory* trj1 = (G4Trajectory*)((*m_trajectoryContainer)[i1]); - G4Trajectory* trj2 = (G4Trajectory*)((*m_trajectoryContainer)[i2]); - G4ThreeVector vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint(0)))->GetPosition(); - G4ThreeVector vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint(0)))->GetPosition(); - // in detector mode the vertex position is stored at the last trajectory point - // not the first one - if ( theMode == TrackingMode::kDetector ) { - // in tracker mode we store the infos about the number of compton and rayleigh - G4int n_points =trj1->GetPointEntries(); - vert1 = ((G4TrajectoryPoint*)(trj1->GetPoint( n_points - 1 )))->GetPosition(); - n_points =trj2->GetPointEntries(); - vert2 = ((G4TrajectoryPoint*)(trj2->GetPoint( n_points - 1 )))->GetPosition(); - } - - G4double dist = (vert1-vert2).mag(); - if (nVerboseLevel > 2) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : distance between gammas vertices : dist (mm) " << dist/mm << G4endl; - } - if (dist/mm < 1E-7) { - if (nVerboseLevel > 1) { - G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; - G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; - } - // we add both photons to the vertex - m_photonIDVec.push_back(trj1->GetTrackID()); - m_photonIDVec.push_back(trj2->GetTrackID()); - // we cancel the 2nd photon from the list to avoid double counting later - photonIndices[j2] = -1; - } - } - } - } - } + FillPhotonIDsForThreePhotons(photonIndices,only2gamma); + if (only2gamma) { + FillPhotonIDsForTwoPhotons(photonIndices); + } } } return m_photonIDVec; From 7f739abfd75ad7cc2688fa003323578769ee96ff Mon Sep 17 00:00:00 2001 From: sparzych Date: Wed, 12 Jul 2023 13:27:37 +0200 Subject: [PATCH 14/15] fixing the photonID saved in output to separate prompt and annihilation gamma --- source/digits_hits/include/GateAnalysis.hh | 3 ++- source/digits_hits/src/GateAnalysis.cc | 13 ++++++++++--- source/digits_hits/src/GateTrajectoryNavigator.cc | 3 ++- 3 files changed, 14 insertions(+), 5 deletions(-) diff --git a/source/digits_hits/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index 741746fbb..9971f3512 100644 --- a/source/digits_hits/include/GateAnalysis.hh +++ b/source/digits_hits/include/GateAnalysis.hh @@ -80,7 +80,8 @@ public: const G4int sourceID, const G4int eventID, const G4int runID, - const G4ThreeVector& sourceVertex + const G4ThreeVector& sourceVertex, + const G4int gammaType ); void CollectCrystalScatterings(std::vector& photon_scatterings, const G4int septalNb, GateCrystalHitsCollection* CHC, const G4int eventID); void SetPhotonIDs(std::vector& photon_scatterings); diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 991139f4b..441d8b5c3 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -230,7 +230,8 @@ void GateAnalysis::UpdateHitDataForAnalysis( const G4int sourceID, const G4int eventID, const G4int runID, - const G4ThreeVector& sourceVertex + const G4ThreeVector& sourceVertex, + const G4int gammaType ) { // fill in values with the branch with C struct G4int trackID = hit->GetTrackID(); @@ -262,7 +263,12 @@ void GateAnalysis::UpdateHitDataForAnalysis( nCrystalCompton = photon_scatterings[index].nCrystalCompton; nCrystalRayleigh = photon_scatterings[index].nCrystalRayleigh; } - + + // if gammaType is PromptPhoton the photonID for output file is set to 0 + if (gammaType==3) { + photonID = 0; + } + // search the primary that originated the track hit->SetSourceID(sourceID); hit->SetSourcePosition(sourceVertex); @@ -292,8 +298,9 @@ void GateAnalysis::CollectCrystalScatterings(std::vector& pho G4cout << "GateAnalysis::RecordEndOfEvent : CrystalHitsCollection: processName : <" << hit->GetProcess() << G4endl; G4cout << "> Particls PDG code : " << hit->GetPDGEncoding() << G4endl; } + G4int gammaType = hit->GetGammaType(); if (hit->GoodForAnalysis()) { - UpdateHitDataForAnalysis(photon_scatterings,hit, septalNb, sourceID, eventID,runID,sourceVertex); + UpdateHitDataForAnalysis(photon_scatterings, hit, septalNb, sourceID, eventID, runID, sourceVertex, gammaType); } } } diff --git a/source/digits_hits/src/GateTrajectoryNavigator.cc b/source/digits_hits/src/GateTrajectoryNavigator.cc index faf2f1d9d..09acfb51b 100644 --- a/source/digits_hits/src/GateTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateTrajectoryNavigator.cc @@ -169,6 +169,7 @@ void GateTrajectoryNavigator::FillPhotonIDsForThreePhotons(std::vector& p // we add both photons to the vertex m_photonIDVec.push_back(trj1->GetTrackID()); m_photonIDVec.push_back(trj2->GetTrackID()); + if ((vert2-vert3).mag()/mm < 1E-7) { m_photonIDVec.push_back(trj3->GetTrackID()); photonIndices[j3] = -1; @@ -312,7 +313,7 @@ G4int GateTrajectoryNavigator::FindPhotonID(G4int trackID) } else if (photonID == photon2ID) { photonID = 2; } else if (photonID == photon3ID) { - photonID = 0; + photonID = 3; } } } From 53a67edc49f2afe6243d475b7b8dcd066ee50ac0 Mon Sep 17 00:00:00 2001 From: sparzych Date: Wed, 12 Jul 2023 13:41:00 +0200 Subject: [PATCH 15/15] fix in GateTrajectoryNavigator to avoid repetition in filling of m_photonIDVec in case of finding 2 photons --- .../digits_hits/src/GateTrajectoryNavigator.cc | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/source/digits_hits/src/GateTrajectoryNavigator.cc b/source/digits_hits/src/GateTrajectoryNavigator.cc index 09acfb51b..9a411f44d 100644 --- a/source/digits_hits/src/GateTrajectoryNavigator.cc +++ b/source/digits_hits/src/GateTrajectoryNavigator.cc @@ -134,6 +134,7 @@ void GateTrajectoryNavigator::FillPhotonIDsForThreePhotons(std::vector& p G4int i3; for (G4int j1=0; j1& p G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common vertex for the two annihilation gammas :"; G4cout << " tracks " << trj1->GetTrackID() << " and " << trj2->GetTrackID() << G4endl; } - // we add both photons to the vertex - m_photonIDVec.push_back(trj1->GetTrackID()); - m_photonIDVec.push_back(trj2->GetTrackID()); if ((vert2-vert3).mag()/mm < 1E-7) { + // we add all three photons to the vertex + m_photonIDVec.push_back(trj1->GetTrackID()); + m_photonIDVec.push_back(trj2->GetTrackID()); m_photonIDVec.push_back(trj3->GetTrackID()); + // we cancel the 3rd photon from the list to avoid double counting later photonIndices[j3] = -1; only2gamma = false; } - // we cancel the 2nd photon from the list to avoid double counting later - photonIndices[j2] = -1; + // we keep in mind that 2nd photon from the list was used to later cancel it to avoid double counting + photonUsed = true; } } } + if(photonUsed) { + // we cancel the 2nd photon from the list to avoid double counting later + photonIndices[j2] = -1; + } } } }