Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4BoldyshevTripletModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 // Author: Sebastien Incerti
 27 //         22 January 2012
 28 //         on base of G4BoldyshevTripletModel (original version)
 29 //         and G4LivermoreRayleighModel (MT version)
 30 
 31 #include "G4BoldyshevTripletModel.hh"
 32 #include "G4PhysicalConstants.hh"
 33 #include "G4SystemOfUnits.hh"
 34 #include "G4Log.hh"
 35 #include "G4Exp.hh"
 36 
 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 38 
 39 using namespace std;
 40 
 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 42 
 43 const G4int G4BoldyshevTripletModel::maxZ;
 44 G4PhysicsFreeVector* G4BoldyshevTripletModel::data[] = {nullptr};
 45 
 46 G4BoldyshevTripletModel::G4BoldyshevTripletModel(const G4ParticleDefinition*, const G4String& nam)
 47   :G4VEmModel(nam),smallEnergy(4.*MeV)
 48 {
 49   fParticleChange = nullptr;
 50   
 51   lowEnergyLimit = 4.0*electron_mass_c2;
 52   momentumThreshold_c = energyThreshold = xb = xn = lowEnergyLimit; 
 53   
 54   verboseLevel= 0;
 55   // Verbosity scale for debugging purposes:
 56   // 0 = nothing 
 57   // 1 = calculation of cross sections, file openings...
 58   // 2 = entering in methods
 59 
 60   if(verboseLevel > 0) 
 61   {
 62     G4cout << "G4BoldyshevTripletModel is constructed " << G4endl;
 63   }
 64 }
 65 
 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 67 
 68 G4BoldyshevTripletModel::~G4BoldyshevTripletModel()
 69 {
 70   if(IsMaster()) {
 71     for(G4int i=0; i<maxZ; ++i) {
 72       if(data[i]) { 
 73   delete data[i];
 74   data[i] = nullptr;
 75       }
 76     }
 77   }
 78 }
 79 
 80 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 81 
 82 void G4BoldyshevTripletModel::Initialise(const G4ParticleDefinition*,
 83                  const G4DataVector&)
 84 {
 85   if (verboseLevel > 1) 
 86   {
 87     G4cout << "Calling Initialise() of G4BoldyshevTripletModel." 
 88      << G4endl
 89      << "Energy range: "
 90      << LowEnergyLimit() / MeV << " MeV - "
 91      << HighEnergyLimit() / GeV << " GeV isMaster: " << IsMaster()
 92      << G4endl;
 93   }
 94   // compute values only once
 95   energyThreshold = 1.1*electron_mass_c2; 
 96   momentumThreshold_c = std::sqrt(energyThreshold * energyThreshold 
 97                                 - electron_mass_c2*electron_mass_c2); 
 98   G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2; 
 99   G4double t = 0.5*G4Log(momentumThreshold_N + 
100        std::sqrt(momentumThreshold_N*momentumThreshold_N + 1.0));
101   //G4cout << 0.5*asinh(momentumThreshold_N) << "  " << t << G4endl;
102   G4double sinht = std::sinh(t);                         
103   G4double cosht = std::cosh(t);  
104   G4double logsinht = G4Log(2.*sinht);                     
105   G4double J1 = 0.5*(t*cosht/sinht - logsinht);
106   G4double J2 = (-2./3.)*logsinht + t*cosht/sinht 
107     + (sinht - t*cosht*cosht*cosht)/(3.*sinht*sinht*sinht);
108 
109   xb = 2.*(J1-J2)/J1; 
110   xn = 1. - xb/6.;
111 
112   if(IsMaster()) 
113   {
114     // Access to elements  
115     const char* path = G4FindDataDir("G4LEDATA");
116 
117     G4ProductionCutsTable* theCoupleTable =
118       G4ProductionCutsTable::GetProductionCutsTable();
119   
120     G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
121   
122     for(G4int i=0; i<numOfCouples; ++i) 
123     {
124       const G4Material* material = 
125         theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
126       const G4ElementVector* theElementVector = material->GetElementVector();
127       std::size_t nelm = material->GetNumberOfElements();
128     
129       for (std::size_t j=0; j<nelm; ++j) 
130       {
131         G4int Z = std::min((*theElementVector)[j]->GetZasInt(), maxZ);
132         if(!data[Z]) { ReadData(Z, path); }
133       }
134     }
135   }
136   if(!fParticleChange) {
137     fParticleChange = GetParticleChangeForGamma();
138   }
139 }
140 
141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
142 
143 G4double 
144 G4BoldyshevTripletModel::MinPrimaryEnergy(const G4Material*,
145             const G4ParticleDefinition*,
146             G4double)
147 {
148   return lowEnergyLimit;
149 }
150 
151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
152 
153 void G4BoldyshevTripletModel::ReadData(size_t Z, const char* path)
154 {
155   if (verboseLevel > 1) 
156   {
157     G4cout << "Calling ReadData() of G4BoldyshevTripletModel" 
158      << G4endl;
159   }
160 
161   if(data[Z]) { return; }
162   
163   const char* datadir = path;
164 
165   if(!datadir) 
166   {
167     datadir = G4FindDataDir("G4LEDATA");
168     if(!datadir) 
169     {
170       G4Exception("G4BoldyshevTripletModel::ReadData()",
171       "em0006",FatalException,
172       "Environment variable G4LEDATA not defined");
173       return;
174     }
175   }
176   
177   data[Z] = new G4PhysicsFreeVector(0,/*spline=*/true);
178   std::ostringstream ost;
179   ost << datadir << "/livermore/tripdata/pp-trip-cs-" << Z <<".dat";
180   std::ifstream fin(ost.str().c_str());
181   
182   if( !fin.is_open()) 
183   {
184     G4ExceptionDescription ed;
185     ed << "G4BoldyshevTripletModel data file <" << ost.str().c_str()
186        << "> is not opened!" << G4endl;
187     G4Exception("G4BoldyshevTripletModel::ReadData()",
188     "em0003",FatalException,
189     ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
190     return;
191   } 
192   
193   else 
194   {
195     
196     if(verboseLevel > 3) { G4cout << "File " << ost.str() 
197        << " is opened by G4BoldyshevTripletModel" << G4endl;}
198     
199     data[Z]->Retrieve(fin, true);
200   } 
201 
202   // Activation of spline interpolation
203   data[Z]->FillSecondDerivatives();
204 }
205 
206 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
207 
208 G4double G4BoldyshevTripletModel::ComputeCrossSectionPerAtom(
209          const G4ParticleDefinition* part,
210    G4double GammaEnergy, G4double Z, G4double, G4double, G4double)
211 {
212   if (verboseLevel > 1) 
213   {
214     G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel" 
215      << G4endl;
216   }
217 
218   if (GammaEnergy < lowEnergyLimit) { return 0.0; } 
219 
220   G4double xs = 0.0;  
221   G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ));
222   G4PhysicsFreeVector* pv = data[intZ];
223 
224   // if element was not initialised
225   // do initialisation safely for MT mode
226   if(!pv) 
227   {
228     InitialiseForElement(part, intZ);
229     pv = data[intZ];
230     if(!pv) { return xs; }
231   }
232   // x-section is taken from the table
233   xs = pv->Value(GammaEnergy); 
234 
235   if(verboseLevel > 1)
236   {
237     G4cout  <<  "*** Triplet conversion xs for Z=" << Z << " at energy E(MeV)=" 
238       << GammaEnergy/MeV <<  "  cs=" << xs/millibarn << " mb" << G4endl;
239   }
240   return xs;
241 }
242 
243 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
244 
245 void G4BoldyshevTripletModel::SampleSecondaries(
246                                  std::vector<G4DynamicParticle*>* fvect,
247          const G4MaterialCutsCouple* /*couple*/,
248          const G4DynamicParticle* aDynamicGamma,
249          G4double, G4double)
250 {
251 
252   // The energies of the secondary particles are sampled using
253   // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26)
254   if (verboseLevel > 1) {
255     G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel" 
256      << G4endl;
257   }
258 
259   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
260   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
261 
262   G4double epsilon;
263 
264   CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
265 
266   // recoil electron thould be 3d particle
267   G4DynamicParticle* particle3 = nullptr;
268   static const G4double costlim = std::cos(4.47*CLHEP::pi/180.);
269   
270   G4double loga, f1_re, greject, cost;
271   G4double cosThetaMax = (energyThreshold - electron_mass_c2 
272      + electron_mass_c2*(energyThreshold + electron_mass_c2)/photonEnergy )
273   /momentumThreshold_c;
274   if (cosThetaMax > 1.) {
275     //G4cout << "G4BoldyshevTripletModel::SampleSecondaries: ERROR cosThetaMax= " 
276     //     << cosThetaMax << G4endl;
277     cosThetaMax = 1.0;
278   }
279 
280   G4double logcostm = G4Log(cosThetaMax);
281   do {
282     cost = G4Exp(logcostm*rndmEngine->flat());
283     G4double are = 1./(14.*cost*cost);
284     G4double bre = (1.-5.*cost*cost)/(2.*cost);
285     loga = G4Log((1.+ cost)/(1.- cost));
286     f1_re = 1. - bre*loga;
287     greject = (cost < costlim) ? are*f1_re : 1.0;
288   } while(greject < rndmEngine->flat());
289       
290   // Calculo de phi - elecron de recoil
291   G4double sint2 = (1. - cost)*(1. + cost);
292   G4double fp = 1. - sint2*loga/(2.*cost) ;
293   G4double rt, phi_re;
294   do {
295     phi_re = twopi*rndmEngine->flat();
296     rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi;
297   } while(rt < rndmEngine->flat());
298 
299   // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
300   G4double S  = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
301   G4double P2 = S - electron_mass_c2*electron_mass_c2;
302 
303   G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2;
304   G4double ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
305       
306   if(ener_re >= energyThreshold) 
307     {
308       G4double electronRKineEnergy = ener_re - electron_mass_c2;
309       G4double sint = std::sqrt(sint2);
310       G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost);
311       electronRDirection.rotateUz(photonDirection);
312       particle3 = new G4DynamicParticle (G4Electron::Electron(),
313            electronRDirection,
314            electronRKineEnergy);
315     }
316   else
317     {
318       // deposito la energia  ener_re - electron_mass_c2
319       // G4cout << "electron de retroceso " << ener_re << G4endl;
320       fParticleChange->ProposeLocalEnergyDeposit(std::max(0.0, ener_re - electron_mass_c2));
321       ener_re = 0.0;
322     }
323   
324   // Depaola (2004) suggested distribution for e+e- energy
325   // VI: very suspect that 1 random number is not enough
326   //     and sampling below is not correct - should be fixed
327   G4double re = rndmEngine->flat();
328   
329   G4double a  = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn);
330   G4double c1 = G4Exp(G4Log((-6. + 12.*re*xn + xb + 2*a)*xb*xb)/3.);
331   epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5;
332   
333   G4double photonEnergy1 = photonEnergy - ener_re ; 
334   // resto al foton la energia del electron de retro.
335   G4double positronTotEnergy = std::max(epsilon*photonEnergy1, electron_mass_c2);
336   G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2);
337   
338   static const G4double a1 = 1.6;
339   static const G4double a2 = 0.5333333333;
340   G4double uu = -G4Log(rndmEngine->flat()*rndmEngine->flat());
341   G4double u = (0.25 > rndmEngine->flat()) ? uu*a1 : uu*a2;
342 
343   G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
344   G4double sinte = std::sin(thetaEle);
345   G4double coste = std::cos(thetaEle);
346 
347   G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
348   G4double sintp = std::sin(thetaPos);
349   G4double costp = std::cos(thetaPos);
350 
351   G4double phi  = twopi * rndmEngine->flat();
352   G4double sinp = std::sin(phi);
353   G4double cosp = std::cos(phi);
354 
355   // Kinematics of the created pair:
356   // the electron and positron are assumed to have a symetric angular
357   // distribution with respect to the Z axis along the parent photon
358 
359   G4double electronKineEnergy = electronTotEnergy - electron_mass_c2;
360 
361   G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
362   electronDirection.rotateUz(photonDirection);
363 
364   G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
365                                                         electronDirection,
366               electronKineEnergy);
367 
368   G4double positronKineEnergy = positronTotEnergy - electron_mass_c2;
369 
370   G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
371   positronDirection.rotateUz(photonDirection);
372 
373   // Create G4DynamicParticle object for the particle2      
374   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
375                                                        positronDirection, positronKineEnergy);
376   // Fill output vector       
377   
378   fvect->push_back(particle1);
379   fvect->push_back(particle2);
380 
381   if(particle3) { fvect->push_back(particle3); }
382   
383   // kill incident photon
384   fParticleChange->SetProposedKineticEnergy(0.);
385   fParticleChange->ProposeTrackStatus(fStopAndKill);
386 }
387 
388 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
389 
390 #include "G4AutoLock.hh"
391 namespace { G4Mutex BoldyshevTripletModelMutex = G4MUTEX_INITIALIZER; }
392 
393 void G4BoldyshevTripletModel::InitialiseForElement(
394      const G4ParticleDefinition*, G4int Z)
395 {
396   G4AutoLock l(&BoldyshevTripletModelMutex);
397   // G4cout << "G4BoldyshevTripletModel::InitialiseForElement Z= " 
398   //    << Z << G4endl;
399   if(!data[Z]) { ReadData(Z); }
400   l.unlock();
401 }
402 
403 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
404