]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 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 -> (DDx) + ISR photon from 3.9 to 4.3 GeV, using CLEO-c data (Brian Lang) | |
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/EvtVPHOtoVISRHi.hh" | |
30 | #include <string> | |
31 | ||
32 | using std::endl; | |
33 | ||
34 | EvtVPHOtoVISRHi::~EvtVPHOtoVISRHi() {} | |
35 | ||
36 | std::string EvtVPHOtoVISRHi::getName(){ | |
37 | ||
38 | return "VPHOTOVISRHI"; | |
39 | ||
40 | } | |
41 | ||
42 | ||
43 | EvtDecayBase* EvtVPHOtoVISRHi::clone(){ | |
44 | ||
45 | return new EvtVPHOtoVISRHi; | |
46 | ||
47 | } | |
48 | ||
49 | void EvtVPHOtoVISRHi::init(){ | |
50 | ||
51 | // check that there are 0 or 1 arguments | |
52 | checkNArg(0,1); | |
53 | ||
54 | // check that there are 2 daughters | |
55 | checkNDaug(2); | |
56 | ||
57 | // check the parent and daughter spins | |
58 | checkSpinParent(EvtSpinType::VECTOR); | |
59 | checkSpinDaughter(0,EvtSpinType::VECTOR); | |
60 | checkSpinDaughter(1,EvtSpinType::PHOTON); | |
61 | } | |
62 | ||
63 | void EvtVPHOtoVISRHi::initProbMax() { | |
64 | ||
65 | setProbMax(20.0); | |
66 | ||
67 | } | |
68 | ||
69 | void EvtVPHOtoVISRHi::decay( EvtParticle *p){ | |
70 | //take photon along z-axis, either forward or backward. | |
71 | //Implement this as generating the photon momentum along | |
72 | //the z-axis uniformly | |
73 | double power=1; | |
74 | if (getNArg()==1) power=getArg(0); | |
75 | // define particle names | |
76 | static EvtId D0=EvtPDL::getId("D0"); | |
77 | static EvtId D0B=EvtPDL::getId("anti-D0"); | |
78 | static EvtId DP=EvtPDL::getId("D+"); | |
79 | static EvtId DM=EvtPDL::getId("D-"); | |
80 | static EvtId DSM=EvtPDL::getId("D_s-"); | |
81 | static EvtId DSP=EvtPDL::getId("D_s+"); | |
82 | static EvtId DSMS=EvtPDL::getId("D_s*-"); | |
83 | static EvtId DSPS=EvtPDL::getId("D_s*+"); | |
84 | static EvtId D0S=EvtPDL::getId("D*0"); | |
85 | static EvtId D0BS=EvtPDL::getId("anti-D*0"); | |
86 | static EvtId DPS=EvtPDL::getId("D*+"); | |
87 | static EvtId DMS=EvtPDL::getId("D*-"); | |
88 | // setup some parameters | |
89 | double w=p->mass(); | |
90 | double s=w*w; | |
91 | double L=2.0*log(w/0.000511); | |
92 | double alpha=1/137.0; | |
93 | double beta=(L-1)*2.0*alpha/EvtConst::pi; | |
94 | // make sure only 2 or 3 body are present | |
95 | assert (p->getDaug(0)->getNDaug() == 2 || p->getDaug(0)->getNDaug() == 3); | |
96 | ||
97 | // determine minimum rest mass of parent | |
98 | double md1 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(0)->getId()); | |
99 | double md2 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(1)->getId()); | |
100 | double minResMass = md1+md2; | |
101 | if (p->getDaug(0)->getNDaug() == 3) { | |
102 | double md3 = EvtPDL::getMeanMass(p->getDaug(0)->getDaug(2)->getId()); | |
103 | minResMass = minResMass + md3; | |
104 | } | |
105 | ||
106 | // calculate the maximum energy of the ISR photon | |
107 | double pgmax=(s-minResMass*minResMass)/(2.0*w); | |
108 | double pgz=0.99*pgmax*exp(log(EvtRandom::Flat(1.0))/(beta*power)); | |
109 | if (EvtRandom::Flat(1.0)<0.5) pgz=-pgz; | |
110 | ||
111 | double k=fabs(pgz); | |
112 | // print of ISR energy | |
113 | // std::cout << "Energy ISR :"<< k <<std::endl; | |
114 | EvtVector4R p4g(k,0.0,0.0,pgz); | |
115 | ||
116 | EvtVector4R p4res=p->getP4Restframe()-p4g; | |
117 | ||
118 | double mres=p4res.mass(); | |
119 | ||
120 | // set masses | |
121 | p->getDaug(0)->init(getDaug(0),p4res); | |
122 | p->getDaug(1)->init(getDaug(1),p4g); | |
123 | ||
124 | ||
125 | // determine XS - langbw | |
126 | // very crude way of determining XS just a simple straight line Approx. | |
127 | // this was determined by eye. | |
128 | // lots of cout statements to make plots to check that things are working as expected | |
129 | double sigma=9.0; | |
130 | if (mres<=3.9) sigma = 0.00001; | |
131 | ||
132 | bool sigmacomputed(false); | |
133 | ||
134 | // DETERMINE XS FOR D*D* | |
135 | if (p->getDaug(0)->getNDaug() == 2 | |
136 | &&((p->getDaug(0)->getDaug(0)->getId()==D0S | |
137 | && p->getDaug(0)->getDaug(1)->getId()==D0BS) | |
138 | ||(p->getDaug(0)->getDaug(0)->getId()==DPS | |
139 | && p->getDaug(0)->getDaug(1)->getId()==DMS))){ | |
140 | if(mres>4.18) { | |
141 | sigma*=5./9.*(1.-1.*sqrt((4.18-mres)*(4.18-mres))/(4.3-4.18)); | |
142 | } | |
143 | else if(mres>4.07 && mres<=4.18) { | |
144 | sigma*=5./9.; | |
145 | } | |
146 | else if (mres<=4.07&&mres>4.03) | |
147 | { | |
148 | sigma*=(5./9. - 1.5/9.*sqrt((4.07-mres)*(4.07-mres))/(4.07-4.03)); | |
149 | } | |
150 | else if (mres<=4.03&& mres>=4.013) | |
151 | { | |
152 | sigma*=(3.5/9. - 3.5/9.*sqrt((4.03-mres)*(4.03-mres))/(4.03-4.013)); | |
153 | } | |
154 | else{ | |
155 | sigma=0.00001; | |
156 | } | |
157 | sigmacomputed = true; | |
158 | // std::cout << "DSDSXS "<<sigma<< " " << mres<<std::endl; | |
159 | } | |
160 | ||
161 | // DETERMINE XS FOR D*D | |
162 | if(p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==D0S | |
163 | && p->getDaug(0)->getDaug(1)->getId()==D0B) | |
164 | ||(p->getDaug(0)->getDaug(0)->getId()==DPS | |
165 | && p->getDaug(0)->getDaug(1)->getId()==DM) | |
166 | ||(p->getDaug(0)->getDaug(0)->getId()==D0BS | |
167 | && p->getDaug(0)->getDaug(1)->getId()==D0) | |
168 | ||(p->getDaug(0)->getDaug(0)->getId()==DMS | |
169 | && p->getDaug(0)->getDaug(1)->getId()==DP)) ) | |
170 | { | |
171 | if(mres>=4.2){ | |
172 | sigma*=1.5/9.; | |
173 | } | |
174 | else if( mres>4.06 && mres<4.2){ | |
175 | sigma*=((1.5/9.+2.5/9.*sqrt((4.2-mres)*(4.2-mres))/(4.2-4.06))); | |
176 | } | |
177 | else if(mres>=4.015 && mres<4.06){ | |
178 | sigma*=((4./9.+3./9.*sqrt((4.06-mres)*(4.06-mres))/(4.06-4.015))); | |
179 | } | |
180 | else if (mres<4.015 && mres>=3.9){ | |
181 | sigma*=((7./9.-7/9.*sqrt((4.015-mres)*(4.015-mres))/(4.015-3.9))); | |
182 | } | |
183 | else { | |
184 | sigma = 0.00001; | |
185 | } | |
186 | sigmacomputed = true; | |
187 | // std::cout << "DSDXS "<<sigma<< " " << mres<<std::endl; | |
188 | } | |
189 | ||
190 | // DETERMINE XS FOR Ds*Ds* | |
191 | if (((p->getDaug(0)->getDaug(0)->getId()==DSPS && p->getDaug(0)->getDaug(1)->getId()==DSMS))) | |
192 | { | |
193 | if(mres>(2.112+2.112)){ | |
194 | sigma=0.4; | |
195 | } | |
196 | else { | |
197 | // sigma=0.4; | |
198 | // sigma = 0 surely below Ds*Ds* threshold? - ponyisi | |
199 | sigma=0.00001; | |
200 | } | |
201 | sigmacomputed = true; | |
202 | // std::cout << "DsSDsSXS "<<sigma<< " " << mres<<std::endl; | |
203 | } | |
204 | ||
205 | // DETERMINE XS FOR Ds*Ds | |
206 | if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==DSPS | |
207 | && p->getDaug(0)->getDaug(1)->getId()==DSM) | |
208 | || (p->getDaug(0)->getDaug(0)->getId()==DSMS | |
209 | && p->getDaug(0)->getDaug(1)->getId()==DSP))) | |
210 | { | |
211 | if(mres>4.26){ | |
212 | sigma=0.05; | |
213 | } | |
214 | else if (mres>4.18 && mres<=4.26){ | |
215 | sigma*=1./9.*(0.05+0.95*sqrt((4.26-mres)*(4.26-mres))/(4.26-4.18)); | |
216 | } | |
217 | else if (mres>4.16 && mres<=4.18){ | |
218 | sigma*=1/9.; | |
219 | } | |
220 | else if (mres<=4.16 && mres>4.08){ | |
221 | sigma*=1/9.*(1-sqrt((4.16-mres)*(4.16-mres))/(4.16-4.08)); | |
222 | } | |
223 | else if (mres<=(4.08)){ | |
224 | sigma=0.00001; | |
225 | } | |
226 | sigmacomputed = true; | |
227 | // std::cout << "DsSDsXS "<<sigma<< " " << mres<<std::endl; | |
228 | } | |
229 | ||
230 | // DETERMINE XS FOR DD | |
231 | if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==D0 | |
232 | && p->getDaug(0)->getDaug(1)->getId()==D0B) | |
233 | ||(p->getDaug(0)->getDaug(0)->getId()==DP | |
234 | && p->getDaug(0)->getDaug(1)->getId()==DM))){ | |
235 | sigma*=0.4/9.; | |
236 | sigmacomputed = true; | |
237 | // std::cout << "DDXS "<<sigma<< " " << mres<<std::endl; | |
238 | } | |
239 | ||
240 | // DETERMINE XS FOR DsDs | |
241 | if (p->getDaug(0)->getNDaug() == 2 && ((p->getDaug(0)->getDaug(0)->getId()==DSP && p->getDaug(0)->getDaug(1)->getId()==DSM))){ | |
242 | sigma*=0.2/9.; | |
243 | sigmacomputed = true; | |
244 | // std::cout << "DsDsXS "<<sigma<< " " << mres<<std::endl; | |
245 | } | |
246 | ||
247 | // DETERMINE XS FOR MULTIBODY | |
248 | if (p->getDaug(0)->getNDaug() == 3){ | |
249 | if(mres>4.03){ | |
250 | sigma*=0.5/9.; | |
251 | } | |
252 | else { | |
253 | sigma=0.00001; | |
254 | } | |
255 | sigmacomputed = true; | |
256 | // std::cout << "DSDpiXS "<<sigma<< " " << mres<<std::endl; | |
257 | } | |
258 | ||
259 | if (! sigmacomputed) { | |
260 | report(ERROR,"EvtGen") << "VPHOTOVISRHI: This model requires daughters to be listed in a particular order." << endl | |
261 | << "The following are acceptable:" << endl | |
262 | << "D0 anti-D0" << endl | |
263 | << "D+ D-" << endl | |
264 | << "D*0 anti-D0" << endl | |
265 | << "anti-D*0 D0" << endl | |
266 | << "D*+ D-" << endl | |
267 | << "D*- D+" << endl | |
268 | << "D*0 anti-D*0" << endl | |
269 | << "D*+ D*-" << endl | |
270 | << "D_s+ D_s-" << endl | |
271 | << "D_s*+ D_s-" << endl | |
272 | << "D_s*- D_s+" << endl | |
273 | << "D_s*+ D_s*-" << endl | |
274 | << "(D* D pi can be in any order)" << endl | |
275 | << "Aborting..." << endl; | |
276 | assert(0); | |
277 | } | |
278 | ||
279 | if(sigma<0) sigma = 0.0; | |
280 | ||
281 | // static double sigmax=sigma; | |
282 | // if (sigma>sigmax){ | |
283 | // sigmax=sigma; | |
284 | // } | |
285 | ||
286 | static int count=0; | |
287 | ||
288 | count++; | |
289 | ||
290 | // if (count%10000==0){ | |
291 | // std::cout << "sigma :"<<sigma<<std::endl; | |
292 | // std::cout << "sigmax:"<<sigmax<<std::endl; | |
293 | // } | |
294 | ||
295 | double norm=sqrt(sigma); | |
296 | ||
297 | // EvtParticle* d=p->getDaug(0); | |
298 | ||
299 | ||
300 | vertex(0,0,0,norm*p->eps(0)*p->epsParent(0).conj()); | |
301 | vertex(1,0,0,norm*p->eps(1)*p->epsParent(0).conj()); | |
302 | vertex(2,0,0,norm*p->eps(2)*p->epsParent(0).conj()); | |
303 | ||
304 | vertex(0,1,0,norm*p->eps(0)*p->epsParent(1).conj()); | |
305 | vertex(1,1,0,norm*p->eps(1)*p->epsParent(1).conj()); | |
306 | vertex(2,1,0,norm*p->eps(2)*p->epsParent(1).conj()); | |
307 | ||
308 | vertex(0,2,0,norm*p->eps(0)*p->epsParent(2).conj()); | |
309 | vertex(1,2,0,norm*p->eps(1)*p->epsParent(2).conj()); | |
310 | vertex(2,2,0,norm*p->eps(2)*p->epsParent(2).conj()); | |
311 | ||
312 | vertex(0,0,1,norm*p->eps(0)*p->epsParent(0).conj()); | |
313 | vertex(1,0,1,norm*p->eps(1)*p->epsParent(0).conj()); | |
314 | vertex(2,0,1,norm*p->eps(2)*p->epsParent(0).conj()); | |
315 | ||
316 | vertex(0,1,1,norm*p->eps(0)*p->epsParent(1).conj()); | |
317 | vertex(1,1,1,norm*p->eps(1)*p->epsParent(1).conj()); | |
318 | vertex(2,1,1,norm*p->eps(2)*p->epsParent(1).conj()); | |
319 | ||
320 | vertex(0,2,1,norm*p->eps(0)*p->epsParent(2).conj()); | |
321 | vertex(1,2,1,norm*p->eps(1)*p->epsParent(2).conj()); | |
322 | vertex(2,2,1,norm*p->eps(2)*p->epsParent(2).conj()); | |
323 | ||
324 | return; | |
325 | } |