Line data Source code
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 0 : EvtVPHOtoVISR::~EvtVPHOtoVISR() {}
33 :
34 : std::string EvtVPHOtoVISR::getName(){
35 :
36 0 : return "VPHOTOVISR";
37 :
38 : }
39 :
40 :
41 : EvtDecayBase* EvtVPHOtoVISR::clone(){
42 :
43 0 : return new EvtVPHOtoVISR;
44 :
45 0 : }
46 :
47 : void EvtVPHOtoVISR::init(){
48 :
49 : // check that there are 0 or 2 arguments
50 0 : checkNArg(0,2);
51 :
52 : // check that there are 2 daughters
53 0 : checkNDaug(2);
54 :
55 : // check the parent and daughter spins
56 0 : checkSpinParent(EvtSpinType::VECTOR);
57 0 : checkSpinDaughter(0,EvtSpinType::VECTOR);
58 0 : checkSpinDaughter(1,EvtSpinType::PHOTON);
59 0 : }
60 :
61 : void EvtVPHOtoVISR::initProbMax() {
62 :
63 : //setProbMax(100000.0);
64 :
65 0 : }
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 0 : double w=p->mass();
74 0 : double s=w*w;
75 :
76 0 : double L=2.0*log(w/0.000511);
77 : double alpha=1/137.0;
78 0 : 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 0 : assert(p->getDaug(0)->getDaug(0)!=0);
83 0 : double md=EvtPDL::getMeanMass(p->getDaug(0)->getDaug(0)->getId());
84 :
85 0 : static double mD0=EvtPDL::getMeanMass(EvtPDL::getId("D0"));
86 0 : static double mDp=EvtPDL::getMeanMass(EvtPDL::getId("D+"));
87 :
88 0 : double pgmax=(s-4.0*md*md)/(2.0*w);
89 :
90 0 : assert(pgmax>0.0);
91 :
92 0 : double pgz=0.99*pgmax*exp(log(EvtRandom::Flat(1.0))/beta);
93 :
94 0 : if (EvtRandom::Flat(1.0)<0.5) pgz=-pgz;
95 :
96 0 : double k=fabs(pgz);
97 :
98 0 : EvtVector4R p4g(k,0.0,0.0,pgz);
99 :
100 0 : EvtVector4R p4res=p->getP4Restframe()-p4g;
101 :
102 0 : double mres=p4res.mass();
103 :
104 0 : double ed=mres/2.0;
105 :
106 0 : assert(ed>md);
107 :
108 0 : 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 0 : p->getDaug(0)->init(getDaug(0),p4res);
114 0 : p->getDaug(1)->init(getDaug(1),p4g);
115 :
116 :
117 0 : double sigma=beta*pow(2/w,beta)*(1+alpha*(1.5*L-2.0+EvtConst::pi*EvtConst::pi/3.0)/EvtConst::pi);
118 :
119 0 : double m=EvtPDL::getMeanMass(p->getDaug(0)->getId());
120 0 : 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 0 : if (ed>mD0) p0=sqrt(ed*ed-mD0*mD0);
126 : double pp=0.0;
127 0 : if (ed>mDp) pp=sqrt(ed*ed-mDp*mDp);
128 :
129 0 : double p0norm=sqrt(0.25*m*m-mD0*mD0);
130 0 : double ppnorm=sqrt(0.25*m*m-mDp*mDp);
131 :
132 : double r0=12.7;
133 : double rp=12.7;
134 :
135 0 : if (getNArg()==2){
136 0 : r0=getArg(0);
137 0 : rp=getArg(1);
138 0 : }
139 :
140 0 : double GammaTot=Gamma*(pp*pp*pp/(1+pp*pp*rp*rp)+p0*p0*p0/(1+p0*p0*r0*r0))/
141 0 : (ppnorm*ppnorm*ppnorm/(1+ppnorm*ppnorm*rp*rp)+
142 0 : p0norm*p0norm*p0norm/(1+p0norm*p0norm*r0*r0));
143 :
144 :
145 0 : sigma*=pd*pd*pd/((mres-m)*(mres-m)+0.25*GammaTot*GammaTot);
146 :
147 0 : assert(sigma>0.0);
148 :
149 0 : static double sigmax=sigma;
150 :
151 0 : if (sigma>sigmax){
152 0 : sigmax=sigma;
153 0 : }
154 :
155 :
156 :
157 : static int count=0;
158 :
159 0 : count++;
160 :
161 : //if (count%10000==0){
162 : // std::cout << "sigma :"<<sigma<<std::endl;
163 : // std::cout << "sigmax:"<<sigmax<<std::endl;
164 : //}
165 :
166 0 : double norm=sqrt(sigma);
167 :
168 0 : vertex(0,0,0,norm*p->eps(0)*p->epsParent(0).conj());
169 0 : vertex(1,0,0,norm*p->eps(1)*p->epsParent(0).conj());
170 0 : vertex(2,0,0,norm*p->eps(2)*p->epsParent(0).conj());
171 :
172 0 : vertex(0,1,0,norm*p->eps(0)*p->epsParent(1).conj());
173 0 : vertex(1,1,0,norm*p->eps(1)*p->epsParent(1).conj());
174 0 : vertex(2,1,0,norm*p->eps(2)*p->epsParent(1).conj());
175 :
176 0 : vertex(0,2,0,norm*p->eps(0)*p->epsParent(2).conj());
177 0 : vertex(1,2,0,norm*p->eps(1)*p->epsParent(2).conj());
178 0 : vertex(2,2,0,norm*p->eps(2)*p->epsParent(2).conj());
179 :
180 0 : vertex(0,0,1,norm*p->eps(0)*p->epsParent(0).conj());
181 0 : vertex(1,0,1,norm*p->eps(1)*p->epsParent(0).conj());
182 0 : vertex(2,0,1,norm*p->eps(2)*p->epsParent(0).conj());
183 :
184 0 : vertex(0,1,1,norm*p->eps(0)*p->epsParent(1).conj());
185 0 : vertex(1,1,1,norm*p->eps(1)*p->epsParent(1).conj());
186 0 : vertex(2,1,1,norm*p->eps(2)*p->epsParent(1).conj());
187 :
188 0 : vertex(0,2,1,norm*p->eps(0)*p->epsParent(2).conj());
189 0 : vertex(1,2,1,norm*p->eps(1)*p->epsParent(2).conj());
190 0 : vertex(2,2,1,norm*p->eps(2)*p->epsParent(2).conj());
191 :
192 : return;
193 0 : }
194 :
|
