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