CMS 3D CMS Logo

PairProductionSimulator.cc
Go to the documentation of this file.
4 
5 #include <cmath>
6 
8 {
9  // Set the minimal photon energy for possible conversion
10  photonEnergy = std::max(0.100,photonEnergyCut);
11 }
12 
13 void
15 {
16 
17  double eGamma = Particle.particle().e();
18 
19  // The photon has enough energy to create a pair
20  if ( eGamma>=photonEnergy ) {
21 
22  // This is a simple version (a la PDG) of a photon conversion generator.
23  // It replaces the buggy GEANT3 -> C++ former version.
24  // Author : Patrick Janot - 7-Jan-2004
25 
26  // Probability to convert is 7/9*(dx/X0)
27  if ( -std::log(random->flatShoot()) <= (7./9.)*radLengths ) {
28 
29  double xe=0;
30  double xm=eMass()/eGamma;
31  double weight = 0.;
32 
33  // Generate electron energy between emass and eGamma-emass
34  do {
35  xe = random->flatShoot()*(1.-2.*xm) + xm;
36  weight = 1. - 4./3.*xe*(1.-xe);
37  } while ( weight < random->flatShoot() );
38 
39  double eElectron = xe * eGamma;
40  double tElectron = eElectron-eMass();
41  double pElectron = std::sqrt(std::max((eElectron+eMass())*tElectron,0.));
42 
43  double ePositron = eGamma-eElectron;
44  double tPositron = ePositron-eMass();
45  double pPositron = std::sqrt((ePositron+eMass())*tPositron);
46 
47  // Generate angles
48  double phi = random->flatShoot()*2.*M_PI;
49  double sphi = std::sin(phi);
50  double cphi = std::cos(phi);
51 
52  double stheta1, stheta2, ctheta1, ctheta2;
53 
54  if ( eElectron > ePositron ) {
55  double theta1 = gbteth(eElectron,eMass(),xe,random)*eMass()/eElectron;
56  stheta1 = std::sin(theta1);
57  ctheta1 = std::cos(theta1);
58  stheta2 = stheta1*pElectron/pPositron;
59  ctheta2 = std::sqrt(std::max(0.,1.0-(stheta2*stheta2)));
60  } else {
61  double theta2 = gbteth(ePositron,eMass(),xe,random)*eMass()/ePositron;
62  stheta2 = std::sin(theta2);
63  ctheta2 = std::cos(theta2);
64  stheta1 = stheta2*pPositron/pElectron;
65  ctheta1 = std::sqrt(std::max(0.,1.0-(stheta1*stheta1)));
66  }
67 
68 
69  double chi = Particle.particle().theta();
70  double psi = Particle.particle().phi();
71  RawParticle::RotationZ rotZ(psi);
72  RawParticle::RotationY rotY(chi);
73 
74  _theUpdatedState.reserve(2);
75  _theUpdatedState.clear();
76 
77  // The electron
78  _theUpdatedState.emplace_back(
80  +11,XYZTLorentzVector(pElectron*stheta1*cphi,
81  pElectron*stheta1*sphi,
82  pElectron*ctheta1,
83  eElectron) ));
84  _theUpdatedState[0].rotate(rotY);
85  _theUpdatedState[0].rotate(rotZ);
86 
87  // The positron
88  _theUpdatedState.emplace_back(
90  -11,XYZTLorentzVector(-pPositron*stheta2*cphi,
91  -pPositron*stheta2*sphi,
92  pPositron*ctheta2,
93  ePositron)));
94  _theUpdatedState[1].rotate(rotY);
95  _theUpdatedState[1].rotate(rotZ);
96 
97  }
98  }
99 }
100 
101 double
102 PairProductionSimulator::gbteth(double ener,double partm,double efrac, RandomEngineAndDistribution const* random)
103 {
104  const double alfa = 0.625;
105 
106  double d = 0.13*(0.8+1.3/theZ())*(100.0+(1.0/ener))*(1.0+efrac);
107  double w1 = 9.0/(9.0+d);
108  double umax = ener*M_PI/partm;
109  double u;
110 
111  do {
112  double beta;
113  if (random->flatShoot()<=w1) beta = alfa;
114  else beta = 3.0*alfa;
115  u = -(std::log(random->flatShoot()*random->flatShoot()))/beta;
116  } while (u>=umax);
117 
118  return u;
119 }
PairProductionSimulator(double photonEnergyCut)
Constructor.
const HepPDT::ParticleDataTable * particleDataTable() const
double photonEnergy
The minimal photon energy for possible conversion.
double flatShoot(double xmin=0.0, double xmax=1.0) const
double gbteth(double ener, double partm, double efrac, RandomEngineAndDistribution const *)
A universal angular distribution - still from GEANT.
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
void compute(ParticlePropagator &Particle, RandomEngineAndDistribution const *) override
Generate an e+e- pair according to the probability that it happens.
Definition: weight.py:1
TRandom random
Definition: MVATrainer.cc:138
RawParticle const & particle() const
The particle being propagated.
std::map< std::string, int, std::less< std::string > > psi
ROOT::Math::RotationZ RotationZ
Definition: RawParticle.h:51
RawParticle makeParticle(HepPDT::ParticleDataTable const *, int id, const math::XYZTLorentzVector &p)
Definition: makeParticle.cc:29
double phi() const
phi of momentum vector
Definition: RawParticle.h:335
double e() const
energy of the momentum
Definition: RawParticle.h:324
T sqrt(T t)
Definition: SSEVec.h:18
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
ROOT::Math::RotationY RotationY
Definition: RawParticle.h:50
double eMass() const
Electron mass in GeV/c2.
#define M_PI
double theta() const
theta of momentum vector
Definition: RawParticle.h:334
std::vector< RawParticle > _theUpdatedState
math::XYZTLorentzVector XYZTLorentzVector
Definition: RawParticle.h:27