Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc (Version 9.6.p4)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id$
 26 //                                                 27 //
 27 // Authors: G.Depaola & F.Longo                    28 // Authors: G.Depaola & F.Longo
 28 //                                                 29 //
 29 // History:                                        30 // History:
 30 // -------                                     <<  31 // --------
 31 //                                             << 
 32 // 05 Apr 2021   J Allison added quantum entan << 
 33 // If the photons have been "tagged" as "quant << 
 34 // G4eplusAnnihilation for annihilation into 2 << 
 35 // here if - and only if - both photons suffer << 
 36 // predictions from Pryce and Ward, Nature No  << 
 37 // Physical Review 73 (1948) p.440. Experiment << 
 38 // entanglement in the MeV regime and its appl << 
 39 // D. Watts, J. Allison et al., Nature Communi << 
 40 // https://doi.org/10.1038/s41467-021-22907-5. << 
 41 //                                             << 
 42 // 02 May 2009   S Incerti as V. Ivanchenko pr     32 // 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreComptonModel.cc
 43 //                                                 33 //
 44 // Cleanup initialisation and generation of se     34 // Cleanup initialisation and generation of secondaries:
 45 //                  - apply internal high-ener     35 //                  - apply internal high-energy limit only in constructor 
 46 //                  - do not apply low-energy      36 //                  - do not apply low-energy limit (default is 0)
 47 //                  - remove GetMeanFreePath m     37 //                  - remove GetMeanFreePath method and table
 48 //                  - added protection against     38 //                  - added protection against numerical problem in energy sampling 
 49 //                  - use G4ElementSelector        39 //                  - use G4ElementSelector
 50                                                    40 
 51 #include "G4LivermorePolarizedComptonModel.hh"     41 #include "G4LivermorePolarizedComptonModel.hh"
 52 #include "G4PhysicalConstants.hh"                  42 #include "G4PhysicalConstants.hh"
 53 #include "G4SystemOfUnits.hh"                      43 #include "G4SystemOfUnits.hh"
 54 #include "G4AutoLock.hh"                       << 
 55 #include "G4Electron.hh"                       << 
 56 #include "G4ParticleChangeForGamma.hh"         << 
 57 #include "G4LossTableManager.hh"               << 
 58 #include "G4VAtomDeexcitation.hh"              << 
 59 #include "G4AtomicShell.hh"                    << 
 60 #include "G4Gamma.hh"                          << 
 61 #include "G4ShellData.hh"                      << 
 62 #include "G4DopplerProfile.hh"                 << 
 63 #include "G4Log.hh"                            << 
 64 #include "G4Exp.hh"                            << 
 65 #include "G4Pow.hh"                            << 
 66 #include "G4LogLogInterpolation.hh"            << 
 67 #include "G4PhysicsModelCatalog.hh"            << 
 68 #include "G4EntanglementAuxInfo.hh"            << 
 69 #include "G4eplusAnnihilationEntanglementClipB << 
 70                                                    44 
 71 //....oooOO0OOooo........oooOO0OOooo........oo     45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 72                                                    46 
 73 using namespace std;                               47 using namespace std;
 74 namespace { G4Mutex LivermorePolarizedComptonM << 
 75                                                << 
 76                                                << 
 77 G4PhysicsFreeVector* G4LivermorePolarizedCompt << 
 78 G4ShellData*       G4LivermorePolarizedCompton << 
 79 G4DopplerProfile*  G4LivermorePolarizedCompton << 
 80 G4CompositeEMDataSet* G4LivermorePolarizedComp << 
 81                                                    48 
 82 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 83                                                    50 
 84 G4LivermorePolarizedComptonModel::G4LivermoreP <<  51 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
 85   :G4VEmModel(nam),isInitialised(false)        <<  52                                              const G4String& nam)
 86 {                                              <<  53   :G4VEmModel(nam),fParticleChange(0),isInitialised(false),
 87   verboseLevel= 1;                             <<  54    meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
                                                   >>  55 {
                                                   >>  56   lowEnergyLimit = 250 * eV; 
                                                   >>  57   highEnergyLimit = 100 * GeV;
                                                   >>  58   //SetLowEnergyLimit(lowEnergyLimit);
                                                   >>  59   SetHighEnergyLimit(highEnergyLimit);
                                                   >>  60   
                                                   >>  61   verboseLevel= 0;
 88   // Verbosity scale:                              62   // Verbosity scale:
 89   // 0 = nothing                                   63   // 0 = nothing 
 90   // 1 = warning for energy non-conservation       64   // 1 = warning for energy non-conservation 
 91   // 2 = details of energy budget                  65   // 2 = details of energy budget
 92   // 3 = calculation of cross sections, file o     66   // 3 = calculation of cross sections, file openings, sampling of atoms
 93   // 4 = entering in methods                       67   // 4 = entering in methods
 94                                                    68 
 95   if( verboseLevel>1 )                         <<  69   if( verboseLevel>0 ) { 
 96     G4cout << "Livermore Polarized Compton is  <<  70   G4cout << "Livermore Polarized Compton is constructed " << G4endl
 97                                                <<  71          << "Energy range: "
 98   //Mark this model as "applicable" for atomic <<  72          << lowEnergyLimit / eV << " eV - "
 99   SetDeexcitationFlag(true);                   <<  73          << highEnergyLimit / GeV << " GeV"
100                                                <<  74          << G4endl;
101   fParticleChange = nullptr;                   <<  75   }
102   fAtomDeexcitation = nullptr;                 << 
103   fEntanglementModelID = G4PhysicsModelCatalog << 
104 }                                                  76 }
105                                                    77 
106 //....oooOO0OOooo........oooOO0OOooo........oo     78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107                                                    79 
108 G4LivermorePolarizedComptonModel::~G4Livermore     80 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
109 {                                                  81 {  
110   if(IsMaster()) {                             <<  82   if (meanFreePathTable)   delete meanFreePathTable;
111     delete shellData;                          <<  83   if (crossSectionHandler) delete crossSectionHandler;
112     shellData = nullptr;                       <<  84   if (scatterFunctionData) delete scatterFunctionData;
113     delete profileData;                        << 
114     profileData = nullptr;                     << 
115     delete scatterFunctionData;                << 
116     scatterFunctionData = nullptr;             << 
117     for(G4int i=0; i<maxZ; ++i) {              << 
118       if(data[i]) {                            << 
119   delete data[i];                              << 
120   data[i] = nullptr;                           << 
121       }                                        << 
122     }                                          << 
123   }                                            << 
124 }                                                  85 }
125                                                    86 
126 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127                                                    88 
128 void G4LivermorePolarizedComptonModel::Initial     89 void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
129                                        const G     90                                        const G4DataVector& cuts)
130 {                                                  91 {
131   if (verboseLevel > 1)                        <<  92   if (verboseLevel > 3)
132     G4cout << "Calling G4LivermorePolarizedCom     93     G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
133                                                    94 
134   // Initialise element selector               <<  95   if (crossSectionHandler)
135   if(IsMaster()) {                             <<  96   {
136     // Access to elements                      <<  97     crossSectionHandler->Clear();
137     const char* path = G4FindDataDir("G4LEDATA <<  98     delete crossSectionHandler;
138                                                << 
139     G4ProductionCutsTable* theCoupleTable =    << 
140       G4ProductionCutsTable::GetProductionCuts << 
141                                                << 
142     G4int numOfCouples = (G4int)theCoupleTable << 
143                                                << 
144     for(G4int i=0; i<numOfCouples; ++i) {      << 
145       const G4Material* material =             << 
146         theCoupleTable->GetMaterialCutsCouple( << 
147       const G4ElementVector* theElementVector  << 
148       std::size_t nelm = material->GetNumberOf << 
149                                                << 
150       for (std::size_t j=0; j<nelm; ++j) {     << 
151         G4int Z = G4lrint((*theElementVector)[ << 
152         if(Z < 1)        { Z = 1; }            << 
153         else if(Z > maxZ){ Z = maxZ; }         << 
154                                                << 
155         if( (!data[Z]) ) { ReadData(Z, path);  << 
156       }                                        << 
157     }                                          << 
158                                                << 
159     // For Doppler broadening                  << 
160     if(!shellData) {                           << 
161       shellData = new G4ShellData();           << 
162       shellData->SetOccupancyData();           << 
163       G4String file = "/doppler/shell-doppler" << 
164       shellData->LoadData(file);               << 
165     }                                          << 
166     if(!profileData) { profileData = new G4Dop << 
167                                                << 
168     // Scattering Function                     << 
169     if(!scatterFunctionData)                   << 
170       {                                        << 
171                                                << 
172   G4VDataSetAlgorithm* scatterInterpolation =  << 
173   G4String scatterFile = "comp/ce-sf-";        << 
174   scatterFunctionData = new G4CompositeEMDataS << 
175   scatterFunctionData->LoadData(scatterFile);  << 
176       }                                        << 
177                                                << 
178     InitialiseElementSelectors(particle, cuts) << 
179   }                                            << 
180                                                << 
181   if (verboseLevel > 2) {                      << 
182     G4cout << "Loaded cross section files" <<  << 
183   }                                                99   }
                                                   >> 100 
                                                   >> 101   // Reading of data files - all materials are read
184                                                   102   
185   if( verboseLevel>1 ) {                       << 103   crossSectionHandler = new G4CrossSectionHandler;
186     G4cout << "G4LivermoreComptonModel is init << 104   crossSectionHandler->Clear();
187      << "Energy range: "                       << 105   G4String crossSectionFile = "comp/ce-cs-";
188      << LowEnergyLimit() / eV << " eV - "      << 106   crossSectionHandler->LoadData(crossSectionFile);
189      << HighEnergyLimit() / GeV << " GeV"      << 107 
190      << G4endl;                                << 108   meanFreePathTable = 0;
                                                   >> 109   meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
                                                   >> 110 
                                                   >> 111   G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
                                                   >> 112   G4String scatterFile = "comp/ce-sf-";
                                                   >> 113   scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
                                                   >> 114   scatterFunctionData->LoadData(scatterFile);
                                                   >> 115 
                                                   >> 116   // For Doppler broadening
                                                   >> 117   shellData.SetOccupancyData();
                                                   >> 118   G4String file = "/doppler/shell-doppler";
                                                   >> 119   shellData.LoadData(file);
                                                   >> 120 
                                                   >> 121   if (verboseLevel > 2) 
                                                   >> 122     G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
                                                   >> 123 
                                                   >> 124   InitialiseElementSelectors(particle,cuts);
                                                   >> 125 
                                                   >> 126   if(  verboseLevel>0 ) { 
                                                   >> 127     G4cout << "Livermore Polarized Compton model is initialized " << G4endl
                                                   >> 128          << "Energy range: "
                                                   >> 129          << LowEnergyLimit() / eV << " eV - "
                                                   >> 130          << HighEnergyLimit() / GeV << " GeV"
                                                   >> 131          << G4endl;
191   }                                               132   }
192   //                                           << 
193   if(isInitialised) { return; }                << 
194                                                   133   
                                                   >> 134   //
                                                   >> 135     
                                                   >> 136   if(isInitialised) return;
195   fParticleChange = GetParticleChangeForGamma(    137   fParticleChange = GetParticleChangeForGamma();
196   fAtomDeexcitation  = G4LossTableManager::Ins << 
197   isInitialised = true;                           138   isInitialised = true;
198 }                                                 139 }
199                                                   140 
200                                                << 
201 void G4LivermorePolarizedComptonModel::Initial << 
202                 G4VEmModel* masterModel)       << 
203 {                                              << 
204   SetElementSelectors(masterModel->GetElementS << 
205 }                                              << 
206                                                << 
207 //....oooOO0OOooo........oooOO0OOooo........oo << 
208                                                << 
209 void G4LivermorePolarizedComptonModel::ReadDat << 
210 {                                              << 
211   if (verboseLevel > 1)                        << 
212     {                                          << 
213       G4cout << "G4LivermorePolarizedComptonMo << 
214        << G4endl;                              << 
215     }                                          << 
216   if(data[Z]) { return; }                      << 
217   const char* datadir = path;                  << 
218   if(!datadir)                                 << 
219     {                                          << 
220       datadir = G4FindDataDir("G4LEDATA");     << 
221       if(!datadir)                             << 
222   {                                            << 
223     G4Exception("G4LivermorePolarizedComptonMo << 
224           "em0006",FatalException,             << 
225           "Environment variable G4LEDATA not d << 
226     return;                                    << 
227   }                                            << 
228     }                                          << 
229                                                << 
230   data[Z] = new G4PhysicsFreeVector();         << 
231                                                << 
232   std::ostringstream ost;                      << 
233   ost << datadir << "/livermore/comp/ce-cs-" < << 
234   std::ifstream fin(ost.str().c_str());        << 
235                                                << 
236   if( !fin.is_open())                          << 
237     {                                          << 
238       G4ExceptionDescription ed;               << 
239       ed << "G4LivermorePolarizedComptonModel  << 
240    << "> is not opened!" << G4endl;            << 
241       G4Exception("G4LivermoreComptonModel::Re << 
242       "em0003",FatalException,                 << 
243       ed,"G4LEDATA version should be G4EMLOW8. << 
244       return;                                  << 
245     } else {                                   << 
246     if(verboseLevel > 3) {                     << 
247       G4cout << "File " << ost.str()           << 
248        << " is opened by G4LivermorePolarizedC << 
249     }                                          << 
250     data[Z]->Retrieve(fin, true);              << 
251     data[Z]->ScaleVector(MeV, MeV*barn);       << 
252   }                                            << 
253   fin.close();                                 << 
254 }                                              << 
255                                                << 
256 //....oooOO0OOooo........oooOO0OOooo........oo    141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
257                                                   142 
258 G4double G4LivermorePolarizedComptonModel::Com    143 G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
259                                        const G    144                                        const G4ParticleDefinition*,
260                                              G    145                                              G4double GammaEnergy,
261                                              G    146                                              G4double Z, G4double,
262                                              G    147                                              G4double, G4double)
263 {                                                 148 {
264   if (verboseLevel > 3)                           149   if (verboseLevel > 3)
265     G4cout << "Calling ComputeCrossSectionPerA    150     G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
266                                                   151 
267   G4double cs = 0.0;                           << 152   if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
268                                                << 
269   if (GammaEnergy < LowEnergyLimit())          << 
270     return 0.0;                                << 
271                                                   153 
272   G4int intZ = G4lrint(Z);                     << 154   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
273   if(intZ < 1 || intZ > maxZ) { return cs; }   << 
274                                                << 
275   G4PhysicsFreeVector* pv = data[intZ];        << 
276                                                << 
277   // if element was not initialised            << 
278   // do initialisation safely for MT mode      << 
279   if(!pv)                                      << 
280     {                                          << 
281       InitialiseForElement(0, intZ);           << 
282       pv = data[intZ];                         << 
283       if(!pv) { return cs; }                   << 
284     }                                          << 
285                                                << 
286   G4int n = G4int(pv->GetVectorLength() - 1);  << 
287   G4double e1 = pv->Energy(0);                 << 
288   G4double e2 = pv->Energy(n);                 << 
289                                                << 
290   if(GammaEnergy <= e1)      { cs = GammaEnerg << 
291   else if(GammaEnergy <= e2) { cs = pv->Value( << 
292   else if(GammaEnergy > e2)  { cs = pv->Value( << 
293                                                << 
294   return cs;                                      155   return cs;
295                                                << 
296 }                                                 156 }
297                                                   157 
298 //....oooOO0OOooo........oooOO0OOooo........oo    158 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299                                                   159 
300 void G4LivermorePolarizedComptonModel::SampleS    160 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
301                 const G4MaterialCutsCouple* co    161                 const G4MaterialCutsCouple* couple,
302                 const G4DynamicParticle* aDyna    162                 const G4DynamicParticle* aDynamicGamma,
303                 G4double,                         163                 G4double,
304                 G4double)                         164                 G4double)
305 {                                                 165 {
306   // The scattered gamma energy is sampled acc    166   // The scattered gamma energy is sampled according to Klein - Nishina formula.
307   // The random number techniques of Butcher &    167   // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
308   // GEANT4 internal units                        168   // GEANT4 internal units
309   //                                              169   //
310   // Note : Effects due to binding of atomic e    170   // Note : Effects due to binding of atomic electrons are negliged.
311                                                   171 
312   if (verboseLevel > 3)                           172   if (verboseLevel > 3)
313     G4cout << "Calling SampleSecondaries() of     173     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
314                                                   174 
315   G4double gammaEnergy0 = aDynamicGamma->GetKi    175   G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
316                                                << 
317   // do nothing below the threshold            << 
318   // should never get here because the XS is z << 
319   if (gammaEnergy0 < LowEnergyLimit())         << 
320     return ;                                   << 
321                                                << 
322   G4ThreeVector gammaPolarization0 = aDynamicG    176   G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
323                                                   177 
324   // Protection: a polarisation parallel to th    178   // Protection: a polarisation parallel to the
325   // direction causes problems;                   179   // direction causes problems;
326   // in that case find a random polarization      180   // in that case find a random polarization
                                                   >> 181 
327   G4ThreeVector gammaDirection0 = aDynamicGamm    182   G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
328                                                   183 
329   // Make sure that the polarization vector is    184   // Make sure that the polarization vector is perpendicular to the
330   // gamma direction. If not                      185   // gamma direction. If not
                                                   >> 186 
331   if(!(gammaPolarization0.isOrthogonal(gammaDi    187   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
332     { // only for testing now                     188     { // only for testing now
333       gammaPolarization0 = GetRandomPolarizati    189       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
334     }                                             190     }
335   else                                            191   else
336     {                                             192     {
337       if ( gammaPolarization0.howOrthogonal(ga    193       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
338   {                                               194   {
339     gammaPolarization0 = GetPerpendicularPolar    195     gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
340   }                                               196   }
341     }                                             197     }
                                                   >> 198 
342   // End of Protection                            199   // End of Protection
343                                                   200 
                                                   >> 201   // Within energy limit?
                                                   >> 202 
                                                   >> 203   if(gammaEnergy0 <= lowEnergyLimit)
                                                   >> 204     {
                                                   >> 205       fParticleChange->ProposeTrackStatus(fStopAndKill);
                                                   >> 206       fParticleChange->SetProposedKineticEnergy(0.);
                                                   >> 207       fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
                                                   >> 208       return;
                                                   >> 209     }
                                                   >> 210 
344   G4double E0_m = gammaEnergy0 / electron_mass    211   G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
345                                                   212 
346   // Select randomly one element in the curren    213   // Select randomly one element in the current material
347   //G4int Z = crossSectionHandler->SelectRando    214   //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
348   const G4ParticleDefinition* particle =  aDyn    215   const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
349   const G4Element* elm = SelectRandomAtom(coup    216   const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
350   G4int Z = (G4int)elm->GetZ();                   217   G4int Z = (G4int)elm->GetZ();
351                                                   218 
352   // Sample the energy and the polarization of    219   // Sample the energy and the polarization of the scattered photon
                                                   >> 220 
353   G4double epsilon, epsilonSq, onecost, sinThe    221   G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
354                                                   222 
355   G4double epsilon0Local = 1./(1. + 2*E0_m);      223   G4double epsilon0Local = 1./(1. + 2*E0_m);
356   G4double epsilon0Sq = epsilon0Local*epsilon0    224   G4double epsilon0Sq = epsilon0Local*epsilon0Local;
357   G4double alpha1   = - G4Log(epsilon0Local);  << 225   G4double alpha1   = - std::log(epsilon0Local);
358   G4double alpha2 = 0.5*(1.- epsilon0Sq);         226   G4double alpha2 = 0.5*(1.- epsilon0Sq);
359                                                   227 
360   G4double wlGamma = h_Planck*c_light/gammaEne    228   G4double wlGamma = h_Planck*c_light/gammaEnergy0;
361   G4double gammaEnergy1;                          229   G4double gammaEnergy1;
362   G4ThreeVector gammaDirection1;                  230   G4ThreeVector gammaDirection1;
363                                                   231 
364   do {                                            232   do {
365     if ( alpha1/(alpha1+alpha2) > G4UniformRan    233     if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
366       {                                           234       {
367   epsilon   = G4Exp(-alpha1*G4UniformRand());  << 235   epsilon   = std::exp(-alpha1*G4UniformRand());  
368   epsilonSq = epsilon*epsilon;                    236   epsilonSq = epsilon*epsilon; 
369       }                                           237       }
370     else                                          238     else 
371       {                                           239       {
372   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4    240   epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
373   epsilon   = std::sqrt(epsilonSq);               241   epsilon   = std::sqrt(epsilonSq);
374       }                                           242       }
375                                                   243 
376     onecost = (1.- epsilon)/(epsilon*E0_m);       244     onecost = (1.- epsilon)/(epsilon*E0_m);
377     sinThetaSqr   = onecost*(2.-onecost);         245     sinThetaSqr   = onecost*(2.-onecost);
378                                                   246 
379     // Protection                                 247     // Protection
380     if (sinThetaSqr > 1.)                         248     if (sinThetaSqr > 1.)
381       {                                           249       {
382   G4cout                                          250   G4cout
383     << " -- Warning -- G4LivermorePolarizedCom    251     << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
384     << "sin(theta)**2 = "                         252     << "sin(theta)**2 = "
385     << sinThetaSqr                                253     << sinThetaSqr
386     << "; set to 1"                               254     << "; set to 1"
387     << G4endl;                                    255     << G4endl;
388   sinThetaSqr = 1.;                               256   sinThetaSqr = 1.;
389       }                                           257       }
390     if (sinThetaSqr < 0.)                         258     if (sinThetaSqr < 0.)
391       {                                           259       {
392   G4cout                                          260   G4cout
393     << " -- Warning -- G4LivermorePolarizedCom    261     << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
394     << "sin(theta)**2 = "                         262     << "sin(theta)**2 = "
395     << sinThetaSqr                                263     << sinThetaSqr
396     << "; set to 0"                               264     << "; set to 0"
397     << G4endl;                                    265     << G4endl;
398   sinThetaSqr = 0.;                               266   sinThetaSqr = 0.;
399       }                                           267       }
400     // End protection                             268     // End protection
401                                                   269 
402     G4double x =  std::sqrt(onecost/2.) / (wlG    270     G4double x =  std::sqrt(onecost/2.) / (wlGamma/cm);;
403     G4double scatteringFunction = scatterFunct    271     G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
404     greject = (1. - epsilon*sinThetaSqr/(1.+ e    272     greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
405                                                   273 
406   } while(greject < G4UniformRand()*Z);           274   } while(greject < G4UniformRand()*Z);
407                                                   275 
408                                                   276 
409   // *****************************************    277   // ****************************************************
410   //    Phi determination                         278   //    Phi determination
411   // *****************************************    279   // ****************************************************
                                                   >> 280 
412   G4double phi = SetPhi(epsilon,sinThetaSqr);     281   G4double phi = SetPhi(epsilon,sinThetaSqr);
413                                                   282 
414   //                                              283   //
415   // scattered gamma angles. ( Z - axis along     284   // scattered gamma angles. ( Z - axis along the parent gamma)
416   //                                              285   //
                                                   >> 286 
417   G4double cosTheta = 1. - onecost;               287   G4double cosTheta = 1. - onecost;
418                                                   288 
419   // Protection                                   289   // Protection
                                                   >> 290 
420   if (cosTheta > 1.)                              291   if (cosTheta > 1.)
421     {                                             292     {
422       G4cout                                      293       G4cout
423   << " -- Warning -- G4LivermorePolarizedCompt    294   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
424   << "cosTheta = "                                295   << "cosTheta = "
425   << cosTheta                                     296   << cosTheta
426   << "; set to 1"                                 297   << "; set to 1"
427   << G4endl;                                      298   << G4endl;
428       cosTheta = 1.;                              299       cosTheta = 1.;
429     }                                             300     }
430   if (cosTheta < -1.)                             301   if (cosTheta < -1.)
431     {                                             302     {
432       G4cout                                      303       G4cout 
433   << " -- Warning -- G4LivermorePolarizedCompt    304   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
434   << "cosTheta = "                                305   << "cosTheta = " 
435   << cosTheta                                     306   << cosTheta
436   << "; set to -1"                                307   << "; set to -1"
437   << G4endl;                                      308   << G4endl;
438       cosTheta = -1.;                             309       cosTheta = -1.;
439     }                                             310     }
440   // End protection                               311   // End protection      
441                                                << 312   
                                                   >> 313   
442   G4double sinTheta = std::sqrt (sinThetaSqr);    314   G4double sinTheta = std::sqrt (sinThetaSqr);
443                                                   315   
444   // Protection                                   316   // Protection
445   if (sinTheta > 1.)                              317   if (sinTheta > 1.)
446     {                                             318     {
447       G4cout                                      319       G4cout 
448   << " -- Warning -- G4LivermorePolarizedCompt    320   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
449   << "sinTheta = "                                321   << "sinTheta = " 
450   << sinTheta                                     322   << sinTheta
451   << "; set to 1"                                 323   << "; set to 1"
452   << G4endl;                                      324   << G4endl;
453       sinTheta = 1.;                              325       sinTheta = 1.;
454     }                                             326     }
455   if (sinTheta < -1.)                             327   if (sinTheta < -1.)
456     {                                             328     {
457       G4cout                                      329       G4cout 
458   << " -- Warning -- G4LivermorePolarizedCompt    330   << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
459   << "sinTheta = "                                331   << "sinTheta = " 
460   << sinTheta                                     332   << sinTheta
461   << "; set to -1"                                333   << "; set to -1" 
462   << G4endl;                                      334   << G4endl;
463       sinTheta = -1.;                             335       sinTheta = -1.;
464     }                                             336     }
465   // End protection                               337   // End protection
466                                                   338   
467   // Check for entanglement and re-sample phi  << 
468                                                << 
469   const auto* auxInfo                          << 
470   = fParticleChange->GetCurrentTrack()->GetAux << 
471   if (auxInfo) {                               << 
472     const auto* entanglementAuxInfo = dynamic_ << 
473     if (entanglementAuxInfo) {                 << 
474       auto* clipBoard = dynamic_cast<G4eplusAn << 
475       (entanglementAuxInfo->GetEntanglementCli << 
476       if (clipBoard) {                         << 
477         // This is an entangled photon from ep << 
478         // If this is the first scatter of the << 
479         // phi on the clipboard.               << 
480         // If this is the first scatter of the << 
481         // phi of the first scatter of the fir << 
482         // theta of the second photon, to samp << 
483         if (clipBoard->IsTrack1Measurement())  << 
484           // Check we have the relevant track. << 
485           // necessary but I want to be sure t << 
486           // entangled system are properly pai << 
487           // Note: the tracking manager pops t << 
488           // will rely on that.  (If not, the  << 
489           // more complicated to ensure we mat << 
490           // So our track 1 is clipboard track << 
491           if (clipBoard->GetTrackB() == fParti << 
492             // This is the first scatter of th << 
493             //            // Debug             << 
494             //            auto* track1 = fPart << 
495             //            G4cout               << 
496             //            << "This is the firs << 
497             //            << "\nTrack: " << tr << 
498             //            << ", Parent: " << t << 
499             //            << ", Name: " << cli << 
500             //            << G4endl;           << 
501             //            // End debug         << 
502             clipBoard->ResetTrack1Measurement( << 
503             // Store cos(theta),phi of first p << 
504             clipBoard->SetComptonCosTheta1(cos << 
505             clipBoard->SetComptonPhi1(phi);    << 
506           }                                    << 
507         } else if (clipBoard->IsTrack2Measurem << 
508           // Check we have the relevant track. << 
509           // Remember our track 2 is clipboard << 
510           if (clipBoard->GetTrackA() == fParti << 
511             // This is the first scatter of th << 
512             //            // Debug             << 
513             //            auto* track2 = fPart << 
514             //            G4cout               << 
515             //            << "This is the firs << 
516             //            << "\nTrack: " << tr << 
517             //            << ", Parent: " << t << 
518             //            << ", Name: " << cli << 
519             //            << G4endl;           << 
520             //            // End debug         << 
521             clipBoard->ResetTrack2Measurement( << 
522                                                << 
523             // Get cos(theta),phi of first pho << 
524             const G4double& cosTheta1 = clipBo << 
525             const G4double& phi1 = clipBoard-> << 
526             // For clarity make aliases for th << 
527             const G4double& cosTheta2 = cosThe << 
528             G4double& phi2 = phi;              << 
529             //            G4cout << "cosTheta1 << 
530             //            G4cout << "cosTheta2 << 
531                                                << 
532             // Re-sample phi                   << 
533             // Draw the difference of azimutha << 
534             // A + B * cos(2*deltaPhi), or rat << 
535             // C = A / (A + |B|) and D = B / ( << 
536             const G4double sin2Theta1 = 1.-cos << 
537             const G4double sin2Theta2 = 1.-cos << 
538                                                << 
539             // Pryce and Ward, Nature No 4065  << 
540             auto* g4Pow = G4Pow::GetInstance() << 
541             const G4double A =                 << 
542             ((g4Pow->powN(1.-cosTheta1,3))+2.) << 
543             ((g4Pow->powN(2.-cosTheta1,3)*g4Po << 
544             const G4double B = -(sin2Theta1*si << 
545             ((g4Pow->powN(2.-cosTheta1,2)*g4Po << 
546                                                << 
547             //    // Snyder et al, Physical Re << 
548             //    // (This is an alternative f << 
549             //    const G4double& k0 = gammaEn << 
550             //    const G4double k1 = k0/(2.-c << 
551             //    const G4double k2 = k0/(2.-c << 
552             //    const G4double gamma1 = k1/k << 
553             //    const G4double gamma2 = k2/k << 
554             //    const G4double A1 = gamma1*g << 
555             //    const G4double B1 = 2.*sin2T << 
556             //    // That's A1 + B1*sin2(delta << 
557             //    const G4double A = A1 + 0.5* << 
558             //    const G4double B = -0.5*B1;  << 
559                                                << 
560             const G4double maxValue = A + std: << 
561             const G4double C = A / maxValue;   << 
562             const G4double D = B / maxValue;   << 
563             //    G4cout << "A,B,C,D: " << A < << 
564                                                << 
565             // Sample delta phi                << 
566             G4double deltaPhi;                 << 
567             const G4int maxCount = 999999;     << 
568             G4int iCount = 0;                  << 
569             for (; iCount < maxCount; ++iCount << 
570               deltaPhi = twopi * G4UniformRand << 
571               if (G4UniformRand() < C + D * co << 
572             }                                  << 
573             if (iCount >= maxCount ) {         << 
574               G4cout << "G4LivermorePolarizedC << 
575               << "Re-sampled delta phi not fou << 
576               << " tries - carrying on anyway. << 
577             }                                  << 
578                                                << 
579             // Thus, the desired second photon << 
580             phi2 = deltaPhi - phi1 + halfpi;   << 
581             // The minus sign is in above stat << 
582             // annihilation photons are in opp << 
583             // are measured in the opposite di << 
584             // halfpi is added for the followi << 
585             // In this function phi is relativ << 
586             // SystemOfRefChange below. We kno << 
587             // the polarisations of the two an << 
588             // to each other, i.e., halfpi dif << 
589             // Furthermore, only sin(phi) and  << 
590             // need to place any range constra << 
591             //            if (phi2 > pi) {     << 
592             //              phi2 -= twopi;     << 
593             //            }                    << 
594             //            if (phi2 < -pi) {    << 
595             //              phi2 += twopi;     << 
596             //            }                    << 
597           }                                    << 
598         }                                      << 
599       }                                        << 
600     }                                          << 
601   }                                            << 
602                                                << 
603   // End of entanglement                       << 
604                                                   339       
605   G4double dirx = sinTheta*std::cos(phi);         340   G4double dirx = sinTheta*std::cos(phi);
606   G4double diry = sinTheta*std::sin(phi);         341   G4double diry = sinTheta*std::sin(phi);
607   G4double dirz = cosTheta ;                      342   G4double dirz = cosTheta ;
608                                                << 343   
                                                   >> 344 
609   // oneCosT , eom                                345   // oneCosT , eom
610                                                   346 
611   // Doppler broadening -  Method based on:       347   // Doppler broadening -  Method based on:
612   // Y. Namito, S. Ban and H. Hirayama,           348   // Y. Namito, S. Ban and H. Hirayama, 
613   // "Implementation of the Doppler Broadening    349   // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 
614   // NIM A 349, pp. 489-494, 1994                 350   // NIM A 349, pp. 489-494, 1994
615                                                   351   
616   // Maximum number of sampling iterations        352   // Maximum number of sampling iterations
617   static G4int maxDopplerIterations = 1000;    << 353 
                                                   >> 354   G4int maxDopplerIterations = 1000;
618   G4double bindingE = 0.;                         355   G4double bindingE = 0.;
619   G4double photonEoriginal = epsilon * gammaEn    356   G4double photonEoriginal = epsilon * gammaEnergy0;
620   G4double photonE = -1.;                         357   G4double photonE = -1.;
621   G4int iteration = 0;                            358   G4int iteration = 0;
622   G4double eMax = gammaEnergy0;                   359   G4double eMax = gammaEnergy0;
623                                                   360 
624   G4int shellIdx = 0;                          << 
625                                                << 
626   if (verboseLevel > 3) {                      << 
627     G4cout << "Started loop to sample broading << 
628   }                                            << 
629                                                << 
630   do                                              361   do
631     {                                             362     {
632       iteration++;                                363       iteration++;
633       // Select shell based on shell occupancy    364       // Select shell based on shell occupancy
634       shellIdx = shellData->SelectRandomShell( << 365       G4int shell = shellData.SelectRandomShell(Z);
635       bindingE = shellData->BindingEnergy(Z,sh << 366       bindingE = shellData.BindingEnergy(Z,shell);
636                                                   367       
637       if (verboseLevel > 3) {                  << 
638   G4cout << "Shell ID= " << shellIdx           << 
639          << " Ebind(keV)= " << bindingE/keV << << 
640       }                                        << 
641       eMax = gammaEnergy0 - bindingE;             368       eMax = gammaEnergy0 - bindingE;
642                                                << 369      
643       // Randomly sample bound electron moment    370       // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
644       G4double pSample = profileData->RandomSe << 371       G4double pSample = profileData.RandomSelectMomentum(Z,shell);
645                                                << 
646       if (verboseLevel > 3) {                  << 
647        G4cout << "pSample= " << pSample << G4e << 
648      }                                         << 
649       // Rescale from atomic units                372       // Rescale from atomic units
650       G4double pDoppler = pSample * fine_struc    373       G4double pDoppler = pSample * fine_structure_const;
651       G4double pDoppler2 = pDoppler * pDoppler    374       G4double pDoppler2 = pDoppler * pDoppler;
652       G4double var2 = 1. + onecost * E0_m;        375       G4double var2 = 1. + onecost * E0_m;
653       G4double var3 = var2*var2 - pDoppler2;      376       G4double var3 = var2*var2 - pDoppler2;
654       G4double var4 = var2 - pDoppler2 * cosTh    377       G4double var4 = var2 - pDoppler2 * cosTheta;
655       G4double var = var4*var4 - var3 + pDoppl    378       G4double var = var4*var4 - var3 + pDoppler2 * var3;
656       if (var > 0.)                               379       if (var > 0.)
657   {                                               380   {
658     G4double varSqrt = std::sqrt(var);            381     G4double varSqrt = std::sqrt(var);        
659     G4double scale = gammaEnergy0 / var3;         382     G4double scale = gammaEnergy0 / var3;  
660           // Random select either root            383           // Random select either root
661     if (G4UniformRand() < 0.5) photonE = (var4    384     if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;               
662     else photonE = (var4 + varSqrt) * scale;      385     else photonE = (var4 + varSqrt) * scale;
663   }                                               386   } 
664       else                                        387       else
665   {                                               388   {
666     photonE = -1.;                                389     photonE = -1.;
667   }                                               390   }
668    } while ( iteration <= maxDopplerIterations    391    } while ( iteration <= maxDopplerIterations && 
669        (photonE < 0. || photonE > eMax || phot    392        (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
670                                                << 393  
671   // End of recalculation of photon energy wit    394   // End of recalculation of photon energy with Doppler broadening
672   // Revert to original if maximum number of i    395   // Revert to original if maximum number of iterations threshold has been reached
673   if (iteration >= maxDopplerIterations)          396   if (iteration >= maxDopplerIterations)
674     {                                             397     {
675       photonE = photonEoriginal;                  398       photonE = photonEoriginal;
676       bindingE = 0.;                              399       bindingE = 0.;
677     }                                             400     }
678                                                   401 
679   gammaEnergy1 = photonE;                         402   gammaEnergy1 = photonE;
680                                                   403  
681   //                                              404   //
682   // update G4VParticleChange for the scattere    405   // update G4VParticleChange for the scattered photon 
683   //                                              406   //
                                                   >> 407 
                                                   >> 408   //  gammaEnergy1 = epsilon*gammaEnergy0;
                                                   >> 409 
                                                   >> 410 
684   // New polarization                             411   // New polarization
                                                   >> 412 
685   G4ThreeVector gammaPolarization1 = SetNewPol    413   G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
686               sinThetaSqr,                        414               sinThetaSqr,
687               phi,                                415               phi,
688               cosTheta);                          416               cosTheta);
689                                                   417 
690   // Set new direction                            418   // Set new direction
691   G4ThreeVector tmpDirection1( dirx,diry,dirz     419   G4ThreeVector tmpDirection1( dirx,diry,dirz );
692   gammaDirection1 = tmpDirection1;                420   gammaDirection1 = tmpDirection1;
693                                                   421 
694   // Change reference frame.                      422   // Change reference frame.
                                                   >> 423 
695   SystemOfRefChange(gammaDirection0,gammaDirec    424   SystemOfRefChange(gammaDirection0,gammaDirection1,
696         gammaPolarization0,gammaPolarization1)    425         gammaPolarization0,gammaPolarization1);
697                                                   426 
698   if (gammaEnergy1 > 0.)                          427   if (gammaEnergy1 > 0.)
699     {                                             428     {
700       fParticleChange->SetProposedKineticEnerg    429       fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
701       fParticleChange->ProposeMomentumDirectio    430       fParticleChange->ProposeMomentumDirection( gammaDirection1 );
702       fParticleChange->ProposePolarization( ga    431       fParticleChange->ProposePolarization( gammaPolarization1 );
703     }                                             432     }
704   else                                            433   else
705     {                                             434     {
706       gammaEnergy1 = 0.;                          435       gammaEnergy1 = 0.;
707       fParticleChange->SetProposedKineticEnerg    436       fParticleChange->SetProposedKineticEnergy(0.) ;
708       fParticleChange->ProposeTrackStatus(fSto    437       fParticleChange->ProposeTrackStatus(fStopAndKill);
709     }                                             438     }
710                                                   439 
711   //                                              440   //
712   // kinematic of the scattered electron          441   // kinematic of the scattered electron
713   //                                              442   //
                                                   >> 443 
714   G4double ElecKineEnergy = gammaEnergy0 - gam    444   G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
715                                                   445 
716   // SI -protection against negative final ene    446   // SI -protection against negative final energy: no e- is created
717   // like in G4LivermoreComptonModel.cc           447   // like in G4LivermoreComptonModel.cc
718   if(ElecKineEnergy < 0.0) {                      448   if(ElecKineEnergy < 0.0) {
719     fParticleChange->ProposeLocalEnergyDeposit    449     fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
720     return;                                       450     return;
721   }                                               451   }
722                                                   452  
                                                   >> 453   // SI - Removed range test
                                                   >> 454   
723   G4double ElecMomentum = std::sqrt(ElecKineEn    455   G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
724                                                   456 
725   G4ThreeVector ElecDirection((gammaEnergy0 *     457   G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
726            gammaEnergy1 * gammaDirection1) * (    458            gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
727                                                   459 
728   G4DynamicParticle* dp =                      << 
729     new G4DynamicParticle (G4Electron::Electro << 
730   fvect->push_back(dp);                        << 
731                                                << 
732   // sample deexcitation                       << 
733   //                                           << 
734   if (verboseLevel > 3) {                      << 
735     G4cout << "Started atomic de-excitation "  << 
736   }                                            << 
737                                                << 
738   if(fAtomDeexcitation && iteration < maxDoppl << 
739     G4int index = couple->GetIndex();          << 
740     if(fAtomDeexcitation->CheckDeexcitationAct << 
741       std::size_t nbefore = fvect->size();     << 
742       G4AtomicShellEnumerator as = G4AtomicShe << 
743       const G4AtomicShell* shell = fAtomDeexci << 
744       fAtomDeexcitation->GenerateParticles(fve << 
745       std::size_t nafter = fvect->size();      << 
746       if(nafter > nbefore) {                   << 
747   for (std::size_t i=nbefore; i<nafter; ++i) { << 
748     //Check if there is enough residual energy << 
749     if (bindingE >= ((*fvect)[i])->GetKineticE << 
750             {                                  << 
751               //Ok, this is a valid secondary: << 
752         bindingE -= ((*fvect)[i])->GetKineticE << 
753             }                                  << 
754     else                                       << 
755             {                                  << 
756         //Invalid secondary: not enough energy << 
757         //Keep its energy in the local deposit << 
758         delete (*fvect)[i];                    << 
759               (*fvect)[i]=0;                   << 
760             }                                  << 
761   }                                            << 
762       }                                        << 
763     }                                          << 
764   }                                            << 
765   //This should never happen                   << 
766   if(bindingE < 0.0)                           << 
767     G4Exception("G4LivermoreComptonModel::Samp << 
768     "em2050",FatalException,"Negative local en << 
769                                                << 
770   fParticleChange->ProposeLocalEnergyDeposit(b    460   fParticleChange->ProposeLocalEnergyDeposit(bindingE);
                                                   >> 461   
                                                   >> 462   G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
                                                   >> 463   fvect->push_back(dp);
771                                                   464 
772 }                                                 465 }
773                                                   466 
774 //....oooOO0OOooo........oooOO0OOooo........oo    467 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
775                                                   468 
776 G4double G4LivermorePolarizedComptonModel::Set    469 G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
777                G4double sinSqrTh)                 470                G4double sinSqrTh)
778 {                                                 471 {
779   G4double rand1;                                 472   G4double rand1;
780   G4double rand2;                                 473   G4double rand2;
781   G4double phiProbability;                        474   G4double phiProbability;
782   G4double phi;                                   475   G4double phi;
783   G4double a, b;                                  476   G4double a, b;
784                                                   477 
785   do                                              478   do
786     {                                             479     {
787       rand1 = G4UniformRand();                    480       rand1 = G4UniformRand();
788       rand2 = G4UniformRand();                    481       rand2 = G4UniformRand();
789       phiProbability=0.;                          482       phiProbability=0.;
790       phi = twopi*rand1;                          483       phi = twopi*rand1;
791                                                   484       
792       a = 2*sinSqrTh;                             485       a = 2*sinSqrTh;
793       b = energyRate + 1/energyRate;              486       b = energyRate + 1/energyRate;
794                                                   487       
795       phiProbability = 1 - (a/b)*(std::cos(phi    488       phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
                                                   >> 489 
                                                   >> 490       
                                                   >> 491  
796     }                                             492     }
797   while ( rand2 > phiProbability );               493   while ( rand2 > phiProbability );
798   return phi;                                     494   return phi;
799 }                                                 495 }
800                                                   496 
801                                                   497 
802 //....oooOO0OOooo........oooOO0OOooo........oo    498 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
803                                                   499 
804 G4ThreeVector G4LivermorePolarizedComptonModel    500 G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
805 {                                                 501 {
806   G4double dx = a.x();                            502   G4double dx = a.x();
807   G4double dy = a.y();                            503   G4double dy = a.y();
808   G4double dz = a.z();                            504   G4double dz = a.z();
809   G4double x = dx < 0.0 ? -dx : dx;               505   G4double x = dx < 0.0 ? -dx : dx;
810   G4double y = dy < 0.0 ? -dy : dy;               506   G4double y = dy < 0.0 ? -dy : dy;
811   G4double z = dz < 0.0 ? -dz : dz;               507   G4double z = dz < 0.0 ? -dz : dz;
812   if (x < y) {                                    508   if (x < y) {
813     return x < z ? G4ThreeVector(-dy,dx,0) : G    509     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
814   }else{                                          510   }else{
815     return y < z ? G4ThreeVector(dz,0,-dx) : G    511     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
816   }                                               512   }
817 }                                                 513 }
818                                                   514 
819 //....oooOO0OOooo........oooOO0OOooo........oo    515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
820                                                   516 
821 G4ThreeVector G4LivermorePolarizedComptonModel    517 G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
822 {                                                 518 {
823   G4ThreeVector d0 = direction0.unit();           519   G4ThreeVector d0 = direction0.unit();
824   G4ThreeVector a1 = SetPerpendicularVector(d0    520   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
825   G4ThreeVector a0 = a1.unit(); // unit vector    521   G4ThreeVector a0 = a1.unit(); // unit vector
826                                                   522 
827   G4double rand1 = G4UniformRand();               523   G4double rand1 = G4UniformRand();
828                                                   524   
829   G4double angle = twopi*rand1; // random pola    525   G4double angle = twopi*rand1; // random polar angle
830   G4ThreeVector b0 = d0.cross(a0); // cross pr    526   G4ThreeVector b0 = d0.cross(a0); // cross product
831                                                   527   
832   G4ThreeVector c;                                528   G4ThreeVector c;
833                                                   529   
834   c.setX(std::cos(angle)*(a0.x())+std::sin(ang    530   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
835   c.setY(std::cos(angle)*(a0.y())+std::sin(ang    531   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
836   c.setZ(std::cos(angle)*(a0.z())+std::sin(ang    532   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
837                                                   533   
838   G4ThreeVector c0 = c.unit();                    534   G4ThreeVector c0 = c.unit();
839                                                   535 
840   return c0;                                      536   return c0;
                                                   >> 537   
841 }                                                 538 }
842                                                   539 
843 //....oooOO0OOooo........oooOO0OOooo........oo    540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
844                                                   541 
845 G4ThreeVector G4LivermorePolarizedComptonModel    542 G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
846 (const G4ThreeVector& gammaDirection, const G4    543 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
847 {                                                 544 {
                                                   >> 545 
848   //                                              546   // 
849   // The polarization of a photon is always pe    547   // The polarization of a photon is always perpendicular to its momentum direction.
850   // Therefore this function removes those vec    548   // Therefore this function removes those vector component of gammaPolarization, which
851   // points in direction of gammaDirection        549   // points in direction of gammaDirection
852   //                                              550   //
853   // Mathematically we search the projection o    551   // Mathematically we search the projection of the vector a on the plane E, where n is the
854   // plains normal vector.                        552   // plains normal vector.
855   // The basic equation can be found in each g    553   // The basic equation can be found in each geometry book (e.g. Bronstein):
856   // p = a - (a o n)/(n o n)*n                    554   // p = a - (a o n)/(n o n)*n
857                                                   555   
858   return gammaPolarization - gammaPolarization    556   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
859 }                                                 557 }
860                                                   558 
861 //....oooOO0OOooo........oooOO0OOooo........oo    559 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
862                                                   560 
863 G4ThreeVector G4LivermorePolarizedComptonModel    561 G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
864                     G4double sinSqrTh,            562                     G4double sinSqrTh, 
865                     G4double phi,                 563                     G4double phi,
866                     G4double costheta)            564                     G4double costheta) 
867 {                                                 565 {
868   G4double rand1;                                 566   G4double rand1;
869   G4double rand2;                                 567   G4double rand2;
870   G4double cosPhi = std::cos(phi);                568   G4double cosPhi = std::cos(phi);
871   G4double sinPhi = std::sin(phi);                569   G4double sinPhi = std::sin(phi);
872   G4double sinTheta = std::sqrt(sinSqrTh);        570   G4double sinTheta = std::sqrt(sinSqrTh);
873   G4double cosSqrPhi = cosPhi*cosPhi;             571   G4double cosSqrPhi = cosPhi*cosPhi;
874   //  G4double cossqrth = 1.-sinSqrTh;            572   //  G4double cossqrth = 1.-sinSqrTh;
875   //  G4double sinsqrphi = sinPhi*sinPhi;         573   //  G4double sinsqrphi = sinPhi*sinPhi;
876   G4double normalisation = std::sqrt(1. - cosS    574   G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
877                                                   575  
                                                   >> 576 
878   // Determination of Theta                       577   // Determination of Theta 
                                                   >> 578   
                                                   >> 579   // ---- MGP ---- Commented out the following 3 lines to avoid compilation 
                                                   >> 580   // warnings (unused variables)
                                                   >> 581   // G4double thetaProbability;
879   G4double theta;                                 582   G4double theta;
                                                   >> 583   // G4double a, b;
                                                   >> 584   // G4double cosTheta;
                                                   >> 585 
                                                   >> 586   /*
                                                   >> 587 
                                                   >> 588   depaola method
                                                   >> 589   
                                                   >> 590   do
                                                   >> 591   {
                                                   >> 592       rand1 = G4UniformRand();
                                                   >> 593       rand2 = G4UniformRand();
                                                   >> 594       thetaProbability=0.;
                                                   >> 595       theta = twopi*rand1;
                                                   >> 596       a = 4*normalisation*normalisation;
                                                   >> 597       b = (epsilon + 1/epsilon) - 2;
                                                   >> 598       thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
                                                   >> 599       cosTheta = std::cos(theta);
                                                   >> 600     }
                                                   >> 601   while ( rand2 > thetaProbability );
                                                   >> 602   
                                                   >> 603   G4double cosBeta = cosTheta;
                                                   >> 604 
                                                   >> 605   */
                                                   >> 606 
880                                                   607 
881   // Dan Xu method (IEEE TNS, 52, 1160 (2005))    608   // Dan Xu method (IEEE TNS, 52, 1160 (2005))
                                                   >> 609 
882   rand1 = G4UniformRand();                        610   rand1 = G4UniformRand();
883   rand2 = G4UniformRand();                        611   rand2 = G4UniformRand();
884                                                   612 
885   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi    613   if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
886     {                                             614     {
887       if (rand2<0.5)                              615       if (rand2<0.5)
888   theta = pi/2.0;                                 616   theta = pi/2.0;
889       else                                        617       else
890   theta = 3.0*pi/2.0;                             618   theta = 3.0*pi/2.0;
891     }                                             619     }
892   else                                            620   else
893     {                                             621     {
894       if (rand2<0.5)                              622       if (rand2<0.5)
895   theta = 0;                                      623   theta = 0;
896       else                                        624       else
897   theta = pi;                                     625   theta = pi;
898     }                                             626     }
899   G4double cosBeta = std::cos(theta);             627   G4double cosBeta = std::cos(theta);
900   G4double sinBeta = std::sqrt(1-cosBeta*cosBe    628   G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
901                                                   629   
902   G4ThreeVector gammaPolarization1;               630   G4ThreeVector gammaPolarization1;
903                                                   631 
904   G4double xParallel = normalisation*cosBeta;     632   G4double xParallel = normalisation*cosBeta;
905   G4double yParallel = -(sinSqrTh*cosPhi*sinPh    633   G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
906   G4double zParallel = -(costheta*sinTheta*cos    634   G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
907   G4double xPerpendicular = 0.;                   635   G4double xPerpendicular = 0.;
908   G4double yPerpendicular = (costheta)*sinBeta    636   G4double yPerpendicular = (costheta)*sinBeta/normalisation;
909   G4double zPerpendicular = -(sinTheta*sinPhi)    637   G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
910                                                   638 
911   G4double xTotal = (xParallel + xPerpendicula    639   G4double xTotal = (xParallel + xPerpendicular);
912   G4double yTotal = (yParallel + yPerpendicula    640   G4double yTotal = (yParallel + yPerpendicular);
913   G4double zTotal = (zParallel + zPerpendicula    641   G4double zTotal = (zParallel + zPerpendicular);
914                                                   642   
915   gammaPolarization1.setX(xTotal);                643   gammaPolarization1.setX(xTotal);
916   gammaPolarization1.setY(yTotal);                644   gammaPolarization1.setY(yTotal);
917   gammaPolarization1.setZ(zTotal);                645   gammaPolarization1.setZ(zTotal);
918                                                   646   
919   return gammaPolarization1;                      647   return gammaPolarization1;
                                                   >> 648 
920 }                                                 649 }
921                                                   650 
922 //....oooOO0OOooo........oooOO0OOooo........oo    651 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
923                                                   652 
924 void G4LivermorePolarizedComptonModel::SystemO    653 void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
925                 G4ThreeVector& direction1,        654                 G4ThreeVector& direction1,
926                 G4ThreeVector& polarization0,     655                 G4ThreeVector& polarization0,
927                 G4ThreeVector& polarization1)     656                 G4ThreeVector& polarization1)
928 {                                                 657 {
929   // direction0 is the original photon directi    658   // direction0 is the original photon direction ---> z
930   // polarization0 is the original photon pola    659   // polarization0 is the original photon polarization ---> x
931   // need to specify y axis in the real refere    660   // need to specify y axis in the real reference frame ---> y 
932   G4ThreeVector Axis_Z0 = direction0.unit();      661   G4ThreeVector Axis_Z0 = direction0.unit();
933   G4ThreeVector Axis_X0 = polarization0.unit()    662   G4ThreeVector Axis_X0 = polarization0.unit();
934   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_    663   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
935                                                   664 
936   G4double direction_x = direction1.getX();       665   G4double direction_x = direction1.getX();
937   G4double direction_y = direction1.getY();       666   G4double direction_y = direction1.getY();
938   G4double direction_z = direction1.getZ();       667   G4double direction_z = direction1.getZ();
939                                                   668   
940   direction1 = (direction_x*Axis_X0 + directio    669   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
941   G4double polarization_x = polarization1.getX    670   G4double polarization_x = polarization1.getX();
942   G4double polarization_y = polarization1.getY    671   G4double polarization_y = polarization1.getY();
943   G4double polarization_z = polarization1.getZ    672   G4double polarization_z = polarization1.getZ();
944                                                   673 
945   polarization1 = (polarization_x*Axis_X0 + po    674   polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
946                                                   675 
947 }                                                 676 }
948                                                   677 
949                                                   678 
950 //....oooOO0OOooo........oooOO0OOooo........oo << 
951 //....oooOO0OOooo........oooOO0OOooo........oo << 
952                                                << 
953 void                                           << 
954 G4LivermorePolarizedComptonModel::InitialiseFo << 
955                 G4int Z)                       << 
956 {                                              << 
957   G4AutoLock l(&LivermorePolarizedComptonModel << 
958   if(!data[Z]) { ReadData(Z); }                << 
959   l.unlock();                                  << 
960 }                                              << 
961                                                   679