Geant4  10.00.p01
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: G4BoldyshevTripletModel.cc 74822 2013-10-22 14:42:13Z gcosmo $
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 
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 
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;
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 
76 {
78 }
79 
80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
81 
82 void
84  const G4DataVector&)
85 {
86  if (verboseLevel > 3)
87  G4cout << "Calling G4BoldyshevTripletModel::Initialise()" << G4endl;
88 
90  {
92  delete crossSectionHandler;
93  }
94 
95  // Read data tables for all materials
96 
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;
120  isInitialised = true;
121 }
122 
123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
124 
125 G4double
127  G4double GammaEnergy,
128  G4double Z, 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  // New Recoil energy calculation
237 
238  G4double momentum_recoil = 2* (electron_mass_c2) * (std::cos(theta_re)/(std::sin(phi_re)*std::sin(phi_re)));
239  G4double ener_recoil = sqrt( momentum_recoil*momentum_recoil + electron_mass_c2*electron_mass_c2);
240  ener_re = ener_recoil;
241 
242  // G4cout << "electron de retroceso " << ener_re << " " << theta_re << " " << phi_re << G4endl;
243 
244  // Recoil electron creation
245  G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
246 
247  G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
248 
249  G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
250  electronRDirection.rotateUz(photonDirection);
251 
253  electronRDirection,
254  electronRKineEnergy);
255  fvect->push_back(particle3);
256 
257  }
258  else
259  {
260  // deposito la energia ener_re - electron_mass_c2
261  // G4cout << "electron de retroceso " << ener_re << G4endl;
262  fParticleChange->ProposeLocalEnergyDeposit(ener_re - electron_mass_c2);
263  }
264 
265  // Depaola (2004) suggested distribution for e+e- energy
266 
267  // G4double t = 0.5*asinh(momentumThreshold_N);
268  G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
269 
270  G4cout << 0.5*asinh(momentumThreshold_N) << " " << t << G4endl;
271 
272  G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
273  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));
274  G4double b = 2.*(J1-J2)/J1;
275 
276  G4double n = 1 - b/6.;
277  G4double re=0.;
278  re = G4UniformRand();
279  G4double a = 0.;
280  G4double b1 = 16. - 3.*b - 36.*b*re*n + 36.*b*pow(re,2.)*pow(n,2.) +
281  6.*pow(b,2.)*re*n;
282  a = pow((b1/b),0.5);
283  G4double c1 = (-6. + 12.*re*n + b + 2*a)*pow(b,2.);
284  epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
285 
286  G4double photonEnergy1 = photonEnergy - ener_re ; // resto al foton la energia del electron de retro.
287  positronTotEnergy = epsilon*photonEnergy1;
288  electronTotEnergy = photonEnergy1 - positronTotEnergy; // temporarly
289 
290  G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy -
291  electron_mass_c2*electron_mass_c2) ;
292  G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy -
293  electron_mass_c2*electron_mass_c2) ;
294 
295  thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
296  thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
297  phi = twopi * G4UniformRand();
298 
299  G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
300  G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
301 
302 
303  // Kinematics of the created pair:
304  // the electron and positron are assumed to have a symetric angular
305  // distribution with respect to the Z axis along the parent photon
306 
307  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
308 
309  // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
310 
311  G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
312  electronDirection.rotateUz(photonDirection);
313 
315  electronDirection,
316  electronKineEnergy);
317 
318  // The e+ is always created (even with kinetic energy = 0) for further annihilation
319  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
320 
321  // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
322 
323  G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
324  positronDirection.rotateUz(photonDirection);
325 
326  // Create G4DynamicParticle object for the particle2
328  positronDirection, positronKineEnergy);
329  // Fill output vector
330 
331 
332  fvect->push_back(particle1);
333  fvect->push_back(particle2);
334 
335 
336 
337 
338  // kill incident photon
341 
342 }
343 
344 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
345 
346 
347 
348 
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static const double MeV
Definition: G4SIunits.hh:193
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
G4BoldyshevTripletModel(const G4ParticleDefinition *p=0, const G4String &nam="BoldyshevTriplet")
const G4double pi
G4double asinh(G4double value)
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
G4ParticleChangeForGamma * fParticleChange
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:683
G4double FindValue(G4int Z, G4double e) const
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
static const double GeV
Definition: G4SIunits.hh:196
const G4int n
static const G4double c1
void Initialise(G4VDataSetAlgorithm *interpolation=0, G4double minE=250 *CLHEP::eV, G4double maxE=100 *CLHEP::GeV, G4int numberOfBins=200, G4double unitE=CLHEP::MeV, G4double unitData=CLHEP::barn, G4int minZ=1, G4int maxZ=99)
static G4Positron * Positron()
Definition: G4Positron.cc:94
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void LoadData(const G4String &dataFile)
G4VCrossSectionHandler * crossSectionHandler
static const G4double b1
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ThreeVector G4ParticleMomentum
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121