Merge pull request #762 from mraedler/acollinearity_effect_correction

Correction of the acollinearity sampling following Toussaint et al.

Author: kochebina

source/geometry/include/GateDetectorConstruction.hh

@@ -159,7 +159,8 @@ public:
   void SetMaterialIoniPotential(G4String n,G4double v){theListOfIonisationPotential[n]=v;}
   G4double GetMaterialIoniPotential(G4String n){ return theListOfIonisationPotential[n];}
 
-
+  void SetMaterialMeanEnergyPerIonPair(G4String n,G4double v){theListOfMeanEnergyPerIonPair[n]=v;}
+  G4double GetMaterialMeanEnergyPerIonPair(G4String n){ return theListOfMeanEnergyPerIonPair[n];}
 
 private :
 
@@ -185,6 +186,7 @@ protected :
   G4bool moveFlag;
 
   std::map<G4String,G4double> theListOfIonisationPotential;
+  std::map<G4String,G4double> theListOfMeanEnergyPerIonPair;
 
 
 private:

source/geometry/include/GateDetectorMessenger.hh

@@ -58,6 +58,7 @@ class GateDetectorMessenger: public G4UImessenger
 
     G4UIcmdWithoutParameter*   pListCreatorsCmd;
     G4UIcmdWithAString*        IoniCmd;
+    G4UIcmdWithAString*        IonPairCmd;
 
     //G4UIcmdWithABool* 	       pEnableAutoUpdateCmd;    
     //G4UIcmdWithABool* 	       pDisableAutoUpdateCmd; 

source/geometry/src/GateDetectorConstruction.cc

@@ -153,6 +153,22 @@ G4VPhysicalVolume* GateDetectorConstruction::Construct()
 
   const G4MaterialTable * theTable = G4Material::GetMaterialTable();
   for(unsigned int i =0;i<(*theTable).size();i++){
+    // Set the default value
+  	G4double defaultMeanEnergyPerIonPair = 5. * eV;
+  	(*theTable)[i]->GetIonisation()->SetMeanEnergyPerIonPair(defaultMeanEnergyPerIonPair);
+
+    // Update, if given
+  	if(theListOfMeanEnergyPerIonPair[(*theTable)[i]->GetName()]){
+  		(*theTable)[i]->GetIonisation()->SetMeanEnergyPerIonPair(theListOfMeanEnergyPerIonPair[(*theTable)[i]->GetName()]);
+  		GateMessage("Physic", 1, " - " << (*theTable)[i]->GetName() << "\t updating the default 'mean energy per ion pair' value of " <<
+					  G4BestUnit(defaultMeanEnergyPerIonPair,"Energy") << "to " <<
+					  G4BestUnit((*theTable)[i]->GetIonisation()->GetMeanEnergyPerIonPair(),"Energy") << Gateendl);
+  	}
+  	else {
+  		GateMessage("Physic", 1, " - " << (*theTable)[i]->GetName() << "\t using the default 'mean energy per ion pair' value of " <<
+					  G4BestUnit((*theTable)[i]->GetIonisation()->GetMeanEnergyPerIonPair(),"Energy") << Gateendl);
+  	}
+
     if(theListOfIonisationPotential[(*theTable)[i]->GetName()]){
       (*theTable)[i]->GetIonisation()->SetMeanExcitationEnergy(theListOfIonisationPotential[(*theTable)[i]->GetName()]);
       GateMessage("Physic", 1, " - " << (*theTable)[i]->GetName() << "\t defaut value: I = " <<

source/geometry/src/GateDetectorMessenger.cc

@@ -171,8 +171,12 @@ GateDetectorMessenger::GateDetectorMessenger(GateDetectorConstruction* GateDet)
   IoniCmd = new G4UIcmdWithAString(cmd,this);
   IoniCmd->SetGuidance("Set the ionisation potential for a material (two parameters 'material' and 'value and unit')");
 
-
-
+  cmd = "/gate/geometry/setMeanEnergyPerIonPair";
+  IonPairCmd = new G4UIcmdWithAString(cmd,this);
+  IonPairCmd->SetGuidance("Set the mean energy per ion pair. "
+                          "This controls the magnitude of the acollinearity effect for electron-positron annihilation. "
+                          "The default value is 5 eV (see https://www.geant4.org/download/release-notes/notes-v10.7.0.html). "
+                          "(two parameters 'material' and 'value with unit')");
 
   GateDistributionListManager::Init();
 }
@@ -245,6 +249,13 @@ void GateDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newValue)
       GetStringAndValueFromCommand(command, newValue, matName, value);
       pDetectorConstruction->SetMaterialIoniPotential(matName,value);
     }
+  else if( command == IonPairCmd )
+    {
+      G4String matName;
+      double value;
+      GetStringAndValueFromCommand(command, newValue, matName, value);
+      pDetectorConstruction->SetMaterialMeanEnergyPerIonPair(matName,value);
+    }
   else
     G4UImessenger::SetNewValue(command,newValue);
 

source/physics/src/GateBackToBack.cc

@@ -52,12 +52,21 @@ void GateBackToBack::GenerateVertex( G4Event* aEvent, G4bool accolinearityFlag)
       if (accoValue == 0.0){
         accoValue = 0.5*pi / 180;
       }
-    
-      G4double dev = CLHEP::RandGauss::shoot( 0.,accoValue / GateConstants::fwhm_to_sigma );
-      G4double Phi1 = ( twopi * G4UniformRand() )/2. ;
-      
-      G4ThreeVector DirectionPhoton( sin( dev ) * cos( Phi1 ),
-                                     sin( dev ) * sin( Phi1 ), cos( dev ) );
+
+      // Adapt the correction suggested by Toussaint et al. (https://doi.org/10.1088/1361-6560/ad70f1)
+
+      // G4double dev = CLHEP::RandGauss::shoot( 0.,accoValue / GateConstants::fwhm_to_sigma );
+      // G4double Phi1 = ( twopi * G4UniformRand() )/2. ;
+      //
+      // G4ThreeVector DirectionPhoton( sin( dev ) * cos( Phi1 ),
+      //                                sin( dev ) * sin( Phi1 ), cos( dev ) );
+
+      G4double phi = CLHEP::RandGauss::shoot(0., accoValue / GateConstants::fwhm_to_sigma);
+      G4double psi = CLHEP::RandGauss::shoot(0., accoValue / GateConstants::fwhm_to_sigma);
+
+      G4double theta = sqrt(pow(phi, 2.0) + pow(psi, 2.0));
+
+      G4ThreeVector DirectionPhoton(sin(theta) * phi / theta, sin(theta) * psi / theta, cos(theta));
 
       // Scale vector to unit length before rotation, re-scale to original length after rotation:
       G4double gammaMom_mag = gammaMom.mag();