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