diff --git a/source/digits_hits/src/GateAnalysis.cc b/source/digits_hits/src/GateAnalysis.cc index 31a15001e..27026d8c9 100644 --- a/source/digits_hits/src/GateAnalysis.cc +++ b/source/digits_hits/src/GateAnalysis.cc @@ -35,6 +35,10 @@ #include "GateToRoot.hh" #include "GateDigitizerMgr.hh" +#include "G4UnitsTable.hh" +#include +#include + //-------------------------------------------------------------------------------------------------- GateAnalysis::GateAnalysis(const G4String& name, GateOutputMgr* outputMgr,DigiMode digiMode) : GateVOutputModule(name,outputMgr,digiMode) @@ -114,6 +118,39 @@ void GateAnalysis::RecordBeginOfEvent(const G4Event* ) //-------------------------------------------------------------------------------------------------- +/////////////////////////////////////////////////////////////////////// + +std::vector argsort(const std::vector& v) { + std::vector idx(v.size()); + std::iota(idx.begin(), idx.end(), 0); + + std::sort(idx.begin(), idx.end(), + [&v](G4int i1, G4int i2) { return v[i1] < v[i2]; }); + + return idx; +} + +std::vector inverse_argsort(const std::vector& unsorted_array) { + std::vector sorted_indices = argsort(unsorted_array); + std::vector original_indices(sorted_indices.size()); + + for (G4int i = 0; i < sorted_indices.size(); ++i) { + original_indices[sorted_indices[i]] = i; + } + + return original_indices; +} +/* +template +void printVector(const std::vector& v) { + for (T ii: v) + std::cout << std::setprecision(16) << ii << ' '; + std::cout << std::endl; +} +*/ +/////////////////////////////////////////////////////////////////////// + + //-------------------------------------------------------------------------------------------------- void GateAnalysis::RecordEndOfEvent(const G4Event* event) { @@ -140,6 +177,148 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) { //OK GND 2022 std::vector CHC_vector = GetOutputMgr()->GetHitCollections(); + + ///////////////////////////////////////////////////////////////// + + //Number of sensitive detectors + //G4cout << CHC_vector.size() << G4endl; + //G4cout << "--" << G4endl; + + //G4bool oneInFirstLayer = false; + //G4bool twoInFirstLayer = false; + //G4double oneKeep, twoKeep; + + // The following bit used to be part of the main loop + // todo: Check if moving this causes some errors + G4int photon1ID = 0; + G4int photon2ID = 0; + + //////////// + // 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"; + + // Print the + //for (G4int ii: photonIDVec) + //std::cout << ii << ' '; + //std::cout << std::endl; + + 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; + + /////////////////////////////////////////////////////////////////// + + + /////////////////////////////////////////////////////////////////// + + + + // Get all the Compton (and Rayleigh) interactions before processing the hit layers individually + std::vector comptonCount1, comptonCount2, rayleighCount1, rayleighCount2; + + G4bool multi_detector = CHC_vector.size() > 1; + + // Only do this if there are multiple layers / geometries + if (multi_detector) { + + std::vector comptonTimes1, comptonTimes2, rayleighTimes1, rayleighTimes2; + + for (size_t i=0; ientries(); + //G4int c1 = 0, c2 = 0, r1 = 0, r2 = 0; + + // Hits loop + for (G4int iHit=0;iHitGetTrackID(); + G4String processName = (*CHC)[iHit]->GetProcess(); + + if (processName.find("ompt") != G4String::npos) { + + if (crystalTrackID == photon1ID) { + comptonTimes1.push_back((*CHC)[iHit]->GetTime()); + //comptonCount1.push_back(c1); + //c1++; + } + if (crystalTrackID == photon2ID) { + comptonTimes2.push_back((*CHC)[iHit]->GetTime()); + //comptonCount2.push_back(c2); + //c2++; + } + } + + if (processName.find("Rayl") != G4String::npos) { + + if (crystalTrackID == photon1ID) { + rayleighTimes1.push_back((*CHC)[iHit]->GetTime()); + //rayleighCount1.push_back(r1); + //r1++; + } + + if (crystalTrackID == photon2ID) { + rayleighTimes2.push_back((*CHC)[iHit]->GetTime()); + //rayleighCount2.push_back(r2); + //r2++; + } + + } + } + } + } + // Get the correct ordering + //comptonCount1 = argsort(comptonTimes1); + //comptonCount2 = argsort(comptonTimes2); + comptonCount1 = inverse_argsort(comptonTimes1); + comptonCount2 = inverse_argsort(comptonTimes2); + + //~ if (comptonCount1.size() > 0) { + //~ printVector(comptonTimes1); + //~ printVector(comptonCount1); + //~ G4cout << G4endl; + //~ } + + //rayleighCount1 = argsort(rayleighTimes1); + //rayleighCount2 = argsort(rayleighTimes2); + rayleighCount1 = inverse_argsort(rayleighTimes1); + rayleighCount2 = inverse_argsort(rayleighTimes2); + } + + G4int ic1=0, ic2=0, ir1=0, ir2=0; + + /////////////////////////////////////////////////////////////////// + for (size_t i=0; iFindPositronTrackID(); - /*if (positronID == 0) - { - if (nVerboseLevel > 0) G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : positronID == 0\n"; - } + //if (positronID == 0) + //{ + //if (nVerboseLevel > 0) G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : positronID == 0\n"; + //} - if (nVerboseLevel > 1) G4cout << "GateAnalysis::RecordEndOfEvent : positronID : " << positronID << Gateendl; - */ + //if (nVerboseLevel > 1) G4cout << "GateAnalysis::RecordEndOfEvent : positronID : " << positronID << Gateendl; + //////////// //search the two gammas @@ -202,6 +385,8 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) if (nVerboseLevel > 0 && photonIDVec.size() > 2) G4cout << "GateAnalysis::RecordEndOfEvent : WARNING : photonID vector size > 2\n"; + //printVector(photonIDVec); + photon1ID = photonIDVec[0]; photon2ID = (photonIDVec.size() >= 2) ? photonIDVec[1] : 0; } @@ -219,6 +404,9 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) << "GateAnalysis::RecordEndOfEvent : photon1ID : " << photon1ID << " photon2ID : " << photon2ID << Gateendl; + */ + ///////////////////////////////////////////////////////// + // analysis of the phantom hits to count the comptons, etc. @@ -381,23 +569,87 @@ void GateAnalysis::RecordEndOfEvent(const G4Event* event) { G4int crystalTrackID = (*CHC)[iHit]->GetTrackID(); G4String processName = (*CHC)[iHit]->GetProcess(); + + ///////////////////////////////////////////////////////////// + // Check the ordering + /* + G4cout << std::setprecision(16) << (*CHC)[iHit]->GetTime()/second << G4endl; // somehow does not show all 16 digits? + if (processName == "Compton") { + std::cout << std::setprecision (16) << (*CHC)[iHit]->GetTime()/second << " " << crystalTrackID << std::endl; + + if ((crystalTrackID == photon1ID) && (i == 1) && oneInFirstLayer){ + if (oneKeep > (*CHC)[iHit]->GetTime()/second){ + std::cout << std::setprecision (16) << oneKeep << std::endl; + std::cout << std::setprecision (16) << (*CHC)[iHit]->GetTime()/second << " " << crystalTrackID << std::endl; + std::cout << std::endl;} + + } + + if ((crystalTrackID == photon2ID) && (i == 1) && twoInFirstLayer){ + if (twoKeep > (*CHC)[iHit]->GetTime()/second){ + std::cout << std::setprecision (16) << twoKeep << std::endl; + std::cout << std::setprecision (16) << (*CHC)[iHit]->GetTime()/second << " " << crystalTrackID << std::endl; + std::cout << std::endl;} + } + + + if ((crystalTrackID == photon1ID) && (i == 0)) oneInFirstLayer = true, oneKeep = (*CHC)[iHit]->GetTime()/second; + if ((crystalTrackID == photon2ID) && (i == 0)) twoInFirstLayer = true, twoKeep = (*CHC)[iHit]->GetTime()/second; + } + */ + // 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++; + if (crystalTrackID == photon1ID) { + if (multi_detector) { + photon1_crystal_compton = comptonCount1[ic1] + 1; + ic1++; + } else { + photon1_crystal_compton++; + } + } + + if (crystalTrackID == photon2ID) { + if (multi_detector) { + photon2_crystal_compton = comptonCount2[ic2] + 1; + ic2++; + } else { + 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++; + if (crystalTrackID == photon1ID) { + if (multi_detector) { + photon1_crystal_Rayleigh = rayleighCount1[ir1] + 1; + ir1++; + } else { + photon1_crystal_Rayleigh++; + } + } + + if (crystalTrackID == photon2ID) { + if (multi_detector) { + photon2_crystal_Rayleigh = rayleighCount2[ir2] + 1; + ir2++; + } else { + photon2_crystal_Rayleigh++; + } + } + } + ///////////////////////////////////////////////////////////// + G4int PDGEncoding = (*CHC)[iHit]->GetPDGEncoding(); if (nVerboseLevel > 2) G4cout << "GateAnalysis::RecordEndOfEvent : HitsCollection: processName : <" << processName diff --git a/source/digits_hits/src/GateToRoot.cc b/source/digits_hits/src/GateToRoot.cc index 269b33879..9a1758991 100644 --- a/source/digits_hits/src/GateToRoot.cc +++ b/source/digits_hits/src/GateToRoot.cc @@ -1572,7 +1572,7 @@ 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.photon2_phantom_Rayleigh = aCRData.photon2_phantom_Rayleigh; theCRData.photon1_phantom_compton = aCRData.photon1_phantom_compton; theCRData.photon2_phantom_compton = aCRData.photon2_phantom_compton; strcpy(theCRData.theComptonVolumeName1, aCRData.theComptonVolumeName1); @@ -1583,7 +1583,7 @@ void GateToRoot::RecordPHData(ComptonRayleighData aCRData) { void GateToRoot::GetPHData(ComptonRayleighData &aCRData) { aCRData.photon1_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; - aCRData.photon2_phantom_Rayleigh = theCRData.photon1_phantom_Rayleigh; + aCRData.photon2_phantom_Rayleigh = theCRData.photon2_phantom_Rayleigh; aCRData.photon1_phantom_compton = theCRData.photon1_phantom_compton; aCRData.photon2_phantom_compton = theCRData.photon2_phantom_compton; strcpy(aCRData.theComptonVolumeName1, theCRData.theComptonVolumeName1); diff --git a/source/geometry/src/GateCPETSystem.cc b/source/geometry/src/GateCPETSystem.cc index 7d9083862..3c7baaa01 100644 --- a/source/geometry/src/GateCPETSystem.cc +++ b/source/geometry/src/GateCPETSystem.cc @@ -41,9 +41,13 @@ GateCPETSystem::GateCPETSystem(const G4String& itsName) // Integrate a coincidence sorter into the digitizer //OK GND 2022 + // MR 2025: For multi-detector geometries: suppress the generation of multiple coincidence sorters called + // "Coincidences" to avoid ambiguity issues GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); - GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); - digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + if (digitizerMgr->m_CoincidenceSortersList.size() == 0) { + GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); + digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } SetOutputIDName((char *)"gantryID",0); SetOutputIDName((char *)"sectorID",1); diff --git a/source/geometry/src/GateCylindricalPETSystem.cc b/source/geometry/src/GateCylindricalPETSystem.cc index a6d2c9f5c..d4f9687a5 100644 --- a/source/geometry/src/GateCylindricalPETSystem.cc +++ b/source/geometry/src/GateCylindricalPETSystem.cc @@ -53,8 +53,28 @@ GateCylindricalPETSystem::GateCylindricalPETSystem(const G4String& itsName) // Integrate a coincidence sorter into the digitizer //OK GND 2022 GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); - GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); - digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + // In case of multiple detector geometries, multiple coincidence sortes are created with the name "Coincidences", + // which creates an ambiguity and eventually crashes with a segmentation violation; using /gate/digitizerMgr/list does + // not list multiple entries with identical names either; therefore we add a check + + std::vector coincidenceSortersList = digitizerMgr->m_CoincidenceSortersList; + + //for (G4int ii = 0; ii < coincidenceSortersList.size(); ++ii) { + // std::cout << std::to_string(ii) + " " << coincidenceSortersList[ii]->GetInputName() << " " << coincidenceSortersList[ii]->GetOutputName() << std::endl; + //} + + // Name options: itsName and GetName() are always "systems/cylindricalPET" in this class; so we go with GetOwnName() + //std::cout << itsName << " " << GetName() << " " << GetOwnName() << std::endl; + + if (coincidenceSortersList.size() == 0) { + GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); + //GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr, "Coincidences_" + GetOwnName()); + digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } + else { + //GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr, "Coincidences_" + GetOwnName()); + //digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } #ifdef GATE_USE_LMF diff --git a/source/geometry/src/GateEcatAccelSystem.cc b/source/geometry/src/GateEcatAccelSystem.cc index ea0b4b962..109e0f275 100644 --- a/source/geometry/src/GateEcatAccelSystem.cc +++ b/source/geometry/src/GateEcatAccelSystem.cc @@ -35,9 +35,13 @@ GateEcatAccelSystem::GateEcatAccelSystem(const G4String& itsName) // Integrate a coincidence sorter into the digitizer //OK GND 2022 + // MR 2025: For multi-detector geometries: suppress the generation of multiple coincidence sorters called + // "Coincidences" to avoid ambiguity issues GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); - GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); - digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + if (digitizerMgr->m_CoincidenceSortersList.size() == 0) { + GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); + digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } // Insert a sinogram maker and a ECAT7 writer into the output manager GateOutputMgr *outputMgr = GateOutputMgr::GetInstance(); diff --git a/source/geometry/src/GateEcatSystem.cc b/source/geometry/src/GateEcatSystem.cc index 3c11063f6..60801200c 100644 --- a/source/geometry/src/GateEcatSystem.cc +++ b/source/geometry/src/GateEcatSystem.cc @@ -35,9 +35,13 @@ GateEcatSystem::GateEcatSystem(const G4String& itsName) // Integrate a coincidence sorter into the digitizer //OK GND 2022 + // MR 2025: For multi-detector geometries: suppress the generation of multiple coincidence sorters called + // "Coincidences" to avoid ambiguity issues GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); - GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); - digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + if (digitizerMgr->m_CoincidenceSortersList.size() == 0) { + GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); + digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } // Insert a sinogram maker and a ECAT7 writer into the output manager GateOutputMgr *outputMgr = GateOutputMgr::GetInstance(); diff --git a/source/geometry/src/GateOPETSystem.cc b/source/geometry/src/GateOPETSystem.cc index 8cdda5afd..1f8c2b4f4 100644 --- a/source/geometry/src/GateOPETSystem.cc +++ b/source/geometry/src/GateOPETSystem.cc @@ -50,9 +50,13 @@ GateOPETSystem::GateOPETSystem(const G4String& itsName) // Integrate a coincidence sorter into the digitizer //OK GND 2022 - GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); + // MR 2025: For multi-detector geometries: suppress the generation of multiple coincidence sorters called + // "Coincidences" to avoid ambiguity issues + GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); + if (digitizerMgr->m_CoincidenceSortersList.size() == 0) { GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } #ifdef GATE_USE_LMF // Insert an LMF output module into the output manager diff --git a/source/geometry/src/GatePETScannerSystem.cc b/source/geometry/src/GatePETScannerSystem.cc index 1f2a09e9e..cc5c5f4ab 100644 --- a/source/geometry/src/GatePETScannerSystem.cc +++ b/source/geometry/src/GatePETScannerSystem.cc @@ -23,7 +23,11 @@ GatePETScannerSystem::GatePETScannerSystem(const G4String& itsName) //G4cout << " Constructeur GatePETScannerSystem \n"; // Integrate a coincidence sorter into the digitizer //OK GND 2022 + // MR 2025: For multi-detector geometries: suppress the generation of multiple coincidence sorters called + // "Coincidences" to avoid ambiguity issues GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance(); - GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); - digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + if (digitizerMgr->m_CoincidenceSortersList.size() == 0) { + GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences"); + digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter); + } }