Merge pull request #700 from Oudihat-Radia/coincidenceSorter

Coincidence sorter :Added checks for global position differences and angle calculations in XY plane to reject invalid coincidences

source/digits_hits/include/GateCoincidenceSorter.hh

@@ -170,6 +170,12 @@ public:
     void SetAcceptancePolicy4CC(const G4String& policy);
 
 
+    G4double GetMinS() const {return m_minS;}
+    G4double GetMaxDeltaZ() const {return m_maxDeltaZ;}
+
+    //! Set the GeometrySelector
+    void SetMinS(G4double val) { m_minS = val;}
+    void SetMaxDeltaZ(G4double val) { m_maxDeltaZ = val;}
 
 
 
@@ -190,7 +196,8 @@ protected:
     acceptance_policy_4CC_t m_acceptance_policy_4CC;    //! <Which is the criteria to accept coincidences in CC sys
     G4bool              m_allDigiOpenCoincGate;        //!< can a digi be part of two coincs?
     G4int               m_depth;                        //!< Depth of system-level for coincidences
-
+    G4double m_minS;  //!< contains the rebirth time.
+    G4double m_maxDeltaZ;  //!< contains the rebirth time.
     G4int coincID_CC;
 
 

source/digits_hits/include/GateCoincidenceSorterMessenger.hh

@@ -69,6 +69,8 @@ private:
     G4UIcmdWithABool            *AllDigiOpenCoincGateCmd;  //!< The UI command "allowMultiples"
     G4UIcmdWithABool            *SetTriggerOnlyByAbsorberCmd;
     G4UIcmdWithABool            *SetEventIDCoincCmd;
+    G4UIcmdWithADoubleAndUnit *minSCmd; //!< set the max time window
+    G4UIcmdWithADoubleAndUnit *maxDeltaZCmd; //!< set the max time window
     GateCoincidenceSorter* m_CoincidenceSorter;
 
 };

source/digits_hits/include/GateCoincidenceTimeDiffSelector.hh

@@ -31,7 +31,9 @@ See LICENSE.md for further details
 class GateCoincidenceTimeDiffSelectorMessenger;
 
 
+
 class GateCoincidenceTimeDiffSelector : public GateVDigitizerModule
+
 {
 public:
 

source/digits_hits/src/GateCoincidenceDigitizerMessenger.cc

@@ -125,6 +125,7 @@ void GateCoincidenceDigitizerMessenger::SetNewValue(G4UIcommand* command,G4Strin
 const G4String& GateCoincidenceDigitizerMessenger::DumpMap()
 {
 
+
    static G4String theList = "deadtime multiplesKiller buffer timeDiffSelector" ;//readout adder energyFraming timeResolution energyResolution spatialResolution efficiency deadtime pileup adderCompton opticaladder noise merger";
 
 
@@ -171,7 +172,6 @@ void GateCoincidenceDigitizerMessenger::DoInsertion(const G4String& childTypeNam
     	  m_CoinDigitizer->AddNewModule(newDM);
       }
 
-
   /*else if (childTypeName=="readout")
   {
 	  newDM = new GateReadout(m_digitizer, DMname);

source/digits_hits/src/GateCoincidenceSorter.cc

@@ -40,6 +40,8 @@ GateCoincidenceSorter::GateCoincidenceSorter(GateDigitizerMgr* itsDigitizerMgr,
     m_offset(0.),
     m_offsetJitter(0.),
     m_minSectorDifference(2),
+	m_minS (-1),
+	m_maxDeltaZ ( -1),
     m_forceMinSecDifferenceToZero(false),	  
     m_multiplesPolicy(kKeepIfAllAreGoods),
     m_allDigiOpenCoincGate(false),
@@ -737,10 +739,37 @@ G4bool GateCoincidenceSorter::IsForbiddenCoincidence(const GateDigi* digi1, cons
       	    G4cout << "[GateCoincidenceSorter::IsForbiddenCoincidence]: coincidence between neighbour blocks --> refused\n";
 	return true;
 	}
+  G4ThreeVector globalPos1 = digi1->GetGlobalPos();
+  G4ThreeVector globalPos2 = digi2->GetGlobalPos();
 
-  return false;
-   }
+  // Check the difference in Z between the two positions
+  if ((m_maxDeltaZ > 0) && (fabs(globalPos2.z() - globalPos1.z()) > m_maxDeltaZ)) {
+      if (nVerboseLevel > 1)
+          G4cout << "[GateCoincidenceSorter::IsForbiddenCoincidence]: difference in Z too large --> refused\n";
+      return true;
+  }
+
+  // Calculate the denominator for distance 's' in the XY plane
+  G4double denom = (globalPos1.y() - globalPos2.y()) * (globalPos1.y() - globalPos2.y()) +
+                   (globalPos2.x() - globalPos1.x()) * (globalPos2.x() - globalPos1.x());
+
+  G4double s = 0.0;
+  if (denom != 0.0) {
+      denom = sqrt(denom);
+      s = (globalPos1.x() * (globalPos1.y() - globalPos2.y()) +
+           globalPos1.y() * (globalPos2.x() - globalPos1.x())) / denom;
+  }
 
+
+  // Check the distance 's' against the maximum threshold
+  if ((m_minS < 0) && (fabs(s) < m_minS)) {
+      if (nVerboseLevel > 1)
+          G4cout << "[GateCoincidenceSorter::IsForbiddenCoincidence]: distance s too large --> refused\n";
+      return true;
+  }
+
+  return false;
+  }
 }
 //------------------------------------------------------------------------------------------------------
 

source/digits_hits/src/GateCoincidenceSorterMessenger.cc

@@ -53,6 +53,16 @@ GateCoincidenceSorterMessenger::GateCoincidenceSorterMessenger(GateCoincidenceSo
   minSectorDiffCmd->SetParameterName("diff",false);
   minSectorDiffCmd->SetRange("diff>=1");
 
+  cmdName = GetDirectoryName() + "setSMin";
+  minSCmd= new G4UIcmdWithADoubleAndUnit(cmdName,this);
+  minSCmd->SetGuidance("Set min S value accepted (<0 --> all is accepted)");
+  minSCmd->SetUnitCategory("Length");
+
+  cmdName = GetDirectoryName() + "setDeltaZMax";
+  maxDeltaZCmd= new G4UIcmdWithADoubleAndUnit(cmdName,this);
+  maxDeltaZCmd->SetGuidance("Set max delta Z value accepted (<0 --> all is accepted)");
+  maxDeltaZCmd->SetUnitCategory("Length");
+
   cmdName = GetDirectoryName()+"forceMinSecDifferenceToZero";
   forceMinSectorDiffCmd = new G4UIcmdWithABool(cmdName,this);
   forceMinSectorDiffCmd->SetGuidance("Force the minimum sector difference for valid coincidences to 0: specsific case for prototype testbench simulations.");
@@ -118,6 +128,8 @@ GateCoincidenceSorterMessenger::~GateCoincidenceSorterMessenger()
     delete SetAcceptancePolicy4CCCmd;
     delete SetEventIDCoincCmd;
     delete forceMinSectorDiffCmd;
+    delete minSCmd;
+    delete maxDeltaZCmd;
 }
 
 
@@ -135,6 +147,10 @@ void GateCoincidenceSorterMessenger::SetNewValue(G4UIcommand* aCommand, G4String
     { m_CoincidenceSorter->SetForcedTo0MinSectorDifference(forceMinSectorDiffCmd->GetNewBoolValue(newValue)); }
   else if( aCommand == minSectorDiffCmd )
     { m_CoincidenceSorter->SetMinSectorDifference(minSectorDiffCmd->GetNewIntValue(newValue)); }
+  else if (aCommand==minSCmd)
+	  m_CoincidenceSorter->SetMinS(minSCmd->GetNewDoubleValue(newValue));
+  else if (aCommand==maxDeltaZCmd)
+	  m_CoincidenceSorter->SetMaxDeltaZ(maxDeltaZCmd->GetNewDoubleValue(newValue));
   else if( aCommand == setDepthCmd )
     { m_CoincidenceSorter->SetDepth(setDepthCmd->GetNewIntValue(newValue)); }
   else if( aCommand == setPresortBufferSizeCmd )