source: trunk/source/processes/electromagnetic/lowenergy/src/G4BoldyshevTripletModel.cc @ 1350

Last change on this file since 1350 was 1350, checked in by garnier, 14 years ago

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