Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PolarizedMollerBhabhaModel.cc
Go to the documentation of this file.
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 // $Id: G4PolarizedMollerBhabhaModel.cc 96114 2016-03-16 18:51:33Z gcosmo $
27 // -------------------------------------------------------------------
28 //
29 // GEANT4 Class file
30 //
31 // File name: G4PolarizedMollerBhabhaModel
32 //
33 // Author: A.Schaelicke on base of Vladimir Ivanchenko code
34 //
35 // Creation date: 10.11.2005
36 //
37 // Modifications:
38 //
39 // 20-08-05, modified interface (A.Schaelicke)
40 //
41 // Class Description:
42 //
43 // Implementation of energy loss and delta-electron production by e+/e-
44 // (including polarization effects)
45 //
46 // -------------------------------------------------------------------
47 //
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50 
52 #include "G4PhysicalConstants.hh"
53 #include "G4Electron.hh"
54 #include "G4Positron.hh"
56 #include "Randomize.hh"
57 
58 #include "G4PolarizationManager.hh"
59 #include "G4PolarizationHelper.hh"
62 
63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64 
66  const G4String& nam)
67  : G4MollerBhabhaModel(p,nam)
68 {
69 
70  // G4cout<<" particle==electron "<<(p==theElectron)<<G4endl;
71  isElectron=(p==theElectron); // necessary due to wrong order in G4MollerBhabhaModel constructor!
72 
73  if (p==nullptr) {
74 
75  }
76  if (!isElectron) {
77  G4cout<<" buildBhabha cross section "<<isElectron<<G4endl;
78  crossSectionCalculator = new G4PolarizedBhabhaCrossSection();
79  } else {
80  G4cout<<" buildMoller cross section "<<isElectron<<G4endl;
81  crossSectionCalculator = new G4PolarizedMollerCrossSection();
82  }
83 }
84 
85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
86 
88 {
89  if (crossSectionCalculator) {
90  delete crossSectionCalculator;
91  }
92 }
93 
94 
95 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96 
98  const G4ParticleDefinition* pd,
99  G4double kinEnergy,
100  G4double cut,
101  G4double emax)
102 {
103  G4double xs =
105  cut,emax);
106 // G4cout<<"calc eIoni xsec "<<xs<<G4endl;
107 // G4cout<<" "<<kinEnergy<<" "<<cut<<" "<<emax<<G4endl;
108  G4double factor=1.;
109  if (xs!=0.) {
110  // G4cout<<"calc asym"<<G4endl;
111  G4double tmax = MaxSecondaryEnergy(pd, kinEnergy);
112  tmax = std::min(emax, tmax);
113 
114  if (std::fabs(cut/emax-1.)<1.e-10) return xs;
115 
116  if(cut < tmax) {
117 
118  G4double xmin = cut/kinEnergy;
119  G4double xmax = tmax/kinEnergy;
120 // G4cout<<"calc asym "<<xmin<<","<<xmax<<G4endl;
121  G4double gam = kinEnergy/electron_mass_c2 + 1.0;
122 
123  G4double crossPol=crossSectionCalculator->
124  TotalXSection(xmin,xmax,gam,
125  theBeamPolarization,
126  theTargetPolarization);
127  G4double crossUnpol=crossSectionCalculator->
128  TotalXSection(xmin,xmax,gam,
131  if (crossUnpol>0.) factor=crossPol/crossUnpol;
132  // G4cout<<" factor="<<factor<<G4endl;
133  }
134  }
135  return xs*factor;
136 }
137 
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139 
140 void G4PolarizedMollerBhabhaModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
141  const G4MaterialCutsCouple* ,
142  const G4DynamicParticle* dp,
143  G4double tmin,
144  G4double maxEnergy)
145 {
146  // *** obtain and save target and beam polarization ***
148 
149  const G4Track * aTrack = fParticleChange->GetCurrentTrack();
150 
151  // obtain polarization of the beam
152  theBeamPolarization = dp->GetPolarization();
153 
154  // obtain polarization of the media
155  G4VPhysicalVolume* aPVolume = aTrack->GetVolume();
156  G4LogicalVolume* aLVolume = aPVolume->GetLogicalVolume();
157  const G4bool targetIsPolarized = polarizationManger->IsPolarized(aLVolume);
158  theTargetPolarization = polarizationManger->GetVolumePolarization(aLVolume);
159 
160  // transfer target polarization in interaction frame
161  if (targetIsPolarized)
162  theTargetPolarization.rotateUz(dp->GetMomentumDirection());
163 
164 
165 
166 
167  G4double tmax = std::min(maxEnergy, MaxSecondaryKinEnergy(dp));
168  if(tmin >= tmax) return;
169  // if(tmin > tmax) tmin = tmax;
170 
171  G4double polL = theBeamPolarization.z()*theTargetPolarization.z();
172  polL=std::fabs(polL);
173  G4double polT = theBeamPolarization.x()*theTargetPolarization.x() +
174  theBeamPolarization.y()*theTargetPolarization.y();
175  polT=std::fabs(polT);
176 
177  G4double kineticEnergy = dp->GetKineticEnergy();
178  G4double energy = kineticEnergy + electron_mass_c2;
179  G4double totalMomentum = std::sqrt(kineticEnergy*(energy + electron_mass_c2));
180  G4double xmin = tmin/kineticEnergy;
181  G4double xmax = tmax/kineticEnergy;
182  G4double gam = energy/electron_mass_c2;
183  G4double gamma2 = gam*gam;
184  G4double gmo = gam - 1.;
185  G4double gmo2 = gmo*gmo;
186  G4double gmo3 = gmo2*gmo;
187  G4double gpo = gam + 1.;
188  G4double gpo2 = gpo*gpo;
189  G4double gpo3 = gpo2*gpo;
190  G4double x, y, q, grej, grej2;
191  G4double z = 0.;
192  G4double xs = 0., phi =0.;
193  G4ThreeVector direction = dp->GetMomentumDirection();
194 
195  //(Polarized) Moller (e-e-) scattering
196  if (isElectron) {
197  // *** dice according to polarized cross section
198  G4double G = ((2.0*gam - 1.0)/gamma2)*(1. - polT - polL*gam);
199  G4double H = (sqr(gam - 1.0)/gamma2)*(1. + polT + polL*((gam + 3.)/(gam - 1.)));
200 
201  y = 1.0 - xmax;
202  grej = 1.0 - G*xmax + xmax*xmax*(H + (1.0 - G*y)/(y*y));
203  grej2 = 1.0 - G*xmin + xmin*xmin*(H + (1.0 - G*y)/(y*y));
204  if (grej2 > grej) grej = grej2;
205  G4double prefM = gamma2*classic_electr_radius*classic_electr_radius/(gmo2*(gam + 1.0));
206  grej *= prefM;
207  do {
208  q = G4UniformRand();
209  x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
210  if (crossSectionCalculator) {
211  crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
212  theTargetPolarization,1);
213  xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
215  z=xs*sqr(x)*4.;
216  if (grej < z) {
217  G4cout<<"WARNING : error in Moller rejection routine! \n"
218  <<" z = "<<z<<" grej="<<grej<<"\n";
219  }
220  } else {
221  G4cout<<"No calculator in Moller scattering"<<G4endl;
222  }
223  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
224  } while(grej * G4UniformRand() > z);
225  //Bhabha (e+e-) scattering
226  } else {
227  // *** dice according to polarized cross section
228  y = xmax*xmax;
229  grej = 0.;
230  grej += y*y*gmo3*(1. + (polL + polT)*(gam + 3.)/gmo);
231  grej += -2.*xmin*xmin*xmin*gam*gmo2*(1. - (polL + polT)*(gam + 3.)/gmo);
232  grej += y*y*gmo*(3.*gamma2 + 6.*gam + 4.)*(1. + (polL*(3.*gam + 1.)*(gamma2 + gam + 1.) + polT*((gam + 2.)*gamma2 + 1.))/(gmo*(3.*gam*(gam + 2.) + 4.)));
233  grej /= gpo3;
234  grej += -xmin*(2.*gamma2 + 4.*gam + 1.)*(1. - gam*(polL*(2.*gam + 1.) + polT)/(2.*gam*(gam + 2.) + 1.))/gpo2;
235  grej += gamma2/(gamma2 - 1.);
237  grej *= prefB;
238 
239  do {
240  q = G4UniformRand();
241  x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
242  if (crossSectionCalculator) {
243  crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
244  theTargetPolarization,1);
245  xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
247  z=xs*sqr(x)*4.;
248  } else {
249  G4cout<<"No calculator in Bhabha scattering"<<G4endl;
250  }
251 
252  if(z > grej) {
253  G4cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<G4endl;
254  G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
255  << "Majorant " << grej << " < "
256  << z << " for x= " << x<<G4endl
257  << " e+e- (Bhabha) scattering"<<" at KinEnergy "<<kineticEnergy<<G4endl;
258  G4cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<G4endl;
259  }
260  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
261  } while(grej * G4UniformRand() > z);
262  }
263  //
264  //
265  // polar asymmetries (due to transverse polarizations)
266  //
267  //
268  if (crossSectionCalculator) {
269  // grej*=1./(sqr(x)*sqr(gamma2-1))*sqr(gam*(1+gam));
270  grej=xs*2.;
271  do {
272  phi = twopi * G4UniformRand() ;
273  crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
274  theTargetPolarization,1);
275  xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
277  if(xs > grej) {
278  if (isElectron){
279  G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
280  << "Majorant " << grej << " < "
281  << xs << " for phi= " << phi<<G4endl
282  << " e-e- (Moller) scattering"<< G4endl
283  <<"PHI DICING"<<G4endl;
284  } else {
285  G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
286  << "Majorant " << grej << " < "
287  << xs << " for phi= " << phi<<G4endl
288  << " e+e- (Bhabha) scattering"<< G4endl
289  <<"PHI DICING"<<G4endl;
290  }
291  }
292  // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
293  } while(grej * G4UniformRand() > xs);
294  }
295 
296  // fix kinematics of delta electron
297  G4double deltaKinEnergy = x * kineticEnergy;
298  G4double deltaMomentum =
299  std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
300  G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
301  (deltaMomentum * totalMomentum);
302  G4double sint = 1.0 - cost*cost;
303  if(sint > 0.0) sint = std::sqrt(sint);
304 
305 
306  G4ThreeVector deltaDirection(-sint*std::cos(phi),-sint*std::sin(phi), cost) ;
307  deltaDirection.rotateUz(direction);
308 
309  // primary change
310  kineticEnergy -= deltaKinEnergy;
312 
313  if(kineticEnergy > DBL_MIN) {
314  G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
315  direction = dir.unit();
317  }
318 
319  // create G4DynamicParticle object for delta ray
320  G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
321  vdp->push_back(delta);
322 
323  // get interaction frame
324  G4ThreeVector nInteractionFrame =
325  G4PolarizationHelper::GetFrame(direction,deltaDirection);
326 
327  if (crossSectionCalculator) {
328  // calculate mean final state polarizations
329 
330  theBeamPolarization.InvRotateAz(nInteractionFrame,direction);
331  theTargetPolarization.InvRotateAz(nInteractionFrame,direction);
332  crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
333  theTargetPolarization,2);
334 
335  // electron/positron
336  fPositronPolarization=crossSectionCalculator->GetPol2();
337  fPositronPolarization.RotateAz(nInteractionFrame,direction);
338 
339  fParticleChange->ProposePolarization(fPositronPolarization);
340 
341  // electron
342  fElectronPolarization=crossSectionCalculator->GetPol3();
343  fElectronPolarization.RotateAz(nInteractionFrame,deltaDirection);
344  delta->SetPolarization(fElectronPolarization.x(),
345  fElectronPolarization.y(),
346  fElectronPolarization.z());
347  }
348  else {
349  fPositronPolarization=G4ThreeVector();
350  fElectronPolarization=G4ThreeVector();
351  }
352 }
353 
354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void ProposePolarization(const G4ThreeVector &dir)
G4double MaxSecondaryKinEnergy(const G4DynamicParticle *dynParticle)
Definition: G4VEmModel.hh:484
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
double x() const
const char * p
Definition: xmltok.h:285
static G4PolarizationManager * GetInstance()
tuple x
Definition: test.py:50
virtual G4double XSection(const G4StokesVector &pol2, const G4StokesVector &pol3)=0
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kinEnergy, G4double cut, G4double emax) override
double z() const
static constexpr double twopi
Definition: G4SIunits.hh:76
virtual G4StokesVector GetPol2()
G4PolarizedMollerBhabhaModel(const G4ParticleDefinition *p=nullptr, const G4String &nam="PolarizedMollerBhabha")
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
G4ParticleChangeForLoss * fParticleChange
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:38
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static G4ThreeVector GetFrame(const G4ThreeVector &, const G4ThreeVector &)
float electron_mass_c2
Definition: hepunit.py:274
void SetPolarization(G4double polX, G4double polY, G4double polZ)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
static const G4double emax
void InvRotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
G4ParticleDefinition * theElectron
G4LogicalVolume * GetLogicalVolume() const
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double ComputeCrossSectionPerElectron(const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy)
Hep3Vector unit() const
double y() const
bool IsPolarized(G4LogicalVolume *lVol) const
const G4ThreeVector & GetPolarization() const
#define DBL_MIN
Definition: templates.hh:75
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4VPhysicalVolume * GetVolume() const
tuple z
Definition: test.py:28
virtual G4StokesVector GetPol3()
#define G4endl
Definition: G4ios.hh:61
const G4Track * GetCurrentTrack() const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
static const G4StokesVector ZERO
const G4ThreeVector & GetVolumePolarization(G4LogicalVolume *lVol) const
virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition *, G4double kinEnergy) final
void RotateAz(G4ThreeVector nInteractionFrame, G4ThreeVector particleDirection)
virtual void Initialize(G4double, G4double, G4double, const G4StokesVector &p0, const G4StokesVector &p1, G4int flag=0)