#include <FastSimulation/MaterialEffects/interface/PairProductionSimulator.h>
Public Member Functions | |
PairProductionSimulator (double photonEnergyCut, const RandomEngine *engine) | |
Constructor. | |
~PairProductionSimulator () | |
Default Destructor. | |
Private Member Functions | |
void | compute (ParticlePropagator &Particle) |
Generate an e+e- pair according to the probability that it happens. | |
double | gbteth (double ener, double partm, double efrac) |
A universal angular distribution - still from GEANT. | |
Private Attributes | |
double | photonEnergy |
The minimal photon energy for possible conversion. |
Definition at line 24 of file PairProductionSimulator.h.
PairProductionSimulator::PairProductionSimulator | ( | double | photonEnergyCut, | |
const RandomEngine * | engine | |||
) |
Constructor.
Definition at line 6 of file PairProductionSimulator.cc.
References max, and photonEnergy.
00007 : 00008 MaterialEffectsSimulator(engine) 00009 { 00010 00011 // Set the minimal photon energy for possible conversion 00012 photonEnergy = std::max(0.100,photonEnergyCut); 00013 00014 }
PairProductionSimulator::~PairProductionSimulator | ( | ) | [inline] |
void PairProductionSimulator::compute | ( | ParticlePropagator & | Particle | ) | [private, virtual] |
Generate an e+e- pair according to the probability that it happens.
Implements MaterialEffectsSimulator.
Definition at line 18 of file PairProductionSimulator.cc.
References MaterialEffectsSimulator::_theUpdatedState, funct::cos(), MaterialEffectsSimulator::eMass(), RandomEngine::flatShoot(), gbteth(), funct::log(), max, phi, photonEnergy, MaterialEffectsSimulator::radLengths, MaterialEffectsSimulator::random, funct::sin(), funct::sqrt(), and weight.
00019 { 00020 00021 double eGamma = Particle.e(); 00022 00023 // The photon has enough energy to create a pair 00024 if ( eGamma>=photonEnergy ) { 00025 00026 // This is a simple version (a la PDG) of a photon conversion generator. 00027 // It replaces the buggy GEANT3 -> C++ former version. 00028 // Author : Patrick Janot - 7-Jan-2004 00029 00030 // Probability to convert is 7/9*(dx/X0) 00031 if ( -std::log(random->flatShoot()) <= (7./9.)*radLengths ) { 00032 00033 double xe=0; 00034 double xm=eMass()/eGamma; 00035 double weight = 0.; 00036 00037 // Generate electron energy between emass and eGamma-emass 00038 do { 00039 xe = random->flatShoot()*(1.-2.*xm) + xm; 00040 weight = 1. - 4./3.*xe*(1.-xe); 00041 } while ( weight < random->flatShoot() ); 00042 00043 double eElectron = xe * eGamma; 00044 double tElectron = eElectron-eMass(); 00045 double pElectron = std::sqrt(std::max((eElectron+eMass())*tElectron,0.)); 00046 00047 double ePositron = eGamma-eElectron; 00048 double tPositron = ePositron-eMass(); 00049 double pPositron = std::sqrt((ePositron+eMass())*tPositron); 00050 00051 // Generate angles 00052 double phi = random->flatShoot()*2.*M_PI; 00053 double sphi = std::sin(phi); 00054 double cphi = std::cos(phi); 00055 00056 double stheta1, stheta2, ctheta1, ctheta2; 00057 00058 if ( eElectron > ePositron ) { 00059 double theta1 = gbteth(eElectron,eMass(),xe)*eMass()/eElectron; 00060 stheta1 = std::sin(theta1); 00061 ctheta1 = std::cos(theta1); 00062 stheta2 = stheta1*pElectron/pPositron; 00063 ctheta2 = std::sqrt(std::max(0.,1.0-(stheta2*stheta2))); 00064 } else { 00065 double theta2 = gbteth(ePositron,eMass(),xe)*eMass()/ePositron; 00066 stheta2 = std::sin(theta2); 00067 ctheta2 = std::cos(theta2); 00068 stheta1 = stheta2*pPositron/pElectron; 00069 ctheta1 = std::sqrt(std::max(0.,1.0-(stheta1*stheta1))); 00070 } 00071 00072 00073 double chi = Particle.theta(); 00074 double psi = Particle.phi(); 00075 RawParticle::RotationZ rotZ(psi); 00076 RawParticle::RotationY rotY(chi); 00077 00078 _theUpdatedState.resize(2,RawParticle()); 00079 00080 // The eletron 00081 _theUpdatedState[0].SetXYZT(pElectron*stheta1*cphi, 00082 pElectron*stheta1*sphi, 00083 pElectron*ctheta1, 00084 eElectron); 00085 _theUpdatedState[0].setID(+11); 00086 _theUpdatedState[0].rotate(rotY); 00087 _theUpdatedState[0].rotate(rotZ); 00088 00089 // The positron 00090 _theUpdatedState[1].SetXYZT(-pPositron*stheta2*cphi, 00091 -pPositron*stheta2*sphi, 00092 pPositron*ctheta2, 00093 ePositron); 00094 _theUpdatedState[1].setID(-11); 00095 _theUpdatedState[1].rotate(rotY); 00096 _theUpdatedState[1].rotate(rotZ); 00097 00098 } 00099 } 00100 }
double PairProductionSimulator::gbteth | ( | double | ener, | |
double | partm, | |||
double | efrac | |||
) | [private] |
A universal angular distribution - still from GEANT.
Definition at line 103 of file PairProductionSimulator.cc.
References DeDxTools::beta(), d, RandomEngine::flatShoot(), funct::log(), MaterialEffectsSimulator::random, MaterialEffectsSimulator::theZ(), and w1.
Referenced by compute().
00104 { 00105 const double alfa = 0.625; 00106 00107 double d = 0.13*(0.8+1.3/theZ())*(100.0+(1.0/ener))*(1.0+efrac); 00108 double w1 = 9.0/(9.0+d); 00109 double umax = ener*M_PI/partm; 00110 double u; 00111 00112 do { 00113 double beta; 00114 if (random->flatShoot()<=w1) beta = alfa; 00115 else beta = 3.0*alfa; 00116 u = -(std::log(random->flatShoot()*random->flatShoot()))/beta; 00117 } while (u>=umax); 00118 00119 return u; 00120 }
double PairProductionSimulator::photonEnergy [private] |
The minimal photon energy for possible conversion.
Definition at line 38 of file PairProductionSimulator.h.
Referenced by compute(), and PairProductionSimulator().