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