BOSS 7.0.8
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPsi3Sdecay.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of models developed at BES collaboration
5// based on the EvtGen framework. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/BesCopyright
9// Copyright (A) 2006 Ping Rong-Gang, Pang Cai-Ying@IHEP
10//
11// Module: EvtPsi3Sdecay.cc
12//
13// Description: Routine to re-select the psi(4040) decay according the .
14// measured xsection at different energy point, see CLEOc measurement:
15// PRD 80, 072001
16// Modification history:
17//
18// Ping R.-G. December, 2010 Module created
19//
20//------------------------------------------------------------------------
21
22#include "EvtGenModels/EvtPsi3Sdecay.hh"
23#include "EvtGenBase/EvtRandom.hh"
24#include "EvtGenBase/EvtParticleFactory.hh"
25#include "EvtGenBase/EvtGenKine.hh"
26#include "EvtGenBase/EvtHelSys.hh"
27#include <string>
28#include <iostream>
29#include <cmath>
30int EvtPsi3Sdecay::psi3Scount=0;
31
32 int EvtPsi3Sdecay::findPoints(){
33 if(Ecms < x[0] ){theLocation =0;} else
34 if(Ecms >=x[nsize-1]) {theLocation = nsize-1;} else{
35 for (int i=0;i<nsize-1;i++){
36 if( x[i] <= Ecms && x[i+1] > Ecms) {theLocation=i;}
37 }
38 }
39 return theLocation;
40}
41
42double EvtPsi3Sdecay::polint(std::vector <double> vy ){
43
44 theLocation = findPoints();
45 double xs;
46 if(theLocation==nsize-1){xs = vy[nsize-1];} else{
47 xs = (vy[theLocation +1 ]- vy[theLocation]) / (x[theLocation+1]-x[theLocation])*(Ecms - x[theLocation])+vy[theLocation];
48 }
49 if(xs <0 ) xs = 0;
50 // for(int i=0;i<vy.size();i++){ std::cout<<"channel "<<i<<" xsection: "<<vy[i]<<std::endl;}
51 // std::cout<<"Ecms, theLocation= "<<Ecms<<" "<<theLocation<<" xs= "<<xs<<std::endl;
52 return xs;
53 }
54
55double EvtPsi3Sdecay::theProb(std::vector<double> myxs,int ich){
56 int Nchannels=myxs.size();
57 //---debuging
58 // std::cout<<"Nchannel= "<<Nchannels<<endl;
59
60 std::vector <double> thexs;
61 thexs.clear();
62 double sumxs=0;
63 for(int j=0;j<Nchannels;j++){
64 sumxs += myxs[j];
65 thexs.push_back(sumxs);
66 //----debugging
67 // std::cout<<"thexs["<<j<<"]= "<<myxs[j]<<std::endl;
68 }
69 if(sumxs == 0) {std::cout<<"EvtPsi3Sdecay::theProb, denominator is 0"<<std::endl;::abort();}
70 return thexs[ich] / sumxs ;
71}
72
73int EvtPsi3Sdecay::findMode(){
74 // mode index: 0) D0D0bar, 1)D+D-; 2)D0D*0bar , 3)D0bar D*0, 4)D*0 D*0 bar, 5)D*+D- 6)D*-D+ 7)D*+ D*- 8) Ds+ Ds-
75 // more modes: 9) D_s^*+ D_s^- ;10) D_s^*- D_s^+; 11) D_s^*+ D_s^*-
76
77 std::string son0,son1,son2;
78 Vson.clear();
79 Vid.clear();
80 for(int i=0;i<Ndaugs;i++){
81 std::string nson=EvtPDL::name(theParent->getDaug(i)->getId());
82 if(nson!="gammaFSR" && nson!="gamma"){ Vson.push_back(nson);Vid.push_back(theParent->getDaug(i)->getId());}
83 }
84 int nh=Vson.size();
85 //debugging
86 //std::cout<<"nh= "<<nh<<" "<<Vson[0]<<" "<<Vson[1]<<" "<<Vson[2]<<std::endl;
87 //theParent->printTree();
88
89 if(nh == 2){
90 //std::cout<<"2 final states: "<<Vson[0]<<" "<<Vson[1]<<std::endl;
91 son0 = Vson[0];son1 = Vson[1];
92 if(son0 == "D0" && son1 == "anti-D0" || son1 == "D0" && son0 == "anti-D0") {return 0;} else
93 if(son0 == "D+" && son1 == "D-" || son1 == "D+" && son0 == "D-" ) {return 1;} else
94 if(son0 == "D0" && son1 == "anti-D*0" || son1 == "D0" && son0 == "anti-D*0") {return 2;} else
95 if(son0 == "anti-D0" && son1 == "D*0" || son1 == "anti-D0" && son0 == "D*0") {return 3;} else
96 if(son0 == "D*0" && son1 == "anti-D*0" || son1 == "D*0" && son0 == "anti-D*0") {return 4;} else
97 if(son0 == "D*+" && son1 == "D-" || son1 == "D*+" && son0 == "D-") {return 5;} else
98 if(son0 == "D*-" && son1 == "D+" || son1 == "D*-" && son0 == "D+") {return 6;} else
99 if(son0 == "D*+" && son1 == "D*-" || son1 == "D*+" && son0 == "D*-") {return 7;} else
100 if(son0 == "D_s+" && son1 == "D_s-" || son1 == "D_s+" && son0 == "D_s-") {return 8;} else
101 if(son0 == "D_s*+" && son1 == "D_s-" || son1 == "D_s*+" && son0 == "D_s-") {return 9;} else
102 if(son0 == "D_s*-" && son1 == "D_s+" || son1 == "D_s*-" && son0 == "D_s+") {return 10;}else
103 if(son0 == "D_s*+" && son1 == "D_s*-" || son1 == "D_s*+" && son0 == "D_s*-") {return 11;}else {goto ErrInfo;}
104 } else if(nh == 3){
105 //std::cout<<"3 final states: "<<Vson[0]<<" "<<Vson[1]<<" "<<Vson[2]<<std::endl;
106 son0 = Vson[0];son1 = Vson[1];son2 = Vson[2];
107 if(son0 == "D*+" && son1 == "D-" && son2 == "pi0" ) {return 12;} else
108 if(son0 == "D*-" && son1 == "D+" && son2 == "pi0" ) {return 13;} else
109 if(son0 == "D*+" && son1 == "anti-D0" && son2 == "pi-" ) {return 14;} else
110 if(son0 == "D*-" && son1 == "D0" && son2 == "pi+" ) {return 15;} else
111 if(son0 == "D+" && son1 == "anti-D*0" && son2 == "pi-" ) {return 16;} else
112 if(son0 == "D-" && son1 == "D*0" && son2 == "pi+" ) {return 17;} else
113 if(son0 == "D*+" && son1 == "D*-" && son2 == "pi0" ) {return 18;} else
114 if(son0 == "anti-D*0" &&son1 == "D*+" && son2 == "pi-" ) {return 19;} else
115 if(son0 == "D*0" && son1 == "D*-" && son2 == "pi+" ) {return 20;} else
116 if(son0 == "D*0" && son1 == "anti-D*0" && son2 == "pi0" ) {return 21;} else
117 if(son0 == "D0" && son1 == "anti-D*0" && son2 == "pi0" ) {return 22;} else
118 if(son0 == "anti-D0" && son1 == "D*0" && son2 == "pi0" ) {return 23;} else {goto ErrInfo;}
119 }
120 ErrInfo:
121 std::cout<<"Not open charm decay"<<std::endl;
122 std::cout<<"final states \"";
123 for(int j=0;j<nh;j++){
124 std::cout<<Vson[j]<<" ";
125 }
126 std::cout<<" \" is not in the Psi3Decay list, see $BESEVTGENROOT/src/EvtGen/EvtGenModels/EvtPsi3Sdecay.hh"<<std::endl;
127 ::abort();
128
129
130}
131
132
133 bool EvtPsi3Sdecay::choseDecay(){ //determing accept or reject a generated decay
134
135 // findPoints();
136 double d0d0bar_xs=polint(d0d0bar);
137 double dpdm_xs = polint(dpdm);
138 double d0dst0bar_xs = polint(d0dst0bar);
139 double d0bardst0_xs = polint(d0bardst0);
140
141 double dst0dst0bar_xs = polint(dst0dst0bar);
142 double dstpdm_xs = polint(dstpdm);
143
144 double dstmdp_xs = polint(dstmdp);
145 double dstpdstm_xs = polint(dstpdstm);
146
147 double dspdsm_xs = polint(dspdsm);
148
149 double dsspdsm_xs = polint(dsspdsm);
150 double dssmdsp_xs = polint(dssmdsp);
151
152 double dsspdssm_xs = polint(dsspdssm);
153 //--- DDpi modes
154 double _xs12 = polint(xs12);
155 double _xs13 = polint(xs13);
156 double _xs14 = polint(xs14);
157 double _xs15 = polint(xs15);
158 double _xs16 = polint(xs16);
159 double _xs17 = polint(xs17);
160 double _xs18 = polint(xs18);
161 double _xs19 = polint(xs19);
162 double _xs20 = polint(xs20);
163 double _xs21 = polint(xs21);
164 double _xs22 = polint(xs22);
165 double _xs23 = polint(xs23);
166
167 int ich = findMode();
168 // std::cout<<"calculated XS "<< d0d0bar_xs<<" "<<dpdm_xs<<" "<<d0dst0bar_xs<<" "<<d0bardst0_xs<< " "<<dst0dst0bar_xs<<" "<<dstpdm_xs<< " "<<dstmdp_xs<<" "<<dstpdstm_xs<<" "<<dspdsm_xs<<std::endl;
169
170
171 double xmtotal=0;
172 for(int i=0;i<Vid.size();i++){
173 xmtotal += EvtPDL::getMinMass(Vid[i]);
174 }
175 double mparent= theParent->getP4().mass();
176 // std::cout<<"mparent= "<<mparent<<", xmtotal= "<<xmtotal<<endl;
177 if (mparent<xmtotal){return false;}
178
179
180 std::vector<double> myxs; myxs.clear();
181 myxs.push_back(d0d0bar_xs); //0
182 myxs.push_back(dpdm_xs);
183 myxs.push_back(d0dst0bar_xs); //2
184 myxs.push_back(d0bardst0_xs);
185 myxs.push_back(dst0dst0bar_xs); //4
186 myxs.push_back(dstpdm_xs);
187 myxs.push_back(dstmdp_xs); //6
188 myxs.push_back(dstpdstm_xs);
189 myxs.push_back(dspdsm_xs); //8
190 myxs.push_back(dsspdsm_xs);
191 myxs.push_back(dssmdsp_xs); //10
192 myxs.push_back(dsspdssm_xs); //11
193
194 myxs.push_back(_xs12); //12
195 myxs.push_back(_xs13); //13
196 myxs.push_back(_xs14); //14
197 myxs.push_back(_xs15); //15
198 myxs.push_back(_xs16); //16
199 myxs.push_back(_xs17); //17
200 myxs.push_back(_xs18); //18
201 myxs.push_back(_xs19); //19
202 myxs.push_back(_xs20); //20
203 myxs.push_back(_xs21); //21
204 myxs.push_back(_xs22); //22
205 myxs.push_back(_xs23); //23
206
207 double Prop0,Prop1;
208 if(ich==0){ Prop0=0;} else
209 {
210 Prop0 = theProb(myxs,ich-1);
211 }
212 Prop1 = theProb(myxs,ich);
213
214 double pm= EvtRandom::Flat(0.,1);
215 bool flag = false;
216 if( Prop0 < pm && pm<= Prop1 ) flag = true;
217
218 //--- debuging
219
220 if(flag) {
221 //std::cout<<"findMode= "<<ich<<std::endl;
222 //for(int i=0;i<myxs.size();i++){ std::cout<<"channel "<<i<<" myxs: "<<myxs[i]<<std::endl;}
223 //std::cout<<"prop0,prop1= "<<Prop0<<" "<<Prop1<<std::endl;
224 }
225
226 //-------------
227
228 return flag;
229}
230//---
231
232EvtParticle* EvtPsi3Sdecay::choseDecay(EvtParticle * par){ //decay par
233 double xm = par->mass();
234 int themode = getDecay(xm);
235 std::vector< EvtId > theid = getVId(themode);
236 int ndaugjs = theid.size();
237 EvtId myId[3];
238 for(int i=0;i<ndaugjs;i++){myId[i]=theid[i];}
239 par->makeDaughters(ndaugjs,myId);
240
241 for(int i=0;i<par->getNDaug();i++){
242 EvtParticle* di=par->getDaug(i);
243 double xmi=EvtPDL::getMinMass(di->getId());
244 di->setMass(xmi);
245 }
246 par->initializePhaseSpace(ndaugjs,myId);
247 _themode = themode;
248 return par;
249}
250
251
252//--
253
254void EvtPsi3Sdecay::PHSPDecay(EvtParticle * par){//decay par
255 v_p4.clear();Vid.clear();
256 double xm = par->mass();
257 EvtVector4R p_ini(xm,0,0,0);
258 EvtParticle* mypar= EvtParticleFactory::particleFactory(par->getId(),p_ini);
259
260 int themode = getDecay(xm);
261 std::vector< EvtId > theid = getVId(themode);
262 int ndaugjs = theid.size();
263 EvtId myId[3];
264 for(int i=0;i<ndaugjs;i++){myId[i]=theid[i];}
265 mypar->makeDaughters(ndaugjs,myId);
266
267 for(int i=0;i<mypar->getNDaug();i++){
268 EvtParticle* di=mypar->getDaug(i);
269 double xmi=EvtPDL::getMinMass(di->getId());
270 di->setMass(xmi);
271 }
272 loop:
273 mypar->initializePhaseSpace(ndaugjs,myId);
274
275 //-- here to do amgular distribution sampling
276 bool pp = (themode == 0||themode == 1 ||themode ==6); //V->PP mode, alpha = -1
277 bool vp = (themode >=2 && themode <=5 || themode >=7 && themode <=9 ); // V->VP mode, alpha = 1
278 bool flag1 = false;
279 double alpha;
280 if(pp){alpha=-1;}else if(vp){alpha=1;} else {alpha=0;}
281 EvtVector4R pp4=par->getP4();
282 EvtVector4R sp4=mypar->getDaug(1)->getP4();
283 flag1=AngSam(pp4,sp4,alpha);
284 if(alpha != 0 && !flag1 ) goto loop;
285 //--
286 _themode = themode;
287 for(int i=0;i<mypar->getNDaug();i++){
288 EvtParticle* di=mypar->getDaug(i);
289 v_p4.push_back( di->getP4() );
290 Vid.push_back(myId[i]);
291 // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
292 }
293
294
295 /*********** same function can be realized
296 double mass[3];
297 EvtVector4R p4[30];
298 for(int i=0;i<ndaugjs;i++){mass[i]=mypar->getDaug(i)->mass();}
299 EvtGenKine::PhaseSpace(ndaugjs,mass,p4,xm);
300 _themode = themode;
301 for(int i=0;i<mypar->getNDaug();i++){
302 v_p4.push_back( p4[i] );
303 Vid.push_back(myId[i]);
304 // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
305 }
306 *************/
307 mypar->deleteTree();
308 return;
309}
310
311
312//--
313
314 int EvtPsi3Sdecay::getDecay(double ecms){ //pick up a decay by the accept-reject sampling method
315 if(ecms<3.97 || ecms >4.66){std::cout<<"EvtPsi3Sdecay::getDecay: You need to set the CMS energy between 3.97~4.66 GeV, but you have ecms= "<<ecms<< " GeV.The lower end can be set at KKMC"<<std::endl; }
316 if(_excflag ==1) return _themode;
317 Ecms = ecms;
318 // findPoints();
319 double d0d0bar_xs=polint(d0d0bar);
320 double dpdm_xs = polint(dpdm);
321 double d0dst0bar_xs = polint(d0dst0bar);
322 double d0bardst0_xs = polint(d0bardst0);
323
324 double dst0dst0bar_xs = polint(dst0dst0bar);
325 double dstpdm_xs = polint(dstpdm);
326
327 double dstmdp_xs = polint(dstmdp);
328 double dstpdstm_xs = polint(dstpdstm);
329
330 double dspdsm_xs = polint(dspdsm);
331
332 double dsspdsm_xs = polint(dsspdsm);
333 double dssmdsp_xs = polint(dssmdsp);
334
335 double dsspdssm_xs = polint(dsspdssm);
336
337 double _xs12 = polint(xs12);
338 double _xs13 = polint(xs13);
339 double _xs14 = polint(xs14);
340 double _xs15 = polint(xs15);
341 double _xs16 = polint(xs16);
342 double _xs17 = polint(xs17);
343 double _xs18 = polint(xs18);
344 double _xs19 = polint(xs19);
345 double _xs20 = polint(xs20);
346 double _xs21 = polint(xs21);
347 double _xs22 = polint(xs22);
348 double _xs23 = polint(xs23);
349
350
351 std::vector<double> myxs; myxs.clear();
352 myxs.push_back(d0d0bar_xs); //0
353 myxs.push_back(dpdm_xs); //1
354 myxs.push_back(d0dst0bar_xs); //2
355 myxs.push_back(d0bardst0_xs); //3
356 myxs.push_back(dst0dst0bar_xs);//4
357 myxs.push_back(dstpdm_xs); //5
358 myxs.push_back(dstmdp_xs); //6
359 myxs.push_back(dstpdstm_xs); //7
360 myxs.push_back(dspdsm_xs); //8
361 myxs.push_back(dsspdsm_xs); //9
362 myxs.push_back(dssmdsp_xs); //10
363 myxs.push_back(dsspdssm_xs); //11
364 myxs.push_back(_xs12); //12
365 myxs.push_back(_xs13); //13
366 myxs.push_back(_xs14); //14
367 myxs.push_back(_xs15); //15
368 myxs.push_back(_xs16); //16
369 myxs.push_back(_xs17); //17
370 myxs.push_back(_xs18); //18
371 myxs.push_back(_xs19); //19
372 myxs.push_back(_xs20); //20
373 myxs.push_back(_xs21); //21
374 myxs.push_back(_xs22); //22
375 myxs.push_back(_xs23); //23
376
377 double mytotxs=0;
378 for(int i=0;i<myxs.size();i++){mytotxs += myxs[i];}
379 if(psi3Scount==0) {std::cout<<"The total xs at "<<ecms<<" is: "<<mytotxs<<" pb"<<std::endl;psi3Scount++;}
380 int niter = 0;
381 loop:
382 int ich = (int)24*EvtRandom::Flat(0.,1);//sampling modes over 24 channel
383
384 niter++;
385 if(niter>1000) {std::cout<<"EvtPsi3Sdecay:getDecay() do loops over 1000 times at Ecm= "<<ecms<<std::endl;abort();}
386
387 double xmtotal=0;
388 std::vector<EvtId> theVid;
389 theVid.clear();
390 theVid = getVId(ich);
391 for(int i=0;i<theVid.size();i++){
392 xmtotal += EvtPDL::getMinMass(theVid[i]);
393 }
394
395 if(Ecms < xmtotal ) {goto loop;}
396
397 double Prop0,Prop1;
398 if(ich==0){ Prop0=0;} else
399 {
400 Prop0 = theProb(myxs,ich-1);
401 }
402 Prop1 = theProb(myxs,ich);
403
404 double pm= EvtRandom::Flat(0.,1);
405
406 if( Prop0 < pm && pm<= Prop1 ) {return ich;}
407 else {goto loop;}
408
409}
410
411//----
412std::vector<EvtId> EvtPsi3Sdecay::getVId(int mode){
413 std::vector<EvtId> theVid;
414 theVid.clear();
415 for(int i=0;i<VmodeSon[mode].size();i++){
416 EvtId theId = EvtPDL::getId(VmodeSon[mode][i]);
417 theVid.push_back(theId);
418 }
419 return theVid;
420}
421
422
423bool EvtPsi3Sdecay::AngSam(EvtVector4R parent_p4cm,EvtVector4R son_p4cm,double alpha){
424 EvtHelSys angles(parent_p4cm,son_p4cm);
425 double theta=angles.getHelAng(1);
426 //double phi =angles.getHelAng(2);
427 //double gamma=0;
428 double costheta=cos(theta); //using helicity angles in parent system
429 double max_alpha;
430 if(alpha>=0) {max_alpha = 1+alpha;}else
431 {max_alpha=1;}
432 double ags = (1+alpha*costheta*costheta)/max_alpha;
433 double rand=EvtRandom::Flat(0.0, 1.0);
434 if(rand <=ags) {return true;}
435 else {return false;}
436}
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
alpha
const double alpha
Definition: FastVertexFit.cxx:4
EvtHelSys::getHelAng
double getHelAng(int i)
Definition: EvtHelSys.cc:54
EvtId
Definition: EvtId.hh:27
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition: EvtPDL.hh:51
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition: EvtParticleFactory.cc:46
EvtParticle::setMass
void setMass(double m)
Definition: EvtParticle.hh:372
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1177
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
EvtParticle::getNDaug
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:127
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1071
EvtParticle::deleteTree
void deleteTree()
Definition: EvtParticle.cc:557
EvtPsi3Sdecay::getDecay
int getDecay(double ecms)
Definition: EvtPsi3Sdecay.cc:314
EvtPsi3Sdecay::theProb
double theProb(std::vector< double > myxs, int ich)
Definition: EvtPsi3Sdecay.cc:55
EvtPsi3Sdecay::PHSPDecay
void PHSPDecay(EvtParticle *par)
Definition: EvtPsi3Sdecay.cc:254
EvtPsi3Sdecay::AngSam
bool AngSam(EvtVector4R parent_p4cm, EvtVector4R son_p4cm, double alpha)
Definition: EvtPsi3Sdecay.cc:423
EvtPsi3Sdecay::getVId
std::vector< EvtId > getVId(int mode)
Definition: EvtPsi3Sdecay.cc:412
EvtPsi3Sdecay::polint
double polint(std::vector< double > points)
Definition: EvtPsi3Sdecay.cc:42
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:74
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cc:39
ecms
const double ecms
Definition: inclkstar.cxx:42
costheta
float costheta
Definition: prepare_data.cxx:28