| 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: G4LivermorePolarizedComptonModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $
>> 27 // GEANT4 tag $Name: geant4-09-02 $
26 // 28 //
27 // Authors: G.Depaola & F.Longo <<
28 // <<
29 // History: <<
30 // ------- <<
31 // <<
32 // 05 Apr 2021 J Allison added quantum entan <<
33 // If the photons have been "tagged" as "quant <<
34 // G4eplusAnnihilation for annihilation into 2 <<
35 // here if - and only if - both photons suffer <<
36 // predictions from Pryce and Ward, Nature No <<
37 // Physical Review 73 (1948) p.440. Experiment <<
38 // entanglement in the MeV regime and its appl <<
39 // D. Watts, J. Allison et al., Nature Communi <<
40 // https://doi.org/10.1038/s41467-021-22907-5. <<
41 // <<
42 // 02 May 2009 S Incerti as V. Ivanchenko pr <<
43 // <<
44 // Cleanup initialisation and generation of se <<
45 // - apply internal high-ener <<
46 // - do not apply low-energy <<
47 // - remove GetMeanFreePath m <<
48 // - added protection against <<
49 // - use G4ElementSelector <<
50 29
51 #include "G4LivermorePolarizedComptonModel.hh" 30 #include "G4LivermorePolarizedComptonModel.hh"
52 #include "G4PhysicalConstants.hh" <<
53 #include "G4SystemOfUnits.hh" <<
54 #include "G4AutoLock.hh" <<
55 #include "G4Electron.hh" <<
56 #include "G4ParticleChangeForGamma.hh" <<
57 #include "G4LossTableManager.hh" <<
58 #include "G4VAtomDeexcitation.hh" <<
59 #include "G4AtomicShell.hh" <<
60 #include "G4Gamma.hh" <<
61 #include "G4ShellData.hh" <<
62 #include "G4DopplerProfile.hh" <<
63 #include "G4Log.hh" <<
64 #include "G4Exp.hh" <<
65 #include "G4Pow.hh" <<
66 #include "G4LogLogInterpolation.hh" <<
67 #include "G4PhysicsModelCatalog.hh" <<
68 #include "G4EntanglementAuxInfo.hh" <<
69 #include "G4eplusAnnihilationEntanglementClipB <<
70 31
71 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72 33
73 using namespace std; 34 using namespace std;
74 namespace { G4Mutex LivermorePolarizedComptonM <<
75 <<
76 <<
77 G4PhysicsFreeVector* G4LivermorePolarizedCompt <<
78 G4ShellData* G4LivermorePolarizedCompton <<
79 G4DopplerProfile* G4LivermorePolarizedCompton <<
80 G4CompositeEMDataSet* G4LivermorePolarizedComp <<
81 35
82 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 37
84 G4LivermorePolarizedComptonModel::G4LivermoreP << 38 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
85 :G4VEmModel(nam),isInitialised(false) << 39 const G4String& nam)
86 { << 40 :G4VEmModel(nam),isInitialised(false)
87 verboseLevel= 1; << 41 {
>> 42 lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
>> 43 highEnergyLimit = 100 * GeV;
>> 44 SetLowEnergyLimit(lowEnergyLimit);
>> 45 SetHighEnergyLimit(highEnergyLimit);
>> 46
>> 47 verboseLevel= 0;
88 // Verbosity scale: 48 // Verbosity scale:
89 // 0 = nothing 49 // 0 = nothing
90 // 1 = warning for energy non-conservation 50 // 1 = warning for energy non-conservation
91 // 2 = details of energy budget 51 // 2 = details of energy budget
92 // 3 = calculation of cross sections, file o 52 // 3 = calculation of cross sections, file openings, sampling of atoms
93 // 4 = entering in methods 53 // 4 = entering in methods
94 54
95 if( verboseLevel>1 ) << 55 G4cout << "Livermore Polarized Compton is constructed " << G4endl
96 G4cout << "Livermore Polarized Compton is << 56 << "Energy range: "
97 << 57 << lowEnergyLimit / keV << " keV - "
98 //Mark this model as "applicable" for atomic << 58 << highEnergyLimit / GeV << " GeV"
99 SetDeexcitationFlag(true); << 59 << G4endl;
100 << 60
101 fParticleChange = nullptr; <<
102 fAtomDeexcitation = nullptr; <<
103 fEntanglementModelID = G4PhysicsModelCatalog <<
104 } 61 }
105 62
106 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 64
108 G4LivermorePolarizedComptonModel::~G4Livermore 65 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
109 { 66 {
110 if(IsMaster()) { << 67 delete meanFreePathTable;
111 delete shellData; << 68 delete crossSectionHandler;
112 shellData = nullptr; << 69 delete scatterFunctionData;
113 delete profileData; <<
114 profileData = nullptr; <<
115 delete scatterFunctionData; <<
116 scatterFunctionData = nullptr; <<
117 for(G4int i=0; i<maxZ; ++i) { <<
118 if(data[i]) { <<
119 delete data[i]; <<
120 data[i] = nullptr; <<
121 } <<
122 } <<
123 } <<
124 } 70 }
125 71
126 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127 73
128 void G4LivermorePolarizedComptonModel::Initial 74 void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
129 const G 75 const G4DataVector& cuts)
130 { 76 {
131 if (verboseLevel > 1) << 77 if (verboseLevel > 3)
132 G4cout << "Calling G4LivermorePolarizedCom 78 G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
133 79
134 // Initialise element selector << 80 InitialiseElementSelectors(particle,cuts);
135 if(IsMaster()) { <<
136 // Access to elements <<
137 const char* path = G4FindDataDir("G4LEDATA <<
138 81
139 G4ProductionCutsTable* theCoupleTable = << 82 // Energy limits
140 G4ProductionCutsTable::GetProductionCuts << 83
141 << 84 if (LowEnergyLimit() < lowEnergyLimit)
142 G4int numOfCouples = (G4int)theCoupleTable << 85 {
143 << 86 G4cout << "G4LivermorePolarizedComptonModel: low energy limit increased from " <<
144 for(G4int i=0; i<numOfCouples; ++i) { << 87 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
145 const G4Material* material = << 88 SetLowEnergyLimit(lowEnergyLimit);
146 theCoupleTable->GetMaterialCutsCouple( <<
147 const G4ElementVector* theElementVector <<
148 std::size_t nelm = material->GetNumberOf <<
149 <<
150 for (std::size_t j=0; j<nelm; ++j) { <<
151 G4int Z = G4lrint((*theElementVector)[ <<
152 if(Z < 1) { Z = 1; } <<
153 else if(Z > maxZ){ Z = maxZ; } <<
154 <<
155 if( (!data[Z]) ) { ReadData(Z, path); <<
156 } <<
157 } 89 }
158 << 90
159 // For Doppler broadening << 91 if (HighEnergyLimit() > highEnergyLimit)
160 if(!shellData) { << 92 {
161 shellData = new G4ShellData(); << 93 G4cout << "G4LivermorePolarizedComptonModel: high energy limit decreased from " <<
162 shellData->SetOccupancyData(); << 94 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
163 G4String file = "/doppler/shell-doppler" << 95 SetHighEnergyLimit(highEnergyLimit);
164 shellData->LoadData(file); <<
165 } 96 }
166 if(!profileData) { profileData = new G4Dop <<
167 97
168 // Scattering Function << 98 // Reading of data files - all materials are read
169 if(!scatterFunctionData) <<
170 { <<
171 <<
172 G4VDataSetAlgorithm* scatterInterpolation = <<
173 G4String scatterFile = "comp/ce-sf-"; <<
174 scatterFunctionData = new G4CompositeEMDataS <<
175 scatterFunctionData->LoadData(scatterFile); <<
176 } <<
177 <<
178 InitialiseElementSelectors(particle, cuts) <<
179 } <<
180 <<
181 if (verboseLevel > 2) { <<
182 G4cout << "Loaded cross section files" << <<
183 } <<
184 <<
185 if( verboseLevel>1 ) { <<
186 G4cout << "G4LivermoreComptonModel is init <<
187 << "Energy range: " <<
188 << LowEnergyLimit() / eV << " eV - " <<
189 << HighEnergyLimit() / GeV << " GeV" <<
190 << G4endl; <<
191 } <<
192 // <<
193 if(isInitialised) { return; } <<
194 99
195 fParticleChange = GetParticleChangeForGamma( << 100 crossSectionHandler = new G4CrossSectionHandler;
196 fAtomDeexcitation = G4LossTableManager::Ins << 101 crossSectionHandler->Clear();
197 isInitialised = true; << 102 G4String crossSectionFile = "comp/ce-cs-";
198 } << 103 crossSectionHandler->LoadData(crossSectionFile);
199 104
>> 105 meanFreePathTable = 0;
>> 106 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
200 107
201 void G4LivermorePolarizedComptonModel::Initial << 108 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
202 G4VEmModel* masterModel) << 109 G4String scatterFile = "comp/ce-sf-";
203 { << 110 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
204 SetElementSelectors(masterModel->GetElementS << 111 scatterFunctionData->LoadData(scatterFile);
205 } <<
206 112
207 //....oooOO0OOooo........oooOO0OOooo........oo << 113 // For Doppler broadening
>> 114 shellData.SetOccupancyData();
>> 115 G4String file = "/doppler/shell-doppler";
>> 116 shellData.LoadData(file);
208 117
209 void G4LivermorePolarizedComptonModel::ReadDat << 118 //
210 { << 119 if (verboseLevel > 2)
211 if (verboseLevel > 1) << 120 G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
212 { << 121
213 G4cout << "G4LivermorePolarizedComptonMo << 122 G4cout << "Livermore Polarized Compton model is initialized " << G4endl
214 << G4endl; << 123 << "Energy range: "
215 } << 124 << LowEnergyLimit() / keV << " keV - "
216 if(data[Z]) { return; } << 125 << HighEnergyLimit() / GeV << " GeV"
217 const char* datadir = path; << 126 << G4endl;
218 if(!datadir) << 127
219 { << 128 //
220 datadir = G4FindDataDir("G4LEDATA"); <<
221 if(!datadir) <<
222 { <<
223 G4Exception("G4LivermorePolarizedComptonMo <<
224 "em0006",FatalException, <<
225 "Environment variable G4LEDATA not d <<
226 return; <<
227 } <<
228 } <<
229 <<
230 data[Z] = new G4PhysicsFreeVector(); <<
231 129
232 std::ostringstream ost; << 130 if(isInitialised) return;
233 ost << datadir << "/livermore/comp/ce-cs-" < << 131
234 std::ifstream fin(ost.str().c_str()); << 132 if(pParticleChange)
235 << 133 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
236 if( !fin.is_open()) << 134 else
237 { << 135 fParticleChange = new G4ParticleChangeForGamma();
238 G4ExceptionDescription ed; << 136
239 ed << "G4LivermorePolarizedComptonModel << 137 isInitialised = true;
240 << "> is not opened!" << G4endl; <<
241 G4Exception("G4LivermoreComptonModel::Re <<
242 "em0003",FatalException, <<
243 ed,"G4LEDATA version should be G4EMLOW8. <<
244 return; <<
245 } else { <<
246 if(verboseLevel > 3) { <<
247 G4cout << "File " << ost.str() <<
248 << " is opened by G4LivermorePolarizedC <<
249 } <<
250 data[Z]->Retrieve(fin, true); <<
251 data[Z]->ScaleVector(MeV, MeV*barn); <<
252 } <<
253 fin.close(); <<
254 } 138 }
255 139
256 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
257 141
258 G4double G4LivermorePolarizedComptonModel::Com 142 G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
259 const G 143 const G4ParticleDefinition*,
260 G 144 G4double GammaEnergy,
261 G 145 G4double Z, G4double,
262 G 146 G4double, G4double)
263 { 147 {
264 if (verboseLevel > 3) 148 if (verboseLevel > 3)
265 G4cout << "Calling ComputeCrossSectionPerA 149 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
266 150
267 G4double cs = 0.0; << 151 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
268 <<
269 if (GammaEnergy < LowEnergyLimit()) <<
270 return 0.0; <<
271 <<
272 G4int intZ = G4lrint(Z); <<
273 if(intZ < 1 || intZ > maxZ) { return cs; } <<
274 <<
275 G4PhysicsFreeVector* pv = data[intZ]; <<
276 <<
277 // if element was not initialised <<
278 // do initialisation safely for MT mode <<
279 if(!pv) <<
280 { <<
281 InitialiseForElement(0, intZ); <<
282 pv = data[intZ]; <<
283 if(!pv) { return cs; } <<
284 } <<
285 <<
286 G4int n = G4int(pv->GetVectorLength() - 1); <<
287 G4double e1 = pv->Energy(0); <<
288 G4double e2 = pv->Energy(n); <<
289 <<
290 if(GammaEnergy <= e1) { cs = GammaEnerg <<
291 else if(GammaEnergy <= e2) { cs = pv->Value( <<
292 else if(GammaEnergy > e2) { cs = pv->Value( <<
293 <<
294 return cs; 152 return cs;
295 <<
296 } 153 }
297 154
298 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299 156
300 void G4LivermorePolarizedComptonModel::SampleS 157 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
301 const G4MaterialCutsCouple* co 158 const G4MaterialCutsCouple* couple,
302 const G4DynamicParticle* aDyna 159 const G4DynamicParticle* aDynamicGamma,
303 G4double, 160 G4double,
304 G4double) 161 G4double)
305 { 162 {
306 // The scattered gamma energy is sampled acc 163 // The scattered gamma energy is sampled according to Klein - Nishina formula.
307 // The random number techniques of Butcher & 164 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
308 // GEANT4 internal units 165 // GEANT4 internal units
309 // 166 //
310 // Note : Effects due to binding of atomic e 167 // Note : Effects due to binding of atomic electrons are negliged.
311 168
312 if (verboseLevel > 3) 169 if (verboseLevel > 3)
313 G4cout << "Calling SampleSecondaries() of 170 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
314 171
315 G4double gammaEnergy0 = aDynamicGamma->GetKi 172 G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
316 <<
317 // do nothing below the threshold <<
318 // should never get here because the XS is z <<
319 if (gammaEnergy0 < LowEnergyLimit()) <<
320 return ; <<
321 <<
322 G4ThreeVector gammaPolarization0 = aDynamicG 173 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
323 174
324 // Protection: a polarisation parallel to th 175 // Protection: a polarisation parallel to the
325 // direction causes problems; 176 // direction causes problems;
326 // in that case find a random polarization 177 // in that case find a random polarization
>> 178
327 G4ThreeVector gammaDirection0 = aDynamicGamm 179 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
328 180
329 // Make sure that the polarization vector is 181 // Make sure that the polarization vector is perpendicular to the
330 // gamma direction. If not 182 // gamma direction. If not
>> 183
331 if(!(gammaPolarization0.isOrthogonal(gammaDi 184 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
332 { // only for testing now 185 { // only for testing now
333 gammaPolarization0 = GetRandomPolarizati 186 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
334 } 187 }
335 else 188 else
336 { 189 {
337 if ( gammaPolarization0.howOrthogonal(ga 190 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
338 { 191 {
339 gammaPolarization0 = GetPerpendicularPolar 192 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
340 } 193 }
341 } 194 }
>> 195
342 // End of Protection 196 // End of Protection
343 197
>> 198 // Within energy limit?
>> 199
>> 200 if(gammaEnergy0 <= lowEnergyLimit)
>> 201 {
>> 202 fParticleChange->ProposeTrackStatus(fStopAndKill);
>> 203 fParticleChange->SetProposedKineticEnergy(0.);
>> 204 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
>> 205 // SI - IS THE FOLLOWING RETURN NECESSARY ?
>> 206 return;
>> 207 }
>> 208
344 G4double E0_m = gammaEnergy0 / electron_mass 209 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
345 210
346 // Select randomly one element in the curren 211 // Select randomly one element in the current material
347 //G4int Z = crossSectionHandler->SelectRando << 212
348 const G4ParticleDefinition* particle = aDyn << 213 G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
349 const G4Element* elm = SelectRandomAtom(coup <<
350 G4int Z = (G4int)elm->GetZ(); <<
351 214
352 // Sample the energy and the polarization of 215 // Sample the energy and the polarization of the scattered photon
>> 216
353 G4double epsilon, epsilonSq, onecost, sinThe 217 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
354 218
355 G4double epsilon0Local = 1./(1. + 2*E0_m); << 219 G4double epsilon0 = 1./(1. + 2*E0_m);
356 G4double epsilon0Sq = epsilon0Local*epsilon0 << 220 G4double epsilon0Sq = epsilon0*epsilon0;
357 G4double alpha1 = - G4Log(epsilon0Local); << 221 G4double alpha1 = - std::log(epsilon0);
358 G4double alpha2 = 0.5*(1.- epsilon0Sq); 222 G4double alpha2 = 0.5*(1.- epsilon0Sq);
359 223
360 G4double wlGamma = h_Planck*c_light/gammaEne 224 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
361 G4double gammaEnergy1; 225 G4double gammaEnergy1;
362 G4ThreeVector gammaDirection1; 226 G4ThreeVector gammaDirection1;
363 227
364 do { 228 do {
365 if ( alpha1/(alpha1+alpha2) > G4UniformRan 229 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
366 { 230 {
367 epsilon = G4Exp(-alpha1*G4UniformRand()); << 231 epsilon = std::exp(-alpha1*G4UniformRand());
368 epsilonSq = epsilon*epsilon; 232 epsilonSq = epsilon*epsilon;
369 } 233 }
370 else 234 else
371 { 235 {
372 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4 236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
373 epsilon = std::sqrt(epsilonSq); 237 epsilon = std::sqrt(epsilonSq);
374 } 238 }
375 239
376 onecost = (1.- epsilon)/(epsilon*E0_m); 240 onecost = (1.- epsilon)/(epsilon*E0_m);
377 sinThetaSqr = onecost*(2.-onecost); 241 sinThetaSqr = onecost*(2.-onecost);
378 242
379 // Protection 243 // Protection
380 if (sinThetaSqr > 1.) 244 if (sinThetaSqr > 1.)
381 { 245 {
382 G4cout 246 G4cout
383 << " -- Warning -- G4LivermorePolarizedCom 247 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
384 << "sin(theta)**2 = " 248 << "sin(theta)**2 = "
385 << sinThetaSqr 249 << sinThetaSqr
386 << "; set to 1" 250 << "; set to 1"
387 << G4endl; 251 << G4endl;
388 sinThetaSqr = 1.; 252 sinThetaSqr = 1.;
389 } 253 }
390 if (sinThetaSqr < 0.) 254 if (sinThetaSqr < 0.)
391 { 255 {
392 G4cout 256 G4cout
393 << " -- Warning -- G4LivermorePolarizedCom 257 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
394 << "sin(theta)**2 = " 258 << "sin(theta)**2 = "
395 << sinThetaSqr 259 << sinThetaSqr
396 << "; set to 0" 260 << "; set to 0"
397 << G4endl; 261 << G4endl;
398 sinThetaSqr = 0.; 262 sinThetaSqr = 0.;
399 } 263 }
400 // End protection 264 // End protection
401 265
402 G4double x = std::sqrt(onecost/2.) / (wlG 266 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
403 G4double scatteringFunction = scatterFunct 267 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
404 greject = (1. - epsilon*sinThetaSqr/(1.+ e 268 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
405 269
406 } while(greject < G4UniformRand()*Z); 270 } while(greject < G4UniformRand()*Z);
407 271
408 272
409 // ***************************************** 273 // ****************************************************
410 // Phi determination 274 // Phi determination
411 // ***************************************** 275 // ****************************************************
>> 276
412 G4double phi = SetPhi(epsilon,sinThetaSqr); 277 G4double phi = SetPhi(epsilon,sinThetaSqr);
413 278
414 // 279 //
415 // scattered gamma angles. ( Z - axis along 280 // scattered gamma angles. ( Z - axis along the parent gamma)
416 // 281 //
>> 282
417 G4double cosTheta = 1. - onecost; 283 G4double cosTheta = 1. - onecost;
418 284
419 // Protection 285 // Protection
>> 286
420 if (cosTheta > 1.) 287 if (cosTheta > 1.)
421 { 288 {
422 G4cout 289 G4cout
423 << " -- Warning -- G4LivermorePolarizedCompt 290 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
424 << "cosTheta = " 291 << "cosTheta = "
425 << cosTheta 292 << cosTheta
426 << "; set to 1" 293 << "; set to 1"
427 << G4endl; 294 << G4endl;
428 cosTheta = 1.; 295 cosTheta = 1.;
429 } 296 }
430 if (cosTheta < -1.) 297 if (cosTheta < -1.)
431 { 298 {
432 G4cout 299 G4cout
433 << " -- Warning -- G4LivermorePolarizedCompt 300 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
434 << "cosTheta = " 301 << "cosTheta = "
435 << cosTheta 302 << cosTheta
436 << "; set to -1" 303 << "; set to -1"
437 << G4endl; 304 << G4endl;
438 cosTheta = -1.; 305 cosTheta = -1.;
439 } 306 }
440 // End protection 307 // End protection
441 << 308
>> 309
442 G4double sinTheta = std::sqrt (sinThetaSqr); 310 G4double sinTheta = std::sqrt (sinThetaSqr);
443 311
444 // Protection 312 // Protection
445 if (sinTheta > 1.) 313 if (sinTheta > 1.)
446 { 314 {
447 G4cout 315 G4cout
448 << " -- Warning -- G4LivermorePolarizedCompt 316 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
449 << "sinTheta = " 317 << "sinTheta = "
450 << sinTheta 318 << sinTheta
451 << "; set to 1" 319 << "; set to 1"
452 << G4endl; 320 << G4endl;
453 sinTheta = 1.; 321 sinTheta = 1.;
454 } 322 }
455 if (sinTheta < -1.) 323 if (sinTheta < -1.)
456 { 324 {
457 G4cout 325 G4cout
458 << " -- Warning -- G4LivermorePolarizedCompt 326 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
459 << "sinTheta = " 327 << "sinTheta = "
460 << sinTheta 328 << sinTheta
461 << "; set to -1" 329 << "; set to -1"
462 << G4endl; 330 << G4endl;
463 sinTheta = -1.; 331 sinTheta = -1.;
464 } 332 }
465 // End protection 333 // End protection
466 334
467 // Check for entanglement and re-sample phi <<
468 <<
469 const auto* auxInfo <<
470 = fParticleChange->GetCurrentTrack()->GetAux <<
471 if (auxInfo) { <<
472 const auto* entanglementAuxInfo = dynamic_ <<
473 if (entanglementAuxInfo) { <<
474 auto* clipBoard = dynamic_cast<G4eplusAn <<
475 (entanglementAuxInfo->GetEntanglementCli <<
476 if (clipBoard) { <<
477 // This is an entangled photon from ep <<
478 // If this is the first scatter of the <<
479 // phi on the clipboard. <<
480 // If this is the first scatter of the <<
481 // phi of the first scatter of the fir <<
482 // theta of the second photon, to samp <<
483 if (clipBoard->IsTrack1Measurement()) <<
484 // Check we have the relevant track. <<
485 // necessary but I want to be sure t <<
486 // entangled system are properly pai <<
487 // Note: the tracking manager pops t <<
488 // will rely on that. (If not, the <<
489 // more complicated to ensure we mat <<
490 // So our track 1 is clipboard track <<
491 if (clipBoard->GetTrackB() == fParti <<
492 // This is the first scatter of th <<
493 // // Debug <<
494 // auto* track1 = fPart <<
495 // G4cout <<
496 // << "This is the firs <<
497 // << "\nTrack: " << tr <<
498 // << ", Parent: " << t <<
499 // << ", Name: " << cli <<
500 // << G4endl; <<
501 // // End debug <<
502 clipBoard->ResetTrack1Measurement( <<
503 // Store cos(theta),phi of first p <<
504 clipBoard->SetComptonCosTheta1(cos <<
505 clipBoard->SetComptonPhi1(phi); <<
506 } <<
507 } else if (clipBoard->IsTrack2Measurem <<
508 // Check we have the relevant track. <<
509 // Remember our track 2 is clipboard <<
510 if (clipBoard->GetTrackA() == fParti <<
511 // This is the first scatter of th <<
512 // // Debug <<
513 // auto* track2 = fPart <<
514 // G4cout <<
515 // << "This is the firs <<
516 // << "\nTrack: " << tr <<
517 // << ", Parent: " << t <<
518 // << ", Name: " << cli <<
519 // << G4endl; <<
520 // // End debug <<
521 clipBoard->ResetTrack2Measurement( <<
522 <<
523 // Get cos(theta),phi of first pho <<
524 const G4double& cosTheta1 = clipBo <<
525 const G4double& phi1 = clipBoard-> <<
526 // For clarity make aliases for th <<
527 const G4double& cosTheta2 = cosThe <<
528 G4double& phi2 = phi; <<
529 // G4cout << "cosTheta1 <<
530 // G4cout << "cosTheta2 <<
531 <<
532 // Re-sample phi <<
533 // Draw the difference of azimutha <<
534 // A + B * cos(2*deltaPhi), or rat <<
535 // C = A / (A + |B|) and D = B / ( <<
536 const G4double sin2Theta1 = 1.-cos <<
537 const G4double sin2Theta2 = 1.-cos <<
538 <<
539 // Pryce and Ward, Nature No 4065 <<
540 auto* g4Pow = G4Pow::GetInstance() <<
541 const G4double A = <<
542 ((g4Pow->powN(1.-cosTheta1,3))+2.) <<
543 ((g4Pow->powN(2.-cosTheta1,3)*g4Po <<
544 const G4double B = -(sin2Theta1*si <<
545 ((g4Pow->powN(2.-cosTheta1,2)*g4Po <<
546 <<
547 // // Snyder et al, Physical Re <<
548 // // (This is an alternative f <<
549 // const G4double& k0 = gammaEn <<
550 // const G4double k1 = k0/(2.-c <<
551 // const G4double k2 = k0/(2.-c <<
552 // const G4double gamma1 = k1/k <<
553 // const G4double gamma2 = k2/k <<
554 // const G4double A1 = gamma1*g <<
555 // const G4double B1 = 2.*sin2T <<
556 // // That's A1 + B1*sin2(delta <<
557 // const G4double A = A1 + 0.5* <<
558 // const G4double B = -0.5*B1; <<
559 <<
560 const G4double maxValue = A + std: <<
561 const G4double C = A / maxValue; <<
562 const G4double D = B / maxValue; <<
563 // G4cout << "A,B,C,D: " << A < <<
564 <<
565 // Sample delta phi <<
566 G4double deltaPhi; <<
567 const G4int maxCount = 999999; <<
568 G4int iCount = 0; <<
569 for (; iCount < maxCount; ++iCount <<
570 deltaPhi = twopi * G4UniformRand <<
571 if (G4UniformRand() < C + D * co <<
572 } <<
573 if (iCount >= maxCount ) { <<
574 G4cout << "G4LivermorePolarizedC <<
575 << "Re-sampled delta phi not fou <<
576 << " tries - carrying on anyway. <<
577 } <<
578 <<
579 // Thus, the desired second photon <<
580 phi2 = deltaPhi - phi1 + halfpi; <<
581 // The minus sign is in above stat <<
582 // annihilation photons are in opp <<
583 // are measured in the opposite di <<
584 // halfpi is added for the followi <<
585 // In this function phi is relativ <<
586 // SystemOfRefChange below. We kno <<
587 // the polarisations of the two an <<
588 // to each other, i.e., halfpi dif <<
589 // Furthermore, only sin(phi) and <<
590 // need to place any range constra <<
591 // if (phi2 > pi) { <<
592 // phi2 -= twopi; <<
593 // } <<
594 // if (phi2 < -pi) { <<
595 // phi2 += twopi; <<
596 // } <<
597 } <<
598 } <<
599 } <<
600 } <<
601 } <<
602 <<
603 // End of entanglement <<
604 335
605 G4double dirx = sinTheta*std::cos(phi); 336 G4double dirx = sinTheta*std::cos(phi);
606 G4double diry = sinTheta*std::sin(phi); 337 G4double diry = sinTheta*std::sin(phi);
607 G4double dirz = cosTheta ; 338 G4double dirz = cosTheta ;
608 << 339
>> 340
609 // oneCosT , eom 341 // oneCosT , eom
610 342
611 // Doppler broadening - Method based on: 343 // Doppler broadening - Method based on:
612 // Y. Namito, S. Ban and H. Hirayama, 344 // Y. Namito, S. Ban and H. Hirayama,
613 // "Implementation of the Doppler Broadening 345 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
614 // NIM A 349, pp. 489-494, 1994 346 // NIM A 349, pp. 489-494, 1994
615 347
616 // Maximum number of sampling iterations 348 // Maximum number of sampling iterations
617 static G4int maxDopplerIterations = 1000; << 349
>> 350 G4int maxDopplerIterations = 1000;
618 G4double bindingE = 0.; 351 G4double bindingE = 0.;
619 G4double photonEoriginal = epsilon * gammaEn 352 G4double photonEoriginal = epsilon * gammaEnergy0;
620 G4double photonE = -1.; 353 G4double photonE = -1.;
621 G4int iteration = 0; 354 G4int iteration = 0;
622 G4double eMax = gammaEnergy0; 355 G4double eMax = gammaEnergy0;
623 356
624 G4int shellIdx = 0; <<
625 <<
626 if (verboseLevel > 3) { <<
627 G4cout << "Started loop to sample broading <<
628 } <<
629 <<
630 do 357 do
631 { 358 {
632 iteration++; 359 iteration++;
633 // Select shell based on shell occupancy 360 // Select shell based on shell occupancy
634 shellIdx = shellData->SelectRandomShell( << 361 G4int shell = shellData.SelectRandomShell(Z);
635 bindingE = shellData->BindingEnergy(Z,sh << 362 bindingE = shellData.BindingEnergy(Z,shell);
636 363
637 if (verboseLevel > 3) { <<
638 G4cout << "Shell ID= " << shellIdx <<
639 << " Ebind(keV)= " << bindingE/keV << <<
640 } <<
641 eMax = gammaEnergy0 - bindingE; 364 eMax = gammaEnergy0 - bindingE;
642 << 365
643 // Randomly sample bound electron moment 366 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
644 G4double pSample = profileData->RandomSe << 367 G4double pSample = profileData.RandomSelectMomentum(Z,shell);
645 <<
646 if (verboseLevel > 3) { <<
647 G4cout << "pSample= " << pSample << G4e <<
648 } <<
649 // Rescale from atomic units 368 // Rescale from atomic units
650 G4double pDoppler = pSample * fine_struc 369 G4double pDoppler = pSample * fine_structure_const;
651 G4double pDoppler2 = pDoppler * pDoppler 370 G4double pDoppler2 = pDoppler * pDoppler;
652 G4double var2 = 1. + onecost * E0_m; 371 G4double var2 = 1. + onecost * E0_m;
653 G4double var3 = var2*var2 - pDoppler2; 372 G4double var3 = var2*var2 - pDoppler2;
654 G4double var4 = var2 - pDoppler2 * cosTh 373 G4double var4 = var2 - pDoppler2 * cosTheta;
655 G4double var = var4*var4 - var3 + pDoppl 374 G4double var = var4*var4 - var3 + pDoppler2 * var3;
656 if (var > 0.) 375 if (var > 0.)
657 { 376 {
658 G4double varSqrt = std::sqrt(var); 377 G4double varSqrt = std::sqrt(var);
659 G4double scale = gammaEnergy0 / var3; 378 G4double scale = gammaEnergy0 / var3;
660 // Random select either root 379 // Random select either root
661 if (G4UniformRand() < 0.5) photonE = (var4 380 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
662 else photonE = (var4 + varSqrt) * scale; 381 else photonE = (var4 + varSqrt) * scale;
663 } 382 }
664 else 383 else
665 { 384 {
666 photonE = -1.; 385 photonE = -1.;
667 } 386 }
668 } while ( iteration <= maxDopplerIterations 387 } while ( iteration <= maxDopplerIterations &&
669 (photonE < 0. || photonE > eMax || phot 388 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
670 << 389
671 // End of recalculation of photon energy wit 390 // End of recalculation of photon energy with Doppler broadening
672 // Revert to original if maximum number of i 391 // Revert to original if maximum number of iterations threshold has been reached
673 if (iteration >= maxDopplerIterations) 392 if (iteration >= maxDopplerIterations)
674 { 393 {
675 photonE = photonEoriginal; 394 photonE = photonEoriginal;
676 bindingE = 0.; 395 bindingE = 0.;
677 } 396 }
678 397
679 gammaEnergy1 = photonE; 398 gammaEnergy1 = photonE;
680 399
681 // 400 //
682 // update G4VParticleChange for the scattere 401 // update G4VParticleChange for the scattered photon
683 // 402 //
>> 403
>> 404 // gammaEnergy1 = epsilon*gammaEnergy0;
>> 405
>> 406
684 // New polarization 407 // New polarization
>> 408
685 G4ThreeVector gammaPolarization1 = SetNewPol 409 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
686 sinThetaSqr, 410 sinThetaSqr,
687 phi, 411 phi,
688 cosTheta); 412 cosTheta);
689 413
690 // Set new direction 414 // Set new direction
691 G4ThreeVector tmpDirection1( dirx,diry,dirz 415 G4ThreeVector tmpDirection1( dirx,diry,dirz );
692 gammaDirection1 = tmpDirection1; 416 gammaDirection1 = tmpDirection1;
693 417
694 // Change reference frame. 418 // Change reference frame.
>> 419
695 SystemOfRefChange(gammaDirection0,gammaDirec 420 SystemOfRefChange(gammaDirection0,gammaDirection1,
696 gammaPolarization0,gammaPolarization1) 421 gammaPolarization0,gammaPolarization1);
697 422
698 if (gammaEnergy1 > 0.) 423 if (gammaEnergy1 > 0.)
699 { 424 {
700 fParticleChange->SetProposedKineticEnerg 425 fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
701 fParticleChange->ProposeMomentumDirectio 426 fParticleChange->ProposeMomentumDirection( gammaDirection1 );
702 fParticleChange->ProposePolarization( ga 427 fParticleChange->ProposePolarization( gammaPolarization1 );
703 } 428 }
704 else 429 else
705 { 430 {
706 gammaEnergy1 = 0.; <<
707 fParticleChange->SetProposedKineticEnerg 431 fParticleChange->SetProposedKineticEnergy(0.) ;
708 fParticleChange->ProposeTrackStatus(fSto 432 fParticleChange->ProposeTrackStatus(fStopAndKill);
709 } 433 }
710 434
711 // 435 //
712 // kinematic of the scattered electron 436 // kinematic of the scattered electron
713 // 437 //
>> 438
714 G4double ElecKineEnergy = gammaEnergy0 - gam 439 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
715 440
716 // SI -protection against negative final ene << 441 // SI - Removed range test
717 // like in G4LivermoreComptonModel.cc << 442
718 if(ElecKineEnergy < 0.0) { <<
719 fParticleChange->ProposeLocalEnergyDeposit <<
720 return; <<
721 } <<
722 <<
723 G4double ElecMomentum = std::sqrt(ElecKineEn 443 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
724 444
725 G4ThreeVector ElecDirection((gammaEnergy0 * 445 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
726 gammaEnergy1 * gammaDirection1) * ( 446 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
727 447
728 G4DynamicParticle* dp = <<
729 new G4DynamicParticle (G4Electron::Electro <<
730 fvect->push_back(dp); <<
731 <<
732 // sample deexcitation <<
733 // <<
734 if (verboseLevel > 3) { <<
735 G4cout << "Started atomic de-excitation " <<
736 } <<
737 <<
738 if(fAtomDeexcitation && iteration < maxDoppl <<
739 G4int index = couple->GetIndex(); <<
740 if(fAtomDeexcitation->CheckDeexcitationAct <<
741 std::size_t nbefore = fvect->size(); <<
742 G4AtomicShellEnumerator as = G4AtomicShe <<
743 const G4AtomicShell* shell = fAtomDeexci <<
744 fAtomDeexcitation->GenerateParticles(fve <<
745 std::size_t nafter = fvect->size(); <<
746 if(nafter > nbefore) { <<
747 for (std::size_t i=nbefore; i<nafter; ++i) { <<
748 //Check if there is enough residual energy <<
749 if (bindingE >= ((*fvect)[i])->GetKineticE <<
750 { <<
751 //Ok, this is a valid secondary: <<
752 bindingE -= ((*fvect)[i])->GetKineticE <<
753 } <<
754 else <<
755 { <<
756 //Invalid secondary: not enough energy <<
757 //Keep its energy in the local deposit <<
758 delete (*fvect)[i]; <<
759 (*fvect)[i]=0; <<
760 } <<
761 } <<
762 } <<
763 } <<
764 } <<
765 //This should never happen <<
766 if(bindingE < 0.0) <<
767 G4Exception("G4LivermoreComptonModel::Samp <<
768 "em2050",FatalException,"Negative local en <<
769 <<
770 fParticleChange->ProposeLocalEnergyDeposit(b 448 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
>> 449
>> 450 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
>> 451 fvect->push_back(dp);
771 452
772 } 453 }
773 454
774 //....oooOO0OOooo........oooOO0OOooo........oo 455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
775 456
776 G4double G4LivermorePolarizedComptonModel::Set 457 G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
777 G4double sinSqrTh) 458 G4double sinSqrTh)
778 { 459 {
779 G4double rand1; 460 G4double rand1;
780 G4double rand2; 461 G4double rand2;
781 G4double phiProbability; 462 G4double phiProbability;
782 G4double phi; 463 G4double phi;
783 G4double a, b; 464 G4double a, b;
784 465
785 do 466 do
786 { 467 {
787 rand1 = G4UniformRand(); 468 rand1 = G4UniformRand();
788 rand2 = G4UniformRand(); 469 rand2 = G4UniformRand();
789 phiProbability=0.; 470 phiProbability=0.;
790 phi = twopi*rand1; 471 phi = twopi*rand1;
791 472
792 a = 2*sinSqrTh; 473 a = 2*sinSqrTh;
793 b = energyRate + 1/energyRate; 474 b = energyRate + 1/energyRate;
794 475
795 phiProbability = 1 - (a/b)*(std::cos(phi 476 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
>> 477
>> 478
>> 479
796 } 480 }
797 while ( rand2 > phiProbability ); 481 while ( rand2 > phiProbability );
798 return phi; 482 return phi;
799 } 483 }
800 484
801 485
802 //....oooOO0OOooo........oooOO0OOooo........oo 486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
803 487
804 G4ThreeVector G4LivermorePolarizedComptonModel 488 G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
805 { 489 {
806 G4double dx = a.x(); 490 G4double dx = a.x();
807 G4double dy = a.y(); 491 G4double dy = a.y();
808 G4double dz = a.z(); 492 G4double dz = a.z();
809 G4double x = dx < 0.0 ? -dx : dx; 493 G4double x = dx < 0.0 ? -dx : dx;
810 G4double y = dy < 0.0 ? -dy : dy; 494 G4double y = dy < 0.0 ? -dy : dy;
811 G4double z = dz < 0.0 ? -dz : dz; 495 G4double z = dz < 0.0 ? -dz : dz;
812 if (x < y) { 496 if (x < y) {
813 return x < z ? G4ThreeVector(-dy,dx,0) : G 497 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
814 }else{ 498 }else{
815 return y < z ? G4ThreeVector(dz,0,-dx) : G 499 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
816 } 500 }
817 } 501 }
818 502
819 //....oooOO0OOooo........oooOO0OOooo........oo 503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
820 504
821 G4ThreeVector G4LivermorePolarizedComptonModel 505 G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
822 { 506 {
823 G4ThreeVector d0 = direction0.unit(); 507 G4ThreeVector d0 = direction0.unit();
824 G4ThreeVector a1 = SetPerpendicularVector(d0 508 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
825 G4ThreeVector a0 = a1.unit(); // unit vector 509 G4ThreeVector a0 = a1.unit(); // unit vector
826 510
827 G4double rand1 = G4UniformRand(); 511 G4double rand1 = G4UniformRand();
828 512
829 G4double angle = twopi*rand1; // random pola 513 G4double angle = twopi*rand1; // random polar angle
830 G4ThreeVector b0 = d0.cross(a0); // cross pr 514 G4ThreeVector b0 = d0.cross(a0); // cross product
831 515
832 G4ThreeVector c; 516 G4ThreeVector c;
833 517
834 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 518 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
835 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 519 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
836 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 520 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
837 521
838 G4ThreeVector c0 = c.unit(); 522 G4ThreeVector c0 = c.unit();
839 523
840 return c0; 524 return c0;
>> 525
841 } 526 }
842 527
843 //....oooOO0OOooo........oooOO0OOooo........oo 528 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
844 529
845 G4ThreeVector G4LivermorePolarizedComptonModel 530 G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
846 (const G4ThreeVector& gammaDirection, const G4 531 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
847 { 532 {
>> 533
848 // 534 //
849 // The polarization of a photon is always pe 535 // The polarization of a photon is always perpendicular to its momentum direction.
850 // Therefore this function removes those vec 536 // Therefore this function removes those vector component of gammaPolarization, which
851 // points in direction of gammaDirection 537 // points in direction of gammaDirection
852 // 538 //
853 // Mathematically we search the projection o 539 // Mathematically we search the projection of the vector a on the plane E, where n is the
854 // plains normal vector. 540 // plains normal vector.
855 // The basic equation can be found in each g 541 // The basic equation can be found in each geometry book (e.g. Bronstein):
856 // p = a - (a o n)/(n o n)*n 542 // p = a - (a o n)/(n o n)*n
857 543
858 return gammaPolarization - gammaPolarization 544 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
859 } 545 }
860 546
861 //....oooOO0OOooo........oooOO0OOooo........oo 547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
862 548
863 G4ThreeVector G4LivermorePolarizedComptonModel 549 G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
864 G4double sinSqrTh, 550 G4double sinSqrTh,
865 G4double phi, 551 G4double phi,
866 G4double costheta) 552 G4double costheta)
867 { 553 {
868 G4double rand1; 554 G4double rand1;
869 G4double rand2; 555 G4double rand2;
870 G4double cosPhi = std::cos(phi); 556 G4double cosPhi = std::cos(phi);
871 G4double sinPhi = std::sin(phi); 557 G4double sinPhi = std::sin(phi);
872 G4double sinTheta = std::sqrt(sinSqrTh); 558 G4double sinTheta = std::sqrt(sinSqrTh);
873 G4double cosSqrPhi = cosPhi*cosPhi; 559 G4double cosSqrPhi = cosPhi*cosPhi;
874 // G4double cossqrth = 1.-sinSqrTh; 560 // G4double cossqrth = 1.-sinSqrTh;
875 // G4double sinsqrphi = sinPhi*sinPhi; 561 // G4double sinsqrphi = sinPhi*sinPhi;
876 G4double normalisation = std::sqrt(1. - cosS 562 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
877 563
>> 564
878 // Determination of Theta 565 // Determination of Theta
>> 566
>> 567 // ---- MGP ---- Commented out the following 3 lines to avoid compilation
>> 568 // warnings (unused variables)
>> 569 // G4double thetaProbability;
879 G4double theta; 570 G4double theta;
>> 571 // G4double a, b;
>> 572 // G4double cosTheta;
>> 573
>> 574 /*
>> 575
>> 576 depaola method
>> 577
>> 578 do
>> 579 {
>> 580 rand1 = G4UniformRand();
>> 581 rand2 = G4UniformRand();
>> 582 thetaProbability=0.;
>> 583 theta = twopi*rand1;
>> 584 a = 4*normalisation*normalisation;
>> 585 b = (epsilon + 1/epsilon) - 2;
>> 586 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
>> 587 cosTheta = std::cos(theta);
>> 588 }
>> 589 while ( rand2 > thetaProbability );
>> 590
>> 591 G4double cosBeta = cosTheta;
>> 592
>> 593 */
>> 594
880 595
881 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) 596 // Dan Xu method (IEEE TNS, 52, 1160 (2005))
>> 597
882 rand1 = G4UniformRand(); 598 rand1 = G4UniformRand();
883 rand2 = G4UniformRand(); 599 rand2 = G4UniformRand();
884 600
885 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi 601 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
886 { 602 {
887 if (rand2<0.5) 603 if (rand2<0.5)
888 theta = pi/2.0; 604 theta = pi/2.0;
889 else 605 else
890 theta = 3.0*pi/2.0; 606 theta = 3.0*pi/2.0;
891 } 607 }
892 else 608 else
893 { 609 {
894 if (rand2<0.5) 610 if (rand2<0.5)
895 theta = 0; 611 theta = 0;
896 else 612 else
897 theta = pi; 613 theta = pi;
898 } 614 }
899 G4double cosBeta = std::cos(theta); 615 G4double cosBeta = std::cos(theta);
900 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 616 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
901 617
902 G4ThreeVector gammaPolarization1; 618 G4ThreeVector gammaPolarization1;
903 619
904 G4double xParallel = normalisation*cosBeta; 620 G4double xParallel = normalisation*cosBeta;
905 G4double yParallel = -(sinSqrTh*cosPhi*sinPh 621 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
906 G4double zParallel = -(costheta*sinTheta*cos 622 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
907 G4double xPerpendicular = 0.; 623 G4double xPerpendicular = 0.;
908 G4double yPerpendicular = (costheta)*sinBeta 624 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
909 G4double zPerpendicular = -(sinTheta*sinPhi) 625 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
910 626
911 G4double xTotal = (xParallel + xPerpendicula 627 G4double xTotal = (xParallel + xPerpendicular);
912 G4double yTotal = (yParallel + yPerpendicula 628 G4double yTotal = (yParallel + yPerpendicular);
913 G4double zTotal = (zParallel + zPerpendicula 629 G4double zTotal = (zParallel + zPerpendicular);
914 630
915 gammaPolarization1.setX(xTotal); 631 gammaPolarization1.setX(xTotal);
916 gammaPolarization1.setY(yTotal); 632 gammaPolarization1.setY(yTotal);
917 gammaPolarization1.setZ(zTotal); 633 gammaPolarization1.setZ(zTotal);
918 634
919 return gammaPolarization1; 635 return gammaPolarization1;
>> 636
920 } 637 }
921 638
922 //....oooOO0OOooo........oooOO0OOooo........oo 639 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
923 640
924 void G4LivermorePolarizedComptonModel::SystemO 641 void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
925 G4ThreeVector& direction1, 642 G4ThreeVector& direction1,
926 G4ThreeVector& polarization0, 643 G4ThreeVector& polarization0,
927 G4ThreeVector& polarization1) 644 G4ThreeVector& polarization1)
928 { 645 {
929 // direction0 is the original photon directi 646 // direction0 is the original photon direction ---> z
930 // polarization0 is the original photon pola 647 // polarization0 is the original photon polarization ---> x
931 // need to specify y axis in the real refere 648 // need to specify y axis in the real reference frame ---> y
932 G4ThreeVector Axis_Z0 = direction0.unit(); 649 G4ThreeVector Axis_Z0 = direction0.unit();
933 G4ThreeVector Axis_X0 = polarization0.unit() 650 G4ThreeVector Axis_X0 = polarization0.unit();
934 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 651 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
935 652
936 G4double direction_x = direction1.getX(); 653 G4double direction_x = direction1.getX();
937 G4double direction_y = direction1.getY(); 654 G4double direction_y = direction1.getY();
938 G4double direction_z = direction1.getZ(); 655 G4double direction_z = direction1.getZ();
939 656
940 direction1 = (direction_x*Axis_X0 + directio 657 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
941 G4double polarization_x = polarization1.getX 658 G4double polarization_x = polarization1.getX();
942 G4double polarization_y = polarization1.getY 659 G4double polarization_y = polarization1.getY();
943 G4double polarization_z = polarization1.getZ 660 G4double polarization_z = polarization1.getZ();
944 661
945 polarization1 = (polarization_x*Axis_X0 + po 662 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
946 663
947 } 664 }
948 665
>> 666 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
949 667
950 //....oooOO0OOooo........oooOO0OOooo........oo << 668 G4double G4LivermorePolarizedComptonModel::GetMeanFreePath(const G4Track& track,
951 //....oooOO0OOooo........oooOO0OOooo........oo << 669 G4double,
952 << 670 G4ForceCondition*)
953 void << 671 {
954 G4LivermorePolarizedComptonModel::InitialiseFo << 672 const G4DynamicParticle* photon = track.GetDynamicParticle();
955 G4int Z) << 673 G4double energy = photon->GetKineticEnergy();
956 { << 674 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
957 G4AutoLock l(&LivermorePolarizedComptonModel << 675 size_t materialIndex = couple->GetIndex();
958 if(!data[Z]) { ReadData(Z); } << 676 G4double meanFreePath;
959 l.unlock(); << 677 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
>> 678 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
>> 679 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
>> 680 return meanFreePath;
960 } 681 }
>> 682
>> 683
961 684