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