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.5.p1)


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