CGEM BOSS 6.6.5.g
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>
30
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
87 if(nh == 2){
88 // std::cout<<"final states: "<<Vson[0]<<" "<<Vson[1]<<std::endl;
89 son0 = Vson[0];son1 = Vson[1];
90 if(son0 == "D0" && son1 == "anti-D0" || son1 == "D0" && son0 == "anti-D0") {return 0;} else
91 if(son0 == "D+" && son1 == "D-" || son1 == "D+" && son0 == "D-" ) {return 1;} else
92 if(son0 == "D0" && son1 == "anti-D*0" || son1 == "D0" && son0 == "anti-D*0") {return 2;} else
93 if(son0 == "anti-D0" && son1 == "D*0" || son1 == "anti-D0" && son0 == "D*0") {return 3;} else
94 if(son0 == "D*0" && son1 == "anti-D*0" || son1 == "D*0" && son0 == "anti-D*0") {return 4;} else
95 if(son0 == "D*+" && son1 == "D-" || son1 == "D*+" && son0 == "D-") {return 5;} else
96 if(son0 == "D*-" && son1 == "D+" || son1 == "D*-" && son0 == "D+") {return 6;} else
97 if(son0 == "D*+" && son1 == "D*-" || son1 == "D*+" && son0 == "D*-") {return 7;} else
98 if(son0 == "D_s+" && son1 == "D_s-" || son1 == "D_s+" && son0 == "D_s-") {return 8;} else
99 if(son0 == "D_s*+" && son1 == "D_s-" || son1 == "D_s*+" && son0 == "D_s-") {return 9;} else
100 if(son0 == "D_s*-" && son1 == "D_s+" || son1 == "D_s*-" && son0 == "D_s+") {return 10;}else
101 if(son0 == "D_s*+" && son1 == "D_s*-" || son1 == "D_s*+" && son0 == "D_s*-") {return 11;}
102 } else if(nh == 3){
103 son0 = Vson[0];son1 = Vson[1];son2 = Vson[2];
104 if(son0 == "D*+" && son1 == "D-" && son2 == "pi0" ) {return 12;} else
105 if(son0 == "D*-" && son1 == "D+" && son2 == "pi0" ) {return 13;} else
106 if(son0 == "D*+" && son1 == "anti-D0" && son2 == "pi-" ) {return 14;} else
107 if(son0 == "D*-" && son1 == "D0" && son2 == "pi+" ) {return 15;} else
108 if(son0 == "D+" && son1 == "anti-D*0" && son2 == "pi-" ) {return 16;} else
109 if(son0 == "D-" && son1 == "D*0" && son2 == "pi+" ) {return 17;} else
110 if(son0 == "D*+" && son1 == "D*-" && son2 == "pi0" ) {return 18;} else
111 if(son0 == "anti-D*0" &&son1 == "D*+" && son2 == "pi-" ) {return 19;} else
112 if(son0 == "D*0" && son1 == "D*-" && son2 == "pi+" ) {return 20;} else
113 if(son0 == "D*0" && son1 == "anti-D*0" && son2 == "pi0" ) {return 21;} else
114 if(son0 == "D0" && son1 == "anti-D*0" && son2 == "pi0" ) {return 22;} else
115 if(son0 == "anti-D0" && son1 == "D*0" && son2 == "pi0" ) {return 23;}
116 }else{
117 std::cout<<"Not open charm decay"<<std::endl;
118 for(int j=0;j<nh;j++){
119 std::cout<<"final states, J= "<<j<<" "<<Vson[j]<<std::endl;
120 }
121 ::abort();
122 }
123
124}
125
126
127 bool EvtPsi3Sdecay::choseDecay(){ //determing accept or reject a generated decay
128
129 // findPoints();
130 double d0d0bar_xs=polint(d0d0bar);
131 double dpdm_xs = polint(dpdm);
132 double d0dst0bar_xs = polint(d0dst0bar);
133 double d0bardst0_xs = polint(d0bardst0);
134
135 double dst0dst0bar_xs = polint(dst0dst0bar);
136 double dstpdm_xs = polint(dstpdm);
137
138 double dstmdp_xs = polint(dstmdp);
139 double dstpdstm_xs = polint(dstpdstm);
140
141 double dspdsm_xs = polint(dspdsm);
142
143 double dsspdsm_xs = polint(dsspdsm);
144 double dssmdsp_xs = polint(dssmdsp);
145
146 double dsspdssm_xs = polint(dsspdssm);
147 //--- DDpi modes
148 double _xs12 = polint(xs12);
149 double _xs13 = polint(xs13);
150 double _xs14 = polint(xs14);
151 double _xs15 = polint(xs15);
152 double _xs16 = polint(xs16);
153 double _xs17 = polint(xs17);
154 double _xs18 = polint(xs18);
155 double _xs19 = polint(xs19);
156 double _xs20 = polint(xs20);
157 double _xs21 = polint(xs21);
158 double _xs22 = polint(xs22);
159 double _xs23 = polint(xs23);
160
161 int ich = findMode();
162 // std::cout<<"calculated XS "<< d0d0bar_xs<<" "<<dpdm_xs<<" "<<d0dst0bar_xs<<" "<<d0bardst0_xs<< " "<<dst0dst0bar_xs<<" "<<dstpdm_xs<< " "<<dstmdp_xs<<" "<<dstpdstm_xs<<" "<<dspdsm_xs<<std::endl;
163
164 double xmtotal=0;
165 for(int i=0;i<Vid.size();i++){
166 xmtotal += EvtPDL::getMinMass(Vid[i]);
167 }
168 double mparent= theParent->getP4().mass();
169 // std::cout<<"mparent= "<<mparent<<", xmtotal= "<<xmtotal<<endl;
170 if (mparent<xmtotal){return false;}
171
172
173 std::vector<double> myxs; myxs.clear();
174 myxs.push_back(d0d0bar_xs); //0
175 myxs.push_back(dpdm_xs);
176 myxs.push_back(d0dst0bar_xs); //2
177 myxs.push_back(d0bardst0_xs);
178 myxs.push_back(dst0dst0bar_xs); //4
179 myxs.push_back(dstpdm_xs);
180 myxs.push_back(dstmdp_xs); //6
181 myxs.push_back(dstpdstm_xs);
182 myxs.push_back(dspdsm_xs); //8
183 myxs.push_back(dsspdsm_xs);
184 myxs.push_back(dssmdsp_xs); //10
185 myxs.push_back(dsspdssm_xs); //11
186
187 myxs.push_back(_xs12); //12
188 myxs.push_back(_xs13); //13
189 myxs.push_back(_xs14); //14
190 myxs.push_back(_xs15); //15
191 myxs.push_back(_xs16); //16
192 myxs.push_back(_xs17); //17
193 myxs.push_back(_xs18); //18
194 myxs.push_back(_xs19); //19
195 myxs.push_back(_xs20); //20
196 myxs.push_back(_xs21); //18
197 myxs.push_back(_xs22); //19
198 myxs.push_back(_xs23); //20
199
200 double Prop0,Prop1;
201 if(ich==0){ Prop0=0;} else
202 {
203 Prop0 = theProb(myxs,ich-1);
204 }
205 Prop1 = theProb(myxs,ich);
206
207 double pm= EvtRandom::Flat(0.,1);
208 bool flag = false;
209 if( Prop0 < pm && pm<= Prop1 ) flag = true;
210
211 //--- debuging
212
213 if(flag) {
214 // std::cout<<"findMode= "<<ich<<std::endl;
215 // for(int i=0;i<myxs.size();i++){ std::cout<<"channel "<<i<<" myxs: "<<myxs[i]<<std::endl;}
216 //std::cout<<"prop0,prop1= "<<Prop0<<" "<<Prop1<<std::endl;
217 }
218
219 //-------------
220
221 return flag;
222}
223//---
224
225EvtParticle* EvtPsi3Sdecay::choseDecay(EvtParticle * par){ //decay par
226 double xm = par->mass();
227 int themode = getDecay(xm);
228 std::vector< EvtId > theid = getVId(themode);
229 int ndaugjs = theid.size();
230 EvtId myId[3];
231 for(int i=0;i<ndaugjs;i++){myId[i]=theid[i];}
232 par->makeDaughters(ndaugjs,myId);
233
234 for(int i=0;i<par->getNDaug();i++){
235 EvtParticle* di=par->getDaug(i);
236 double xmi=EvtPDL::getMinMass(di->getId());
237 di->setMass(xmi);
238 }
239 par->initializePhaseSpace(ndaugjs,myId);
240 _themode = themode;
241 return par;
242}
243
244
245//--
246
247void EvtPsi3Sdecay::PHSPDecay(EvtParticle * par){//decay par
248 v_p4.clear();Vid.clear();
249 double xm = par->mass();
250 EvtVector4R p_ini(xm,0,0,0);
251 EvtParticle* mypar= EvtParticleFactory::particleFactory(par->getId(),p_ini);
252
253 int themode = getDecay(xm);
254 std::vector< EvtId > theid = getVId(themode);
255 int ndaugjs = theid.size();
256 EvtId myId[3];
257 for(int i=0;i<ndaugjs;i++){myId[i]=theid[i];}
258 mypar->makeDaughters(ndaugjs,myId);
259
260 for(int i=0;i<mypar->getNDaug();i++){
261 EvtParticle* di=mypar->getDaug(i);
262 double xmi=EvtPDL::getMinMass(di->getId());
263 di->setMass(xmi);
264 }
265 loop:
266 mypar->initializePhaseSpace(ndaugjs,myId);
267
268 //-- here to do amgular distribution sampling
269 bool pp = (themode == 0||themode == 1 ||themode ==6); //V->PP mode, alpha = -1
270 bool vp = (themode >=2 && themode <=5 || themode >=7 && themode <=9 ); // V->VP mode, alpha = 1
271 bool flag1 = false;
272 double alpha;
273 if(pp){alpha=-1;}else if(vp){alpha=1;} else {alpha=0;}
274 EvtVector4R pp4=par->getP4();
275 EvtVector4R sp4=mypar->getDaug(1)->getP4();
276 flag1=AngSam(pp4,sp4,alpha);
277 if(alpha != 0 && !flag1 ) goto loop;
278 //--
279 _themode = themode;
280 for(int i=0;i<mypar->getNDaug();i++){
281 EvtParticle* di=mypar->getDaug(i);
282 v_p4.push_back( di->getP4() );
283 Vid.push_back(myId[i]);
284 // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
285 }
286
287
288 /*********** same function can be realized
289 double mass[3];
290 EvtVector4R p4[30];
291 for(int i=0;i<ndaugjs;i++){mass[i]=mypar->getDaug(i)->mass();}
292 EvtGenKine::PhaseSpace(ndaugjs,mass,p4,xm);
293 _themode = themode;
294 for(int i=0;i<mypar->getNDaug();i++){
295 v_p4.push_back( p4[i] );
296 Vid.push_back(myId[i]);
297 // std::cout<<"Vid "<<EvtPDL::name(Vid[i])<<" p4: "<<v_p4[i]<<std::endl;
298 }
299 *************/
300 mypar->deleteTree();
301 return;
302}
303
304
305//--
306
307 int EvtPsi3Sdecay::getDecay(double ecms){ //pick up a decay by the accept-reject sampling method
308 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; }
309 if(_excflag ==1) return _themode;
310 Ecms = ecms;
311 // findPoints();
312 double d0d0bar_xs=polint(d0d0bar);
313 double dpdm_xs = polint(dpdm);
314 double d0dst0bar_xs = polint(d0dst0bar);
315 double d0bardst0_xs = polint(d0bardst0);
316
317 double dst0dst0bar_xs = polint(dst0dst0bar);
318 double dstpdm_xs = polint(dstpdm);
319
320 double dstmdp_xs = polint(dstmdp);
321 double dstpdstm_xs = polint(dstpdstm);
322
323 double dspdsm_xs = polint(dspdsm);
324
325 double dsspdsm_xs = polint(dsspdsm);
326 double dssmdsp_xs = polint(dssmdsp);
327
328 double dsspdssm_xs = polint(dsspdssm);
329
330 double _xs12 = polint(xs12);
331 double _xs13 = polint(xs13);
332 double _xs14 = polint(xs14);
333 double _xs15 = polint(xs15);
334 double _xs16 = polint(xs16);
335 double _xs17 = polint(xs17);
336 double _xs18 = polint(xs18);
337 double _xs19 = polint(xs19);
338 double _xs20 = polint(xs20);
339 double _xs21 = polint(xs21);
340 double _xs22 = polint(xs22);
341 double _xs23 = polint(xs23);
342
343
344 std::vector<double> myxs; myxs.clear();
345 myxs.push_back(d0d0bar_xs); //0
346 myxs.push_back(dpdm_xs); //1
347 myxs.push_back(d0dst0bar_xs); //2
348 myxs.push_back(d0bardst0_xs); //3
349 myxs.push_back(dst0dst0bar_xs);//4
350 myxs.push_back(dstpdm_xs); //5
351 myxs.push_back(dstmdp_xs); //6
352 myxs.push_back(dstpdstm_xs); //7
353 myxs.push_back(dspdsm_xs); //8
354 myxs.push_back(dsspdsm_xs); //9
355 myxs.push_back(dssmdsp_xs); //10
356 myxs.push_back(dsspdssm_xs); //11
357 myxs.push_back(_xs12); //12
358 myxs.push_back(_xs13); //13
359 myxs.push_back(_xs14); //14
360 myxs.push_back(_xs15); //15
361 myxs.push_back(_xs16); //16
362 myxs.push_back(_xs17); //17
363 myxs.push_back(_xs18); //18
364 myxs.push_back(_xs19); //19
365 myxs.push_back(_xs20); //20
366 myxs.push_back(_xs21); //21
367 myxs.push_back(_xs22); //22
368 myxs.push_back(_xs23); //23
369
370 int niter = 0;
371 loop:
372 int ich = (int)24*EvtRandom::Flat(0.,1);//sampling modes over 24 channel
373
374 niter++;
375 if(niter>1000) {std::cout<<"EvtPsi3Sdecay:getDecay() do loops over 1000 times at Ecm= "<<ecms<<std::endl;abort();}
376
377 double xmtotal=0;
378 std::vector<EvtId> theVid;
379 theVid.clear();
380 theVid = getVId(ich);
381 for(int i=0;i<theVid.size();i++){
382 xmtotal += EvtPDL::getMinMass(theVid[i]);
383 }
384 if(Ecms < xmtotal ) {goto loop;}
385
386 double Prop0,Prop1;
387 if(ich==0){ Prop0=0;} else
388 {
389 Prop0 = theProb(myxs,ich-1);
390 }
391 Prop1 = theProb(myxs,ich);
392
393 double pm= EvtRandom::Flat(0.,1);
394
395 if( Prop0 < pm && pm<= Prop1 ) {return ich;}
396 else {goto loop;}
397
398}
399
400//----
401std::vector<EvtId> EvtPsi3Sdecay::getVId(int mode){
402 std::vector<EvtId> theVid;
403 theVid.clear();
404 for(int i=0;i<VmodeSon[mode].size();i++){
405 EvtId theId = EvtPDL::getId(VmodeSon[mode][i]);
406 theVid.push_back(theId);
407 }
408 return theVid;
409}
410
411
412bool EvtPsi3Sdecay::AngSam(EvtVector4R parent_p4cm,EvtVector4R son_p4cm,double alpha){
413 EvtHelSys angles(parent_p4cm,son_p4cm);
414 double theta=angles.getHelAng(1);
415 //double phi =angles.getHelAng(2);
416 //double gamma=0;
417 double costheta=cos(theta); //using helicity angles in parent system
418 double max_alpha;
419 if(alpha>=0) {max_alpha = 1+alpha;}else
420 {max_alpha=1;}
421 double ags = (1+alpha*costheta*costheta)/max_alpha;
422 double rand=EvtRandom::Flat(0.0, 1.0);
423 if(rand <=ags) {return true;}
424 else {return false;}
425}
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:358
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1175
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:1069
EvtParticle::deleteTree
void deleteTree()
Definition: EvtParticle.cc:555
EvtPsi3Sdecay::getDecay
int getDecay(double ecms)
Definition: EvtPsi3Sdecay.cc:307
EvtPsi3Sdecay::theProb
double theProb(std::vector< double > myxs, int ich)
Definition: EvtPsi3Sdecay.cc:55
EvtPsi3Sdecay::PHSPDecay
void PHSPDecay(EvtParticle *par)
Definition: EvtPsi3Sdecay.cc:247
EvtPsi3Sdecay::AngSam
bool AngSam(EvtVector4R parent_p4cm, EvtVector4R son_p4cm, double alpha)
Definition: EvtPsi3Sdecay.cc:412
EvtPsi3Sdecay::getVId
std::vector< EvtId > getVId(int mode)
Definition: EvtPsi3Sdecay.cc:401
EvtPsi3Sdecay::polint
double polint(std::vector< double > points)
Definition: EvtPsi3Sdecay.cc:42
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:73
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cc:39
ecms
const double ecms
Definition: inclkstar.cxx:42