Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedRayleighModel.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/G4LivermorePolarizedRayleighModel.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4LivermorePolarizedRayleighModel.cc (Version 10.7)


  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 //                                                 26 //
 27 // Author: Sebastien Incerti                       27 // Author: Sebastien Incerti
 28 //         30 October 2008                         28 //         30 October 2008
 29 //         on base of G4LowEnergyPolarizedRayl     29 //         on base of G4LowEnergyPolarizedRayleigh developed by R. Capra
 30 //                                                 30 //
 31 // History:                                        31 // History:
 32 // --------                                        32 // --------
 33 // 02 May 2009   S Incerti as V. Ivanchenko pr     33 // 02 May 2009   S Incerti as V. Ivanchenko proposed in G4LivermoreRayleighModel.cc
 34 //                                                 34 //
 35 // Cleanup initialisation and generation of se     35 // Cleanup initialisation and generation of secondaries:
 36 //                  - apply internal high-ener     36 //                  - apply internal high-energy limit only in constructor 
 37 //                  - do not apply low-energy      37 //                  - do not apply low-energy limit (default is 0)
 38 //                  - remove GetMeanFreePath m     38 //                  - remove GetMeanFreePath method and table
 39 //                  - remove initialisation of     39 //                  - remove initialisation of element selector 
 40 //                  - use G4ElementSelector        40 //                  - use G4ElementSelector
 41                                                    41 
 42 #include "G4LivermorePolarizedRayleighModel.hh     42 #include "G4LivermorePolarizedRayleighModel.hh"
 43 #include "G4PhysicalConstants.hh"                  43 #include "G4PhysicalConstants.hh"
 44 #include "G4SystemOfUnits.hh"                      44 #include "G4SystemOfUnits.hh"
 45 #include "G4LogLogInterpolation.hh"                45 #include "G4LogLogInterpolation.hh"
 46 #include "G4CompositeEMDataSet.hh"                 46 #include "G4CompositeEMDataSet.hh"
 47 #include "G4AutoLock.hh"                       << 
 48                                                    47 
 49 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 50                                                    49 
 51 using namespace std;                               50 using namespace std;
 52 namespace { G4Mutex LivermorePolarizedRayleigh << 
 53                                                    51 
 54 //....oooOO0OOooo........oooOO0OOooo........oo     52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 55                                                    53 
 56 G4PhysicsFreeVector* G4LivermorePolarizedRayle <<  54 G4int G4LivermorePolarizedRayleighModel::maxZ = 100;
 57 G4PhysicsFreeVector* G4LivermorePolarizedRayle <<  55 G4LPhysicsFreeVector* G4LivermorePolarizedRayleighModel::dataCS[] = {0};
                                                   >>  56 G4VEMDataSet* G4LivermorePolarizedRayleighModel::formFactorData = 0;
 58                                                    57 
 59 G4LivermorePolarizedRayleighModel::G4Livermore     58 G4LivermorePolarizedRayleighModel::G4LivermorePolarizedRayleighModel(const G4ParticleDefinition*,
 60                    const G4String& nam)            59                    const G4String& nam)
 61   :G4VEmModel(nam),fParticleChange(nullptr),is <<  60   :G4VEmModel(nam),fParticleChange(0),isInitialised(false)
 62 {                                                  61 {
 63   lowEnergyLimit = 250 * CLHEP::eV;            <<  62   fParticleChange =0;
                                                   >>  63   lowEnergyLimit = 250 * eV; 
                                                   >>  64   //SetLowEnergyLimit(lowEnergyLimit);
                                                   >>  65   //SetHighEnergyLimit(highEnergyLimit);
 64   //                                               66   //
 65   verboseLevel= 0;                                 67   verboseLevel= 0;
 66   // Verbosity scale:                              68   // Verbosity scale:
 67   // 0 = nothing                                   69   // 0 = nothing 
 68   // 1 = warning for energy non-conservation       70   // 1 = warning for energy non-conservation 
 69   // 2 = details of energy budget                  71   // 2 = details of energy budget
 70   // 3 = calculation of cross sections, file o     72   // 3 = calculation of cross sections, file openings, sampling of atoms
 71   // 4 = entering in methods                       73   // 4 = entering in methods
 72                                                    74 
 73   if(verboseLevel > 0) {                           75   if(verboseLevel > 0) {
 74     G4cout << "Livermore Polarized Rayleigh is     76     G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
 75            << "Energy range: " << LowEnergyLim <<  77          << "Energy range: "
 76      << HighEnergyLimit() / CLHEP::GeV << " Ge <<  78      << LowEnergyLimit() / eV << " eV - "
 77            << G4endl;                          <<  79      << HighEnergyLimit() / GeV << " GeV"
                                                   >>  80          << G4endl;
 78   }                                                81   }
 79 }                                                  82 }
 80                                                    83 
 81 //....oooOO0OOooo........oooOO0OOooo........oo     84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 82                                                    85 
 83 G4LivermorePolarizedRayleighModel::~G4Livermor     86 G4LivermorePolarizedRayleighModel::~G4LivermorePolarizedRayleighModel()
 84 {                                                  87 {  
 85   if(IsMaster()) {                             <<  88  if(IsMaster()) {
 86     for(G4int i=0; i<maxZ; ++i) {              <<  89    for(G4int i=0; i<maxZ; ++i) {
 87       if(dataCS[i]) {                          <<  90      if(dataCS[i]) { 
 88   delete dataCS[i];                            <<  91        delete dataCS[i];
 89   dataCS[i] = nullptr;                         <<  92        dataCS[i] = 0;
 90   delete formFactorData[i];                    <<  93      }
 91   formFactorData[i] = nullptr;                 <<  94    }
 92       }                                        <<  95    delete formFactorData;
 93     }                                          <<  96    formFactorData = 0; 
 94   }                                            <<  97    
                                                   >>  98  }
 95 }                                                  99 }
 96                                                   100 
 97 //....oooOO0OOooo........oooOO0OOooo........oo    101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 98                                                   102 
 99 void G4LivermorePolarizedRayleighModel::Initia    103 void G4LivermorePolarizedRayleighModel::Initialise(const G4ParticleDefinition* particle,
100                                        const G    104                                        const G4DataVector& cuts)
101 {                                                 105 {
102   // Rayleigh process:                      Th << 106 // Rayleigh process:                      The Quantum Theory of Radiation
103   //                                        W. << 107 //                                        W. Heitler,       Oxford at the Clarendon Press, Oxford (1954)                                                 
104   // Scattering function:                   A  << 108 // Scattering function:                   A simple model of photon transport
105   //                                        D. << 109 //                                        D.E. Cullen,      Nucl. Instr. Meth. in Phys. Res. B 101 (1995) 499-510                                       
106   // Polarization of the outcoming photon:  Be << 110 // Polarization of the outcoming photon:  Beam test of a prototype detector array for the PoGO astronomical hard X-ray/soft gamma-ray polarimeter
107   //                                        X- << 111 //                                        T. Mizuno et al., Nucl. Instr. Meth. in Phys. Res. A 540 (2005) 158-168                                        
108   //                                        T. << 112 
109   if (verboseLevel > 3)                           113   if (verboseLevel > 3)
110     G4cout << "Calling G4LivermorePolarizedRay    114     G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl;
111                                                   115 
                                                   >> 116 
112   if(IsMaster()) {                                117   if(IsMaster()) {
113                                                   118     
                                                   >> 119     // Form Factor 
                                                   >> 120     
                                                   >> 121     G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
                                                   >> 122     G4String formFactorFile = "rayl/re-ff-";
                                                   >> 123     formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
                                                   >> 124     formFactorData->LoadData(formFactorFile);
                                                   >> 125     
114     // Initialise element selector                126     // Initialise element selector
115     InitialiseElementSelectors(particle, cuts)    127     InitialiseElementSelectors(particle, cuts);
116                                                   128     
117     // Access to elements                         129     // Access to elements
118     const char* path = G4FindDataDir("G4LEDATA << 130     char* path = std::getenv("G4LEDATA");
119     auto elmTable = G4Element::GetElementTable << 131     G4ProductionCutsTable* theCoupleTable =
120     for (auto const & elm : *elmTable) {       << 132       G4ProductionCutsTable::GetProductionCutsTable();
121       G4int Z = std::min(elm->GetZasInt(), max << 133     G4int numOfCouples = theCoupleTable->GetTableSize();
122       if( nullptr == dataCS[Z] ) { ReadData(Z, << 134    
123     }                                          << 135      for(G4int i=0; i<numOfCouples; ++i) 
                                                   >> 136        {
                                                   >> 137    const G4MaterialCutsCouple* couple = 
                                                   >> 138      theCoupleTable->GetMaterialCutsCouple(i);
                                                   >> 139    const G4Material* material = couple->GetMaterial();
                                                   >> 140    const G4ElementVector* theElementVector = material->GetElementVector();
                                                   >> 141    G4int nelm = material->GetNumberOfElements();
                                                   >> 142    
                                                   >> 143    for (G4int j=0; j<nelm; ++j) 
                                                   >> 144      {
                                                   >> 145        G4int Z = G4lrint((*theElementVector)[j]->GetZ());
                                                   >> 146        if(Z < 1)          { Z = 1; }
                                                   >> 147        else if(Z > maxZ)  { Z = maxZ; }
                                                   >> 148        if( (!dataCS[Z]) ) { ReadData(Z, path); }
                                                   >> 149      }
                                                   >> 150        }
124   }                                               151   }
125                                                   152   
126   if(isInitialised) { return; }                   153   if(isInitialised) { return; }
127   fParticleChange = GetParticleChangeForGamma(    154   fParticleChange = GetParticleChangeForGamma();
128   isInitialised = true;                           155   isInitialised = true;
                                                   >> 156   
129 }                                                 157 }
130                                                   158 
131                                                   159 
132 //....oooOO0OOooo........oooOO0OOooo........oo    160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
133                                                   161 
134 void G4LivermorePolarizedRayleighModel::Initia << 162 void G4LivermorePolarizedRayleighModel::InitialiseLocal(const G4ParticleDefinition*,
135                 const G4ParticleDefinition*, G << 163                 G4VEmModel* masterModel)
136 {                                                 164 {
137   SetElementSelectors(masterModel->GetElementS    165   SetElementSelectors(masterModel->GetElementSelectors());
138 }                                                 166 }
139                                                   167 
140 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141                                                   169 
142 void G4LivermorePolarizedRayleighModel::ReadDa << 170 void G4LivermorePolarizedRayleighModel::ReadData(size_t Z, const char* path)
143 {                                                 171 {
144   if (verboseLevel > 1) {                      << 172   if (verboseLevel > 1) 
145     G4cout << "Calling ReadData() of G4Livermo << 173     {
146      << G4endl;                                << 174       G4cout << "Calling ReadData() of G4LivermoreRayleighModel" 
147   }                                            << 175        << G4endl;
                                                   >> 176     }
148                                                   177   
149   if(nullptr != dataCS[Z]) { return; }         << 178   if(dataCS[Z]) { return; }
150                                                   179   
151   const char* datadir = path;                     180   const char* datadir = path;
152                                                   181   
153   if(nullptr == datadir)                       << 182   if(!datadir) 
154     {                                             183     {
155       datadir = G4FindDataDir("G4LEDATA");     << 184       datadir = std::getenv("G4LEDATA");
156       if(nullptr == datadir)                   << 185       if(!datadir) 
157   {                                               186   {
158     G4Exception("G4LivermoreRayleighModelModel    187     G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
159           FatalException,                         188           FatalException,
160           "Environment variable G4LEDATA not d    189           "Environment variable G4LEDATA not defined");
161     return;                                       190     return;
162   }                                               191   }
163     }                                             192     }
164   dataCS[Z] = new G4PhysicsFreeVector();       << 193   
165   formFactorData[Z] = new G4PhysicsFreeVector( << 194   //
                                                   >> 195   
                                                   >> 196   dataCS[Z] = new G4LPhysicsFreeVector();
                                                   >> 197   
                                                   >> 198   // Activation of spline interpolation
                                                   >> 199   //dataCS[Z] ->SetSpline(true);
166                                                   200   
167   std::ostringstream ostCS;                       201   std::ostringstream ostCS;
168   ostCS << datadir << "/livermore/rayl/re-cs-"    202   ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
169   std::ifstream finCS(ostCS.str().c_str());       203   std::ifstream finCS(ostCS.str().c_str());
170                                                   204   
171   if( !finCS .is_open() ) {                    << 205   if( !finCS .is_open() ) 
172     G4ExceptionDescription ed;                 << 206     {
173     ed << "G4LivermorePolarizedRayleighModel d << 207      G4ExceptionDescription ed;
174        << "> is not opened!" << G4endl;        << 208      ed << "G4LivermorePolarizedRayleighModel data file <" << ostCS.str().c_str()
175     G4Exception("G4LivermorePolarizedRayleighM << 209         << "> is not opened!" << G4endl;
176     FatalException,                            << 210      G4Exception("G4LivermorePolarizedRayleighModel::ReadData()","em0003",FatalException,
177     ed,"G4LEDATA version should be G4EMLOW8.0  << 211      ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
178     return;                                    << 212      return;
179   } else {                                     << 213    } 
180     if(verboseLevel > 3) {                     << 214    else 
181       G4cout << "File " << ostCS.str()         << 215    {
182        << " is opened by G4LivermoreRayleighMo << 216      if(verboseLevel > 3) { 
183     }                                          << 217        G4cout << "File " << ostCS.str() 
184     dataCS[Z]->Retrieve(finCS, true);          << 218         << " is opened by G4LivermoreRayleighModel" << G4endl;
185   }                                            << 219      }
186                                                << 220      dataCS[Z]->Retrieve(finCS, true);
187   std::ostringstream ostFF;                    << 221    } 
188   ostFF << datadir << "/livermore/rayl/re-ff-" << 222  }
189   std::ifstream finFF(ostFF.str().c_str());    << 
190                                                << 
191   if( !finFF.is_open() ) {                     << 
192     G4ExceptionDescription ed;                 << 
193     ed << "G4LivermorePolarizedRayleighModel d << 
194        << "> is not opened!" << G4endl;        << 
195     G4Exception("G4LivermorePolarizedRayleighM << 
196     FatalException,                            << 
197     ed,"G4LEDATA version should be G4EMLOW8.0  << 
198     return;                                    << 
199   } else {                                     << 
200     if(verboseLevel > 3) {                     << 
201       G4cout << "File " << ostFF.str()         << 
202                << " is opened by G4LivermoreRa << 
203     }                                          << 
204     formFactorData[Z]->Retrieve(finFF, true);  << 
205   }                                            << 
206 }                                              << 
207                                                   223  
208 //....oooOO0OOooo........oooOO0OOooo........oo << 224  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209                                                   225  
210 G4double G4LivermorePolarizedRayleighModel::Co << 226  G4double G4LivermorePolarizedRayleighModel::ComputeCrossSectionPerAtom(
211                                         const     227                                         const G4ParticleDefinition*,
212           G4double GammaEnergy,                   228           G4double GammaEnergy,
213           G4double Z, G4double,                   229           G4double Z, G4double,
214           G4double, G4double)                     230           G4double, G4double)
215 {                                              << 231  {
216   if (verboseLevel > 1) {                      << 232    if (verboseLevel > 1) 
217     G4cout << "G4LivermoreRayleighModel::Compu << 233      {
218      << G4endl;                                << 234        G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()" 
219   }                                            << 235         << G4endl;
                                                   >> 236      }
220                                                   237  
221   if(GammaEnergy < lowEnergyLimit) { return 0. << 238    if(GammaEnergy < lowEnergyLimit) { return 0.0; }
                                                   >> 239    
                                                   >> 240    G4double xs = 0.0;
222                                                   241    
223   G4double xs = 0.0;                           << 242    G4int intZ = G4lrint(Z);
224                                                   243    
225   G4int intZ = G4lrint(Z);                     << 244    if(intZ < 1 || intZ > maxZ) { return xs; }
226   if(intZ < 1 || intZ > maxZ) { return xs; }   << 
227                                                   245    
228   G4PhysicsFreeVector* pv = dataCS[intZ];      << 246    G4LPhysicsFreeVector* pv = dataCS[intZ];
229                                                   247  
230   // if element was not initialised            << 248    // if element was not initialised
231   // do initialisation safely for MT mode      << 249    // do initialisation safely for MT mode
232   if(nullptr == pv) {                          << 250    if(!pv) { 
233     InitialiseForElement(0, intZ);             << 251      InitialiseForElement(0, intZ);
234     pv = dataCS[intZ];                         << 252      pv = dataCS[intZ];
235     if(nullptr == pv) { return xs; }           << 253      if(!pv) { return xs; }
236   }                                            << 254    }
                                                   >> 255  
                                                   >> 256    G4int n = pv->GetVectorLength() - 1;
                                                   >> 257    G4double e = GammaEnergy/MeV;
                                                   >> 258    if(e >= pv->Energy(n)) {
                                                   >> 259      xs = (*pv)[n]/(e*e);  
                                                   >> 260    } else if(e >= pv->Energy(0)) {
                                                   >> 261      xs = pv->Value(e)/(e*e);  
                                                   >> 262    }
237                                                   263  
238   G4int n = G4int(pv->GetVectorLength() - 1);  << 264    /*   if(verboseLevel > 0)
239   G4double e = GammaEnergy/MeV;                << 265   {
240   if(e >= pv->Energy(n)) {                     << 266   G4cout  <<  "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" 
241     xs = (*pv)[n]/(e*e);                       << 267      << e << G4endl;
242   } else if(e >= pv->Energy(0)) {              << 268        G4cout  <<  "  cs (Geant4 internal unit)=" << xs << G4endl;
243     xs = pv->Value(e)/(e*e);                   << 269      G4cout  <<  "    -> first E*E*cs value in CS data file (iu) =" << (*pv)[0] 
244   }                                            << 270      << G4endl;
245   return xs;                                   << 271        G4cout  <<  "    -> last  E*E*cs value in CS data file (iu) =" << (*pv)[n] 
246 }                                              << 272      << G4endl;
                                                   >> 273        G4cout  <<  "*********************************************************" 
                                                   >> 274      << G4endl;
                                                   >> 275        }
                                                   >> 276    */
                                                   >> 277    
                                                   >> 278    return xs;
                                                   >> 279  }
247                                                   280 
248 //....oooOO0OOooo........oooOO0OOooo........oo    281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
249                                                   282 
250 void G4LivermorePolarizedRayleighModel::Sample << 283 void G4LivermorePolarizedRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
251                 std::vector<G4DynamicParticle* << 284                 const G4MaterialCutsCouple* couple,
252     const G4MaterialCutsCouple* couple,        << 285                 const G4DynamicParticle* aDynamicGamma,
253     const G4DynamicParticle* aDynamicGamma,    << 286                 G4double,
254     G4double, G4double)                        << 287                 G4double)
255 {                                                 288 {
256   if (verboseLevel > 3)                           289   if (verboseLevel > 3)
257     G4cout << "Calling SampleSecondaries() of     290     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedRayleighModel" << G4endl;
258                                                   291 
259   G4double photonEnergy0 = aDynamicGamma->GetK    292   G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
260                                                   293   
261   if (photonEnergy0 <= lowEnergyLimit)            294   if (photonEnergy0 <= lowEnergyLimit)
262   {                                               295   {
263     fParticleChange->ProposeTrackStatus(fStopA << 296       fParticleChange->ProposeTrackStatus(fStopAndKill);
264     fParticleChange->SetProposedKineticEnergy( << 297       fParticleChange->SetProposedKineticEnergy(0.);
265     fParticleChange->ProposeLocalEnergyDeposit << 298       fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
266     return;                                    << 299       return ;
267   }                                               300   }
268                                                   301 
269   G4ParticleMomentum photonDirection0 = aDynam    302   G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
270                                                   303 
271   // Select randomly one element in the curren    304   // Select randomly one element in the current material
                                                   >> 305   // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
272   const G4ParticleDefinition* particle =  aDyn    306   const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
273   const G4Element* elm = SelectRandomAtom(coup    307   const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
274   G4int Z = elm->GetZasInt();                  << 308   G4int Z = (G4int)elm->GetZ();
275                                                   309 
276   G4double outcomingPhotonCosTheta = GenerateC    310   G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z);
277   G4double outcomingPhotonPhi = GeneratePhi(ou    311   G4double outcomingPhotonPhi = GeneratePhi(outcomingPhotonCosTheta);
278   G4double beta = GeneratePolarizationAngle(); << 312   G4double beta=GeneratePolarizationAngle();
279                                                   313  
280   // incomingPhoton reference frame:              314   // incomingPhoton reference frame:
281   // z = versor parallel to the incomingPhoton    315   // z = versor parallel to the incomingPhotonDirection
282   // x = versor parallel to the incomingPhoton    316   // x = versor parallel to the incomingPhotonPolarization
283   // y = defined as z^x                           317   // y = defined as z^x
284                                                   318  
285   // outgoingPhoton reference frame:              319   // outgoingPhoton reference frame:
286   // z' = versor parallel to the outgoingPhoto    320   // z' = versor parallel to the outgoingPhotonDirection
287   // x' = defined as x-x*z'z' normalized          321   // x' = defined as x-x*z'z' normalized
288   // y' = defined as z'^x'                     << 322   // y' = defined as z'^x'
                                                   >> 323  
289   G4ThreeVector z(aDynamicGamma->GetMomentumDi    324   G4ThreeVector z(aDynamicGamma->GetMomentumDirection().unit()); 
290   G4ThreeVector x(GetPhotonPolarization(*aDyna    325   G4ThreeVector x(GetPhotonPolarization(*aDynamicGamma));
291   G4ThreeVector y(z.cross(x));                    326   G4ThreeVector y(z.cross(x));
292                                                   327  
293   // z' = std::cos(phi)*std::sin(theta)        << 328   // z' = std::cos(phi)*std::sin(theta) x + std::sin(phi)*std::sin(theta) y + std::cos(theta) z
294   // x + std::sin(phi)*std::sin(theta) y + std << 
295   G4double xDir;                                  329   G4double xDir;
296   G4double yDir;                                  330   G4double yDir;
297   G4double zDir;                                  331   G4double zDir;
298   zDir=outcomingPhotonCosTheta;                   332   zDir=outcomingPhotonCosTheta;
299   xDir=std::sqrt(1-outcomingPhotonCosTheta*out    333   xDir=std::sqrt(1-outcomingPhotonCosTheta*outcomingPhotonCosTheta);
300   yDir=xDir;                                      334   yDir=xDir;
301   xDir*=std::cos(outcomingPhotonPhi);             335   xDir*=std::cos(outcomingPhotonPhi);
302   yDir*=std::sin(outcomingPhotonPhi);             336   yDir*=std::sin(outcomingPhotonPhi);
303                                                   337  
304   G4ThreeVector zPrime((xDir*x + yDir*y + zDir    338   G4ThreeVector zPrime((xDir*x + yDir*y + zDir*z).unit());
305   G4ThreeVector xPrime(x.perpPart(zPrime).unit    339   G4ThreeVector xPrime(x.perpPart(zPrime).unit());
306   G4ThreeVector yPrime(zPrime.cross(xPrime));     340   G4ThreeVector yPrime(zPrime.cross(xPrime));
307                                                   341  
308   // outgoingPhotonPolarization is directed as << 342   // outgoingPhotonPolarization is directed as x' std::cos(beta) + y' std::sin(beta)
309   // x' std::cos(beta) + y' std::sin(beta)     << 
310   G4ThreeVector outcomingPhotonPolarization(xP    343   G4ThreeVector outcomingPhotonPolarization(xPrime*std::cos(beta) + yPrime*std::sin(beta));
311                                                   344  
312   fParticleChange->ProposeMomentumDirection(zP    345   fParticleChange->ProposeMomentumDirection(zPrime);
313   fParticleChange->ProposePolarization(outcomi    346   fParticleChange->ProposePolarization(outcomingPhotonPolarization);
314   fParticleChange->SetProposedKineticEnergy(ph    347   fParticleChange->SetProposedKineticEnergy(photonEnergy0); 
                                                   >> 348 
315 }                                                 349 }
316                                                   350 
317 //....oooOO0OOooo........oooOO0OOooo........oo    351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
318                                                   352 
319 G4double G4LivermorePolarizedRayleighModel::Ge << 353 G4double G4LivermorePolarizedRayleighModel::GenerateCosTheta(G4double incomingPhotonEnergy, G4int zAtom) const
320 {                                                 354 {
321   //  d sigma                                     355   //  d sigma                                                                    k0
322   // --------- =  r0^2 * pi * F^2(x, Z) * ( 2     356   // --------- =  r0^2 * pi * F^2(x, Z) * ( 2 - sin^2 theta) * std::sin (theta), x = ---- std::sin(theta/2)
323   //  d theta                                     357   //  d theta                                                                    hc
324                                                   358  
325   //  d sigma                                     359   //  d sigma                                             k0          1 - y
326   // --------- = r0^2 * pi * F^2(x, Z) * ( 1 +    360   // --------- = r0^2 * pi * F^2(x, Z) * ( 1 + y^2), x = ---- std::sqrt ( ------- ), y = std::cos(theta)
327   //    d y                                       361   //    d y                                               hc            2
328                                                   362 
329   //              Z                               363   //              Z
330   // F(x, Z) ~ --------                           364   // F(x, Z) ~ --------
331   //            a + bx                            365   //            a + bx
332   //                                              366   //
333   // The time to exit from the outer loop grow    367   // The time to exit from the outer loop grows as ~ k0
334   // On pcgeant2 the time is ~ 1 s for k0 ~ 1     368   // On pcgeant2 the time is ~ 1 s for k0 ~ 1 MeV on the oxygen element. A 100 GeV
335   // event will take ~ 10 hours.                  369   // event will take ~ 10 hours.
336   //                                              370   //
337   // On the avarage the inner loop does 1.5 it    371   // On the avarage the inner loop does 1.5 iterations before exiting
338   const G4double xxfact = CLHEP::cm/(CLHEP::h_ << 372  
339   const G4double xFactor = incomingPhotonEnerg << 373   const G4double xFactor = (incomingPhotonEnergy*cm)/(h_Planck*c_light);
                                                   >> 374   //const G4VEMDataSet * formFactorData = GetScatterFunctionData();
340                                                   375 
341   G4double cosTheta;                              376   G4double cosTheta;
342   G4double fCosTheta;                             377   G4double fCosTheta;
343   G4double x;                                     378   G4double x;
344   G4double fValue;                                379   G4double fValue;
345                                                   380 
346   if (incomingPhotonEnergy > 5.*CLHEP::MeV)    << 381   if (incomingPhotonEnergy > 5.*MeV)
347   {                                               382   {
348     cosTheta = 1.;                                383     cosTheta = 1.;
349   }                                               384   }
350   else                                            385   else
351   {                                               386   {
352     do                                            387     do
353     {                                             388     {
354       do                                          389       do
355   {                                               390   {
356     cosTheta = 2.*G4UniformRand()-1.;             391     cosTheta = 2.*G4UniformRand()-1.;
357     fCosTheta = (1.+cosTheta*cosTheta)/2.;        392     fCosTheta = (1.+cosTheta*cosTheta)/2.;
358   }                                               393   }
359       while (fCosTheta < G4UniformRand());        394       while (fCosTheta < G4UniformRand());
360                                                   395   
361       x = xFactor*std::sqrt((1.-cosTheta)/2.);    396       x = xFactor*std::sqrt((1.-cosTheta)/2.);
362                                                   397   
363       if (x > 1.e+005)                            398       if (x > 1.e+005)
364   fValue = formFactorData[Z]->Value(x);        << 399   fValue = formFactorData->FindValue(x, zAtom-1);
365       else                                        400       else
366   fValue = formFactorData[Z]->Value(0.);       << 401   fValue = formFactorData->FindValue(0., zAtom-1);
367                                                   402    
368       fValue /= Z;                             << 403       fValue/=zAtom;
369       fValue *= fValue;                        << 404       fValue*=fValue;
370     }                                             405     }
371     while(fValue < G4UniformRand());              406     while(fValue < G4UniformRand());
372   }                                               407   }
373                                                   408 
374   return cosTheta;                                409   return cosTheta;
375 }                                                 410 }
376                                                   411 
377 //....oooOO0OOooo........oooOO0OOooo........oo    412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
378                                                   413 
379 G4double G4LivermorePolarizedRayleighModel::Ge    414 G4double G4LivermorePolarizedRayleighModel::GeneratePhi(G4double cosTheta) const
380 {                                                 415 {
381   //  d sigma                                     416   //  d sigma
382   // --------- = alpha * ( 1 - sin^2 (theta) *    417   // --------- = alpha * ( 1 - sin^2 (theta) * cos^2 (phi) )
383   //   d phi                                      418   //   d phi
384                                                   419  
385   // On the average the loop takes no more tha    420   // On the average the loop takes no more than 2 iterations before exiting 
386                                                   421 
387   G4double phi;                                   422   G4double phi;
388   G4double cosPhi;                                423   G4double cosPhi;
389   G4double phiProbability;                        424   G4double phiProbability;
390   G4double sin2Theta;                             425   G4double sin2Theta;
391                                                   426  
392   sin2Theta=1.-cosTheta*cosTheta;                 427   sin2Theta=1.-cosTheta*cosTheta;
393                                                   428  
394   do                                              429   do
395     {                                             430     {
396       phi = CLHEP::twopi * G4UniformRand();    << 431       phi = twopi * G4UniformRand();
397       cosPhi = std::cos(phi);                     432       cosPhi = std::cos(phi);
398       phiProbability= 1. - sin2Theta*cosPhi*co    433       phiProbability= 1. - sin2Theta*cosPhi*cosPhi;
399     }                                             434     }
400   while (phiProbability < G4UniformRand());       435   while (phiProbability < G4UniformRand());
401                                                   436  
402   return phi;                                     437   return phi;
403 }                                                 438 }
404                                                   439 
405 //....oooOO0OOooo........oooOO0OOooo........oo    440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
406                                                   441 
407 G4double G4LivermorePolarizedRayleighModel::Ge    442 G4double G4LivermorePolarizedRayleighModel::GeneratePolarizationAngle(void) const
408 {                                                 443 {
409   // Rayleigh polarization is always on the x'    444   // Rayleigh polarization is always on the x' direction
                                                   >> 445 
410   return 0;                                       446   return 0;
411 }                                                 447 }
412                                                   448 
413 //....oooOO0OOooo........oooOO0OOooo........oo    449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
414                                                   450 
415 G4ThreeVector G4LivermorePolarizedRayleighMode    451 G4ThreeVector G4LivermorePolarizedRayleighModel::GetPhotonPolarization(const G4DynamicParticle&  photon)
416 {                                                 452 {
417   // From G4VLowEnergyDiscretePhotonProcess.cc << 453 
                                                   >> 454 // SI - From G4VLowEnergyDiscretePhotonProcess.cc
                                                   >> 455  
418   G4ThreeVector photonMomentumDirection;          456   G4ThreeVector photonMomentumDirection;
419   G4ThreeVector photonPolarization;               457   G4ThreeVector photonPolarization;
420                                                   458 
421   photonPolarization = photon.GetPolarization(    459   photonPolarization = photon.GetPolarization(); 
422   photonMomentumDirection = photon.GetMomentum    460   photonMomentumDirection = photon.GetMomentumDirection();
423                                                   461 
424   if ((!photonPolarization.isOrthogonal(photon    462   if ((!photonPolarization.isOrthogonal(photonMomentumDirection, 1e-6)) || photonPolarization.mag()==0.)
425     {                                             463     {
426       // if |photonPolarization|==0. or |photo    464       // if |photonPolarization|==0. or |photonPolarization * photonDirection0| > 1e-6 * |photonPolarization ^ photonDirection0|
427       // then polarization is choosen randomly << 465       // then polarization is choosen randomly.
                                                   >> 466   
428       G4ThreeVector e1(photonMomentumDirection    467       G4ThreeVector e1(photonMomentumDirection.orthogonal().unit());
429       G4ThreeVector e2(photonMomentumDirection    468       G4ThreeVector e2(photonMomentumDirection.cross(e1).unit());
430                                                   469   
431       G4double angle(G4UniformRand() * CLHEP:: << 470       G4double angle(G4UniformRand() * twopi);
432                                                   471   
433       e1*=std::cos(angle);                        472       e1*=std::cos(angle);
434       e2*=std::sin(angle);                        473       e2*=std::sin(angle);
435                                                   474   
436       photonPolarization=e1+e2;                   475       photonPolarization=e1+e2;
437     }                                             476     }
438   else if (photonPolarization.howOrthogonal(ph    477   else if (photonPolarization.howOrthogonal(photonMomentumDirection) != 0.)
439     {                                             478     {
440       // if |photonPolarization * photonDirect    479       // if |photonPolarization * photonDirection0| != 0.
441       // then polarization is made orthonormal << 480       // then polarization is made orthonormal;
                                                   >> 481   
442       photonPolarization=photonPolarization.pe    482       photonPolarization=photonPolarization.perpPart(photonMomentumDirection);
443     }                                             483     }
444                                                   484  
445   return photonPolarization.unit();               485   return photonPolarization.unit();
446 }                                                 486 }
447                                                   487 
448 //....oooOO0OOooo........oooOO0OOooo........oo    488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
449                                                   489  
450 void G4LivermorePolarizedRayleighModel::Initia << 490  #include "G4AutoLock.hh"
451                 const G4ParticleDefinition*, G << 491  namespace { G4Mutex LivermorePolarizedRayleighModelMutex = G4MUTEX_INITIALIZER; }
452 {                                              << 492  
453   G4AutoLock l(&LivermorePolarizedRayleighMode << 493 void  G4LivermorePolarizedRayleighModel::InitialiseForElement(const G4ParticleDefinition*, 
454   if(nullptr == dataCS[Z]) { ReadData(Z); }    << 494                   G4int Z)
455   l.unlock();                                  << 495  {
456 }                                              << 496    G4AutoLock l(&LivermorePolarizedRayleighModelMutex);
                                                   >> 497    //  G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= " 
                                                   >> 498    //   << Z << G4endl;
                                                   >> 499    if(!dataCS[Z]) { ReadData(Z); }
                                                   >> 500    l.unlock();
                                                   >> 501  }
457                                                   502