ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedAnnihilationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PolarizedAnnihilationModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 //
32 // File name: G4PolarizedAnnihilationModel
33 //
34 // Author: Andreas Schaelicke
35 //
36 // Creation date: 01.05.2005
37 //
38 // Modifications:
39 // 18-07-06 use newly calculated cross sections (P. Starovoitov)
40 // 21-08-06 update interface (A. Schaelicke)
41 // 17-11-06 add protection agaist e+ zero energy PostStep (V.Ivanchenko)
42 // 10-07-07 copied Initialise() method from G4eeToTwoGammaModel to provide a
43 // local ParticleChangeForGamma object and reduce overhead
44 // in SampleSecondaries() (A. Schaelicke)
45 //
46 //
47 // Class Description:
48 //
49 // Implementation of polarized gamma Annihilation scattering on free electron
50 //
51 
52 // -------------------------------------------------------------------
53 #include "G4PolarizedAnnihilationModel.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4PolarizationManager.hh"
56 #include "G4PolarizationHelper.hh"
57 #include "G4StokesVector.hh"
58 #include "G4PolarizedAnnihilationCrossSection.hh"
59 #include "G4ParticleChangeForGamma.hh"
60 #include "G4TrackStatus.hh"
61 #include "G4Gamma.hh"
62 
63 G4PolarizedAnnihilationModel::G4PolarizedAnnihilationModel(const G4ParticleDefinition* p,
64  const G4String& nam)
65  : G4eeToTwoGammaModel(p,nam),
66  crossSectionCalculator(nullptr),
67  verboseLevel(0),
68  gParticleChange(nullptr)
69 {
70  crossSectionCalculator = new G4PolarizedAnnihilationCrossSection();
71 }
72 
73 G4PolarizedAnnihilationModel::~G4PolarizedAnnihilationModel()
74 {
75  delete crossSectionCalculator;
76 }
77 
78 
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80 
81 void G4PolarizedAnnihilationModel::Initialise(const G4ParticleDefinition* part,
82  const G4DataVector& dv)
83 {
84  G4eeToTwoGammaModel::Initialise(part, dv);
85  if(gParticleChange) { return; }
86  gParticleChange = GetParticleChangeForGamma();
87 }
88 
89 G4double
90 G4PolarizedAnnihilationModel::ComputeCrossSectionPerElectron(G4double kinEnergy)
91 {
92  // cross section from base model
93  G4double xs = G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(kinEnergy);
94 
95  G4double polzz = theBeamPolarization.z()*theTargetPolarization.z();
96  G4double poltt = theBeamPolarization.x()*theTargetPolarization.x()
97  + theBeamPolarization.y()*theTargetPolarization.y();
98  if (polzz!=0 || poltt!=0) {
99  G4double xval,lasym,tasym;
100  ComputeAsymmetriesPerElectron(kinEnergy,xval,lasym,tasym);
101  xs*=(1.+polzz*lasym+poltt*tasym);
102  }
103 
104  return xs;
105 }
106 
107 void G4PolarizedAnnihilationModel::ComputeAsymmetriesPerElectron(G4double ene,
108  G4double & valueX,
109  G4double & valueA,
110  G4double & valueT)
111 {
112  // *** calculate asymmetries
113  G4double gam = 1. + ene/electron_mass_c2;
114  G4double xs0=crossSectionCalculator->TotalXSection(0.,1.,gam,
115  G4StokesVector::ZERO,
116  G4StokesVector::ZERO);
117  G4double xsA=crossSectionCalculator->TotalXSection(0.,1.,gam,
118  G4StokesVector::P3,
119  G4StokesVector::P3);
120  G4double xsT1=crossSectionCalculator->TotalXSection(0.,1.,gam,
121  G4StokesVector::P1,
122  G4StokesVector::P1);
123  G4double xsT2=crossSectionCalculator->TotalXSection(0.,1.,gam,
124  G4StokesVector::P2,
125  G4StokesVector::P2);
126  G4double xsT=0.5*(xsT1+xsT2);
127 
128  valueX=xs0;
129  valueA=xsA/xs0-1.;
130  valueT=xsT/xs0-1.;
131  // G4cout<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
132  if ( (valueA < -1) || (1 < valueA)) {
133  G4cout<< " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n";
134  G4cout<< " something wrong in total cross section calculation (valueA)\n";
135  G4cout<< " LONG: "<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
136  }
137  if ( (valueT < -1) || (1 < valueT)) {
138  G4cout<< " ERROR PolarizedAnnihilationPS::ComputeAsymmetries \n";
139  G4cout<< " something wrong in total cross section calculation (valueT)\n";
140  G4cout<< " TRAN: "<<valueX<<"\t"<<valueA<<"\t"<<valueT<<" energy = "<<gam<<G4endl;
141  }
142 }
143 
144 void G4PolarizedAnnihilationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
145  const G4MaterialCutsCouple*,
146  const G4DynamicParticle* dp,
147  G4double, G4double)
148 {
149  const G4Track * aTrack = gParticleChange->GetCurrentTrack();
150 
151  // kill primary
152  gParticleChange->SetProposedKineticEnergy(0.);
153  gParticleChange->ProposeTrackStatus(fStopAndKill);
154 
155  // V.Ivanchenko add protection against zero kin energy
156  G4double PositKinEnergy = dp->GetKineticEnergy();
157 
158  if(PositKinEnergy == 0.0) {
159 
160  G4double cosTeta = 2.*G4UniformRand()-1.;
161  G4double sinTeta = std::sqrt((1.0 - cosTeta)*(1.0 + cosTeta));
162  G4double phi = twopi * G4UniformRand();
163  G4ThreeVector dir(sinTeta*std::cos(phi), sinTeta*std::sin(phi), cosTeta);
164  fvect->push_back( new G4DynamicParticle(G4Gamma::Gamma(), dir, electron_mass_c2));
165  fvect->push_back( new G4DynamicParticle(G4Gamma::Gamma(),-dir, electron_mass_c2));
166  return;
167  }
168 
169  // *** obtain and save target and beam polarization ***
170  G4PolarizationManager * polarizationManager = G4PolarizationManager::GetInstance();
171 
172  // obtain polarization of the beam
173  theBeamPolarization = aTrack->GetPolarization();
174 
175  // obtain polarization of the media
176  G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
177  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
178  const G4bool targetIsPolarized = polarizationManager->IsPolarized(aLVolume);
179  theTargetPolarization = polarizationManager->GetVolumePolarization(aLVolume);
180 
181  if (verboseLevel >= 1) {
182  G4cout << "G4PolarizedComptonModel::SampleSecondaries in "
183  << aLVolume->GetName() << G4endl;
184  }
185 
186  // transfer target electron polarization in frame of positron
187  if (targetIsPolarized)
188  theTargetPolarization.rotateUz(dp->GetMomentumDirection());
189 
190  G4ParticleMomentum PositDirection = dp->GetMomentumDirection();
191 
192  // polar asymmetry:
193  G4double polarization = theBeamPolarization.p3()*theTargetPolarization.p3();
194 
195  G4double gamam1 = PositKinEnergy/electron_mass_c2;
196  G4double gama = gamam1+1. , gamap1 = gamam1+2.;
197  G4double sqgrate = std::sqrt(gamam1/gamap1)/2. , sqg2m1 = std::sqrt(gamam1*gamap1);
198 
199  // limits of the energy sampling
200  G4double epsilmin = 0.5 - sqgrate , epsilmax = 0.5 + sqgrate;
201  G4double epsilqot = epsilmax/epsilmin;
202 
203  //
204  // sample the energy rate of the created gammas
205  // note: for polarized partices, the actual dicing strategy
206  // will depend on the energy, and the degree of polarization !!
207  //
208  G4double epsil;
209  G4double gmax=1. + std::fabs(polarization); // crude estimate
210 
211  //G4bool check_range=true;
212 
213  crossSectionCalculator->Initialize(epsilmin, gama, 0., theBeamPolarization, theTargetPolarization);
214  if (crossSectionCalculator->DiceEpsilon()<0) {
215  G4cout<<"ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
216  <<"epsilmin DiceRoutine not appropriate ! "<<crossSectionCalculator->DiceEpsilon()<<G4endl;
217  //check_range=false;
218  }
219 
220  crossSectionCalculator->Initialize(epsilmax, gama, 0., theBeamPolarization, theTargetPolarization);
221  if (crossSectionCalculator->DiceEpsilon()<0) {
222  G4cout<<"ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
223  <<"epsilmax DiceRoutine not appropriate ! "<<crossSectionCalculator->DiceEpsilon()<<G4endl;
224  //check_range=false;
225  }
226 
227  G4int ncount=0;
228  G4double trejectmax=0.;
229  G4double treject;
230 
231 
232  do {
233  //
234  epsil = epsilmin*std::pow(epsilqot,G4UniformRand());
235 
236  crossSectionCalculator->Initialize(epsil, gama, 0., theBeamPolarization, theTargetPolarization,1);
237 
238  treject = crossSectionCalculator->DiceEpsilon();
239  treject*=epsil;
240 
241  if (treject>gmax || treject<0.)
242  G4cout<<"ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
243  <<" eps ("<<epsil<<") rejection does not work properly: "<<treject<<G4endl;
244  ++ncount;
245  if (treject>trejectmax) trejectmax=treject;
246  if (ncount>1000) {
247  G4cout<<"WARNING in PolarizedAnnihilationPS::PostStepDoIt\n"
248  <<"eps dicing very inefficient ="<<trejectmax/gmax
249  <<", "<<treject/gmax<<". For secondary energy = "<<epsil<<" "<<ncount<<G4endl;
250  break;
251  }
252 
253  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
254  } while( treject < gmax*G4UniformRand() );
255 
256  //
257  // scattered Gamma angles. ( Z - axis along the parent positron)
258  //
259 
260  G4double cost = (epsil*gamap1-1.)/(epsil*sqg2m1);
261  G4double sint = std::sqrt((1.+cost)*(1.-cost));
262  G4double phi = 0.;
263  G4double beamTrans = std::sqrt(sqr(theBeamPolarization.p1()) + sqr(theBeamPolarization.p2()));
264  G4double targetTrans = std::sqrt(sqr(theTargetPolarization.p1()) + sqr(theTargetPolarization.p2()));
265 
266  // G4cout<<"phi dicing START"<<G4endl;
267  do{
268  phi = twopi * G4UniformRand();
269  crossSectionCalculator->Initialize(epsil, gama, 0., theBeamPolarization, theTargetPolarization,2);
270 
271  G4double gdiced =crossSectionCalculator->getVar(0);
272  gdiced += crossSectionCalculator->getVar(3)*theBeamPolarization.p3()*theTargetPolarization.p3();
273  gdiced += 1.*(std::fabs(crossSectionCalculator->getVar(1))
274  + std::fabs(crossSectionCalculator->getVar(2)))*beamTrans*targetTrans;
275  gdiced += 1.*std::fabs(crossSectionCalculator->getVar(4))
276  *(std::fabs(theBeamPolarization.p3())*targetTrans + std::fabs(theTargetPolarization.p3())*beamTrans);
277 
278  G4double gdist = crossSectionCalculator->getVar(0);
279  gdist += crossSectionCalculator->getVar(3)*theBeamPolarization.p3()*theTargetPolarization.p3();
280  gdist += crossSectionCalculator->getVar(1)*(std::cos(phi)*theBeamPolarization.p1()
281  + std::sin(phi)*theBeamPolarization.p2())
282  *(std::cos(phi)*theTargetPolarization.p1()
283  + std::sin(phi)*theTargetPolarization.p2());
284  gdist += crossSectionCalculator->getVar(2)*(std::cos(phi)*theBeamPolarization.p2()
285  - std::sin(phi)*theBeamPolarization.p1())
286  *(std::cos(phi)*theTargetPolarization.p2()
287  - std::sin(phi)*theTargetPolarization.p1());
288  gdist += crossSectionCalculator->getVar(4)
289  *(std::cos(phi)*theBeamPolarization.p3()*theTargetPolarization.p1()
290  + std::cos(phi)*theBeamPolarization.p1()*theTargetPolarization.p3()
291  + std::sin(phi)*theBeamPolarization.p3()*theTargetPolarization.p2()
292  + std::sin(phi)*theBeamPolarization.p2()*theTargetPolarization.p3());
293 
294  treject = gdist/gdiced;
295  //G4cout<<" treject = "<<treject<<" at phi = "<<phi<<G4endl;
296  if (treject>1.+1.e-10 || treject<0){
297  G4cout<<"!!!ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
298  <<" phi rejection does not work properly: "<<treject<<G4endl;
299  G4cout<<" gdiced = "<<gdiced<<G4endl;
300  G4cout<<" gdist = "<<gdist<<G4endl;
301  G4cout<<" epsil = "<<epsil<<G4endl;
302  }
303 
304  if (treject<1.e-3) {
305  G4cout<<"!!!ERROR in PolarizedAnnihilationPS::PostStepDoIt\n"
306  <<" phi rejection does not work properly: "<<treject<<"\n";
307  G4cout<<" gdiced="<<gdiced<<" gdist="<<gdist<<"\n";
308  G4cout<<" epsil = "<<epsil<<G4endl;
309  }
310 
311  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
312  } while( treject < G4UniformRand() );
313  // G4cout<<"phi dicing END"<<G4endl;
314 
315  G4double dirx = sint*std::cos(phi) , diry = sint*std::sin(phi) , dirz = cost;
316 
317  //
318  // kinematic of the created pair
319  //
320  G4double TotalAvailableEnergy = PositKinEnergy + 2*electron_mass_c2;
321  G4double Phot1Energy = epsil*TotalAvailableEnergy;
322  G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy;
323 
324  // *** prepare calculation of polarization transfer ***
325  G4ThreeVector Phot1Direction (dirx, diry, dirz);
326 
327  // get interaction frame
328  G4ThreeVector nInteractionFrame =
329  G4PolarizationHelper::GetFrame(PositDirection,Phot1Direction);
330 
331  // define proper in-plane and out-of-plane component of initial spins
332  theBeamPolarization.InvRotateAz(nInteractionFrame,PositDirection);
333  theTargetPolarization.InvRotateAz(nInteractionFrame,PositDirection);
334 
335  // calculate spin transfere matrix
336 
337  crossSectionCalculator->Initialize(epsil,gama,phi,theBeamPolarization,theTargetPolarization,2);
338 
339  // **********************************************************************
340 
341  Phot1Direction.rotateUz(PositDirection);
342  // create G4DynamicParticle object for the particle1
343  G4DynamicParticle* aParticle1= new G4DynamicParticle (G4Gamma::Gamma(),
344  Phot1Direction, Phot1Energy);
345  finalGamma1Polarization=crossSectionCalculator->GetPol2();
346  G4double n1=finalGamma1Polarization.mag2();
347  if (n1>1) {
348  G4cout<<"ERROR: PolarizedAnnihilation Polarization Vector at epsil = "
349  <<epsil<<" is too large!!! \n"
350  <<"annihi pol1= "<<finalGamma1Polarization<<", ("<<n1<<")\n";
351  finalGamma1Polarization*=1./std::sqrt(n1);
352  }
353 
354  // define polarization of first final state photon
355  finalGamma1Polarization.SetPhoton();
356  finalGamma1Polarization.RotateAz(nInteractionFrame,Phot1Direction);
357  aParticle1->SetPolarization(finalGamma1Polarization.p1(),
358  finalGamma1Polarization.p2(),
359  finalGamma1Polarization.p3());
360 
361  fvect->push_back(aParticle1);
362 
363 
364  // **********************************************************************
365 
366  G4double Eratio= Phot1Energy/Phot2Energy;
367  G4double PositP= std::sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2));
368  G4ThreeVector Phot2Direction (-dirx*Eratio, -diry*Eratio,
369  (PositP-dirz*Phot1Energy)/Phot2Energy);
370  Phot2Direction.rotateUz(PositDirection);
371  // create G4DynamicParticle object for the particle2
372  G4DynamicParticle* aParticle2= new G4DynamicParticle (G4Gamma::Gamma(),
373  Phot2Direction, Phot2Energy);
374 
375  // define polarization of second final state photon
376  finalGamma2Polarization=crossSectionCalculator->GetPol3();
377  G4double n2=finalGamma2Polarization.mag2();
378  if (n2>1) {
379  G4cout<<"ERROR: PolarizedAnnihilation Polarization Vector at epsil = "<<epsil<<" is too large!!! \n";
380  G4cout<<"annihi pol2= "<<finalGamma2Polarization<<", ("<<n2<<")\n";
381 
382  finalGamma2Polarization*=1./std::sqrt(n2);
383  }
384  finalGamma2Polarization.SetPhoton();
385  finalGamma2Polarization.RotateAz(nInteractionFrame,Phot2Direction);
386  aParticle2->SetPolarization(finalGamma2Polarization.p1(),
387  finalGamma2Polarization.p2(),
388  finalGamma2Polarization.p3());
389 
390  fvect->push_back(aParticle2);
391 }