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/include/GateAnalysis.hh b/source/digits_hits/include/GateAnalysis.hh index b2f8133a3..9971f3512 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; @@ -26,6 +27,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(); @@ -57,6 +68,23 @@ 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); + 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, + const G4int gammaType + ); + 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/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/include/GateToRoot.hh b/source/digits_hits/include/GateToRoot.hh index 530ca17d5..ba237d43c 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/include/GateTrajectoryNavigator.hh b/source/digits_hits/include/GateTrajectoryNavigator.hh index b062fce25..4737995fa 100644 --- a/source/digits_hits/include/GateTrajectoryNavigator.hh +++ b/source/digits_hits/include/GateTrajectoryNavigator.hh @@ -27,43 +27,46 @@ 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: + + std::vector GetPhotonIndices(); + void FillPhotonIDsForThreePhotons(std::vector& photonIndices, bool& only2gamma); + void FillPhotonIDsForTwoPhotons(std::vector& photonIndices); 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/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;} //----------------------------------------------------------------------------- diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 7e0a53715..441d8b5c3 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) @@ -110,356 +111,270 @@ 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::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::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, + const G4int gammaType +) { + // 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; + } + + // 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); + 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; + } + G4int gammaType = hit->GetGammaType(); + if (hit->GoodForAnalysis()) { + UpdateHitDataForAnalysis(photon_scatterings, hit, septalNb, sourceID, eventID, runID, sourceVertex, gammaType); + } + } +} + +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) { + const G4int eventID = event->GetEventID(); 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(); - //G4cout << "GateAnalysis::EventID et RunID : " < 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(); - //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 rootID = 0; - G4int primaryID = 0; - - G4int photon1_phantom_compton = 0; - G4int photon2_phantom_compton = 0; - - G4int photon1_crystal_compton = 0; - G4int photon2_crystal_compton = 0; - - G4int photon1_phantom_Rayleigh = 0; - G4int photon2_phantom_Rayleigh = 0; - - G4int photon1_crystal_Rayleigh = 0; - G4int photon2_crystal_Rayleigh = 0; - - 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(); - 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 - { - // 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; - } - - if (photon1ID == 0) - { - if (nVerboseLevel > 0) G4cout - << "GateAnalysis::RecordEndOfEvent : WARNING : photon1ID == 0\n"; - } - if (photon2ID == 0) { - if (nVerboseLevel > 1) G4cout - << "GateAnalysis::RecordEndOfEvent : WARNING : photon2ID == 0\n"; - } - 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 theRayleighVolumeName("NULL"); - G4String theRayleighVolumeName1("NULL"); - G4String theRayleighVolumeName2("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; - } - } - - // 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; - } - } - } // 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 - { // G4cout << " GateAnalysis eventID "<GetModule("root")); - ComptonRayleighData aCRData; - aCRData.photon1_phantom_Rayleigh = photon1_phantom_Rayleigh; - aCRData.photon2_phantom_Rayleigh = photon2_phantom_Rayleigh; - aCRData.photon1_phantom_compton = photon1_phantom_compton; - aCRData.photon2_phantom_compton = photon2_phantom_compton; - strcpy(aCRData.theComptonVolumeName1 , theComptonVolumeName1.c_str() ); - strcpy(aCRData.theComptonVolumeName2 , theComptonVolumeName2.c_str() ); - strcpy(aCRData.theRayleighVolumeName1 , theRayleighVolumeName1.c_str() ); - strcpy(aCRData.theRayleighVolumeName2 , theRayleighVolumeName2.c_str() ); - gateToRoot->RecordPHData( aCRData ); - // return; - } - - 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; - photon1_phantom_compton += aCRData.photon1_phantom_compton; - photon2_phantom_compton += aCRData.photon2_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; - theRayleighVolumeName1 = aCRData.theRayleighVolumeName1; - theRayleighVolumeName2 = aCRData.theRayleighVolumeName2; - - } - - - ///////// - // 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 || processName.find("Rayleigh") != G4String::npos) { - if (processName.find("ompt") != G4String::npos) - { - - if (crystalTrackID == photon1ID) photon1_crystal_compton++; - if (crystalTrackID == photon2ID) photon2_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++; - } - - 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 parentID = (*CHC)[iHit]->GetParentID(); - - G4int photonID = 0; - G4int nPhantomCompton = 0; - G4int nCrystalCompton = 0; - - 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 - 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; - } - - // 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) + } else { + GateCrystalHitsCollection* CHC = GetOutputMgr()->GetCrystalHitCollection(); + if (CHC) { + std::vector photon_scatterings(3); + 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) { + return; + } + CollectCrystalScatterings(photon_scatterings,septalNb,CHC,eventID); + } // end if (CHC) + } // end if (!trajectoryContainer) } // end function //-------------------------------------------------------------------------------------------------- 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 diff --git a/source/digits_hits/src/GateToRoot.cc b/source/digits_hits/src/GateToRoot.cc index 15a7976aa..c67d8cff6 100644 --- a/source/digits_hits/src/GateToRoot.cc +++ b/source/digits_hits/src/GateToRoot.cc @@ -69,26 +69,34 @@ 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); + 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; + 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); - return *this; + 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; + 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"); //////////////////////// @@ -621,15 +633,16 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) { 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; - 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.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.; @@ -639,12 +652,16 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) { 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; - 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.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("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("RunID").c_str(), &m_RSRunID); // rewind counters @@ -1432,24 +1455,31 @@ 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; - strcpy(theCRData.theComptonVolumeName1, aCRData.theComptonVolumeName1); - strcpy(theCRData.theComptonVolumeName2, aCRData.theComptonVolumeName2); - strcpy(theCRData.theRayleighVolumeName1, aCRData.theRayleighVolumeName1); - strcpy(theCRData.theRayleighVolumeName2, aCRData.theRayleighVolumeName2); + 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; - strcpy(aCRData.theComptonVolumeName1, theCRData.theComptonVolumeName1); - strcpy(aCRData.theComptonVolumeName2, theCRData.theComptonVolumeName2); - strcpy(aCRData.theRayleighVolumeName1, theCRData.theRayleighVolumeName1); - strcpy(aCRData.theRayleighVolumeName2, theCRData.theRayleighVolumeName2); - + 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 115899d58..9a411f44d 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,173 @@ 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; } +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; + } + + 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 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; + } + } + } +} - if (nVerboseLevel > 2) - G4cout << "GateTrajectoryNavigator::FindAnnihilationGammasTrackID\n"; - - std::vector photonIndices; - photonIndices.clear(); +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; + } 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); - } - } 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); - - } - } - } - + 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 @@ -181,73 +263,13 @@ 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; - 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; - } - } - } + bool only2gamma = true; + FillPhotonIDsForThreePhotons(photonIndices,only2gamma); + if (only2gamma) { + FillPhotonIDsForTwoPhotons(photonIndices); } } } - - /*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; } @@ -255,49 +277,50 @@ 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; 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))) { - found = false; - G4int n_trajectories = m_trajectoryContainer->entries(); - for (G4int iTrj=0; (iTrjGetTrackID()) { - photonID = trj->GetParentID(); - found = true; - } - } + 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; + } + } } 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 = 3; } - } } @@ -310,10 +333,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 @@ -322,12 +344,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); } @@ -338,11 +359,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; } @@ -352,8 +374,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; }