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