Geant4: G4BoldyshevTripletModel Class Reference

G4BoldyshevTripletModel::G4BoldyshevTripletModel	(	const G4ParticleDefinition *	p = `0`,
		const G4String &	nam = `"BoldyshevTriplet"`
	)

Definition at line 47 of file G4BoldyshevTripletModel.cc.

References G4cout, G4endl, and G4VEmModel::SetHighEnergyLimit().

00049   :G4VEmModel(nam),fParticleChange(0),smallEnergy(4.*MeV),isInitialised(false),
00050    crossSectionHandler(0),meanFreePathTable(0)
00051 {
00052   lowEnergyLimit = 4.0*electron_mass_c2;
00053   highEnergyLimit = 100 * GeV;
00054   SetHighEnergyLimit(highEnergyLimit);
00055          
00056   verboseLevel= 0;
00057   // Verbosity scale:
00058   // 0 = nothing 
00059   // 1 = warning for energy non-conservation 
00060   // 2 = details of energy budget
00061   // 3 = calculation of cross sections, file openings, sampling of atoms
00062   // 4 = entering in methods
00063 
00064   if(verboseLevel > 0) {
00065     G4cout << "Triplet Gamma conversion is constructed " << G4endl
00066            << "Energy range: "
00067            << lowEnergyLimit / MeV << " MeV - "
00068            << highEnergyLimit / GeV << " GeV"
00069            << G4endl;
00070   }
00071 }

G4BoldyshevTripletModel::~G4BoldyshevTripletModel ( ) [virtual]

Definition at line 75 of file G4BoldyshevTripletModel.cc.

00076 {  
00077   if (crossSectionHandler) delete crossSectionHandler;
00078 }

G4double G4BoldyshevTripletModel::ComputeCrossSectionPerAtom	(	const G4ParticleDefinition *	,
		G4double	kinEnergy,
		G4double	Z,
		G4double	A = `0`,
		G4double	cut = `0`,
		G4double	emax = `DBL_MAX`
	)			`[virtual]`

Reimplemented from G4VEmModel.

Definition at line 126 of file G4BoldyshevTripletModel.cc.

References G4VCrossSectionHandler::FindValue(), G4cout, and G4endl.

00130 {
00131   if (verboseLevel > 3) {
00132     G4cout << "Calling ComputeCrossSectionPerAtom() of G4BoldyshevTripletModel" 
00133            << G4endl;
00134   }
00135   if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0;
00136 
00137   G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
00138   return cs;
00139 }

G4double G4BoldyshevTripletModel::GetMeanFreePath	(	const G4Track &	aTrack,
		G4double	previousStepSize,
		G4ForceCondition *	condition
	)			`[protected]`

void G4BoldyshevTripletModel::Initialise	(	const G4ParticleDefinition *	,
		const G4DataVector &
	)			`[virtual]`

Implements G4VEmModel.

Definition at line 83 of file G4BoldyshevTripletModel.cc.

References G4VCrossSectionHandler::Clear(), fParticleChange, G4cout, G4endl, G4VEmModel::GetParticleChangeForGamma(), G4VEmModel::HighEnergyLimit(), G4VCrossSectionHandler::Initialise(), G4VCrossSectionHandler::LoadData(), and G4VEmModel::LowEnergyLimit().

00085 {
00086   if (verboseLevel > 3)
00087     G4cout << "Calling G4BoldyshevTripletModel::Initialise()" << G4endl;
00088 
00089   if (crossSectionHandler)
00090   {
00091     crossSectionHandler->Clear();
00092     delete crossSectionHandler;
00093   }
00094 
00095   // Read data tables for all materials
00096   
00097   crossSectionHandler = new G4CrossSectionHandler();
00098   crossSectionHandler->Initialise(0,lowEnergyLimit,100.*GeV,400);
00099   G4String crossSectionFile = "tripdata/pp-trip-cs-"; // here only pair in electron field cs should be used
00100   crossSectionHandler->LoadData(crossSectionFile);
00101 
00102   //
00103   
00104   if (verboseLevel > 0) {
00105     G4cout << "Loaded cross section files for Livermore GammaConversion" << G4endl;
00106     G4cout << "To obtain the total cross section this should be used only " << G4endl 
00107            << "in connection with G4NuclearGammaConversion " << G4endl;
00108   }
00109 
00110   if (verboseLevel > 0) { 
00111     G4cout << "Livermore Electron Gamma Conversion model is initialized " << G4endl
00112            << "Energy range: "
00113            << LowEnergyLimit() / MeV << " MeV - "
00114            << HighEnergyLimit() / GeV << " GeV"
00115            << G4endl;
00116   }
00117 
00118   if(isInitialised) return;
00119   fParticleChange = GetParticleChangeForGamma();
00120   isInitialised = true;
00121 }

void G4BoldyshevTripletModel::SampleSecondaries	(	std::vector< G4DynamicParticle * > *	,
		const G4MaterialCutsCouple *	,
		const G4DynamicParticle *	,
		G4double	tmin,
		G4double	maxEnergy
	)			`[virtual]`

Implements G4VEmModel.

Definition at line 143 of file G4BoldyshevTripletModel.cc.

References G4Electron::Electron(), fParticleChange, fStopAndKill, G4cout, G4endl, G4UniformRand, G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), CLHEP::detail::n, G4INCL::Math::pi, G4Positron::Positron(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4VParticleChange::ProposeTrackStatus(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().

00148 {
00149 
00150 // The energies of the secondary particles are sampled using
00151 // a modified Wheeler-Lamb model (see PhysRevD 7 (1973), 26) 
00152 
00153   if (verboseLevel > 3)
00154     G4cout << "Calling SampleSecondaries() of G4BoldyshevTripletModel" << G4endl;
00155 
00156   G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00157   G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00158 
00159   G4double epsilon ;
00160   G4double p0 = electron_mass_c2; 
00161   
00162   G4double positronTotEnergy, electronTotEnergy, thetaEle, thetaPos;
00163   G4double ener_re=0., theta_re, phi_re, phi;
00164 
00165   // Calculo de theta - elecron de recoil
00166   
00167   G4double energyThreshold = sqrt(2.)*electron_mass_c2; // -> momentumThreshold_N = 1
00168   energyThreshold = 1.1*electron_mass_c2;
00169   // G4cout << energyThreshold << G4endl; 
00170  
00171   G4double momentumThreshold_c = sqrt(energyThreshold * energyThreshold - electron_mass_c2*electron_mass_c2); // momentun in MeV/c unit
00172   G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2; // momentun in mc unit
00173   
00174   // Calculation of recoil electron production 
00175   
00176   G4double SigmaTot = (28./9.) * std::log ( 2.* photonEnergy / electron_mass_c2 ) - 218. / 27. ; 
00177   G4double X_0 = 2. * ( sqrt(momentumThreshold_N*momentumThreshold_N + 1) -1 );
00178   G4double SigmaQ = (82./27. - (14./9.) * log (X_0) + 4./15.*X_0 - 0.0348 * X_0 * X_0); 
00179   G4double recoilProb = G4UniformRand();
00180   //G4cout << "SIGMA TOT " << SigmaTot <<  " " << "SigmaQ " << SigmaQ << " " << SigmaQ/SigmaTot << " " << recoilProb << G4endl;
00181   
00182   if (recoilProb >= SigmaQ/SigmaTot) // create electron recoil 
00183     { 
00184       
00185       G4double cosThetaMax = (  ( energyThreshold - electron_mass_c2 ) / (momentumThreshold_c) + electron_mass_c2*
00186                                 ( energyThreshold + electron_mass_c2 ) / (photonEnergy*momentumThreshold_c) );
00187       
00188       if (cosThetaMax > 1) G4cout << "ERRORE " << G4endl;
00189       
00190       G4double r1;
00191       G4double r2;
00192       G4double are, bre, loga, f1_re, greject, cost;
00193       
00194       do {
00195         r1 = G4UniformRand();
00196         r2 = G4UniformRand();
00197         //      cost = (pow(4./enern,0.5*r1)) ;
00198         cost = pow(cosThetaMax,r1);
00199         theta_re = acos(cost);
00200         are = 1./(14.*cost*cost);
00201         bre = (1.-5.*cost*cost)/(2.*cost);
00202         loga = log((1.+ cost)/(1.- cost));
00203         f1_re = 1. - bre*loga;
00204         
00205         if ( theta_re >= 4.47*CLHEP::pi/180.)
00206           {
00207             greject = are*f1_re;
00208           } else {
00209           greject = 1. ;
00210         }
00211       } while(greject < r2);
00212       
00213       // Calculo de phi - elecron de recoil
00214       
00215       G4double r3, r4, rt;
00216       
00217       do {
00218         
00219         r3 = G4UniformRand();
00220         r4 = G4UniformRand();
00221         phi_re = twopi*r3 ;
00222         G4double sint2 = 1. - cost*cost ;
00223         G4double fp = 1. - sint2*loga/(2.*cost) ;
00224         rt = (1.-cos(2.*phi_re)*fp/f1_re)/(2.*pi) ;
00225         
00226       } while(rt < r4);
00227       
00228       // Calculo de la energia - elecron de recoil - relacion momento maximo <-> angulo
00229       
00230       G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
00231       G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2
00232         + (S - electron_mass_c2*electron_mass_c2)
00233         *(S - electron_mass_c2*electron_mass_c2)*sin(theta_re)*sin(theta_re);
00234       ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
00235       
00236       //      G4cout << "electron de retroceso " << ener_re << " " << theta_re << " " << phi_re << G4endl;
00237       
00238       // Recoil electron creation
00239       G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
00240       
00241       G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
00242       
00243       G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
00244       electronRDirection.rotateUz(photonDirection);
00245       
00246       G4DynamicParticle* particle3 = new G4DynamicParticle (G4Electron::Electron(),
00247                                                             electronRDirection,
00248                                                             electronRKineEnergy);
00249       fvect->push_back(particle3);          
00250       
00251     }
00252   else
00253     {
00254       // deposito la energia  ener_re - electron_mass_c2
00255       // G4cout << "electron de retroceso " << ener_re << G4endl;
00256       fParticleChange->ProposeLocalEnergyDeposit(ener_re - electron_mass_c2);
00257     }
00258   
00259   // Depaola (2004) suggested distribution for e+e- energy
00260   
00261   //  G4double t = 0.5*asinh(momentumThreshold_N);
00262   G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
00263   
00264   G4cout << 0.5*asinh(momentumThreshold_N) << "  " << t << G4endl;
00265 
00266   G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
00267   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));
00268   G4double b = 2.*(J1-J2)/J1;
00269   
00270   G4double n = 1 - b/6.;
00271   G4double re=0.;
00272   re = G4UniformRand();
00273   G4double a = 0.;
00274   G4double b1 =  16. - 3.*b - 36.*b*re*n + 36.*b*pow(re,2.)*pow(n,2.) + 
00275     6.*pow(b,2.)*re*n;
00276   a = pow((b1/b),0.5);
00277   G4double c1 = (-6. + 12.*re*n + b + 2*a)*pow(b,2.);
00278   epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
00279   
00280   G4double photonEnergy1 = photonEnergy - ener_re ; // resto al foton la energia del electron de retro.
00281   positronTotEnergy = epsilon*photonEnergy1;
00282   electronTotEnergy = photonEnergy1 - positronTotEnergy; // temporarly
00283   
00284   G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy - 
00285                             electron_mass_c2*electron_mass_c2) ;
00286   G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy - 
00287                             electron_mass_c2*electron_mass_c2) ;
00288   
00289   thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
00290   thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
00291   phi  = twopi * G4UniformRand();
00292   
00293   G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
00294   G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
00295   
00296   
00297   // Kinematics of the created pair:
00298   // the electron and positron are assumed to have a symetric angular 
00299   // distribution with respect to the Z axis along the parent photon
00300   
00301   G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
00302   
00303   // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
00304   
00305   G4ThreeVector electronDirection (dxEle, dyEle, dzEle);
00306   electronDirection.rotateUz(photonDirection);
00307   
00308   G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
00309                                                         electronDirection,
00310                                                         electronKineEnergy);
00311   
00312   // The e+ is always created (even with kinetic energy = 0) for further annihilation
00313   G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
00314   
00315   // SI - The range test has been removed wrt original G4LowEnergyGammaconversion class
00316   
00317   G4ThreeVector positronDirection (dxPos, dyPos, dzPos);
00318   positronDirection.rotateUz(photonDirection);   
00319   
00320   // Create G4DynamicParticle object for the particle2 
00321   G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
00322                                                        positronDirection, positronKineEnergy);
00323   // Fill output vector
00324   
00325   
00326   fvect->push_back(particle1);
00327   fvect->push_back(particle2);
00328   
00329 
00330   
00331   
00332   // kill incident photon
00333   fParticleChange->SetProposedKineticEnergy(0.);
00334   fParticleChange->ProposeTrackStatus(fStopAndKill);   
00335   
00336 }

G4ParticleChangeForGamma* G4BoldyshevTripletModel::fParticleChange [protected]

Definition at line 70 of file G4BoldyshevTripletModel.hh.

Referenced by Initialise(), and SampleSecondaries().


Public Member Functions
	G4BoldyshevTripletModel (const G4ParticleDefinition *p=0, const G4String &nam="BoldyshevTriplet")
virtual	~G4BoldyshevTripletModel ()
virtual void	Initialise (const G4ParticleDefinition *, const G4DataVector &)
virtual G4double	ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
virtual void	SampleSecondaries (std::vector< G4DynamicParticle * > , const G4MaterialCutsCouple , const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
Protected Member Functions
G4double	GetMeanFreePath (const G4Track &aTrack, G4double previousStepSize, G4ForceCondition *condition)
Protected Attributes
G4ParticleChangeForGamma *	fParticleChange

G4BoldyshevTripletModel Class Reference

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation