Clean up of the code, working for XYZ case

OpenGATE · Sep 18, 2024 · c6c1587 · c6c1587
1 parent 6c581e1
commit c6c1587
Show file tree

Hide file tree

Showing 5 changed files with 129 additions and 69 deletions.
diff --git a/source/digits_hits/include/GateDigi.hh b/source/digits_hits/include/GateDigi.hh
@@ -125,7 +125,7 @@ public:
       inline G4double GetScannerRotAngle() const                   	      	{ return m_scannerRotAngle; }
 
       inline void  SetOutputVolumeID(const GateOutputVolumeID& outputVolumeID)        	{ m_outputVolumeID = outputVolumeID; }
-      inline void  SetOutputVolumeID(const G4int outputVolumeID,G4int depth)        	{ m_outputVolumeID[depth] = outputVolumeID; }
+      inline void  SetOutputVolumeID(const G4int copyNum,G4int depth)        	{ m_outputVolumeID[depth] = copyNum; }
       inline const GateOutputVolumeID& GetOutputVolumeID()  const             	      	{ return m_outputVolumeID; }
       inline G4int GetComponentID(size_t depth) const    { return (m_outputVolumeID.size()>depth) ? m_outputVolumeID[depth] : -1; }
 

diff --git a/source/digits_hits/include/GateDiscretizerModule.hh b/source/digits_hits/include/GateDiscretizerModule.hh
@@ -67,7 +67,8 @@ public:
      void SetResolutionZ(G4double val)   { m_resolutionZ = val;  }
 
 
-     void SetVirtualIDs(int nBinsX, int nBinsY,int nBinsZ, G4ThreeVector& pos);
+     void SetVirtualIDs(int nBinsX, int nBinsY,int nBinsZ,double pitchX,double pitchY,double pitchZ, G4ThreeVector& pos);
+     void SetVirtualID(int nBins,double pitch, G4double pos, int depth);
      void DescribeMyself(size_t );
 
 protected:

diff --git a/source/digits_hits/src/GateDiscretizerModule.cc b/source/digits_hits/src/GateDiscretizerModule.cc
@@ -30,6 +30,7 @@
 
 #include "GateDigitizerMgr.hh"
 
+#include "G4Box.hh"
 #include "G4SystemOfUnits.hh"
 #include "G4EventManager.hh"
 #include "G4Event.hh"
@@ -84,31 +85,36 @@ void GateDiscretizerModule::Digitize()
 {
 
 
+//From here
 
 	GateSpatialResolution* digi_SpatialResolution;
 	digi_SpatialResolution =  dynamic_cast<GateSpatialResolution*>(m_digitizer->FindDigitizerModule("digitizerMgr/"
 	    		  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  	  +m_digitizer->GetSD()->GetName()
 																						  +"/SinglesDigitizer/"
 																						  + m_digitizer->GetName()
-																						  + "/SpatialResolution"));
+																						  + "/spatialResolution"));
 
 
 
 	G4double resolutionX;
 	G4double resolutionY;
 	G4double resolutionZ;
 
+	G4int nBinsX;
+	G4int nBinsY;
+	G4int nBinsZ;
 
+	G4double pitchX;
+	G4double pitchY;
+	G4double pitchZ;
 
 
 
-
-
-
-
+	//TODO bulletproof this!
 
 	if (digi_SpatialResolution->GetFWHM())
 	{
+
 		if(m_nameAxis.find('X') != std::string::npos)
 			{resolutionX = digi_SpatialResolution->GetFWHM();}
 
@@ -121,6 +127,7 @@ void GateDiscretizerModule::Digitize()
 
 
 	else if(digi_SpatialResolution->GetFWHMx() || digi_SpatialResolution->GetFWHMy()|| digi_SpatialResolution->GetFWHMz()){
+
 		if(digi_SpatialResolution->GetFWHMx() && m_nameAxis.find('X') != std::string::npos)
 				{resolutionX = digi_SpatialResolution->GetFWHMx();}
 
@@ -131,25 +138,12 @@ void GateDiscretizerModule::Digitize()
 				{resolutionZ = digi_SpatialResolution->GetFWHMz();}
 	}
 
-	else if (m_nameAxis == "XYZ" && (m_resolution==0) || (m_resolutionX==0 || m_resolutionY==0 ||m_resolutionZ==0))
-	{
-
-		GateError("***ERROR*** Discretization in XYZ has been selected but at least one axis has resolution = 0. Please ensure that all axis have a non zero value of resolution.");
-	}
 
 
-	if (m_resolution != 0)
+	else if (m_nameAxis == "XYZ" && (m_resolution==0) || (m_resolutionX==0 || m_resolutionY==0 ||m_resolutionZ==0))
 	{
-		resolutionX = m_resolution;
-		resolutionY = m_resolution;
-		resolutionZ = m_resolution;
-	}
 
-	else
-	{
-		resolutionX = m_resolutionX;
-		resolutionY = m_resolutionY;
-		resolutionZ = m_resolutionZ;
+		GateError("***ERROR*** Discretization in XYZ has been selected but at least one axis has resolution = 0. Please ensure that all axis have a non zero value of resolution.");
 	}
 
 
@@ -167,14 +161,15 @@ void GateDiscretizerModule::Digitize()
 				 "Please, add them to your geometry macro in /gate/systems/cylindricalPET/XXX/attach    YYY. Abort.\n");
 	}
 
+
+
 	m_systemDepth = m_system->GetTreeDepth();
 	//TODO I need to create another parameter asking if spatial resolution is provided by hand.
 	//I need to get system, system depth. Check that the size for the last 4 solids is the same (Monolithic crystal)
 	//Need to access spatial resolution
 
 
 
-
 	G4String digitizerName = m_digitizer->m_digitizerName;
 	G4String outputCollName = m_digitizer-> GetOutputName();
 
@@ -195,6 +190,10 @@ void GateDiscretizerModule::Digitize()
 
 
 
+	G4double xLength,yLength,zLength;
+
+
+
   if (IDC)
      {
 	  G4int n_digi = IDC->entries();
@@ -205,45 +204,62 @@ void GateDiscretizerModule::Digitize()
 		  inputDigi=(*IDC)[i];
 		  // ***** the following part of the code to adapt
 		  /// *** This part is from ProcessPulseList
+		  if (inputDigi->GetVolumeID().size()){
+		  G4Box* box = dynamic_cast<G4Box*>(inputDigi->GetVolumeID().GetBottomCreator()->GetLogicalVolume()->GetSolid());
+	        xLength = 2*box->GetXHalfLength();
+	        yLength = 2*box->GetYHalfLength();
+	        zLength = 2*box->GetZHalfLength();
+		  }
+		  else std::cout<<"Well that's another problem here..."<< std::endl;
+
+
+
+
+
 
-		  ////// ** This part is from ProcessOnePulse
-		     for (iter=OutputDigiCollectionVector->begin(); iter!= OutputDigiCollectionVector->end() ; ++iter)
-		     {
-		    	 if ( (*iter)->GetVolumeID()   == inputDigi->GetVolumeID() )
-		    	 {
-		    		continue;
-		    		if(m_nameAxis=="XYZ"){continue;}
-		    		/* if(m_parameter=="numAxis"){
-		                 DummyMethod1(inputDigi);
-		    		 }
-		    		 else{
-		    			 DummyMethod2( inputDigi );
-		    		 }
-					*/
 		      if (nVerboseLevel>1)
 		      {
 		    	 G4cout << "Merged previous pulse for volume " << inputDigi->GetVolumeID()
 		 		 << " with new pulse of energy " << G4BestUnit(inputDigi->GetEnergy(),"Energy") <<".\n"
 		 		 << "Resulting pulse is: \n"
 		 		 << **iter << Gateendl << Gateendl ;
 		      }
-		 	  break;
-		       }
-		     }
 
-		     if ( iter == OutputDigiCollectionVector->end() )
-		     {
+
 		       m_outputDigi = new GateDigi(*inputDigi);
 		       m_outputDigi->SetEnergyIniTrack(-1);
 		       m_outputDigi->SetEnergyFin(-1);
+
+
+
+
+		       if(m_nameAxis=="XYZ"){
+
+		       			    		pitchX = calculatePitch(xLength,resolutionX);
+		       			    		nBinsX = int(xLength/pitchX);
+		       			    		pitchY = calculatePitch(yLength,resolutionY);
+		       			    		nBinsY = int(yLength/pitchY);
+		       			    		pitchZ = calculatePitch(zLength,resolutionZ);
+		       			    		nBinsZ = int(zLength/pitchZ);
+
+
+		       			    		G4ThreeVector localPos = inputDigi->GetLocalPos();
+		       			    		//std::cout<<"Local position = " << localPos <<" nBinsX "<<nBinsX<<" nBinsY "<<nBinsY<<" nBinsZ "<<nBinsZ<<std::endl;
+		       			    		SetVirtualIDs(nBinsX,nBinsY,nBinsZ,pitchX,pitchY,pitchZ,localPos);
+
+		       		    			}
+
+		       std::cout<<"m_outputdigi "<<m_outputDigi->GetOutputVolumeID()<<std::endl;
+
+
 		       if (nVerboseLevel>1)
 		 	  G4cout << "Created new pulse for volume " << inputDigi->GetVolumeID() << ".\n"
 		 		 << "Resulting pulse is: \n"
 		 		 << *m_outputDigi << Gateendl << Gateendl ;
 		      /// !!!!!! The following line should be kept !!!! -> inserts the outputdigi to collection
 		       m_OutputDigiCollection->insert(m_outputDigi);
 
-		     }
+
 		 ////// ** End of the part from ProcessOnePulse
 
 
@@ -286,10 +302,13 @@ void GateDiscretizerModule::SetNameAxis(const G4String &param)
 
 
 
-
+/*
 G4double GateDiscretizerModule::calculatePitch(G4double crystal_size, G4double spatial_resolution) {
-         // Target pitch size
+
+
+	// Target pitch size
          double target_pitch = spatial_resolution / 2.0;
+         std::cout<<"Spatial Resolution: "<<spatial_resolution<<std::endl;
          double best_pitch = 0;
          double min_diff = std::numeric_limits<double>::max();
 
@@ -302,51 +321,91 @@ G4double GateDiscretizerModule::calculatePitch(G4double crystal_size, G4double s
              }
          }
 
+         std::cout<<"Best pitch is: "<<best_pitch<<" weird since the condition is to be smaller than: "<< min_diff<<std::endl;
+
          return best_pitch;
      }
+*/
 
+//TODO: This function needs much optimising
+G4double GateDiscretizerModule::calculatePitch(G4double crystal_size, G4double spatial_resolution) {
+    // Target pitch size
+    double target_pitch = spatial_resolution / 2.0;
+    double best_pitch = 0;
+    double min_diff = 999;
+    int num_pitches = 1;
+
+    while (true) {
+        double pitch = crystal_size / num_pitches;
+
+        // Check if pitch meets both conditions
+        if (pitch <= target_pitch && std::fmod(crystal_size, pitch) == 0) {
+            double diff = std::fabs(pitch - target_pitch);
+            if (diff <= min_diff) {
+                min_diff = diff;
+                best_pitch = pitch;
+            }
+        }
+
+        // Break if pitch is smaller than a threshold, to prevent infinite loop
+        if ((best_pitch <= target_pitch)&&(best_pitch!=0)) {
+            break;
+        }
+        else if (num_pitches>1000){
+        	std::cout<<"ERROR! No pitch found in 1000 iterations!"<<std::endl;
+        	break;
+        }
+
+        ++num_pitches;
+    }
 
 
+    return best_pitch;
+}
+
+
+void GateDiscretizerModule::SetVirtualIDs( int nBinsX, int nBinsY,int nBinsZ,double pitchX,double pitchY,double pitchZ, G4ThreeVector& pos ){
 
-void GateDiscretizerModule::SetVirtualIDs( int nBinsX, int nBinsY,int nBinsZ, G4ThreeVector& pos ){
-     int crystalID;
-     int subModuleID;
-     int layerID;
 
      int binX,binY,binZ;
 
-	 binX = (pos.X/nBinsX+nBinsX/2);
-	 binY = (pos.Y/nBinsY+nBinsY/2);
-	 binZ = (pos.Z/nBinsY+nBinsY/2);
+
+
+	 binX = static_cast<int>(pos.getX()/pitchX)+static_cast<int>(nBinsX/2);
+	 binY = static_cast<int>(pos.getY()/pitchY)+static_cast<int>(nBinsY/2);
+	 binZ = static_cast<int>(pos.getZ()/pitchZ)+static_cast<int>(nBinsY/2);
+
 
 	 //Change the OutputVolumeID at depths 3,4,5
 	 //(SubmoduleID=binX, crystalID = binY and layerID = binZ)
-	 m_outputDigi->SetOutputVolumeID(3,binX);
-	 m_outputDigi->SetOutputVolumeID(4,binY);
-	 m_outputDigi->SetOutputVolumeID(5,binZ);
+	 m_outputDigi->SetOutputVolumeID(binX,3);
+	 m_outputDigi->SetOutputVolumeID(binY,4);
+	 m_outputDigi->SetOutputVolumeID(binZ,5);
+
+
  }
 
 
 
 
-/*
-void GateDiscretizerModule::DummyMethod1(GateDigi *right)
-{
-	//to copy all variables that are not changed
-	m_outputDigi=right;
 
+void GateDiscretizerModule::SetVirtualID( int nBins, double pitch, G4double pos , int depth){
 
-}
 
-void GateDiscretizerModule::DummyMethod2(GateDigi *right)
-{
-	//to copy all variables that are not changed
-	m_outputDigi=right;
+
+     int bin;
+
+
+	 bin = static_cast<int>(pos/pitch)+static_cast<int>(nBins/2);
+
+
+	 m_outputDigi->SetOutputVolumeID(bin,depth);
+
+ }
+
 
 
-}
 
-*/
 void GateDiscretizerModule::DescribeMyself(size_t indent )
 {
 	if(m_resolution)

diff --git a/source/digits_hits/src/GateDummyDigitizerModule.cc b/source/digits_hits/src/GateDummyDigitizerModule.cc
@@ -46,7 +46,7 @@
    	  	  inputPulse -> inputDigi
 	      outputPulse -> m_outputDigi + correct the first declaration (as in this Dummy module)
 	      outputPulseList.push_back(outputPulse) ->  m_OutputDigiCollection->insert(m_outputDigi);
-	10) Add YourDigitizerModule to GateSinglesDigitizer.cc
+	10) Add YourDigitizerModule to GateSinglesDigitizer
 			- #include "YourDigitizerModule.hh"
 			- in DumpMap() method in
 				static G4String theList = " ...."

diff --git a/source/digits_hits/src/GateSinglesDigitizerMessenger.cc b/source/digits_hits/src/GateSinglesDigitizerMessenger.cc
@@ -127,7 +127,7 @@ const G4String& GateSinglesDigitizerMessenger::DumpMap()
 {
 
 
-   static G4String theList = "readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger intrinsicResolution buffer crosstalk doIModel timeDelay clustering adderComptPhotIdeal gridDiscretizator multipleRejection";
+   static G4String theList = "readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger intrinsicResolution buffer crosstalk doIModel timeDelay clustering adderComptPhotIdeal gridDiscretizator multipleRejection discretizer";
    return theList;
 }
-Original file line number
+Diff line change
@@ Expand Up / @@ -127,7 +127,7 @@ const G4String& GateSinglesDigitizerMessenger::DumpMap() @@
     {
-       static G4String theList = "readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger intrinsicResolution buffer crosstalk doIModel timeDelay clustering adderComptPhotIdeal gridDiscretizator multipleRejection";
+       static G4String theList = "readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger intrinsicResolution buffer crosstalk doIModel timeDelay clustering adderComptPhotIdeal gridDiscretizator multipleRejection discretizer";
        return theList;
     }
@@ Expand Down @@