Geant4 Cross Reference

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


  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: G4LivermorePolarizedGammaConversionModel.hh,v 1.2 2010-11-23 16:42:15 flongo Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
 26 //                                                 28 //
 27 // Authors: G.Depaola & F.Longo                    29 // Authors: G.Depaola & F.Longo
 28 //                                                 30 //
 29 //                                                 31 //
 30                                                    32 
 31 #include "G4LivermorePolarizedGammaConversionM     33 #include "G4LivermorePolarizedGammaConversionModel.hh"
 32 #include "G4PhysicalConstants.hh"                  34 #include "G4PhysicalConstants.hh"
 33 #include "G4SystemOfUnits.hh"                      35 #include "G4SystemOfUnits.hh"
 34 #include "G4Electron.hh"                       << 
 35 #include "G4Positron.hh"                       << 
 36 #include "G4ParticleChangeForGamma.hh"         << 
 37 #include "G4Log.hh"                            << 
 38 #include "G4AutoLock.hh"                       << 
 39 #include "G4Exp.hh"                            << 
 40 #include "G4ProductionCutsTable.hh"            << 
 41                                                    36 
 42 //....oooOO0OOooo........oooOO0OOooo........oo     37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 43                                                    38 
 44 using namespace std;                               39 using namespace std;
 45 namespace { G4Mutex LivermorePolarizedGammaCon << 
 46                                                << 
 47 G4PhysicsFreeVector* G4LivermorePolarizedGamma << 
 48                                                    40 
 49 //....oooOO0OOooo........oooOO0OOooo........oo     41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 50                                                    42 
 51 G4LivermorePolarizedGammaConversionModel::G4Li     43 G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel(
 52    const G4ParticleDefinition*, const G4String     44    const G4ParticleDefinition*, const G4String& nam)
 53   :G4VEmModel(nam), smallEnergy(2.*MeV), isIni <<  45   :G4VEmModel(nam),fParticleChange(0),
                                                   >>  46    isInitialised(false),meanFreePathTable(0),crossSectionHandler(0)
 54 {                                                  47 {
 55   fParticleChange = nullptr;                   << 
 56   lowEnergyLimit = 2*electron_mass_c2;             48   lowEnergyLimit = 2*electron_mass_c2;
 57                                                <<  49   highEnergyLimit = 100 * GeV;
                                                   >>  50   SetLowEnergyLimit(lowEnergyLimit);
                                                   >>  51   SetHighEnergyLimit(highEnergyLimit);
                                                   >>  52   smallEnergy = 2.*MeV;
                                                   >>  53 
 58   Phi=0.;                                          54   Phi=0.;
 59   Psi=0.;                                          55   Psi=0.;
 60                                                    56   
 61   verboseLevel= 0;                                 57   verboseLevel= 0;
 62   // Verbosity scale:                              58   // Verbosity scale:
 63   // 0 = nothing                                   59   // 0 = nothing 
 64   // 1 = calculation of cross sections, file o <<  60   // 1 = warning for energy non-conservation 
 65   // 2 = entering in methods                   <<  61   // 2 = details of energy budget
 66                                                <<  62   // 3 = calculation of cross sections, file openings, samping of atoms
                                                   >>  63   // 4 = entering in methods
                                                   >>  64 
 67   if(verboseLevel > 0) {                           65   if(verboseLevel > 0) {
 68     G4cout << "Livermore Polarized GammaConver <<  66     G4cout << "Livermore Polarized GammaConversion is constructed " << G4endl
                                                   >>  67      << "Energy range: "
                                                   >>  68      << lowEnergyLimit / keV << " keV - "
                                                   >>  69      << highEnergyLimit / GeV << " GeV"
 69      << G4endl;                                    70      << G4endl;
 70   }                                                71   }
 71                                                <<  72 
                                                   >>  73   crossSectionHandler = new G4CrossSectionHandler();
                                                   >>  74   crossSectionHandler->Initialise(0,lowEnergyLimit,highEnergyLimit,400);
 72 }                                                  75 }
 73                                                    76 
 74 //....oooOO0OOooo........oooOO0OOooo........oo     77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 75                                                    78 
 76 G4LivermorePolarizedGammaConversionModel::~G4L     79 G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel()
 77 {                                                  80 {  
 78   if(IsMaster()) {                             <<  81   delete crossSectionHandler;
 79     for(G4int i=0; i<maxZ; ++i) {              << 
 80       if(data[i]) {                            << 
 81   delete data[i];                              << 
 82   data[i] = nullptr;                           << 
 83       }                                        << 
 84     }                                          << 
 85   }                                            << 
 86 }                                                  82 }
 87                                                    83 
                                                   >>  84 
                                                   >>  85 
 88 //....oooOO0OOooo........oooOO0OOooo........oo     86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 89                                                    87 
 90 void G4LivermorePolarizedGammaConversionModel:     88 void G4LivermorePolarizedGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
 91                                        const G     89                                        const G4DataVector& cuts)
 92 {                                                  90 {
 93   if (verboseLevel > 1)                        <<  91   if (verboseLevel > 3)
 94     {                                          <<  92     G4cout << "Calling G4LivermorePolarizedGammaConversionModel::Initialise()" << G4endl;
 95       G4cout << "Calling1 G4LivermorePolarized << 
 96        << G4endl                               << 
 97              << "Energy range: "               << 
 98        << LowEnergyLimit() / MeV << " MeV - "  << 
 99              << HighEnergyLimit() / GeV << " G << 
100              << G4endl;                        << 
101     }                                          << 
102                                                << 
103   if(IsMaster())                               << 
104     {                                          << 
105       // Initialise element selector           << 
106       InitialiseElementSelectors(particle, cut << 
107                                                << 
108       // Access to elements                    << 
109       const char* path = G4FindDataDir("G4LEDA << 
110                                                << 
111       G4ProductionCutsTable* theCoupleTable =  << 
112   G4ProductionCutsTable::GetProductionCutsTabl << 
113                                                << 
114       G4int numOfCouples = (G4int)theCoupleTab << 
115                                                << 
116       for(G4int i=0; i<numOfCouples; ++i)      << 
117   {                                            << 
118     const G4Material* material =               << 
119       theCoupleTable->GetMaterialCutsCouple(i) << 
120     const G4ElementVector* theElementVector =  << 
121     std::size_t nelm = material->GetNumberOfEl << 
122                                                << 
123     for (std::size_t j=0; j<nelm; ++j)         << 
124       {                                        << 
125         G4int Z = (G4int)(*theElementVector)[j << 
126         if(Z < 1)          { Z = 1; }          << 
127         else if(Z > maxZ)  { Z = maxZ; }       << 
128         if(!data[Z]) { ReadData(Z, path); }    << 
129       }                                        << 
130   }                                            << 
131     }                                          << 
132   if(isInitialised) { return; }                << 
133   fParticleChange = GetParticleChangeForGamma( << 
134   isInitialised = true;                        << 
135 }                                              << 
136                                                << 
137 //....oooOO0OOooo........oooOO0OOooo........oo << 
138                                                    93 
139 void G4LivermorePolarizedGammaConversionModel: <<  94   if (crossSectionHandler)
140      const G4ParticleDefinition*, G4VEmModel*  <<  95   {
141 {                                              <<  96     crossSectionHandler->Clear();
142   SetElementSelectors(masterModel->GetElementS <<  97     delete crossSectionHandler;
143 }                                              <<  98   }
144                                                    99 
145 //....oooOO0OOooo........oooOO0OOooo........oo << 100   // Energy limits
                                                   >> 101   /*
                                                   >> 102   // V.Ivanchenko: this was meanless check  
                                                   >> 103   if (LowEnergyLimit() < lowEnergyLimit)
                                                   >> 104     {
                                                   >> 105       G4cout << "G4LivermorePolarizedGammaConversionModel: low energy limit increased from " << 
                                                   >> 106   LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
                                                   >> 107       //      SetLowEnergyLimit(lowEnergyLimit);
                                                   >> 108     }
                                                   >> 109   */
                                                   >> 110   if (HighEnergyLimit() > highEnergyLimit)
                                                   >> 111     {
                                                   >> 112       G4cout << "G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " << 
                                                   >> 113   HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
                                                   >> 114       // V.Ivanchenko: this is forbidden 
                                                   >> 115      // SetHighEnergyLimit(highEnergyLimit);
                                                   >> 116     }
                                                   >> 117 
                                                   >> 118   // Reading of data files - all materials are read
                                                   >> 119   
                                                   >> 120   crossSectionHandler = new G4CrossSectionHandler;
                                                   >> 121   crossSectionHandler->Clear();
                                                   >> 122   G4String crossSectionFile = "pair/pp-cs-";
                                                   >> 123   crossSectionHandler->LoadData(crossSectionFile);
146                                                   124 
147 G4double G4LivermorePolarizedGammaConversionMo << 125   //
148            const G4ParticleDefinition*, G4doub << 126   if (verboseLevel > 2) {
149 {                                              << 127     G4cout << "Loaded cross section files for Livermore Polarized GammaConversion model" 
150   return lowEnergyLimit;                       << 128      << G4endl;
151 }                                              << 129   }
                                                   >> 130   InitialiseElementSelectors(particle,cuts);
152                                                   131 
153 //....oooOO0OOooo........oooOO0OOooo........oo << 132   if(verboseLevel > 0) {
                                                   >> 133     G4cout << "Livermore Polarized GammaConversion model is initialized " << G4endl
                                                   >> 134      << "Energy range: "
                                                   >> 135      << LowEnergyLimit() / keV << " keV - "
                                                   >> 136      << HighEnergyLimit() / GeV << " GeV"
                                                   >> 137      << G4endl;
                                                   >> 138   }
154                                                   139 
155 void G4LivermorePolarizedGammaConversionModel: << 140   //    
156 {                                              << 141   if(!isInitialised) { 
157   if (verboseLevel > 1)                        << 142     isInitialised = true;
158     {                                          << 143     fParticleChange = GetParticleChangeForGamma();
159       G4cout << "Calling ReadData() of G4Liver << 144   }
160        << G4endl;                              << 
161     }                                          << 
162                                                << 
163   if(data[Z]) { return; }                      << 
164                                                << 
165   const char* datadir = path;                  << 
166                                                << 
167   if(!datadir)                                 << 
168     {                                          << 
169       datadir = G4FindDataDir("G4LEDATA");     << 
170       if(!datadir)                             << 
171   {                                            << 
172     G4Exception("G4LivermorePolarizedGammaConv << 
173           "em0006",FatalException,             << 
174           "Environment variable G4LEDATA not d << 
175     return;                                    << 
176   }                                            << 
177     }                                          << 
178   //                                           << 
179   data[Z] = new G4PhysicsFreeVector(0,/*spline << 
180   //                                           << 
181   std::ostringstream ost;                      << 
182   ost << datadir << "/livermore/pair/pp-cs-" < << 
183   std::ifstream fin(ost.str().c_str());        << 
184                                                << 
185   if( !fin.is_open())                          << 
186     {                                          << 
187       G4ExceptionDescription ed;               << 
188       ed << "G4LivermorePolarizedGammaConversi << 
189    << "> is not opened!" << G4endl;            << 
190       G4Exception("G4LivermorePolarizedGammaCo << 
191       "em0003",FatalException,                 << 
192       ed,"G4LEDATA version should be G4EMLOW6. << 
193       return;                                  << 
194     }                                          << 
195   else                                         << 
196     {                                          << 
197                                                << 
198       if(verboseLevel > 3) { G4cout << "File " << 
199             << " is opened by G4LivermorePolar << 
200                                                << 
201       data[Z]->Retrieve(fin, true);            << 
202     }                                          << 
203                                                << 
204   // Activation of spline interpolation        << 
205   data[Z]->FillSecondDerivatives();            << 
206                                                << 
207 }                                                 145 }
208                                                   146 
209 //....oooOO0OOooo........oooOO0OOooo........oo    147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210                                                   148 
211 G4double G4LivermorePolarizedGammaConversionMo    149 G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom(
212                                        const G    150                                        const G4ParticleDefinition*,
213                G4double GammaEnergy,              151                G4double GammaEnergy,
214                G4double Z, G4double,              152                G4double Z, G4double,
215                G4double, G4double)                153                G4double, G4double)
216 {                                                 154 {
217   if (verboseLevel > 1) {                      << 155   if (verboseLevel > 3) {
218     G4cout << "G4LivermorePolarizedGammaConver    156     G4cout << "G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom()" 
219      << G4endl;                                   157      << G4endl;
220   }                                               158   }
221   if (GammaEnergy < lowEnergyLimit) { return 0 << 159   if(Z < 0.9 || GammaEnergy <= lowEnergyLimit) { return 0.0; }
222                                                << 160   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
223   G4double xs = 0.0;                           << 161   return cs;
224                                                << 
225   G4int intZ=G4int(Z);                         << 
226                                                << 
227   if(intZ < 1 || intZ > maxZ) { return xs; }   << 
228                                                << 
229   G4PhysicsFreeVector* pv = data[intZ];        << 
230                                                << 
231   // if element was not initialised            << 
232   // do initialisation safely for MT mode      << 
233   if(!pv)                                      << 
234     {                                          << 
235       InitialiseForElement(0, intZ);           << 
236       pv = data[intZ];                         << 
237       if(!pv) { return xs; }                   << 
238     }                                          << 
239   // x-section is taken from the table         << 
240   xs = pv->Value(GammaEnergy);                 << 
241                                                << 
242   if(verboseLevel > 0)                         << 
243     {                                          << 
244       G4int n = G4int(pv->GetVectorLength() -  << 
245       G4cout  <<  "****** DEBUG: tcs value for << 
246         << GammaEnergy/MeV << G4endl;          << 
247       G4cout  <<  "  cs (Geant4 internal unit) << 
248       G4cout  <<  "    -> first cs value in EA << 
249       G4cout  <<  "    -> last  cs value in EA << 
250       G4cout  <<  "*************************** << 
251     }                                          << 
252                                                << 
253   return xs;                                   << 
254 }                                                 162 }
255                                                   163 
256 //....oooOO0OOooo........oooOO0OOooo........oo    164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
257                                                   165 
258 void                                              166 void 
259 G4LivermorePolarizedGammaConversionModel::Samp    167 G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
260                   const G4MaterialCutsCouple*     168                   const G4MaterialCutsCouple* couple,
261                   const G4DynamicParticle* aDy    169                   const G4DynamicParticle* aDynamicGamma,
262                   G4double,                       170                   G4double,
263                   G4double)                       171                   G4double)
264 {                                                 172 {
265                                                   173 
266   // Fluorescence generated according to:         174   // Fluorescence generated according to:
267   // J. Stepanek ,"A program to determine the     175   // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
268   // subshell ionisation by a particle or due     176   // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
269   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)       177   // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
                                                   >> 178 
270   if (verboseLevel > 3)                           179   if (verboseLevel > 3)
271     G4cout << "Calling SampleSecondaries() of     180     G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
272                                                   181 
273   G4double photonEnergy = aDynamicGamma->GetKi    182   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
                                                   >> 183   // Within energy limit?
274                                                   184 
275   if(photonEnergy <= lowEnergyLimit)              185   if(photonEnergy <= lowEnergyLimit)
276     {                                             186     {
277       fParticleChange->ProposeTrackStatus(fSto    187       fParticleChange->ProposeTrackStatus(fStopAndKill);
278       fParticleChange->SetProposedKineticEnerg    188       fParticleChange->SetProposedKineticEnergy(0.);
279       return;                                     189       return;
280     }                                             190     }
281                                                   191 
                                                   >> 192 
282   G4ThreeVector gammaPolarization0 = aDynamicG    193   G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
283   G4ThreeVector gammaDirection0 = aDynamicGamm    194   G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
284                                                   195 
285   // Make sure that the polarization vector is    196   // Make sure that the polarization vector is perpendicular to the
286   // gamma direction. If not                      197   // gamma direction. If not
                                                   >> 198 
287   if(!(gammaPolarization0.isOrthogonal(gammaDi    199   if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
288     { // only for testing now                     200     { // only for testing now
289       gammaPolarization0 = GetRandomPolarizati    201       gammaPolarization0 = GetRandomPolarization(gammaDirection0);
290     }                                             202     }
291   else                                            203   else
292     {                                             204     {
293       if ( gammaPolarization0.howOrthogonal(ga    205       if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
294   {                                               206   {
295     gammaPolarization0 = GetPerpendicularPolar    207     gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
296   }                                               208   }
297     }                                             209     }
298                                                   210 
299   // End of Protection                            211   // End of Protection
300                                                   212 
                                                   >> 213 
301   G4double epsilon ;                              214   G4double epsilon ;
302   G4double epsilon0Local = electron_mass_c2 /     215   G4double epsilon0Local = electron_mass_c2 / photonEnergy ;
303                                                   216 
304   // Do it fast if photon energy < 2. MeV         217   // Do it fast if photon energy < 2. MeV
305                                                   218 
306   if (photonEnergy < smallEnergy )                219   if (photonEnergy < smallEnergy )
307     {                                             220     {
308       epsilon = epsilon0Local + (0.5 - epsilon    221       epsilon = epsilon0Local + (0.5 - epsilon0Local) * G4UniformRand();
309     }                                             222     }
310   else                                            223   else
311     {                                             224     {
312       // Select randomly one element in the cu << 225 
                                                   >> 226  // Select randomly one element in the current material
                                                   >> 227 
                                                   >> 228       //     G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
                                                   >> 229       //const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
                                                   >> 230 
313       const G4ParticleDefinition* particle =      231       const G4ParticleDefinition* particle =  aDynamicGamma->GetDefinition();
314       const G4Element* element = SelectRandomA    232       const G4Element* element = SelectRandomAtom(couple,particle,photonEnergy);
315                                                << 233 
316       if (element == nullptr)                  << 234       /*
                                                   >> 235       if (element == 0)
317         {                                         236         {
318           G4cout << "G4LivermorePolarizedGamma << 237           G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - element = 0" << G4endl;
319     return;                                    << 
320         }                                         238         }
                                                   >> 239       */
321                                                   240       
                                                   >> 241       G4IonisParamElm* ionisation = element->GetIonisation();
322                                                   242       
323       G4IonisParamElm* ionisation = element->G << 243       /*
324       if (ionisation == nullptr)               << 244       if (ionisation == 0)
325         {                                         245         {
326           G4cout << "G4LivermorePolarizedGamma << 246           G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - ionisation = 0" << G4endl;
327     return;                                    << 
328         }                                         247         }
329                                                << 248       */
                                                   >> 249 
                                                   >> 250 
330       // Extract Coulomb factor for this Eleme    251       // Extract Coulomb factor for this Element
                                                   >> 252 
331       G4double fZ = 8. * (ionisation->GetlogZ3    253       G4double fZ = 8. * (ionisation->GetlogZ3());
332       if (photonEnergy > 50. * MeV) fZ += 8. *    254       if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
333                                                   255 
334       // Limits of the screening variable         256       // Limits of the screening variable
335       G4double screenFactor = 136. * epsilon0L    257       G4double screenFactor = 136. * epsilon0Local / (element->GetIonisation()->GetZ3()) ;
336       G4double screenMax = G4Exp ((42.24 - fZ) << 258       G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
337       G4double screenMin = std::min(4.*screenF    259       G4double screenMin = std::min(4.*screenFactor,screenMax) ;
338                                                   260 
339       // Limits of the energy sampling            261       // Limits of the energy sampling
340       G4double epsilon1 = 0.5 - 0.5 * sqrt(1.     262       G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
341       G4double epsilonMin = std::max(epsilon0L    263       G4double epsilonMin = std::max(epsilon0Local,epsilon1);
342       G4double epsilonRange = 0.5 - epsilonMin    264       G4double epsilonRange = 0.5 - epsilonMin ;
343                                                   265 
344       // Sample the energy rate of the created    266       // Sample the energy rate of the created electron (or positron)
345       G4double screen;                            267       G4double screen;
346       G4double gReject ;                          268       G4double gReject ;
347                                                   269 
348       G4double f10 = ScreenFunction1(screenMin    270       G4double f10 = ScreenFunction1(screenMin) - fZ;
349       G4double f20 = ScreenFunction2(screenMin    271       G4double f20 = ScreenFunction2(screenMin) - fZ;
350       G4double normF1 = std::max(f10 * epsilon    272       G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
351       G4double normF2 = std::max(1.5 * f20,0.)    273       G4double normF2 = std::max(1.5 * f20,0.);
352                                                   274 
353       do {                                        275       do {
354         if (normF1 / (normF1 + normF2) > G4Uni    276         if (normF1 / (normF1 + normF2) > G4UniformRand() )
355           {                                       277           {
356             epsilon = 0.5 - epsilonRange * pow    278             epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ;
357             screen = screenFactor / (epsilon *    279             screen = screenFactor / (epsilon * (1. - epsilon));
358             gReject = (ScreenFunction1(screen)    280             gReject = (ScreenFunction1(screen) - fZ) / f10 ;
359           }                                       281           }
360         else                                      282         else
361           {                                       283           {
362             epsilon = epsilonMin + epsilonRang    284             epsilon = epsilonMin + epsilonRange * G4UniformRand();
363             screen = screenFactor / (epsilon *    285             screen = screenFactor / (epsilon * (1 - epsilon));
364             gReject = (ScreenFunction2(screen)    286             gReject = (ScreenFunction2(screen) - fZ) / f20 ;
                                                   >> 287 
                                                   >> 288 
365     }                                             289     }
366       } while ( gReject < G4UniformRand() );      290       } while ( gReject < G4UniformRand() );
                                                   >> 291 
367     }   //  End of epsilon sampling               292     }   //  End of epsilon sampling
368                                                   293   
369   // Fix charges randomly                         294   // Fix charges randomly
                                                   >> 295   
370   G4double electronTotEnergy;                     296   G4double electronTotEnergy;
371   G4double positronTotEnergy;                     297   G4double positronTotEnergy;
372                                                   298 
373   if (G4UniformRand() > 0.5)                   << 299 
                                                   >> 300   if (G4int(2*G4UniformRand()))  
374     {                                             301     {
375       electronTotEnergy = (1. - epsilon) * pho    302       electronTotEnergy = (1. - epsilon) * photonEnergy;
376       positronTotEnergy = epsilon * photonEner    303       positronTotEnergy = epsilon * photonEnergy;
377     }                                             304     }
378   else                                            305   else
379     {                                             306     {
380       positronTotEnergy = (1. - epsilon) * pho    307       positronTotEnergy = (1. - epsilon) * photonEnergy;
381       electronTotEnergy = epsilon * photonEner    308       electronTotEnergy = epsilon * photonEnergy;
382     }                                             309     }
383                                                << 310 
384   // Scattered electron (positron) angles. ( Z    311   // Scattered electron (positron) angles. ( Z - axis along the parent photon)
385   // Universal distribution suggested by L. Ur    312   // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
386   // derived from Tsai distribution (Rev. Mod.    313   // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
                                                   >> 314 
                                                   >> 315 /*
                                                   >> 316   G4double u;
                                                   >> 317   const G4double a1 = 0.625;
                                                   >> 318   G4double a2 = 3. * a1;
                                                   >> 319 
                                                   >> 320   if (0.25 > G4UniformRand())
                                                   >> 321     {
                                                   >> 322       u = - log(G4UniformRand() * G4UniformRand()) / a1 ;
                                                   >> 323     }
                                                   >> 324   else
                                                   >> 325     {
                                                   >> 326       u = - log(G4UniformRand() * G4UniformRand()) / a2 ;
                                                   >> 327     }
                                                   >> 328 */
                                                   >> 329 
387   G4double Ene = electronTotEnergy/electron_ma    330   G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy
388                                                   331 
389   G4double cosTheta = 0.;                         332   G4double cosTheta = 0.;
390   G4double sinTheta = 0.;                         333   G4double sinTheta = 0.;
391                                                   334 
392   SetTheta(&cosTheta,&sinTheta,Ene);              335   SetTheta(&cosTheta,&sinTheta,Ene);
                                                   >> 336 
                                                   >> 337   //  G4double theta = u * electron_mass_c2 / photonEnergy ;
                                                   >> 338   //  G4double phi  = twopi * G4UniformRand() ;
                                                   >> 339 
393   G4double phi,psi=0.;                            340   G4double phi,psi=0.;
394                                                   341 
395   //corrected e+ e- angular angular distributi    342   //corrected e+ e- angular angular distribution //preliminary!
                                                   >> 343 
                                                   >> 344   //  if(photonEnergy>50*MeV)
                                                   >> 345   // {
396   phi = SetPhi(photonEnergy);                     346   phi = SetPhi(photonEnergy);
397   psi = SetPsi(photonEnergy,phi);                 347   psi = SetPsi(photonEnergy,phi);
                                                   >> 348   //  }
                                                   >> 349   //else
                                                   >> 350   // {
                                                   >> 351   //psi = G4UniformRand()*2.*pi;
                                                   >> 352   //phi = pi; // coplanar
                                                   >> 353   // }
                                                   >> 354 
398   Psi = psi;                                      355   Psi = psi;
399   Phi = phi;                                      356   Phi = phi;
                                                   >> 357   //G4cout << "PHI " << phi << G4endl;
                                                   >> 358   //G4cout << "PSI " << psi << G4endl;
400                                                   359 
401   G4double phie, phip;                         << 360   G4double phie = psi; //azimuthal angle for the electron
402   G4double choice, choice2;                    << 
403   choice = G4UniformRand();                    << 
404   choice2 = G4UniformRand();                   << 
405                                                   361 
406   if (choice2 <= 0.5)                          << 
407     {                                          << 
408       // do nothing                            << 
409       //  phi = phi;                           << 
410     }                                          << 
411   else                                         << 
412     {                                          << 
413       phi = -phi;                              << 
414     }                                          << 
415                                                << 
416   if (choice <= 0.5)                           << 
417     {                                          << 
418       phie = psi; //azimuthal angle for the el << 
419       phip = phie+phi; //azimuthal angle for t << 
420     }                                          << 
421   else                                         << 
422     {                                          << 
423       // opzione 1 phie / phip equivalenti     << 
424       phip = psi; //azimuthal angle for the po << 
425       phie = phip + phi; //azimuthal angle for << 
426     }                                          << 
427                                                << 
428                                                << 
429   // Electron Kinematics                       << 
430   G4double dirX = sinTheta*cos(phie);             362   G4double dirX = sinTheta*cos(phie);
431   G4double dirY = sinTheta*sin(phie);             363   G4double dirY = sinTheta*sin(phie);
432   G4double dirZ = cosTheta;                       364   G4double dirZ = cosTheta;
433   G4ThreeVector electronDirection(dirX,dirY,di    365   G4ThreeVector electronDirection(dirX,dirY,dirZ);
434                                                << 
435   // Kinematics of the created pair:              366   // Kinematics of the created pair:
436   // the electron and positron are assumed to     367   // the electron and positron are assumed to have a symetric angular
437   // distribution with respect to the Z axis a    368   // distribution with respect to the Z axis along the parent photon
438                                                   369 
                                                   >> 370   //G4double localEnergyDeposit = 0. ;
                                                   >> 371 
439   G4double electronKineEnergy = std::max(0.,el    372   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
440                                                   373 
441   SystemOfRefChange(gammaDirection0,electronDi    374   SystemOfRefChange(gammaDirection0,electronDirection,
442         gammaPolarization0);                      375         gammaPolarization0);
443                                                   376 
444   G4DynamicParticle* particle1 = new G4Dynamic    377   G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
445               electronDirection,                  378               electronDirection,
446               electronKineEnergy);                379               electronKineEnergy);
447                                                   380 
448   // The e+ is always created (even with kinet    381   // The e+ is always created (even with kinetic energy = 0) for further annihilation
                                                   >> 382 
449   Ene = positronTotEnergy/electron_mass_c2; //    383   Ene = positronTotEnergy/electron_mass_c2; // Normalized energy
450                                                   384 
451   cosTheta = 0.;                                  385   cosTheta = 0.;
452   sinTheta = 0.;                                  386   sinTheta = 0.;
453                                                   387 
454   SetTheta(&cosTheta,&sinTheta,Ene);              388   SetTheta(&cosTheta,&sinTheta,Ene);
                                                   >> 389   G4double phip = phie+phi; //azimuthal angle for the positron
455                                                   390 
456   // Positron Kinematics                       << 
457   dirX = sinTheta*cos(phip);                      391   dirX = sinTheta*cos(phip);
458   dirY = sinTheta*sin(phip);                      392   dirY = sinTheta*sin(phip);
459   dirZ = cosTheta;                                393   dirZ = cosTheta;
460   G4ThreeVector positronDirection(dirX,dirY,di    394   G4ThreeVector positronDirection(dirX,dirY,dirZ);
461                                                   395 
462   G4double positronKineEnergy = std::max(0.,po    396   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
463   SystemOfRefChange(gammaDirection0,positronDi    397   SystemOfRefChange(gammaDirection0,positronDirection,
464         gammaPolarization0);                      398         gammaPolarization0);
465                                                   399 
466   // Create G4DynamicParticle object for the p    400   // Create G4DynamicParticle object for the particle2
467   G4DynamicParticle* particle2 = new G4Dynamic    401   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
468                                                   402                                                        positronDirection, positronKineEnergy);
                                                   >> 403 
                                                   >> 404 
469   fvect->push_back(particle1);                    405   fvect->push_back(particle1);
470   fvect->push_back(particle2);                    406   fvect->push_back(particle2);
471                                                   407 
                                                   >> 408 
                                                   >> 409 
472   // Kill the incident photon                     410   // Kill the incident photon
                                                   >> 411 
                                                   >> 412 
                                                   >> 413 
                                                   >> 414   // Create lists of pointers to DynamicParticles (photons and electrons)
                                                   >> 415   // (Is the electron vector necessary? To be checked)
                                                   >> 416   //  std::vector<G4DynamicParticle*>* photonVector = 0;
                                                   >> 417   //std::vector<G4DynamicParticle*> electronVector;
                                                   >> 418 
                                                   >> 419   fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
473   fParticleChange->SetProposedKineticEnergy(0.    420   fParticleChange->SetProposedKineticEnergy(0.);
474   fParticleChange->ProposeTrackStatus(fStopAnd    421   fParticleChange->ProposeTrackStatus(fStopAndKill);
                                                   >> 422 
475 }                                                 423 }
476                                                   424 
477 //....oooOO0OOooo........oooOO0OOooo........oo    425 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
478                                                   426 
479 G4double G4LivermorePolarizedGammaConversionMo    427 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1(G4double screenVariable)
480 {                                                 428 {
481   // Compute the value of the screening functi    429   // Compute the value of the screening function 3*phi1 - phi2
                                                   >> 430 
482   G4double value;                                 431   G4double value;
                                                   >> 432 
483   if (screenVariable > 1.)                        433   if (screenVariable > 1.)
484     value = 42.24 - 8.368 * log(screenVariable    434     value = 42.24 - 8.368 * log(screenVariable + 0.952);
485   else                                            435   else
486     value = 42.392 - screenVariable * (7.796 -    436     value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
487                                                   437 
488   return value;                                   438   return value;
489 }                                                 439 }
490                                                   440 
491                                                   441 
492                                                   442 
493 G4double G4LivermorePolarizedGammaConversionMo    443 G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2(G4double screenVariable)
494 {                                                 444 {
495   // Compute the value of the screening functi    445   // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
                                                   >> 446 
496   G4double value;                                 447   G4double value;
497                                                   448 
498   if (screenVariable > 1.)                        449   if (screenVariable > 1.)
499     value = 42.24 - 8.368 * log(screenVariable    450     value = 42.24 - 8.368 * log(screenVariable + 0.952);
500   else                                            451   else
501     value = 41.405 - screenVariable * (5.828 -    452     value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
502                                                   453 
503   return value;                                   454   return value;
504 }                                                 455 }
505                                                   456 
506                                                   457 
507 void G4LivermorePolarizedGammaConversionModel:    458 void G4LivermorePolarizedGammaConversionModel::SetTheta(G4double* p_cosTheta, G4double* p_sinTheta, G4double Energy)
508 {                                                 459 {
                                                   >> 460 
509   // to avoid computational errors since Theta    461   // to avoid computational errors since Theta could be very small
510   // Energy in Normalized Units (!)               462   // Energy in Normalized Units (!)
511                                                   463 
512   G4double Momentum = sqrt(Energy*Energy -1);     464   G4double Momentum = sqrt(Energy*Energy -1);
513   G4double Rand = G4UniformRand();                465   G4double Rand = G4UniformRand();
514                                                   466 
515   *p_cosTheta = (Energy*((2*Rand)- 1) + Moment    467   *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
516   *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momen    468   *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
517 }                                                 469 }
518                                                   470 
519                                                   471 
520                                                   472 
521 G4double G4LivermorePolarizedGammaConversionMo    473 G4double G4LivermorePolarizedGammaConversionModel::SetPhi(G4double Energy)
522 {                                                 474 {
                                                   >> 475 
                                                   >> 476 
523   G4double value = 0.;                            477   G4double value = 0.;
524   G4double Ene = Energy/MeV;                      478   G4double Ene = Energy/MeV;
525                                                   479 
526   G4double pl[4];                                 480   G4double pl[4];
                                                   >> 481 
                                                   >> 482 
527   G4double pt[2];                                 483   G4double pt[2];
528   G4double xi = 0;                                484   G4double xi = 0;
529   G4double xe = 0.;                               485   G4double xe = 0.;
530   G4double n1=0.;                                 486   G4double n1=0.;
531   G4double n2=0.;                                 487   G4double n2=0.;
532                                                   488 
                                                   >> 489 
533   if (Ene>=50.)                                   490   if (Ene>=50.)
534     {                                             491     {
535       const G4double ay0=5.6, by0=18.6, aa0=2.    492       const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
536       const G4double aw = 0.0151, bw = 10.7, c    493       const G4double aw = 0.0151, bw = 10.7, cw = -410.;
537                                                   494 
538       const G4double axc = 3.1455, bxc = -1.11    495       const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
539                                                   496 
540       pl[0] = Fln(ay0,by0,Ene);                   497       pl[0] = Fln(ay0,by0,Ene);
541       pl[1] = aa0 + ba0*(Ene);                    498       pl[1] = aa0 + ba0*(Ene);
542       pl[2] = Poli(aw,bw,cw,Ene);                 499       pl[2] = Poli(aw,bw,cw,Ene);
543       pl[3] = Poli(axc,bxc,cxc,Ene);              500       pl[3] = Poli(axc,bxc,cxc,Ene);
544                                                   501 
545       const G4double abf = 3.1216, bbf = 2.68;    502       const G4double abf = 3.1216, bbf = 2.68;
546       pt[0] = -1.4;                               503       pt[0] = -1.4;
547       pt[1] = abf + bbf/Ene;                      504       pt[1] = abf + bbf/Ene;
548                                                   505 
                                                   >> 506 
                                                   >> 507 
                                                   >> 508       //G4cout << "PL > 50. "<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
                                                   >> 509 
549       xi = 3.0;                                   510       xi = 3.0;
550       xe = Encu(pl,pt,xi);                        511       xe = Encu(pl,pt,xi);
                                                   >> 512       //G4cout << "ENCU "<< xe << G4endl;
551       n1 = Fintlor(pl,pi) - Fintlor(pl,xe);       513       n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
552       n2 = Finttan(pt,xe) - Finttan(pt,0.);       514       n2 = Finttan(pt,xe) - Finttan(pt,0.);
553     }                                             515     }
554   else                                            516   else
555     {                                             517     {
556       const G4double ay0=0.144, by0=0.11;         518       const G4double ay0=0.144, by0=0.11;
557       const G4double aa0=2.7, ba0 = 2.74;         519       const G4double aa0=2.7, ba0 = 2.74;
558       const G4double aw = 0.21, bw = 10.8, cw     520       const G4double aw = 0.21, bw = 10.8, cw = -58.;
559       const G4double axc = 3.17, bxc = -0.87,     521       const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
560                                                   522 
561       pl[0] = Fln(ay0, by0, Ene);                 523       pl[0] = Fln(ay0, by0, Ene);
562       pl[1] = Fln(aa0, ba0, Ene);                 524       pl[1] = Fln(aa0, ba0, Ene);
563       pl[2] = Poli(aw,bw,cw,Ene);                 525       pl[2] = Poli(aw,bw,cw,Ene);
564       pl[3] = Poli(axc,bxc,cxc,Ene);              526       pl[3] = Poli(axc,bxc,cxc,Ene);
565                                                   527 
                                                   >> 528       //G4cout << "PL < 50."<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
                                                   >> 529       //G4cout << "ENCU "<< xe << G4endl;
566       n1 = Fintlor(pl,pi) - Fintlor(pl,xe);       530       n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
                                                   >> 531 
567     }                                             532     }
568                                                   533 
569                                                   534 
570   G4double n=0.;                                  535   G4double n=0.;
571   n = n1+n2;                                      536   n = n1+n2;
572                                                   537 
573   G4double c1 = 0.;                               538   G4double c1 = 0.;
574   c1 = Glor(pl, xe);                              539   c1 = Glor(pl, xe);
575                                                   540 
                                                   >> 541 /*
                                                   >> 542   G4double xm = 0.;
                                                   >> 543   xm = Flor(pl,pl[3])*Glor(pl,pl[3]);
                                                   >> 544 */
                                                   >> 545 
576   G4double r1,r2,r3;                              546   G4double r1,r2,r3;
577   G4double xco=0.;                                547   G4double xco=0.;
578                                                   548 
579   if (Ene>=50.)                                   549   if (Ene>=50.)
580     {                                             550     {
581       r1= G4UniformRand();                        551       r1= G4UniformRand();
582       if( r1>=n2/n)                               552       if( r1>=n2/n)
583         {                                         553         {
584           do                                      554           do
585       {                                           555       {
586               r2 = G4UniformRand();               556               r2 = G4UniformRand();
587               value = Finvlor(pl,xe,r2);          557               value = Finvlor(pl,xe,r2);
588               xco = Glor(pl,value)/c1;            558               xco = Glor(pl,value)/c1;
589               r3 = G4UniformRand();               559               r3 = G4UniformRand();
590             } while(r3>=xco);                     560             } while(r3>=xco);
591         }                                         561         }
592       else                                        562       else
593         {                                         563         {
594           value = Finvtan(pt,n,r1);               564           value = Finvtan(pt,n,r1);
595         }                                         565         }
596     }                                             566     }
597   else                                            567   else
598     {                                             568     {
599       do                                          569       do
600         {                                         570         {
601           r2 = G4UniformRand();                   571           r2 = G4UniformRand();
602           value = Finvlor(pl,xe,r2);              572           value = Finvlor(pl,xe,r2);
603           xco = Glor(pl,value)/c1;                573           xco = Glor(pl,value)/c1;
604           r3 = G4UniformRand();                   574           r3 = G4UniformRand();
605         } while(r3>=xco);                         575         } while(r3>=xco);
606     }                                             576     }
                                                   >> 577 
                                                   >> 578   //  G4cout << "PHI = " <<value <<  G4endl;
607   return value;                                   579   return value;
608 }                                                 580 }
609                                                << 
610 //....oooOO0OOooo........oooOO0OOooo........oo << 
611                                                << 
612 G4double G4LivermorePolarizedGammaConversionMo    581 G4double G4LivermorePolarizedGammaConversionModel::SetPsi(G4double Energy, G4double PhiLocal)
613 {                                                 582 {
                                                   >> 583 
614   G4double value = 0.;                            584   G4double value = 0.;
615   G4double Ene = Energy/MeV;                      585   G4double Ene = Energy/MeV;
616                                                   586 
617   G4double p0l[4];                                587   G4double p0l[4];
618   G4double ppml[4];                               588   G4double ppml[4];
619   G4double p0t[2];                                589   G4double p0t[2];
620   G4double ppmt[2];                               590   G4double ppmt[2];
621                                                   591 
622   G4double xi = 0.;                               592   G4double xi = 0.;
623   G4double xe0 = 0.;                              593   G4double xe0 = 0.;
624   G4double xepm = 0.;                             594   G4double xepm = 0.;
625                                                   595 
626   if (Ene>=50.)                                   596   if (Ene>=50.)
627     {                                             597     {
628       const G4double ay00 = 3.4, by00 = 9.8, a    598       const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
629       const G4double aw0 = 0.014, bw0 = 9.7, c    599       const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
630       const G4double axc0 = 3.1423, bxc0 = -2.    600       const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
631       const G4double ay0p = 1.53, by0p = 3.2,     601       const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
632       const G4double awp = 6.9E-3, bwp = 12.6,    602       const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
633       const G4double axcp = 2.8E-3,bxcp = -3.1    603       const G4double axcp = 2.8E-3,bxcp = -3.133;
634       const G4double abf0 = 3.1213, bbf0 = 2.6    604       const G4double abf0 = 3.1213, bbf0 = 2.61;
635       const G4double abfpm = 3.1231, bbfpm = 2    605       const G4double abfpm = 3.1231, bbfpm = 2.84;
636                                                   606 
637       p0l[0] = Fln(ay00, by00, Ene);              607       p0l[0] = Fln(ay00, by00, Ene);
638       p0l[1] = Fln(aa00, ba00, Ene);              608       p0l[1] = Fln(aa00, ba00, Ene);
639       p0l[2] = Poli(aw0, bw0, cw0, Ene);          609       p0l[2] = Poli(aw0, bw0, cw0, Ene);
640       p0l[3] = Poli(axc0, bxc0, cxc0, Ene);       610       p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
641                                                   611 
642       ppml[0] = Fln(ay0p, by0p, Ene);             612       ppml[0] = Fln(ay0p, by0p, Ene);
643       ppml[1] = aap + bap*(Ene);                  613       ppml[1] = aap + bap*(Ene);
644       ppml[2] = Poli(awp, bwp, cwp, Ene);         614       ppml[2] = Poli(awp, bwp, cwp, Ene);
645       ppml[3] = Fln(axcp,bxcp,Ene);               615       ppml[3] = Fln(axcp,bxcp,Ene);
646                                                   616 
647       p0t[0] = -0.81;                             617       p0t[0] = -0.81;
648       p0t[1] = abf0 + bbf0/Ene;                   618       p0t[1] = abf0 + bbf0/Ene;
649       ppmt[0] = -0.6;                             619       ppmt[0] = -0.6;
650       ppmt[1] = abfpm + bbfpm/Ene;                620       ppmt[1] = abfpm + bbfpm/Ene;
651                                                   621 
                                                   >> 622       //G4cout << "P0L > 50"<< p0l[0] << " " << p0l[1] << " " << p0l[2] << " " <<p0l[3] << " " << G4endl;
                                                   >> 623       //G4cout << "PPML > 50"<< ppml[0] << " " << ppml[1] << " " << ppml[2] << " " <<ppml[3] << " " << G4endl;
                                                   >> 624 
652       xi = 3.0;                                   625       xi = 3.0;
653       xe0 = Encu(p0l, p0t, xi);                   626       xe0 = Encu(p0l, p0t, xi);
                                                   >> 627       //G4cout << "ENCU1 "<< xe0 << G4endl;
654       xepm = Encu(ppml, ppmt, xi);                628       xepm = Encu(ppml, ppmt, xi);
                                                   >> 629       //G4cout << "ENCU2 "<< xepm << G4endl;
655     }                                             630     }
656   else                                            631   else
657     {                                             632     {
658       const G4double ay00 = 2.82, by00 = 6.35;    633       const G4double ay00 = 2.82, by00 = 6.35;
659       const G4double aa00 = -1.75, ba00 = 0.25    634       const G4double aa00 = -1.75, ba00 = 0.25;
660                                                   635 
661       const G4double aw0 = 0.028, bw0 = 5., cw    636       const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
662       const G4double axc0 = 3.14213, bxc0 = -2    637       const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
663       const G4double ay0p = 1.56, by0p = 3.6;     638       const G4double ay0p = 1.56, by0p = 3.6;
664       const G4double aap = 0.86, bap = 8.3E-3;    639       const G4double aap = 0.86, bap = 8.3E-3;
665       const G4double awp = 0.022, bwp = 7.4, c    640       const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
666       const G4double xcp = 3.1486;                641       const G4double xcp = 3.1486;
667                                                   642 
668       p0l[0] = Fln(ay00, by00, Ene);              643       p0l[0] = Fln(ay00, by00, Ene);
669       p0l[1] = aa00+pow(Ene, ba00);               644       p0l[1] = aa00+pow(Ene, ba00);
670       p0l[2] = Poli(aw0, bw0, cw0, Ene);          645       p0l[2] = Poli(aw0, bw0, cw0, Ene);
671       p0l[3] = Poli(axc0, bxc0, cxc0, Ene);       646       p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
672       ppml[0] = Fln(ay0p, by0p, Ene);             647       ppml[0] = Fln(ay0p, by0p, Ene);
673       ppml[1] = aap + bap*(Ene);                  648       ppml[1] = aap + bap*(Ene);
674       ppml[2] = Poli(awp, bwp, cwp, Ene);         649       ppml[2] = Poli(awp, bwp, cwp, Ene);
675       ppml[3] = xcp;                              650       ppml[3] = xcp;
                                                   >> 651 
676     }                                             652     }
677                                                   653 
678   G4double a,b=0.;                                654   G4double a,b=0.;
679                                                   655 
680   if (Ene>=50.)                                   656   if (Ene>=50.)
681     {                                             657     {
682       if (PhiLocal>xepm)                          658       if (PhiLocal>xepm)
683   {                                               659   {
684           b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(    660           b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
685         }                                         661         }
686       else                                        662       else
687         {                                         663         {
688           b = Ftan(ppmt,PhiLocal);                664           b = Ftan(ppmt,PhiLocal);
689         }                                         665         }
690       if (PhiLocal>xe0)                           666       if (PhiLocal>xe0)
691         {                                         667         {
692           a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l    668           a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
693         }                                         669         }
694       else                                        670       else
695         {                                         671         {
696           a = Ftan(p0t,PhiLocal);                 672           a = Ftan(p0t,PhiLocal);
697         }                                         673         }
698     }                                             674     }
699   else                                            675   else
700     {                                             676     {
701       b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml    677       b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,PhiLocal));
702       a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi    678       a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,PhiLocal));
703     }                                             679     }
704   G4double nr =0.;                                680   G4double nr =0.;
705                                                   681 
706   if (b>a)                                        682   if (b>a)
707     {                                             683     {
708       nr = 1./b;                                  684       nr = 1./b;
709     }                                             685     }
710   else                                            686   else
711     {                                             687     {
712       nr = 1./a;                                  688       nr = 1./a;
713     }                                             689     }
714                                                   690 
715   G4double r1,r2=0.;                              691   G4double r1,r2=0.;
716   G4double r3 =-1.;                               692   G4double r3 =-1.;
717   do                                              693   do
718     {                                             694     {
719       r1 = G4UniformRand();                       695       r1 = G4UniformRand();
720       r2 = G4UniformRand();                       696       r2 = G4UniformRand();
721       //value = r2*pi;                         << 697       value = r2*pi;
722       value = 2.*r2*pi;                        << 
723       r3 = nr*(a*cos(value)*cos(value) + b*sin    698       r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
724     }while(r1>r3);                                699     }while(r1>r3);
725                                                   700 
726   return value;                                   701   return value;
727 }                                                 702 }
728                                                   703 
729 //....oooOO0OOooo........oooOO0OOooo........oo << 
730                                                   704 
731 G4double G4LivermorePolarizedGammaConversionMo    705 G4double G4LivermorePolarizedGammaConversionModel::Poli
732 (G4double a, G4double b, G4double c, G4double     706 (G4double a, G4double b, G4double c, G4double x)
733 {                                                 707 {
734   G4double value=0.;                              708   G4double value=0.;
735   if(x>0.)                                        709   if(x>0.)
736     {                                             710     {
737       value =(a + b/x + c/(x*x*x));               711       value =(a + b/x + c/(x*x*x));
738     }                                             712     }
739   else                                            713   else
740     {                                             714     {
741       //G4cout << "ERROR in Poli! " << G4endl;    715       //G4cout << "ERROR in Poli! " << G4endl;
742     }                                             716     }
743   return value;                                   717   return value;
744 }                                                 718 }
745                                                << 
746 //....oooOO0OOooo........oooOO0OOooo........oo << 
747                                                << 
748 G4double G4LivermorePolarizedGammaConversionMo    719 G4double G4LivermorePolarizedGammaConversionModel::Fln
749 (G4double a, G4double b, G4double x)              720 (G4double a, G4double b, G4double x)
750 {                                                 721 {
751   G4double value=0.;                              722   G4double value=0.;
752   if(x>0.)                                        723   if(x>0.)
753     {                                             724     {
754       value =(a*log(x)-b);                        725       value =(a*log(x)-b);
755     }                                             726     }
756   else                                            727   else
757     {                                             728     {
758       //G4cout << "ERROR in Fln! " << G4endl;     729       //G4cout << "ERROR in Fln! " << G4endl;
759     }                                             730     }
760   return value;                                   731   return value;
761 }                                                 732 }
762                                                   733 
763 //....oooOO0OOooo........oooOO0OOooo........oo << 
764                                                   734 
765 G4double G4LivermorePolarizedGammaConversionMo    735 G4double G4LivermorePolarizedGammaConversionModel::Encu
766 (G4double* p_p1, G4double* p_p2, G4double x0)     736 (G4double* p_p1, G4double* p_p2, G4double x0)
767 {                                                 737 {
768   G4int i=0;                                      738   G4int i=0;
769   G4double fx = 1.;                               739   G4double fx = 1.;
770   G4double x = x0;                                740   G4double x = x0;
771   const G4double xmax = 3.0;                      741   const G4double xmax = 3.0;
772                                                   742 
773   for(i=0; i<100; ++i)                            743   for(i=0; i<100; ++i)
774     {                                             744     {
775       fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p    745       fx = (Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x))/
776   (Fdlor(p_p1,x) - Fdtan(p_p2,x));                746   (Fdlor(p_p1,x) - Fdtan(p_p2,x));
777       x -= fx;                                    747       x -= fx;
778       if(x > xmax) { return xmax; }               748       if(x > xmax) { return xmax; }
                                                   >> 749       //      x -= (Flor(p_p1, x)*Glor(p_p1,x) - Ftan(p_p2, x))/
                                                   >> 750       //  (Fdlor(p_p1,x) - Fdtan(p_p2,x));
                                                   >> 751       // fx = Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x);
                                                   >> 752       // G4cout << std::fabs(fx) << " " << i << " " << x << "dentro ENCU " << G4endl;
779       if(std::fabs(fx) <= x*1.0e-6) { break; }    753       if(std::fabs(fx) <= x*1.0e-6) { break; }
780     }                                             754     } 
781                                                   755 
782   if(x < 0.0) { x = 0.0; }                        756   if(x < 0.0) { x = 0.0; }
783   return x;                                       757   return x;
784 }                                                 758 }
785                                                   759 
786 //....oooOO0OOooo........oooOO0OOooo........oo << 
787                                                   760 
788 G4double G4LivermorePolarizedGammaConversionMo    761 G4double G4LivermorePolarizedGammaConversionModel::Flor(G4double* p_p1, G4double x)
789 {                                                 762 {
790   G4double value =0.;                             763   G4double value =0.;
                                                   >> 764   // G4double y0 = p_p1[0];
                                                   >> 765   // G4double A = p_p1[1];
791   G4double w = p_p1[2];                           766   G4double w = p_p1[2];
792   G4double xc = p_p1[3];                          767   G4double xc = p_p1[3];
793                                                   768 
794   value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));       769   value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
795   return value;                                   770   return value;
796 }                                                 771 }
797                                                   772 
798 //....oooOO0OOooo........oooOO0OOooo........oo << 
799                                                   773 
800 G4double G4LivermorePolarizedGammaConversionMo    774 G4double G4LivermorePolarizedGammaConversionModel::Glor(G4double* p_p1, G4double x)
801 {                                                 775 {
802   G4double value =0.;                             776   G4double value =0.;
803   G4double y0 = p_p1[0];                          777   G4double y0 = p_p1[0];
804   G4double A = p_p1[1];                           778   G4double A = p_p1[1];
805   G4double w = p_p1[2];                           779   G4double w = p_p1[2];
806   G4double xc = p_p1[3];                          780   G4double xc = p_p1[3];
807                                                   781 
808   value = (y0 *pi*(w*w +  4.*(x-xc)*(x-xc)) +     782   value = (y0 *pi*(w*w +  4.*(x-xc)*(x-xc)) + 2.*A*w);
809   return value;                                   783   return value;
810 }                                                 784 }
811                                                   785 
812 //....oooOO0OOooo........oooOO0OOooo........oo << 
813                                                   786 
814 G4double G4LivermorePolarizedGammaConversionMo    787 G4double G4LivermorePolarizedGammaConversionModel::Fdlor(G4double* p_p1, G4double x)
815 {                                                 788 {
816   G4double value =0.;                             789   G4double value =0.;
                                                   >> 790   //G4double y0 = p_p1[0];
817   G4double A = p_p1[1];                           791   G4double A = p_p1[1];
818   G4double w = p_p1[2];                           792   G4double w = p_p1[2];
819   G4double xc = p_p1[3];                          793   G4double xc = p_p1[3];
820                                                   794 
821   value = (-16.*A*w*(x-xc))/                      795   value = (-16.*A*w*(x-xc))/
822     (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*    796     (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
823   return value;                                   797   return value;
824 }                                                 798 }
825                                                   799 
826 //....oooOO0OOooo........oooOO0OOooo........oo << 
827                                                   800 
828 G4double G4LivermorePolarizedGammaConversionMo    801 G4double G4LivermorePolarizedGammaConversionModel::Fintlor(G4double* p_p1, G4double x)
829 {                                                 802 {
830   G4double value =0.;                             803   G4double value =0.;
831   G4double y0 = p_p1[0];                          804   G4double y0 = p_p1[0];
832   G4double A = p_p1[1];                           805   G4double A = p_p1[1];
833   G4double w = p_p1[2];                           806   G4double w = p_p1[2];
834   G4double xc = p_p1[3];                          807   G4double xc = p_p1[3];
835                                                   808 
836   value = y0*x + A*atan( 2*(x-xc)/w) / pi;        809   value = y0*x + A*atan( 2*(x-xc)/w) / pi;
837   return value;                                   810   return value;
838 }                                                 811 }
839                                                   812 
840                                                   813 
841 G4double G4LivermorePolarizedGammaConversionMo    814 G4double G4LivermorePolarizedGammaConversionModel::Finvlor(G4double* p_p1, G4double x, G4double r)
842 {                                                 815 {
843   G4double value = 0.;                            816   G4double value = 0.;
844   G4double nor = 0.;                              817   G4double nor = 0.;
                                                   >> 818   //G4double y0 = p_p1[0];
                                                   >> 819   //  G4double A = p_p1[1];
845   G4double w = p_p1[2];                           820   G4double w = p_p1[2];
846   G4double xc = p_p1[3];                          821   G4double xc = p_p1[3];
847                                                   822 
848   nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.    823   nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w);
849   value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(    824   value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w));
850                                                   825 
851   return value;                                   826   return value;
852 }                                                 827 }
853                                                   828 
854 //....oooOO0OOooo........oooOO0OOooo........oo << 
855                                                   829 
856 G4double G4LivermorePolarizedGammaConversionMo    830 G4double G4LivermorePolarizedGammaConversionModel::Ftan(G4double* p_p1, G4double x)
857 {                                                 831 {
858   G4double value =0.;                             832   G4double value =0.;
859   G4double a = p_p1[0];                           833   G4double a = p_p1[0];
860   G4double b = p_p1[1];                           834   G4double b = p_p1[1];
861                                                   835 
862   value = a /(x-b);                               836   value = a /(x-b);
863   return value;                                   837   return value;
864 }                                                 838 }
865                                                   839 
866 //....oooOO0OOooo........oooOO0OOooo........oo << 
867                                                   840 
868 G4double G4LivermorePolarizedGammaConversionMo    841 G4double G4LivermorePolarizedGammaConversionModel::Fdtan(G4double* p_p1, G4double x)
869 {                                                 842 {
870   G4double value =0.;                             843   G4double value =0.;
871   G4double a = p_p1[0];                           844   G4double a = p_p1[0];
872   G4double b = p_p1[1];                           845   G4double b = p_p1[1];
873                                                   846 
874   value = -1.*a / ((x-b)*(x-b));                  847   value = -1.*a / ((x-b)*(x-b));
875   return value;                                   848   return value;
876 }                                                 849 }
877                                                   850 
878 //....oooOO0OOooo........oooOO0OOooo........oo << 
879                                                   851 
880 G4double G4LivermorePolarizedGammaConversionMo    852 G4double G4LivermorePolarizedGammaConversionModel::Finttan(G4double* p_p1, G4double x)
881 {                                                 853 {
882   G4double value =0.;                             854   G4double value =0.;
883   G4double a = p_p1[0];                           855   G4double a = p_p1[0];
884   G4double b = p_p1[1];                           856   G4double b = p_p1[1];
885                                                   857 
                                                   >> 858 
886   value = a*log(b-x);                             859   value = a*log(b-x);
887   return value;                                   860   return value;
888 }                                                 861 }
889                                                   862 
890 //....oooOO0OOooo........oooOO0OOooo........oo << 
891                                                   863 
892 G4double G4LivermorePolarizedGammaConversionMo    864 G4double G4LivermorePolarizedGammaConversionModel::Finvtan(G4double* p_p1, G4double cnor, G4double r)
893 {                                                 865 {
894   G4double value =0.;                             866   G4double value =0.;
895   G4double a = p_p1[0];                           867   G4double a = p_p1[0];
896   G4double b = p_p1[1];                           868   G4double b = p_p1[1];
897                                                   869 
898   value = b*(1-G4Exp(r*cnor/a));               << 870   value = b*(1-exp(r*cnor/a));
899                                                   871 
900   return value;                                   872   return value;
901 }                                                 873 }
902                                                   874 
                                                   >> 875 
                                                   >> 876 
                                                   >> 877 
903 //....oooOO0OOooo........oooOO0OOooo........oo    878 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
904                                                   879 
905 G4ThreeVector G4LivermorePolarizedGammaConvers    880 G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector(G4ThreeVector& a)
906 {                                                 881 {
907   G4double dx = a.x();                            882   G4double dx = a.x();
908   G4double dy = a.y();                            883   G4double dy = a.y();
909   G4double dz = a.z();                            884   G4double dz = a.z();
910   G4double x = dx < 0.0 ? -dx : dx;               885   G4double x = dx < 0.0 ? -dx : dx;
911   G4double y = dy < 0.0 ? -dy : dy;               886   G4double y = dy < 0.0 ? -dy : dy;
912   G4double z = dz < 0.0 ? -dz : dz;               887   G4double z = dz < 0.0 ? -dz : dz;
913   if (x < y) {                                    888   if (x < y) {
914     return x < z ? G4ThreeVector(-dy,dx,0) : G    889     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
915   }else{                                          890   }else{
916     return y < z ? G4ThreeVector(dz,0,-dx) : G    891     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
917   }                                               892   }
918 }                                                 893 }
919                                                   894 
920 //....oooOO0OOooo........oooOO0OOooo........oo    895 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
921                                                   896 
922 G4ThreeVector G4LivermorePolarizedGammaConvers    897 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization(G4ThreeVector& direction0)
923 {                                                 898 {
924   G4ThreeVector d0 = direction0.unit();           899   G4ThreeVector d0 = direction0.unit();
925   G4ThreeVector a1 = SetPerpendicularVector(d0    900   G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
926   G4ThreeVector a0 = a1.unit(); // unit vector    901   G4ThreeVector a0 = a1.unit(); // unit vector
927                                                   902 
928   G4double rand1 = G4UniformRand();               903   G4double rand1 = G4UniformRand();
929                                                   904   
930   G4double angle = twopi*rand1; // random pola    905   G4double angle = twopi*rand1; // random polar angle
931   G4ThreeVector b0 = d0.cross(a0); // cross pr    906   G4ThreeVector b0 = d0.cross(a0); // cross product
932                                                   907   
933   G4ThreeVector c;                                908   G4ThreeVector c;
934                                                   909   
935   c.setX(std::cos(angle)*(a0.x())+std::sin(ang    910   c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
936   c.setY(std::cos(angle)*(a0.y())+std::sin(ang    911   c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
937   c.setZ(std::cos(angle)*(a0.z())+std::sin(ang    912   c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
938                                                   913   
939   G4ThreeVector c0 = c.unit();                    914   G4ThreeVector c0 = c.unit();
940                                                   915 
941   return c0;                                   << 916   return c0;
                                                   >> 917   
942 }                                                 918 }
943                                                   919 
944 //....oooOO0OOooo........oooOO0OOooo........oo    920 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
945                                                   921 
946 G4ThreeVector G4LivermorePolarizedGammaConvers    922 G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization
947 (const G4ThreeVector& gammaDirection, const G4    923 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
948 {                                                 924 {
                                                   >> 925 
949   //                                              926   // 
950   // The polarization of a photon is always pe    927   // The polarization of a photon is always perpendicular to its momentum direction.
951   // Therefore this function removes those vec    928   // Therefore this function removes those vector component of gammaPolarization, which
952   // points in direction of gammaDirection        929   // points in direction of gammaDirection
953   //                                              930   //
954   // Mathematically we search the projection o    931   // Mathematically we search the projection of the vector a on the plane E, where n is the
955   // plains normal vector.                        932   // plains normal vector.
956   // The basic equation can be found in each g    933   // The basic equation can be found in each geometry book (e.g. Bronstein):
957   // p = a - (a o n)/(n o n)*n                    934   // p = a - (a o n)/(n o n)*n
958                                                   935   
959   return gammaPolarization - gammaPolarization    936   return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
960 }                                                 937 }
961                                                   938 
962 //....oooOO0OOooo........oooOO0OOooo........oo    939 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
963                                                   940 
                                                   >> 941 
964 void G4LivermorePolarizedGammaConversionModel:    942 void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange
965     (G4ThreeVector& direction0,G4ThreeVector&     943     (G4ThreeVector& direction0,G4ThreeVector& direction1,
966      G4ThreeVector& polarization0)                944      G4ThreeVector& polarization0)
967 {                                                 945 {
968   // direction0 is the original photon directi    946   // direction0 is the original photon direction ---> z
969   // polarization0 is the original photon pola    947   // polarization0 is the original photon polarization ---> x
970   // need to specify y axis in the real refere    948   // need to specify y axis in the real reference frame ---> y 
971   G4ThreeVector Axis_Z0 = direction0.unit();      949   G4ThreeVector Axis_Z0 = direction0.unit();
972   G4ThreeVector Axis_X0 = polarization0.unit()    950   G4ThreeVector Axis_X0 = polarization0.unit();
973   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_    951   G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
974                                                   952   
975   G4double direction_x = direction1.getX();       953   G4double direction_x = direction1.getX();
976   G4double direction_y = direction1.getY();       954   G4double direction_y = direction1.getY();
977   G4double direction_z = direction1.getZ();       955   G4double direction_z = direction1.getZ();
978                                                   956   
979   direction1 = (direction_x*Axis_X0 + directio << 957   direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 +  direction_z*Axis_Z0).unit();
                                                   >> 958   
980 }                                                 959 }
981                                                   960 
982 //....oooOO0OOooo........oooOO0OOooo........oo << 
983                                                   961 
984 void G4LivermorePolarizedGammaConversionModel: << 962 
985                       const G4ParticleDefiniti << 963 
986                       G4int Z)                 << 
987 {                                              << 
988   G4AutoLock l(&LivermorePolarizedGammaConvers << 
989   if(!data[Z]) { ReadData(Z); }                << 
990   l.unlock();                                  << 
991 }                                              << 
992                                                   964