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