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