V9.4.1 cluster

source/digits_hits/include/GateToRoot.hh

@@ -370,7 +370,14 @@ private:
     G4ThreeVector m_positronGenerationPos;
     G4ThreeVector m_positronAnnihilPos;
 
-    G4double dxg1, dyg1, dzg1, dxg2, dyg2, dzg2;
+    G4double dxg1, dyg1, dzg1, dxg2, dyg2, dzg2, dxg3, dyg3, dzg3;
+
+   // new third photon
+    double photon3_source_x;
+    double photon3_source_y;
+    double photon3_source_z;
+    double photon3_energy;
+   // new third photon
 
     G4int saveRndm;
 
@@ -439,7 +446,7 @@ private:
     G4ThreeVector m_ionDecayPos_copy;
     G4ThreeVector m_positronGenerationPos_copy;
     G4ThreeVector m_positronAnnihilPos_copy;
-    G4double dxg1_copy, dyg1_copy, dzg1_copy, dxg2_copy, dyg2_copy, dzg2_copy;
+    G4double dxg1_copy, dyg1_copy, dzg1_copy, dxg2_copy, dyg2_copy, dzg2_copy, dxg3_copy, dyg3_copy, dzg3_copy;
     G4int m_RSEventID, m_RSRunID;
 
     /// ROOT FILE & DATAS FOR TRACKS INFOS
@@ -474,6 +481,7 @@ private:
     G4String m_particleName;
     Char_t m_volumeName[40], m_processName[40], m_parentparticleName[40];
     /* PY Descourt 08/09/2009 */
+
 };
 //--------------------------------------------------------------------------
 

source/digits_hits/include/GateTrajectoryNavigator.hh

@@ -21,6 +21,12 @@ class G4Trajectory;
 class GateTrajectoryNavigator
 {
 public:
+// new thied photon
+G4ThreeVector GetPhotonInitialPosition(G4int photonTrackID);
+G4ThreeVector GetPhotonInitialDirection(G4int photonTrackID);
+G4double  	GetPhotonInitialEnergy(G4int photonTrackID);
+
+// new third photon
 
   GateTrajectoryNavigator();
 
@@ -74,5 +80,4 @@ private:
   void               FillPhotonIDsForTwoPhotons(std::vector<G4int>& photonIndices);
 };
 
-
 #endif

source/digits_hits/src/GateToRoot.cc

@@ -113,6 +113,16 @@ GateToRoot::GateToRoot(const G4String &name, GateOutputMgr *outputMgr, DigiMode
         // the VOutputModule pure virtual method GiveNameOfFile()
 	, m_rootMessenger(0)
 {
+// new third photon 
+photon3_source_x = 0.0;
+photon3_source_y = 0.0;
+photon3_source_z = 0.0;
+photon3_energy   = 0.0;
+
+dxg3 = 0.0;
+dyg3 = 0.0;
+dzg3 = 0.0;
+// new third phoiton
 
 	//G4cout<<"GateToRoot::GateToRoot "<<G4endl;
     /*
@@ -525,6 +535,9 @@ 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("dxg3").c_str(), &dxg3, "dxg3/D");
+        m_RecStepTree->Branch(G4String("dyg3").c_str(), &dyg3, "dyg3/D");
+        m_RecStepTree->Branch(G4String("dzg3").c_str(), &dzg3, "dzg3/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");
@@ -740,6 +753,9 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) {
     dxg2 = 0.;
     dyg2 = 0.;
     dzg2 = 0.;
+	dxg3 = 0.;
+    dyg3 = 0.;
+    dzg3 = 0.;
 
     theCRData.photon1_phantom_Rayleigh = 0;
     theCRData.photon2_phantom_Rayleigh = 0;
@@ -775,6 +791,9 @@ void GateToRoot::RecordBeginOfEvent(const G4Event *evt) {
             dxg2 = dxg2_copy;
             dyg2 = dyg2_copy;
             dzg2 = dzg2_copy;
+			dxg3 = dxg3_copy;
+            dyg3 = dyg3_copy;
+            dzg3 = dzg3_copy;
             fSkipRecStepData = 0;
             //G4cout << "GateToRoot::RecordBeginOfEvent \n";
             //PrintRecStep();
@@ -878,28 +897,89 @@ void GateToRoot::RecordEndOfEvent(const G4Event *event) {
 						G4cout
 								<< "GateToRoot:  ROOT: Cannot find histo" << hist_name << Gateendl;
 				}
+				//new edited 
+auto photonIDs = m_trajectoryNavigator->FindAnnihilationGammasTrackID();
+G4int photon3ID = (photonIDs.size() >= 3) ? photonIDs[2] : 0;
+
+				// new 3 rd photon
+				if (photon3ID > 0) {
+					// Source position
+					auto pos3 = m_trajectoryNavigator->GetPhotonInitialPosition(photon3ID);
+					photon3_source_x = pos3.x();
+					photon3_source_y = pos3.y();
+					photon3_source_z = pos3.z();
+					// Energy
+					photon3_energy = m_trajectoryNavigator->GetPhotonInitialEnergy(photon3ID);
+					// Direction
+					auto dir3 = m_trajectoryNavigator->GetPhotonInitialDirection(photon3ID);
+					dxg3 = dir3.x();
+					dyg3 = dir3.y();
+					dzg3 = dir3.z();
+				}
+                // new third 
 
 				// Histo of acolinearity angle distribution
 
-				G4double dev = (dxg1 * dxg2 + dyg1 * dyg2 + dzg1 * dzg2) /
+				G4double dev12 = (dxg1 * dxg2 + dyg1 * dyg2 + dzg1 * dzg2) /
 							   ((sqrt(dxg1 * dxg1 + dyg1 * dyg1 + dzg1 * dzg1)) *
 								(sqrt(dxg2 * dxg2 + dyg2 * dyg2 + dzg2 * dzg2)));
-				if (dzg1 > dzg2) { dev = rad2deg(acos(-dev)); }
-				else { dev = rad2deg(acos(dev)) - 180; }
 
-				if (std::isnan(dev)) dev = 0.;
+				G4double dev13 = (dxg1 * dxg3 + dyg1 * dyg3 + dzg1 * dzg3) /
+							   ((sqrt(dxg1 * dxg1 + dyg1 * dyg1 + dzg1 * dzg1)) *
+								(sqrt(dxg3 * dxg3 + dyg3 * dyg3 + dzg3 * dzg3)));
+
+				G4double dev23 = (dxg2 * dxg3 + dyg2 * dyg3 + dzg2 * dzg3) /
+							   ((sqrt(dxg2 * dxg2 + dyg2 * dyg2 + dzg2 * dzg2)) *
+								(sqrt(dxg3 * dxg3 + dyg3 * dyg3 + dzg3 * dzg3)));
+
+				if (dzg1 > dzg2) { dev12 = rad2deg(acos(-dev12)); }
+				else { dev12 = rad2deg(acos(dev12)) - 180; }
+
+				if (std::isnan(dev12)) dev12 = 0.;
 
 				// G4cout<< " dev = " << dev << Gateendl;
 
 				hist_name = "Acolinea_Angle_Distribution_deg";
 				hist = NULL;
 				if ((hist = (TH1F *) m_working_root_directory->GetList()->FindObject(hist_name)) != NULL) {
-					hist->Fill(dev);
+					hist->Fill(dev12);
+				} else {
+					//if (nVerboseLevel > 0)
+					G4cout << "GateToRoot:  ROOT: Cannot find histo " << hist_name << Gateendl;
+				}
+
+				if (dzg1 > dzg3) { dev13 = rad2deg(acos(-dev13)); }
+				else { dev13 = rad2deg(acos(dev13)) - 180; }
+
+				if (std::isnan(dev13)) dev13 = 0.;
+
+				// G4cout<< " dev13 = " << dev13 << Gateendl;
+
+				hist_name = "Acolinea_Angle_Distribution_deg";
+				hist = NULL;
+				if ((hist = (TH1F *) m_working_root_directory->GetList()->FindObject(hist_name)) != NULL) {
+					hist->Fill(dev13);
 				} else {
 					//if (nVerboseLevel > 0)
 					G4cout << "GateToRoot:  ROOT: Cannot find histo " << hist_name << Gateendl;
 				}
 
+				if (dzg2 > dzg3) { dev23 = rad2deg(acos(-dev23)); }
+				else { dev23 = rad2deg(acos(dev23)) - 180; }
+
+				if (std::isnan(dev23)) dev23 = 0.;
+
+				// G4cout<< " dev23 = " << dev23 << Gateendl;
+
+				hist_name = "Acolinea_Angle_Distribution_deg";
+				hist = NULL;
+				if ((hist = (TH1F *) m_working_root_directory->GetList()->FindObject(hist_name)) != NULL) {
+					hist->Fill(dev23);
+				} else {
+					//if (nVerboseLevel > 0)
+					G4cout << "GateToRoot:  ROOT: Cannot find histo " << hist_name << Gateendl;
+				}
+
 				TNtuple *ntuple;
 				G4String ntuple_name = "Gate";
 				if ((ntuple = (TNtuple *) m_working_root_directory->GetList()->FindObject(ntuple_name)) == NULL) {
@@ -1150,6 +1230,12 @@ void GateToRoot::RecordStepWithVolume(const GateVVolume *, const G4Step *aStep)
             } else {
                 procName = "";
             }
+			if (aStep->GetTrack()->GetTrackID() == 1 && procName == "annihil") {
+                dxg3 = momentumDirection.x();
+                dyg3 = momentumDirection.y();
+                dzg3 = momentumDirection.z();
+            }
+
             if (aStep->GetTrack()->GetTrackID() == 2 && procName == "annihil") {
 
                 dxg1 = momentumDirection.x();
@@ -1377,7 +1463,10 @@ void GateToRoot::CloseTracksRootFile() {
         dxg2_copy = dxg2;
         dyg2_copy = dyg2;
         dzg2_copy = dzg2;
-
+		dxg3_copy = dxg3;
+        dyg3_copy = dyg3;
+        dzg3_copy = dzg3;
+
         m_positronKinEnergy_copy = m_positronKinEnergy;
         m_ionDecayPos_copy = m_ionDecayPos;
         m_positronGenerationPos_copy = m_positronGenerationPos;
@@ -1408,6 +1497,9 @@ void GateToRoot::PrintRecStep() {
     G4cout << "dxg2 = " << dxg2 << Gateendl;
     G4cout << "dyg2 = " << dyg2 << Gateendl;
     G4cout << "dzg2 = " << dzg2 << Gateendl;
+	G4cout << "dxg3 = " << dxg3 << Gateendl;
+    G4cout << "dyg3 = " << dyg3 << Gateendl;
+    G4cout << "dzg3 = " << dzg3 << Gateendl;
     G4cout << "photon1_phantom_Rayleigh = " << theCRData.photon2_phantom_Rayleigh << Gateendl;
     G4cout << "photon2_phantom_Rayleigh = " << theCRData.photon2_phantom_Rayleigh << Gateendl;
     G4cout << "photon3_phantom_Rayleigh = " << theCRData.photon3_phantom_Rayleigh << Gateendl;
@@ -1431,6 +1523,9 @@ void GateToRoot::PrintRecStep() {
     G4cout << "dxg2 = " << dxg2_copy << Gateendl;
     G4cout << "dyg2 = " << dyg2_copy << Gateendl;
     G4cout << "dzg2 = " << dzg2_copy << Gateendl;
+	G4cout << "dxg3 = " << dxg3_copy << Gateendl;
+    G4cout << "dyg3 = " << dyg3_copy << Gateendl;
+    G4cout << "dzg3 = " << dzg3_copy << Gateendl;
     G4cout << "photon1_phantom_Rayleigh = " << theCRData_copy.photon2_phantom_Rayleigh << Gateendl;
     G4cout << "photon2_phantom_Rayleigh = " << theCRData_copy.photon2_phantom_Rayleigh << Gateendl;
     G4cout << "photon3_phantom_Rayleigh = " << theCRData_copy.photon3_phantom_Rayleigh << Gateendl;
@@ -1579,6 +1674,9 @@ void GateToRoot::OpenTracksFile() {
     m_RecStepTree->SetBranchAddress(G4String("dxg2").c_str(), &dxg2);
     m_RecStepTree->SetBranchAddress(G4String("dyg2").c_str(), &dyg2);
     m_RecStepTree->SetBranchAddress(G4String("dzg2").c_str(), &dzg2);
+	m_RecStepTree->SetBranchAddress(G4String("dxg3").c_str(), &dxg3);
+	m_RecStepTree->SetBranchAddress(G4String("dyg3").c_str(), &dyg3);
+	m_RecStepTree->SetBranchAddress(G4String("dyg3").c_str(), &dyg3);
     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);

source/digits_hits/src/GateTrajectoryNavigator.cc

@@ -189,7 +189,7 @@ void GateTrajectoryNavigator::FillPhotonIDsForThreePhotons(std::vector<G4int>& p
                       "gammas vertices : dist (mm) "
                    << dist12 / mm << Gateendl;
 
-          if (dist12 / mm < 1E-7) {
+            if (dist12 / mm < 1E-7) {
             if (nVerboseLevel > 1) {
               G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common "
                         "vertex for the two annihilation gammas :"
@@ -198,7 +198,7 @@ void GateTrajectoryNavigator::FillPhotonIDsForThreePhotons(std::vector<G4int>& p
             }
 
             // check if third one also from same vertex
-            if ((vert2 - vert3).mag() / mm < 1E-7) {
+             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());
@@ -260,7 +260,7 @@ void GateTrajectoryNavigator::FillPhotonIDsForTwoPhotons(std::vector<G4int>& pho
                     "gammas vertices : dist (mm) "
                  << dist / mm << Gateendl;
 
-        if (dist / mm < 1E-7) {
+          if (dist / mm < 1E-7) {
           if (nVerboseLevel > 1) {
             G4cout << "[GateTrajectoryNavigator::FindAnnihilationGammasTrackID] : Found common "
                       "vertex for the two annihilation gammas :"
@@ -425,4 +425,38 @@ void GateTrajectoryNavigator::Initialize()
   m_photonIDVec.clear();
   m_positronTrackID = -1;
 }
+  // capturing the third photon
+G4ThreeVector GateTrajectoryNavigator::GetPhotonInitialPosition(G4int photonTrackID)
+{
+	for (size_t i = 0; i < m_trajectoryContainer->entries(); i++) {
+    	G4Trajectory* trj = (G4Trajectory*)((*m_trajectoryContainer)[i]);
+    	if (trj->GetTrackID() == photonTrackID) {
+        	return trj->GetPoint(0)->GetPosition();
+    	}
+	}
+	return G4ThreeVector(0,0,0);
+}
+
+G4ThreeVector GateTrajectoryNavigator::GetPhotonInitialDirection(G4int photonTrackID)
+{
+	for (size_t i = 0; i < m_trajectoryContainer->entries(); i++) {
+    	G4Trajectory* trj = (G4Trajectory*)((*m_trajectoryContainer)[i]);
+    	if (trj->GetTrackID() == photonTrackID) {
+        	return trj->GetInitialMomentum().unit();
+    	}
+	}
+	return G4ThreeVector(0,0,0);
+}
+
+G4double GateTrajectoryNavigator::GetPhotonInitialEnergy(G4int photonTrackID)
+{
+	for (size_t i = 0; i < m_trajectoryContainer->entries(); i++) {
+    	G4Trajectory* trj = (G4Trajectory*)((*m_trajectoryContainer)[i]);
+    	if (trj->GetTrackID() == photonTrackID) {
+        	return trj->GetInitialKineticEnergy();
+    	}
+	}
+	return 0.0;
+}
+  // Capturing the third photon
 

source/geometry/src/GateCylindricalPETSystem.cc

@@ -52,9 +52,9 @@ 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);
+  GateDigitizerMgr* digitizerMgr = GateDigitizerMgr::GetInstance();
+  GateCoincidenceSorter* coincidenceSorter = new GateCoincidenceSorter(digitizerMgr,"Coincidences");
+  digitizerMgr->AddNewCoincidenceSorter(coincidenceSorter);
 
 #ifdef GATE_USE_LMF