source: trunk/source/processes/electromagnetic/polarisation/src/G4PolarizedMollerBhabhaModel.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 12.9 KB
RevLine 
[819]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,v 1.4 2007/05/23 08:52:20 vnivanch Exp $
[1007]27// GEANT4 tag $Name: geant4-09-02 $
[819]28// -------------------------------------------------------------------
29//
30// GEANT4 Class file
31//
32// File name: G4PolarizedMollerBhabhaModel
33//
34// Author: A.Schaelicke on base of Vladimir Ivanchenko code
35//
36// Creation date: 10.11.2005
37//
38// Modifications:
39//
40// 20-08-05, modified interface (A.Schaelicke)
41//
42// Class Description:
43//
44// Implementation of energy loss and delta-electron production by e+/e-
45// (including polarization effects)
46//
47// -------------------------------------------------------------------
48//
49//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
50//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51
52#include "G4PolarizedMollerBhabhaModel.hh"
53#include "G4Electron.hh"
54#include "G4Positron.hh"
55#include "G4ParticleChangeForLoss.hh"
56#include "Randomize.hh"
57
58#include "G4PolarizationManager.hh"
59#include "G4PolarizationHelper.hh"
60#include "G4PolarizedBhabhaCrossSection.hh"
61#include "G4PolarizedMollerCrossSection.hh"
62
63//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
64
65G4PolarizedMollerBhabhaModel::G4PolarizedMollerBhabhaModel(const G4ParticleDefinition* p,
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==0) {
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
87G4PolarizedMollerBhabhaModel::~G4PolarizedMollerBhabhaModel()
88{
89 if (crossSectionCalculator) {
90 delete crossSectionCalculator;
91 }
92}
93
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
96
97G4double G4PolarizedMollerBhabhaModel::ComputeCrossSectionPerElectron(
98 const G4ParticleDefinition* pd,
99 G4double kinEnergy,
100 G4double cut,
101 G4double emax)
102{
103 G4double xs =
104 G4MollerBhabhaModel::ComputeCrossSectionPerElectron(pd,kinEnergy,
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,
129 G4StokesVector::ZERO,
130 G4StokesVector::ZERO);
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
140void 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 ***
147 G4PolarizationManager * polarizationManger = G4PolarizationManager::GetInstance();
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,
214 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 } while(grej * G4UniformRand() > z);
224 //Bhabha (e+e-) scattering
225 } else {
226 // *** dice according to polarized cross section
227 y = xmax*xmax;
228 grej = 0.;
229 grej += y*y*gmo3*(1. + (polL + polT)*(gam + 3.)/gmo);
230 grej += -2.*xmin*xmin*xmin*gam*gmo2*(1. - (polL + polT)*(gam + 3.)/gmo);
231 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.)));
232 grej /= gpo3;
233 grej += -xmin*(2.*gamma2 + 4.*gam + 1.)*(1. - gam*(polL*(2.*gam + 1.) + polT)/(2.*gam*(gam + 2.) + 1.))/gpo2;
234 grej += gamma2/(gamma2 - 1.);
235 G4double prefB = classic_electr_radius*classic_electr_radius/(gam - 1.0);
236 grej *= prefB;
237
238 do {
239 q = G4UniformRand();
240 x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
241 if (crossSectionCalculator) {
242 crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
243 theTargetPolarization,1);
244 xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
245 G4StokesVector::ZERO);
246 z=xs*sqr(x)*4.;
247 } else {
248 G4cout<<"No calculator in Bhabha scattering"<<G4endl;
249 }
250
251 if(z > grej) {
252 G4cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<G4endl;
253 G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
254 << "Majorant " << grej << " < "
255 << z << " for x= " << x<<G4endl
256 << " e+e- (Bhabha) scattering"<<" at KinEnergy "<<kineticEnergy<<G4endl;
257 G4cout<<"&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&"<<G4endl;
258 }
259 } while(grej * G4UniformRand() > z);
260 }
261 //
262 //
263 // polar asymmetries (due to transverse polarizations)
264 //
265 //
266 if (crossSectionCalculator) {
267 // grej*=1./(sqr(x)*sqr(gamma2-1))*sqr(gam*(1+gam));
268 grej=xs*2.;
269 do {
270 phi = twopi * G4UniformRand() ;
271 crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
272 theTargetPolarization,1);
273 xs=crossSectionCalculator->XSection(G4StokesVector::ZERO,
274 G4StokesVector::ZERO);
275 if(xs > grej) {
276 if (isElectron){
277 G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
278 << "Majorant " << grej << " < "
279 << xs << " for phi= " << phi<<G4endl
280 << " e-e- (Moller) scattering"<< G4endl
281 <<"PHI DICING"<<G4endl;
282 } else {
283 G4cout << "G4PolarizedMollerBhabhaModel::SampleSecondaries Warning! "<<G4endl
284 << "Majorant " << grej << " < "
285 << xs << " for phi= " << phi<<G4endl
286 << " e+e- (Bhabha) scattering"<< G4endl
287 <<"PHI DICING"<<G4endl;
288 }
289 }
290 } while(grej * G4UniformRand() > xs);
291 }
292
293 // fix kinematics of delta electron
294 G4double deltaKinEnergy = x * kineticEnergy;
295 G4double deltaMomentum =
296 std::sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
297 G4double cost = deltaKinEnergy * (energy + electron_mass_c2) /
298 (deltaMomentum * totalMomentum);
299 G4double sint = 1.0 - cost*cost;
300 if(sint > 0.0) sint = std::sqrt(sint);
301
302
303 G4ThreeVector deltaDirection(-sint*std::cos(phi),-sint*std::sin(phi), cost) ;
304 deltaDirection.rotateUz(direction);
305
306 // primary change
307 kineticEnergy -= deltaKinEnergy;
308 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
309
310 if(kineticEnergy > DBL_MIN) {
311 G4ThreeVector dir = totalMomentum*direction - deltaMomentum*deltaDirection;
312 direction = dir.unit();
313 fParticleChange->SetProposedMomentumDirection(direction);
314 }
315
316 // create G4DynamicParticle object for delta ray
317 G4DynamicParticle* delta = new G4DynamicParticle(theElectron,deltaDirection,deltaKinEnergy);
318 vdp->push_back(delta);
319
320 // get interaction frame
321 G4ThreeVector nInteractionFrame =
322 G4PolarizationHelper::GetFrame(direction,deltaDirection);
323
324 if (crossSectionCalculator) {
325 // calculate mean final state polarizations
326
327 theBeamPolarization.InvRotateAz(nInteractionFrame,direction);
328 theTargetPolarization.InvRotateAz(nInteractionFrame,direction);
329 crossSectionCalculator->Initialize(x,gam,phi,theBeamPolarization,
330 theTargetPolarization,2);
331
332 // electron/positron
333 fPositronPolarization=crossSectionCalculator->GetPol2();
334 fPositronPolarization.RotateAz(nInteractionFrame,direction);
335
336 fParticleChange->ProposePolarization(fPositronPolarization);
337
338 // electron
339 fElectronPolarization=crossSectionCalculator->GetPol3();
340 fElectronPolarization.RotateAz(nInteractionFrame,deltaDirection);
341 delta->SetPolarization(fElectronPolarization.x(),
342 fElectronPolarization.y(),
343 fElectronPolarization.z());
344 }
345 else {
346 fPositronPolarization=G4ThreeVector();
347 fElectronPolarization=G4ThreeVector();
348 }
349}
350
351//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.