#include <G4eeToHadronsModel.hh>
Inheritance diagram for G4eeToHadronsModel:
Public Member Functions | |
G4eeToHadronsModel (G4Vee2hadrons *, G4int ver=0, const G4String &nam="eeToHadrons") | |
virtual | ~G4eeToHadronsModel () |
virtual void | Initialise (const G4ParticleDefinition *, const G4DataVector &) |
virtual G4double | CrossSectionPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) |
virtual G4double | ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kineticEnergy, G4double Z, G4double A, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) |
virtual G4double | ComputeCrossSectionPerElectron (const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX) |
virtual void | SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double maxEnergy=DBL_MAX) |
G4DynamicParticle * | GenerateCMPhoton (G4double) |
G4double | PeakEnergy () const |
Definition at line 59 of file G4eeToHadronsModel.hh.
G4eeToHadronsModel::G4eeToHadronsModel | ( | G4Vee2hadrons * | , | |
G4int | ver = 0 , |
|||
const G4String & | nam = "eeToHadrons" | |||
) |
Definition at line 68 of file G4eeToHadronsModel.cc.
References G4Gamma::Gamma(), G4VEmModel::HighEnergyLimit(), and G4VEmModel::LowEnergyLimit().
00070 : G4VEmModel(nam), 00071 model(mod), 00072 crossPerElectron(0), 00073 crossBornPerElectron(0), 00074 isInitialised(false), 00075 nbins(100), 00076 verbose(ver) 00077 { 00078 theGamma = G4Gamma::Gamma(); 00079 highKinEnergy = HighEnergyLimit(); 00080 lowKinEnergy = LowEnergyLimit(); 00081 emin = lowKinEnergy; 00082 emax = highKinEnergy; 00083 peakKinEnergy = highKinEnergy; 00084 epeak = emax; 00085 }
G4eeToHadronsModel::~G4eeToHadronsModel | ( | ) | [virtual] |
Definition at line 89 of file G4eeToHadronsModel.cc.
00090 { 00091 delete model; 00092 delete crossPerElectron; 00093 delete crossBornPerElectron; 00094 }
G4double G4eeToHadronsModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kineticEnergy, | |||
G4double | Z, | |||
G4double | A, | |||
G4double | cutEnergy = 0.0 , |
|||
G4double | maxEnergy = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 195 of file G4eeToHadronsModel.cc.
References ComputeCrossSectionPerElectron().
00200 { 00201 return Z*ComputeCrossSectionPerElectron(p, kineticEnergy); 00202 }
G4double G4eeToHadronsModel::ComputeCrossSectionPerElectron | ( | const G4ParticleDefinition * | , | |
G4double | kineticEnergy, | |||
G4double | cutEnergy = 0.0 , |
|||
G4double | maxEnergy = DBL_MAX | |||
) | [virtual] |
Definition at line 206 of file G4eeToHadronsModel.cc.
References G4PhysicsVector::GetValue().
Referenced by ComputeCrossSectionPerAtom(), and CrossSectionPerVolume().
00210 { 00211 G4double cross = 0.0; 00212 if(crossPerElectron) { 00213 G4bool b; 00214 G4double e = 2.0*electron_mass_c2* 00215 sqrt(1.0 + 0.5*kineticEnergy/electron_mass_c2); 00216 cross = crossPerElectron->GetValue(e, b); 00217 } 00218 // G4cout << "e= " << kineticEnergy << " cross= " << cross << G4endl; 00219 return cross; 00220 }
G4double G4eeToHadronsModel::CrossSectionPerVolume | ( | const G4Material * | , | |
const G4ParticleDefinition * | , | |||
G4double | kineticEnergy, | |||
G4double | cutEnergy, | |||
G4double | maxEnergy | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 183 of file G4eeToHadronsModel.cc.
References ComputeCrossSectionPerElectron(), and G4Material::GetElectronDensity().
00188 { 00189 return mat->GetElectronDensity()* 00190 ComputeCrossSectionPerElectron(p, kineticEnergy); 00191 }
G4DynamicParticle * G4eeToHadronsModel::GenerateCMPhoton | ( | G4double | ) |
Definition at line 303 of file G4eeToHadronsModel.cc.
References G4cout, G4endl, G4UniformRand, G4PhysicsVector::GetValue(), G4INCL::Math::pi, and G4InuclParticleNames::s0.
Referenced by SampleSecondaries().
00304 { 00305 G4bool b; 00306 G4double x; 00307 G4DynamicParticle* gamma = 0; 00308 G4double L = 2.0*log(e/electron_mass_c2); 00309 G4double bt = 2.0*fine_structure_const*(L - 1.)/pi; 00310 G4double btm1= bt - 1.0; 00311 G4double del = 1. + fine_structure_const*(1.5*L + pi*pi/3. -2.)/pi; 00312 00313 G4double s0 = crossBornPerElectron->GetValue(e, b); 00314 G4double de = (emax - emin)/(G4double)nbins; 00315 G4double x0 = min(de,e - emin)/e; 00316 G4double ds = crossBornPerElectron->GetValue(e, b) 00317 *(del*pow(x0,bt) - bt*(x0 - 0.25*x0*x0)); 00318 G4double e1 = e*(1. - x0); 00319 00320 if(e1 < emax && s0*G4UniformRand()<ds) { 00321 x = x0*pow(G4UniformRand(),1./bt); 00322 } else { 00323 00324 x = 1. - e1/e; 00325 G4double s1 = crossBornPerElectron->GetValue(e1, b); 00326 G4double w1 = bt*(del*pow(x,btm1) - 1.0 + 0.5*x); 00327 G4double grej = s1*w1; 00328 G4double f; 00329 // G4cout << "e= " << e/GeV << " epeak= " << epeak/GeV 00330 // << " s1= " << s1 << " w1= " << w1 00331 // << " grej= " << grej << G4endl; 00332 // Above emax cross section is 0 00333 if(e1 > emax) { 00334 x = 1. - emax/e; 00335 G4double s2 = crossBornPerElectron->GetValue(emax, b); 00336 G4double w2 = bt*(del*pow(x,btm1) - 1.0 + 0.5*x); 00337 grej = s2*w2; 00338 // G4cout << "emax= " << emax << " s2= " << s2 << " w2= " << w2 00339 // << " grej= " << grej << G4endl; 00340 } 00341 00342 if(e1 > epeak) { 00343 x = 1. - epeak/e; 00344 G4double s2 = crossBornPerElectron->GetValue(epeak, b); 00345 G4double w2 = bt*(del*pow(x,btm1) - 1.0 + 0.5*x); 00346 grej = max(grej,s2*w2); 00347 //G4cout << "epeak= " << epeak << " s2= " << s2 << " w2= " << w2 00348 // << " grej= " << grej << G4endl; 00349 } 00350 G4double xmin = 1. - e1/e; 00351 if(e1 > emax) xmin = 1. - emax/e; 00352 G4double xmax = 1. - emin/e; 00353 do { 00354 x = xmin + G4UniformRand()*(xmax - xmin); 00355 G4double s2 = crossBornPerElectron->GetValue((1.0 - x)*e, b); 00356 G4double w2 = bt*(del*pow(x,btm1) - 1.0 + 0.5*x); 00357 //G4cout << "x= " << x << " xmin= " << xmin << " xmax= " << xmax 00358 // << " s2= " << s2 << " w2= " << w2 00359 // << G4endl; 00360 f = s2*w2; 00361 if(f > grej) { 00362 G4cout << "G4DynamicParticle* G4eeToHadronsModel:WARNING " 00363 << f << " > " << grej << " majorant is`small!" 00364 << G4endl; 00365 } 00366 } while (f < grej*G4UniformRand()); 00367 } 00368 00369 G4ThreeVector dir(0.0,0.0,1.0); 00370 gamma = new G4DynamicParticle(theGamma,dir,x*e); 00371 return gamma; 00372 }
void G4eeToHadronsModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Definition at line 98 of file G4eeToHadronsModel.cc.
References G4Vee2hadrons::ComputeCrossSection(), G4cout, G4endl, G4PhysicsVector::GetLowEdgeEnergy(), G4PhysicsVector::GetValue(), G4PhysicsVector::GetVectorLength(), G4Vee2hadrons::HighEnergy(), G4VEmModel::HighEnergyLimit(), G4Vee2hadrons::LowEnergy(), G4VEmModel::LowEnergyLimit(), G4Vee2hadrons::PeakEnergy(), G4Vee2hadrons::PhysicsVector(), G4PhysicsVector::PutValue(), G4VEmModel::SetHighEnergyLimit(), G4Vee2hadrons::SetLowEnergy(), and G4VEmModel::SetLowEnergyLimit().
00100 { 00101 if(isInitialised) { return; } 00102 isInitialised = true; 00103 00104 // Lab system 00105 highKinEnergy = HighEnergyLimit(); 00106 lowKinEnergy = LowEnergyLimit(); 00107 00108 // CM system 00109 emin = model->LowEnergy(); 00110 emax = model->HighEnergy(); 00111 00112 G4double emin0 = 00113 2.0*electron_mass_c2*sqrt(1.0 + 0.5*lowKinEnergy/electron_mass_c2); 00114 G4double emax0 = 00115 2.0*electron_mass_c2*sqrt(1.0 + 0.5*highKinEnergy/electron_mass_c2); 00116 00117 // recompute low energy 00118 if(emin0 > emax) { 00119 emin0 = emax; 00120 model->SetLowEnergy(emin0); 00121 } 00122 if(emin > emin0) { 00123 emin0 = emin; 00124 lowKinEnergy = 0.5*emin*emin/electron_mass_c2 - 2.0*electron_mass_c2; 00125 SetLowEnergyLimit(lowKinEnergy); 00126 } 00127 00128 // recompute high energy 00129 if(emax < emax0) { 00130 emax0 = emax; 00131 highKinEnergy = 0.5*emax*emax/electron_mass_c2 - 2.0*electron_mass_c2; 00132 SetHighEnergyLimit(highKinEnergy); 00133 } 00134 00135 // peak energy 00136 epeak = std::min(model->PeakEnergy(), emax); 00137 peakKinEnergy = 0.5*epeak*epeak/electron_mass_c2 - 2.0*electron_mass_c2; 00138 00139 if(verbose>0) { 00140 G4cout << "G4eeToHadronsModel::Initialise: " << G4endl; 00141 G4cout << "LabSystem: emin(GeV)= " << lowKinEnergy/GeV 00142 << " epeak(GeV)= " << peakKinEnergy/GeV 00143 << " emax(GeV)= " << highKinEnergy/GeV 00144 << G4endl; 00145 G4cout << "SM System: emin(MeV)= " << emin/MeV 00146 << " epeak(MeV)= " << epeak/MeV 00147 << " emax(MeV)= " << emax/MeV 00148 << G4endl; 00149 } 00150 00151 if(lowKinEnergy < peakKinEnergy) { 00152 crossBornPerElectron = model->PhysicsVector(emin, emax); 00153 crossPerElectron = model->PhysicsVector(emin, emax); 00154 nbins = crossPerElectron->GetVectorLength(); 00155 for(G4int i=0; i<nbins; i++) { 00156 G4double e = crossPerElectron->GetLowEdgeEnergy(i); 00157 G4double cs = model->ComputeCrossSection(e); 00158 crossBornPerElectron->PutValue(i, cs); 00159 } 00160 ComputeCMCrossSectionPerElectron(); 00161 } 00162 if(verbose>1) { 00163 G4cout << "G4eeToHadronsModel: Cross secsions per electron" 00164 << " nbins= " << nbins 00165 << " emin(MeV)= " << emin/MeV 00166 << " emax(MeV)= " << emax/MeV 00167 << G4endl; 00168 G4bool b; 00169 for(G4int i=0; i<nbins; i++) { 00170 G4double e = crossPerElectron->GetLowEdgeEnergy(i); 00171 G4double s1 = crossPerElectron->GetValue(e, b); 00172 G4double s2 = crossBornPerElectron->GetValue(e, b); 00173 G4cout << "E(MeV)= " << e/MeV 00174 << " cross(nb)= " << s1/nanobarn 00175 << " crossBorn(nb)= " << s2/nanobarn 00176 << G4endl; 00177 } 00178 } 00179 }
G4double G4eeToHadronsModel::PeakEnergy | ( | ) | const [inline] |
void G4eeToHadronsModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin = 0.0 , |
|||
G4double | maxEnergy = DBL_MAX | |||
) | [virtual] |
Implements G4VEmModel.
Definition at line 224 of file G4eeToHadronsModel.cc.
References GenerateCMPhoton(), G4DynamicParticle::Get4Momentum(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMomentumDirection(), G4Vee2hadrons::SampleSecondaries(), and G4DynamicParticle::Set4Momentum().
00229 { 00230 if(crossPerElectron) { 00231 G4double t = dParticle->GetKineticEnergy(); 00232 G4double e = 2.0*electron_mass_c2*sqrt(1.0 + 0.5*t/electron_mass_c2); 00233 G4LorentzVector inlv = dParticle->Get4Momentum(); 00234 G4ThreeVector inBoost = inlv.boostVector(); 00235 if(e > emin) { 00236 G4DynamicParticle* gamma = GenerateCMPhoton(e); 00237 G4LorentzVector gLv = gamma->Get4Momentum(); 00238 G4LorentzVector lv(0.0,0.0,0.0,e); 00239 lv -= gLv; 00240 G4double mass = lv.m(); 00241 G4ThreeVector boost = lv.boostVector(); 00242 const G4ThreeVector dir = gamma->GetMomentumDirection(); 00243 model->SampleSecondaries(newp, mass, dir); 00244 G4int np = newp->size(); 00245 for(G4int j=0; j<np; j++) { 00246 G4DynamicParticle* dp = (*newp)[j]; 00247 G4LorentzVector v = dp->Get4Momentum(); 00248 v.boost(boost); 00249 v.boost(inBoost); 00250 dp->Set4Momentum(v); 00251 } 00252 gLv.boost(inBoost); 00253 gamma->Set4Momentum(gLv); 00254 newp->push_back(gamma); 00255 } 00256 } 00257 }