90 std::string son0, son1, son2;
93 for (
int i = 0; i < Ndaugs; i++ )
95 std::string nson =
EvtPDL::name( theParent->getDaug( i )->getId() );
96 if ( nson !=
"gammaFSR" && nson !=
"gamma" )
98 Vson.push_back( nson );
99 Vid.push_back( theParent->getDaug( i )->getId() );
102 int nh = Vson.size();
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" )
116 else if ( son0 ==
"anti-D0" && son1 ==
"D*0" || son1 ==
"anti-D0" && son0 ==
"D*0" )
118 else if ( son0 ==
"D*0" && son1 ==
"anti-D*0" || son1 ==
"D*0" && son0 ==
"anti-D*0" )
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-" )
125 else if ( son0 ==
"D_s*+" && son1 ==
"D_s-" || son1 ==
"D_s*+" && son0 ==
"D_s-" )
127 else if ( son0 ==
"D_s*-" && son1 ==
"D_s+" || son1 ==
"D_s*-" && son0 ==
"D_s+" )
129 else if ( son0 ==
"D_s*+" && son1 ==
"D_s*-" || son1 ==
"D_s*+" && son0 ==
"D_s*-" )
131 else {
goto ErrInfo; }
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; }
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"
166 double d0d0bar_xs =
polint( d0d0bar );
167 double dpdm_xs =
polint( dpdm );
168 double d0dst0bar_xs =
polint( d0dst0bar );
169 double d0bardst0_xs =
polint( d0bardst0 );
171 double dst0dst0bar_xs =
polint( dst0dst0bar );
172 double dstpdm_xs =
polint( dstpdm );
174 double dstmdp_xs =
polint( dstmdp );
175 double dstpdstm_xs =
polint( dstpdstm );
177 double dspdsm_xs =
polint( dspdsm );
179 double dsspdsm_xs =
polint( dsspdsm );
180 double dssmdsp_xs =
polint( dssmdsp );
182 double dsspdssm_xs =
polint( dsspdssm );
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 );
204 double mparent = theParent->getP4().mass();
206 if ( mparent < xmtotal ) {
return false; }
208 std::vector<double> myxs;
210 myxs.push_back( d0d0bar_xs );
211 myxs.push_back( dpdm_xs );
212 myxs.push_back( d0dst0bar_xs );
213 myxs.push_back( d0bardst0_xs );
214 myxs.push_back( dst0dst0bar_xs );
215 myxs.push_back( dstpdm_xs );
216 myxs.push_back( dstmdp_xs );
217 myxs.push_back( dstpdstm_xs );
218 myxs.push_back( dspdsm_xs );
219 myxs.push_back( dsspdsm_xs );
220 myxs.push_back( dssmdsp_xs );
221 myxs.push_back( dsspdssm_xs );
223 myxs.push_back( _xs12 );
224 myxs.push_back( _xs13 );
225 myxs.push_back( _xs14 );
226 myxs.push_back( _xs15 );
227 myxs.push_back( _xs16 );
228 myxs.push_back( _xs17 );
229 myxs.push_back( _xs18 );
230 myxs.push_back( _xs19 );
231 myxs.push_back( _xs20 );
232 myxs.push_back( _xs21 );
233 myxs.push_back( _xs22 );
234 myxs.push_back( _xs23 );
237 if ( ich == 0 ) { Prop0 = 0; }
238 else { Prop0 =
theProb( myxs, ich - 1 ); }
243 if ( Prop0 < pm && pm <= Prop1 )
flag =
true;
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;
354 if ( _excflag == 1 )
return _themode;
357 double d0d0bar_xs =
polint( d0d0bar );
358 double dpdm_xs =
polint( dpdm );
359 double d0dst0bar_xs =
polint( d0dst0bar );
360 double d0bardst0_xs =
polint( d0bardst0 );
362 double dst0dst0bar_xs =
polint( dst0dst0bar );
363 double dstpdm_xs =
polint( dstpdm );
365 double dstmdp_xs =
polint( dstmdp );
366 double dstpdstm_xs =
polint( dstpdstm );
368 double dspdsm_xs =
polint( dspdsm );
370 double dsspdsm_xs =
polint( dsspdsm );
371 double dssmdsp_xs =
polint( dssmdsp );
373 double dsspdssm_xs =
polint( dsspdssm );
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 );
388 std::vector<double> myxs;
390 myxs.push_back( d0d0bar_xs );
391 myxs.push_back( dpdm_xs );
392 myxs.push_back( d0dst0bar_xs );
393 myxs.push_back( d0bardst0_xs );
394 myxs.push_back( dst0dst0bar_xs );
395 myxs.push_back( dstpdm_xs );
396 myxs.push_back( dstmdp_xs );
397 myxs.push_back( dstpdstm_xs );
398 myxs.push_back( dspdsm_xs );
399 myxs.push_back( dsspdsm_xs );
400 myxs.push_back( dssmdsp_xs );
401 myxs.push_back( dsspdssm_xs );
402 myxs.push_back( _xs12 );
403 myxs.push_back( _xs13 );
404 myxs.push_back( _xs14 );
405 myxs.push_back( _xs15 );
406 myxs.push_back( _xs16 );
407 myxs.push_back( _xs17 );
408 myxs.push_back( _xs18 );
409 myxs.push_back( _xs19 );
410 myxs.push_back( _xs20 );
411 myxs.push_back( _xs21 );
412 myxs.push_back( _xs22 );
413 myxs.push_back( _xs23 );
416 for (
int i = 0; i < myxs.size(); i++ ) { mytotxs += myxs[i]; }
417 if ( psi3Scount == 0 )
419 std::cout <<
"The total xs at " <<
ecms <<
" is: " << mytotxs <<
" pb" << std::endl;
429 std::cout <<
"EvtPsi3Sdecay:getDecay() do loops over 1000 times at Ecm= " <<
ecms
435 std::vector<EvtId> theVid;
438 for (
int i = 0; i < theVid.size(); i++ ) { xmtotal +=
EvtPDL::getMinMass( theVid[i] ); }
440 if ( Ecms < xmtotal ) {
goto loop; }
443 if ( ich == 0 ) { Prop0 = 0; }
444 else { Prop0 =
theProb( myxs, ich - 1 ); }
449 if ( Prop0 < pm && pm <= Prop1 ) {
return ich; }
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
void makeDaughters(int ndaug, EvtId *id)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)