]>
Commit | Line | Data |
---|---|---|
1 | //-------------------------------------------------------------------------- | |
2 | // | |
3 | // Environment: | |
4 | // This software is part of the EvtGen package developed jointly | |
5 | // for the BaBar and CLEO collaborations. If you use all or part | |
6 | // of it, please give an appropriate acknowledgement. | |
7 | // | |
8 | // Copyright Information: See EvtGen/COPYRIGHT | |
9 | // Copyright (C) 2004 Cornell | |
10 | // | |
11 | // Module: EvtVPHOtoVISR.cc | |
12 | // | |
13 | // Description: Routine to decay vpho -> vector ISR photon | |
14 | // | |
15 | // Modification history: | |
16 | // | |
17 | // Ryd March 20, 2004 Module created | |
18 | // | |
19 | //------------------------------------------------------------------------ | |
20 | // | |
21 | #include <stdlib.h> | |
22 | #include "EvtGenBase/EvtParticle.hh" | |
23 | #include "EvtGenBase/EvtGenKine.hh" | |
24 | #include "EvtGenBase/EvtPDL.hh" | |
25 | #include "EvtGenBase/EvtVector4C.hh" | |
26 | #include "EvtGenBase/EvtVector4R.hh" | |
27 | #include "EvtGenBase/EvtReport.hh" | |
28 | #include "EvtGenBase/EvtRandom.hh" | |
29 | #include "EvtGenModels/EvtVPHOtoVISR.hh" | |
30 | #include <string> | |
31 | ||
32 | EvtVPHOtoVISR::~EvtVPHOtoVISR() {} | |
33 | ||
34 | std::string EvtVPHOtoVISR::getName(){ | |
35 | ||
36 | return "VPHOTOVISR"; | |
37 | ||
38 | } | |
39 | ||
40 | ||
41 | EvtDecayBase* EvtVPHOtoVISR::clone(){ | |
42 | ||
43 | return new EvtVPHOtoVISR; | |
44 | ||
45 | } | |
46 | ||
47 | void EvtVPHOtoVISR::init(){ | |
48 | ||
49 | // check that there are 0 or 2 arguments | |
50 | checkNArg(0,2); | |
51 | ||
52 | // check that there are 2 daughters | |
53 | checkNDaug(2); | |
54 | ||
55 | // check the parent and daughter spins | |
56 | checkSpinParent(EvtSpinType::VECTOR); | |
57 | checkSpinDaughter(0,EvtSpinType::VECTOR); | |
58 | checkSpinDaughter(1,EvtSpinType::PHOTON); | |
59 | } | |
60 | ||
61 | void EvtVPHOtoVISR::initProbMax() { | |
62 | ||
63 | //setProbMax(100000.0); | |
64 | ||
65 | } | |
66 | ||
67 | void EvtVPHOtoVISR::decay( EvtParticle *p){ | |
68 | ||
69 | //take photon along z-axis, either forward or backward. | |
70 | //Implement this as generating the photon momentum along | |
71 | //the z-axis uniformly | |
72 | ||
73 | double w=p->mass(); | |
74 | double s=w*w; | |
75 | ||
76 | double L=2.0*log(w/0.000511); | |
77 | double alpha=1/137.0; | |
78 | double beta=(L-1)*2.0*alpha/EvtConst::pi; | |
79 | ||
80 | //This uses the fact that there is a daughter of the | |
81 | //psi(3770) | |
82 | assert(p->getDaug(0)->getDaug(0)!=0); | |
83 | double md=EvtPDL::getMeanMass(p->getDaug(0)->getDaug(0)->getId()); | |
84 | ||
85 | static double mD0=EvtPDL::getMeanMass(EvtPDL::getId("D0")); | |
86 | static double mDp=EvtPDL::getMeanMass(EvtPDL::getId("D+")); | |
87 | ||
88 | double pgmax=(s-4.0*md*md)/(2.0*w); | |
89 | ||
90 | assert(pgmax>0.0); | |
91 | ||
92 | double pgz=0.99*pgmax*exp(log(EvtRandom::Flat(1.0))/beta); | |
93 | ||
94 | if (EvtRandom::Flat(1.0)<0.5) pgz=-pgz; | |
95 | ||
96 | double k=fabs(pgz); | |
97 | ||
98 | EvtVector4R p4g(k,0.0,0.0,pgz); | |
99 | ||
100 | EvtVector4R p4res=p->getP4Restframe()-p4g; | |
101 | ||
102 | double mres=p4res.mass(); | |
103 | ||
104 | double ed=mres/2.0; | |
105 | ||
106 | assert(ed>md); | |
107 | ||
108 | double pd=sqrt(ed*ed-md*md); | |
109 | ||
110 | ||
111 | //std::cout << "k, mres, w, md, ed, pd:"<<k<<" "<<mres<<" "<<w<<" "<<md<<" "<<ed<<" "<<pd<<std::endl; | |
112 | ||
113 | p->getDaug(0)->init(getDaug(0),p4res); | |
114 | p->getDaug(1)->init(getDaug(1),p4g); | |
115 | ||
116 | ||
117 | double sigma=beta*pow(2/w,beta)*(1+alpha*(1.5*L-2.0+EvtConst::pi*EvtConst::pi/3.0)/EvtConst::pi); | |
118 | ||
119 | double m=EvtPDL::getMeanMass(p->getDaug(0)->getId()); | |
120 | double Gamma=EvtPDL::getWidth(p->getDaug(0)->getId()); | |
121 | ||
122 | //mres is the energy of the psi(3770) | |
123 | ||
124 | double p0=0.0; | |
125 | if (ed>mD0) p0=sqrt(ed*ed-mD0*mD0); | |
126 | double pp=0.0; | |
127 | if (ed>mDp) pp=sqrt(ed*ed-mDp*mDp); | |
128 | ||
129 | double p0norm=sqrt(0.25*m*m-mD0*mD0); | |
130 | double ppnorm=sqrt(0.25*m*m-mDp*mDp); | |
131 | ||
132 | double r0=12.7; | |
133 | double rp=12.7; | |
134 | ||
135 | if (getNArg()==2){ | |
136 | r0=getArg(0); | |
137 | rp=getArg(1); | |
138 | } | |
139 | ||
140 | double GammaTot=Gamma*(pp*pp*pp/(1+pp*pp*rp*rp)+p0*p0*p0/(1+p0*p0*r0*r0))/ | |
141 | (ppnorm*ppnorm*ppnorm/(1+ppnorm*ppnorm*rp*rp)+ | |
142 | p0norm*p0norm*p0norm/(1+p0norm*p0norm*r0*r0)); | |
143 | ||
144 | ||
145 | sigma*=pd*pd*pd/((mres-m)*(mres-m)+0.25*GammaTot*GammaTot); | |
146 | ||
147 | assert(sigma>0.0); | |
148 | ||
149 | static double sigmax=sigma; | |
150 | ||
151 | if (sigma>sigmax){ | |
152 | sigmax=sigma; | |
153 | } | |
154 | ||
155 | ||
156 | ||
157 | static int count=0; | |
158 | ||
159 | count++; | |
160 | ||
161 | //if (count%10000==0){ | |
162 | // std::cout << "sigma :"<<sigma<<std::endl; | |
163 | // std::cout << "sigmax:"<<sigmax<<std::endl; | |
164 | //} | |
165 | ||
166 | double norm=sqrt(sigma); | |
167 | ||
168 | vertex(0,0,0,norm*p->eps(0)*p->epsParent(0).conj()); | |
169 | vertex(1,0,0,norm*p->eps(1)*p->epsParent(0).conj()); | |
170 | vertex(2,0,0,norm*p->eps(2)*p->epsParent(0).conj()); | |
171 | ||
172 | vertex(0,1,0,norm*p->eps(0)*p->epsParent(1).conj()); | |
173 | vertex(1,1,0,norm*p->eps(1)*p->epsParent(1).conj()); | |
174 | vertex(2,1,0,norm*p->eps(2)*p->epsParent(1).conj()); | |
175 | ||
176 | vertex(0,2,0,norm*p->eps(0)*p->epsParent(2).conj()); | |
177 | vertex(1,2,0,norm*p->eps(1)*p->epsParent(2).conj()); | |
178 | vertex(2,2,0,norm*p->eps(2)*p->epsParent(2).conj()); | |
179 | ||
180 | vertex(0,0,1,norm*p->eps(0)*p->epsParent(0).conj()); | |
181 | vertex(1,0,1,norm*p->eps(1)*p->epsParent(0).conj()); | |
182 | vertex(2,0,1,norm*p->eps(2)*p->epsParent(0).conj()); | |
183 | ||
184 | vertex(0,1,1,norm*p->eps(0)*p->epsParent(1).conj()); | |
185 | vertex(1,1,1,norm*p->eps(1)*p->epsParent(1).conj()); | |
186 | vertex(2,1,1,norm*p->eps(2)*p->epsParent(1).conj()); | |
187 | ||
188 | vertex(0,2,1,norm*p->eps(0)*p->epsParent(2).conj()); | |
189 | vertex(1,2,1,norm*p->eps(1)*p->epsParent(2).conj()); | |
190 | vertex(2,2,1,norm*p->eps(2)*p->epsParent(2).conj()); | |
191 | ||
192 | return; | |
193 | } | |
194 |