Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4BoldyshevTripletModel.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$
27 // GEANT4 tag $Name: $
28 //
29 //
30 // Author: Gerardo Depaola & Francesco Longo
31 //
32 // History:
33 // --------
34 // 23-06-2010 First implementation as model
35 
36 
37 #include "G4BoldyshevTripletModel.hh"
38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh"
40 
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 
43 using namespace std;
44 
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 
47 G4BoldyshevTripletModel::G4BoldyshevTripletModel(const G4ParticleDefinition*,
48  const G4String& nam)
49  :G4VEmModel(nam),fParticleChange(0),smallEnergy(4.*MeV),isInitialised(false),
50  crossSectionHandler(0),meanFreePathTable(0)
51 {
52  lowEnergyLimit = 4.0*electron_mass_c2;
53  highEnergyLimit = 100 * GeV;
54  SetHighEnergyLimit(highEnergyLimit);
55 
56  verboseLevel= 0;
57  // Verbosity scale:
58  // 0 = nothing
59  // 1 = warning for energy non-conservation
60  // 2 = details of energy budget
61  // 3 = calculation of cross sections, file openings, sampling of atoms
62  // 4 = entering in methods
63 
64  if(verboseLevel > 0) {
65  G4cout << "Triplet Gamma conversion is constructed " << G4endl
66  << "Energy range: "
67  << lowEnergyLimit / MeV << " MeV - "
68  << highEnergyLimit / GeV << " GeV"
69  << G4endl;
70  }
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
75 G4BoldyshevTripletModel::~G4BoldyshevTripletModel()
76 {
77  if (crossSectionHandler) delete crossSectionHandler;
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
82 void
83 G4BoldyshevTripletModel::Initialise(const G4ParticleDefinition*,
84  const G4DataVector&)
85 {
86  if (verboseLevel > 3)
87  G4cout << "Calling G4BoldyshevTripletModel::Initialise()" << G4endl;
88 
89  if (crossSectionHandler)
90  {
91  crossSectionHandler->Clear();
92  delete crossSectionHandler;
93  }
94 
95  // Read data tables for all materials
96 
97  crossSectionHandler = new G4CrossSectionHandler();
98  crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
99  G4String crossSectionFile = "tripdata/pp-trip-cs-"; // here only pair in electron field cs should be used
100  crossSectionHandler->LoadData(crossSectionFile);
101 
102  //
103 
104  if (verboseLevel > 0) {
105  G4cout << "Loaded cross section files for Livermore GammaConversion" << G4endl;
106  G4cout << "To obtain the total cross section this should be used only " << G4endl
107  << "in connection with G4NuclearGammaConversion " << G4endl;
108  }
109 
110  if (verboseLevel > 0) {
111  G4cout << "Livermore Electron Gamma Conversion model is initialized " << G4endl
112  << "Energy range: "
113  << LowEnergyLimit() / MeV << " MeV - "
114  << HighEnergyLimit() / GeV << " GeV"
115  << G4endl;
116  }
117 
118  if(isInitialised) return;
119  fParticleChange = GetParticleChangeForGamma();
120  isInitialised = true;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
125 G4double
126 G4BoldyshevTripletModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
127  G4double GammaEnergy,
128  G4double Z, G4double,
129  G4double, G4double)
130 {
131  if (verboseLevel > 3) {
132  G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel"
133  << G4endl;
134  }
135  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
136 
137  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
138  return cs;
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 
143 void G4BoldyshevTripletModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
144  const G4MaterialCutsCouple* ,
145  const G4DynamicParticle* aDynamicGamma,
146  G4double,
147  G4double)
148 {
149 
150 // The energies of the secondary particles are sampled using
151 // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26)
152 
153  if (verboseLevel > 3)
154  G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel" << G4endl;
155 
156  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
157  G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
158 
159  G4double epsilon ;
160  G4double p0 = electron_mass_c2;
161 
162  G4double positronTotEnergy, electronTotEnergy, thetaEle, thetaPos;
163  G4double ener_re=0., theta_re, phi_re, phi;
164 
165  // Calculo de theta - elecron de recoil
166 
167  G4double energyThreshold = sqrt(2.)*electron_mass_c2; // -> momentumThreshold_N = 1
168  energyThreshold = 1.1*electron_mass_c2;
169  // G4cout << energyThreshold << G4endl;
170 
171  G4double momentumThreshold_c = sqrt(energyThreshold * energyThreshold - electron_mass_c2*electron_mass_c2); // momentun in MeV/c unit
172  G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2; // momentun in mc unit
173 
174  // Calculation of recoil electron production
175 
176  G4double SigmaTot = (28./9.) * std::log ( 2.* photonEnergy / electron_mass_c2 ) - 218. / 27. ;
177  G4double X_0 = 2. * ( sqrt(momentumThreshold_N*momentumThreshold_N + 1) -1 );
178  G4double SigmaQ = (82./27. - (14./9.) * log (X_0) + 4./15.*X_0 - 0.0348 * X_0 * X_0);
179  G4double recoilProb = G4UniformRand();
180  //G4cout << "SIGMA TOT " << SigmaTot << " " << "SigmaQ " << SigmaQ << " " << SigmaQ/SigmaTot << " " << recoilProb << G4endl;
181 
182  if (recoilProb >= SigmaQ/SigmaTot) // create electron recoil
183  {
184 
185  G4double cosThetaMax = ( ( energyThreshold - electron_mass_c2 ) / (momentumThreshold_c) + electron_mass_c2*
186  ( energyThreshold + electron_mass_c2 ) / (photonEnergy*momentumThreshold_c) );
187 
188  if (cosThetaMax > 1) G4cout << "ERRORE " << G4endl;
189 
190  G4double r1;
191  G4double r2;
192  G4double are, bre, loga, f1_re, greject, cost;
193 
194  do {
195  r1 = G4UniformRand();
196  r2 = G4UniformRand();
197  // cost = (pow(4./enern,0.5*r1)) ;
198  cost = pow(cosThetaMax,r1);
199  theta_re = acos(cost);
200  are = 1./(14.*cost*cost);
201  bre = (1.-5.*cost*cost)/(2.*cost);
202  loga = log((1.+ cost)/(1.- cost));
203  f1_re = 1. - bre*loga;
204 
205  if ( theta_re >= 4.47*CLHEP::pi/180.)
206  {
207  greject = are*f1_re;
208  } else {
209  greject = 1. ;
210  }
211  } while(greject < r2);
212 
213  // Calculo de phi - elecron de recoil
214 
215  G4double r3, r4, rt;
216 
217  do {
218 
219  r3 = G4UniformRand();
220  r4 = G4UniformRand();
221  phi_re = twopi*r3 ;
222  G4double sint2 = 1. - cost*cost ;
223  G4double fp = 1. - sint2*loga/(2.*cost) ;
224  rt = (1.-cos(2.*phi_re)*fp/f1_re)/(2.*pi) ;
225 
226  } while(rt < r4);
227 
228  // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
229 
230  G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
231  G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2
232  + (S - electron_mass_c2*electron_mass_c2)
233  *(S - electron_mass_c2*electron_mass_c2)*sin(theta_re)*sin(theta_re);
234  ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
235 
236  // G4cout << "electron de retroceso " << ener_re << " " << theta_re << " " << phi_re << G4endl;
237 
238  // Recoil electron creation
239  G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
240 
241  G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
242 
243  G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
244  electronRDirection.rotateUz(photonDirection);
245 
246  G4DynamicParticle* particle3 = new G4DynamicParticle (G4Electron::Electron(),
247  electronRDirection,
248  electronRKineEnergy);
249  fvect->push_back(particle3);
250 
251  }
252  else
253  {
254  // deposito la energia ener_re - electron_mass_c2
255  // G4cout << "electron de retroceso " << ener_re << G4endl;
256  fParticleChange->ProposeLocalEnergyDeposit(ener_re - electron_mass_c2);
257  }
258 
259  // Depaola (2004) suggested distribution for e+e- energy
260 
261  // G4double t = 0.5*asinh(momentumThreshold_N);
262  G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
263 
264  G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
265 
266  G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
267  G4double J2 = (-2./3.)*log(2.*sinh(t)) + t*cosh(t)/sinh(t) + (sinh(t)-t*pow(cosh(t),3))/(3.*pow(sinh(t),3));
268  G4double b = 2.*(J1-J2)/J1;
269 
270  G4double n = 1 - b/6.;
271  G4double re=0.;
272  re = G4UniformRand();
273  G4double a = 0.;
274  G4double b1 = 16. - 3.*b - 36.*b*re*n + 36.*b*pow(re,2.)*pow(n,2.) +
275  6.*pow(b,2.)*re*n;
276  a = pow((b1/b),0.5);
277  G4double c1 = (-6. + 12.*re*n + b + 2*a)*pow(b,2.);
278  epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
279 
280  G4double photonEnergy1 = photonEnergy - ener_re ; // resto al foton la energia del electron de retro.
281  positronTotEnergy = epsilon*photonEnergy1;
282  electronTotEnergy = photonEnergy1 - positronTotEnergy; // temporarly
283 
284  G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy -
285  electron_mass_c2*electron_mass_c2) ;
286  G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy -
287  electron_mass_c2*electron_mass_c2) ;
288 
289  thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
290  thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
291  phi = twopi * G4UniformRand();
292 
293  G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
294  G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
295 
296 
297  // Kinematics of the created pair:
298  // the electron and positron are assumed to have a symetric angular
299  // distribution with respect to the Z axis along the parent photon
300 
301  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
302 
303  // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
304 
305  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
306  electronDirection.rotateUz(photonDirection);
307 
308  G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
309  electronDirection,
310  electronKineEnergy);
311 
312  // The e+ is always created (even with kinetic energy = 0) for further annihilation
313  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
314 
315  // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
316 
317  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
318  positronDirection.rotateUz(photonDirection);
319 
320  // Create G4DynamicParticle object for the particle2
321  G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
322  positronDirection, positronKineEnergy);
323  // Fill output vector
324 
325 
326  fvect->push_back(particle1);
327  fvect->push_back(particle2);
328 
329 
330 
331 
332  // kill incident photon
333  fParticleChange->SetProposedKineticEnergy(0.);
334  fParticleChange->ProposeTrackStatus(fStopAndKill);
335 
336 }
337 
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339 
340 
341 
342 
Generated on Sat May 25 2013 14:33:33 for Geant4 by   doxygen 1.8.4