| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // -------------------------------------------
27 //
28 // Geant4 Class file
29 //
30 // File name: G4PolarizedAnnihilationModel
31 //
32 // Author: Andreas Schaelicke
33 //
34 // Class Description:
35 // Implementation of polarized gamma Annihil
36
37 #include "G4PolarizedAnnihilationModel.hh"
38
39 #include "G4Gamma.hh"
40 #include "G4ParticleChangeForGamma.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4PolarizationHelper.hh"
43 #include "G4PolarizationManager.hh"
44 #include "G4PolarizedAnnihilationXS.hh"
45 #include "G4StokesVector.hh"
46 #include "G4TrackStatus.hh"
47
48 G4PolarizedAnnihilationModel::G4PolarizedAnnih
49 const G4ParticleDefinition* p, const G4Strin
50 : G4eeToTwoGammaModel(p, nam)
51 , fCrossSectionCalculator(nullptr)
52 , fParticleChange(nullptr)
53 , fVerboseLevel(0)
54 {
55 fCrossSectionCalculator = new G4PolarizedAn
56 fBeamPolarization = G4StokesVector::Z
57 fTargetPolarization = G4StokesVector::Z
58 fFinalGamma1Polarization = G4StokesVector::Z
59 fFinalGamma2Polarization = G4StokesVector::Z
60 }
61
62 //....oooOO0OOooo........oooOO0OOooo........oo
63 G4PolarizedAnnihilationModel::~G4PolarizedAnni
64 {
65 delete fCrossSectionCalculator;
66 }
67
68 //....oooOO0OOooo........oooOO0OOooo........oo
69 void G4PolarizedAnnihilationModel::Initialise(
70
71 {
72 G4eeToTwoGammaModel::Initialise(part, dv);
73 if(fParticleChange)
74 {
75 return;
76 }
77 fParticleChange = GetParticleChangeForGamma(
78 }
79
80 //....oooOO0OOooo........oooOO0OOooo........oo
81 G4double G4PolarizedAnnihilationModel::Compute
82 G4double kinEnergy)
83 {
84 // cross section from base model
85 G4double xs = G4eeToTwoGammaModel::ComputeCr
86
87 G4double polzz = fBeamPolarization.z() * fTa
88 G4double poltt = fBeamPolarization.x() * fTa
89 fBeamPolarization.y() * fTa
90 if(polzz != 0 || poltt != 0)
91 {
92 G4double xval, lasym, tasym;
93 ComputeAsymmetriesPerElectron(kinEnergy, x
94 xs *= (1. + polzz * lasym + poltt * tasym)
95 }
96
97 return xs;
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101 void G4PolarizedAnnihilationModel::ComputeAsym
102 G4double ene, G4double& valueX, G4double& va
103 {
104 // *** calculate asymmetries
105 G4double gam = 1. + ene / electron_mass_c2;
106 G4double xs0 = fCrossSectionCalculator->Tota
107 0., 1., gam, G4StokesVector::ZERO, G4Stoke
108 G4double xsA = fCrossSectionCalculator->Tota
109 0., 1., gam, G4StokesVector::P3, G4StokesV
110 G4double xsT1 = fCrossSectionCalculator->Tot
111 0., 1., gam, G4StokesVector::P1, G4StokesV
112 G4double xsT2 = fCrossSectionCalculator->Tot
113 0., 1., gam, G4StokesVector::P2, G4StokesV
114 G4double xsT = 0.5 * (xsT1 + xsT2);
115
116 valueX = xs0;
117 valueA = xsA / xs0 - 1.;
118 valueT = xsT / xs0 - 1.;
119
120 if((valueA < -1) || (1 < valueA))
121 {
122 G4ExceptionDescription ed;
123 ed << " ERROR PolarizedAnnihilationPS::Com
124 ed << " something wrong in total cross sec
125 ed << " LONG: " << valueX << "\t" << value
126 << " energy = " << gam << G4endl;
127 G4Exception("G4PolarizedAnnihilationModel:
128 "pol004", JustWarning, ed);
129 }
130 if((valueT < -1) || (1 < valueT))
131 {
132 G4ExceptionDescription ed;
133 ed << " ERROR PolarizedAnnihilationPS::Com
134 ed << " something wrong in total cross sec
135 ed << " TRAN: " << valueX << "\t" << value
136 << " energy = " << gam << G4endl;
137 G4Exception("G4PolarizedAnnihilationModel:
138 "pol005", JustWarning, ed);
139 }
140 }
141
142 void G4PolarizedAnnihilationModel::SampleSecon
143 std::vector<G4DynamicParticle*>* fvect, cons
144 const G4DynamicParticle* dp, G4double, G4dou
145 {
146 const G4Track* aTrack = fParticleChange->Get
147
148 // kill primary
149 fParticleChange->SetProposedKineticEnergy(0.
150 fParticleChange->ProposeTrackStatus(fStopAnd
151
152 // V.Ivanchenko add protection against zero
153 G4double PositKinEnergy = dp->GetKineticEner
154
155 if(PositKinEnergy == 0.0)
156 {
157 G4double cosTeta = 2. * G4UniformRand() -
158 G4double sinTeta = std::sqrt((1.0 - cosTet
159 G4double phi = twopi * G4UniformRand()
160 G4ThreeVector dir(sinTeta * std::cos(phi),
161 cosTeta);
162 fvect->push_back(
163 new G4DynamicParticle(G4Gamma::Gamma(),
164 fvect->push_back(
165 new G4DynamicParticle(G4Gamma::Gamma(),
166 return;
167 }
168
169 // *** obtain and save target and beam polar
170 G4PolarizationManager* polarizationManager =
171 G4PolarizationManager::GetInstance();
172
173 // obtain polarization of the beam
174 fBeamPolarization = G4StokesVector(aTrack->G
175
176 // obtain polarization of the media
177 G4VPhysicalVolume* aPVolume = aTrack->Get
178 G4LogicalVolume* aLVolume = aPVolume->G
179 const G4bool targetIsPolarized = polarizatio
180 fTargetPolarization = polarizationManager->G
181
182 if(fVerboseLevel >= 1)
183 {
184 G4cout << "G4PolarizedComptonModel::Sample
185 << aLVolume->GetName() << G4endl;
186 }
187
188 // transfer target electron polarization in
189 if(targetIsPolarized)
190 fTargetPolarization.rotateUz(dp->GetMoment
191
192 G4ParticleMomentum PositDirection = dp->GetM
193
194 // polar asymmetry:
195 G4double polarization = fBeamPolarization.p3
196
197 G4double gamam1 = PositKinEnergy / electron_
198 G4double gama = gamam1 + 1., gamap1 = gamam1
199 G4double sqgrate = std::sqrt(gamam1 / gamap1
200 sqg2m1 = std::sqrt(gamam1 * gamap1
201
202 // limits of the energy sampling
203 G4double epsilmin = 0.5 - sqgrate, epsilmax
204 G4double epsilqot = epsilmax / epsilmin;
205
206 // sample the energy rate of the created gam
207 // note: for polarized partices, the actual
208 // will depend on the energy, and the
209 G4double epsil;
210 G4double gmax = 1. + std::fabs(polarization)
211
212 fCrossSectionCalculator->Initialize(epsilmin
213 fTargetP
214 if(fCrossSectionCalculator->DiceEpsilon() <
215 {
216 G4ExceptionDescription ed;
217 ed << "ERROR in PolarizedAnnihilationPS::P
218 << "epsilmin DiceRoutine not appropriat
219 << fCrossSectionCalculator->DiceEpsilon
220 G4Exception("G4PolarizedAnnihilationModel:
221 JustWarning, ed);
222 }
223
224 fCrossSectionCalculator->Initialize(epsilmax
225 fTargetP
226 if(fCrossSectionCalculator->DiceEpsilon() <
227 {
228 G4ExceptionDescription ed;
229 ed << "ERROR in PolarizedAnnihilationPS::P
230 << "epsilmax DiceRoutine not appropriat
231 << fCrossSectionCalculator->DiceEpsilon
232 G4Exception("G4PolarizedAnnihilationModel:
233 JustWarning, ed);
234 }
235
236 G4int ncount = 0;
237 G4double trejectmax = 0.;
238 G4double treject;
239
240 do
241 {
242 epsil = epsilmin * std::pow(epsilqot, G4Un
243
244 fCrossSectionCalculator->Initialize(epsil,
245 fTarge
246
247 treject = fCrossSectionCalculator->DiceEps
248 treject *= epsil;
249
250 if(treject > gmax || treject < 0.)
251 {
252 G4ExceptionDescription ed;
253 ed << "ERROR in PolarizedAnnihilationPS:
254 << " eps (" << epsil
255 << ") rejection does not work properl
256 G4Exception("G4PolarizedAnnihilationMode
257 JustWarning, ed);
258 }
259 ++ncount;
260 if(treject > trejectmax)
261 trejectmax = treject;
262 if(ncount > 1000)
263 {
264 G4ExceptionDescription ed;
265 ed << "WARNING in PolarizedAnnihilation
266 << "eps dicing very inefficient =" <<
267 << treject / gmax << ". For secondar
268 << ncount << G4endl;
269 G4Exception("G4PolarizedAnnihilationMode
270 JustWarning, ed);
271 break;
272 }
273
274 // Loop checking, 03-Aug-2015, Vladimir Iv
275 } while(treject < gmax * G4UniformRand());
276
277 // scattered Gamma angles. ( Z - axis along
278 G4double cost = (epsil * gamap1 - 1.) / (eps
279 G4double sint = std::sqrt((1. + cost) * (1.
280 G4double phi = 0.;
281 G4double beamTrans =
282 std::sqrt(sqr(fBeamPolarization.p1()) + sq
283 G4double targetTrans =
284 std::sqrt(sqr(fTargetPolarization.p1()) +
285
286 do
287 {
288 phi = twopi * G4UniformRand();
289 fCrossSectionCalculator->Initialize(epsil,
290 fTarge
291
292 G4double gdiced = fCrossSectionCalculator-
293 gdiced += fCrossSectionCalculator->getVar(
294 fTargetPolarization.p3();
295 gdiced += 1. *
296 (std::fabs(fCrossSectionCalculat
297 std::fabs(fCrossSectionCalculat
298 beamTrans * targetTrans;
299 gdiced += 1. * std::fabs(fCrossSectionCalc
300 (std::fabs(fBeamPolarization.p3(
301 std::fabs(fTargetPolarization.p
302
303 G4double gdist = fCrossSectionCalculator->
304 gdist += fCrossSectionCalculator->getVar(3
305 fTargetPolarization.p3();
306 gdist += fCrossSectionCalculator->getVar(1
307 (std::cos(phi) * fBeamPolarizatio
308 std::sin(phi) * fBeamPolarizatio
309 (std::cos(phi) * fTargetPolarizat
310 std::sin(phi) * fTargetPolarizat
311 gdist += fCrossSectionCalculator->getVar(2
312 (std::cos(phi) * fBeamPolarizatio
313 std::sin(phi) * fBeamPolarizatio
314 (std::cos(phi) * fTargetPolarizat
315 std::sin(phi) * fTargetPolarizat
316 gdist +=
317 fCrossSectionCalculator->getVar(4) *
318 (std::cos(phi) * fBeamPolarization.p3()
319 std::cos(phi) * fBeamPolarization.p1()
320 std::sin(phi) * fBeamPolarization.p3()
321 std::sin(phi) * fBeamPolarization.p2()
322
323 treject = gdist / gdiced;
324 if(treject > 1. + 1.e-10 || treject < 0)
325 {
326 G4ExceptionDescription ed;
327 ed << "!!!ERROR in PolarizedAnnihilation
328 << " phi rejection does not work prop
329 G4cout << " gdiced = " << gdiced << G4en
330 G4cout << " gdist = " << gdist << G4endl
331 G4cout << " epsil = " << epsil << G4endl
332 G4Exception("G4PolarizedAnnihilationMode
333 JustWarning, ed);
334 }
335
336 if(treject < 1.e-3)
337 {
338 G4ExceptionDescription ed;
339 ed << "!!!ERROR in PolarizedAnnihilation
340 << " phi rejection does not work prop
341 G4cout << " gdiced=" << gdiced << " gd
342 G4cout << " epsil = " << epsil << G4endl
343 G4Exception("G4PolarizedAnnihilationMode
344 JustWarning, ed);
345 }
346
347 // Loop checking, 03-Aug-2015, Vladimir Iv
348 } while(treject < G4UniformRand());
349
350 G4double dirx = sint * std::cos(phi);
351 G4double diry = sint * std::sin(phi);
352 G4double dirz = cost;
353
354 // kinematic of the created pair
355 G4double TotalAvailableEnergy = PositKinEner
356 G4double Phot1Energy = epsil * Tota
357 G4double Phot2Energy = (1. - epsil)
358
359 // *** prepare calculation of polarization t
360 G4ThreeVector Phot1Direction(dirx, diry, dir
361
362 // get interaction frame
363 G4ThreeVector nInteractionFrame =
364 G4PolarizationHelper::GetFrame(PositDirect
365
366 // define proper in-plane and out-of-plane c
367 fBeamPolarization.InvRotateAz(nInteractionFr
368 fTargetPolarization.InvRotateAz(nInteraction
369
370 // calculate spin transfere matrix
371
372 fCrossSectionCalculator->Initialize(epsil, g
373 fTargetP
374
375 Phot1Direction.rotateUz(PositDirection);
376 // create G4DynamicParticle object for the p
377 G4DynamicParticle* aParticle1 =
378 new G4DynamicParticle(G4Gamma::Gamma(), Ph
379 fFinalGamma1Polarization = fCrossSectionCalc
380 G4double n1 = fFinalGamma1Polar
381 if(n1 > 1.)
382 {
383 G4ExceptionDescription ed;
384 ed << "ERROR: PolarizedAnnihilation Polari
385 << epsil << " is too large!!! \n"
386 << "annihi pol1= " << fFinalGamma1Polar
387 fFinalGamma1Polarization *= 1. / std::sqrt
388 G4Exception("G4PolarizedAnnihilationModel:
389 JustWarning, ed);
390 }
391
392 // define polarization of first final state
393 fFinalGamma1Polarization.SetPhoton();
394 fFinalGamma1Polarization.RotateAz(nInteracti
395 aParticle1->SetPolarization(fFinalGamma1Pola
396 fFinalGamma1Pola
397 fFinalGamma1Pola
398
399 fvect->push_back(aParticle1);
400
401 // *****************************************
402
403 G4double Eratio = Phot1Energy / Phot2Energy;
404 G4double PositP =
405 std::sqrt(PositKinEnergy * (PositKinEnergy
406 G4ThreeVector Phot2Direction(-dirx * Eratio,
407 (PositP - dirz
408 Phot2Direction.rotateUz(PositDirection);
409 // create G4DynamicParticle object for the p
410 G4DynamicParticle* aParticle2 =
411 new G4DynamicParticle(G4Gamma::Gamma(), Ph
412
413 // define polarization of second final state
414 fFinalGamma2Polarization = fCrossSectionCalc
415 G4double n2 = fFinalGamma2Polar
416 if(n2 > 1.)
417 {
418 G4ExceptionDescription ed;
419 ed << "ERROR: PolarizedAnnihilation Polari
420 << epsil << " is too large!!! \n";
421 ed << "annihi pol2= " << fFinalGamma2Polar
422
423 G4Exception("G4PolarizedAnnihilationModel:
424 JustWarning, ed);
425 fFinalGamma2Polarization *= 1. / std::sqrt
426 }
427 fFinalGamma2Polarization.SetPhoton();
428 fFinalGamma2Polarization.RotateAz(nInteracti
429 aParticle2->SetPolarization(fFinalGamma2Pola
430 fFinalGamma2Pola
431 fFinalGamma2Pola
432
433 fvect->push_back(aParticle2);
434 }
435