BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtConExc.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP, Wang Weiping @USTC
9//
10// Module: EvtXsection.hh
11//
12// Description: To define cross section for the continuum exclusive process
13// Experimental cross section taken from PRD73,012005, PRD76,092006, also
14// see a review: Rev. Mod. Phys. 83,1545
15// Modification history:
16//
17// Ping R.-G. Nov., 2012 Module created
18// Wang Weiping: Nov., 2017 added the mode==85: e+e- --> D_s*+D_s*-
19// Wang Weiping: Nov., 2017 fixed final states of mode==28,33,92
20// Wang Weiping: Nov., 2017 fixed the angle sampling of VP states
21//
22/*******************--- mode definition: also see EvtXsection.cc
23mode==0: ppbar
24mode==1: nnbar
25mode==2: Lambda0 anti-Lambda0
26mode==3: Sigma0 anti-Sigma0
27mode==4: Lambda0 anti-Sigma0
28mode==5: Sigma0 anti-Lambda0
29mode==6: pi+ pi-
30mode==7: pi+ pi- pi0
31mode==8: K+K- pi0
32mode==9: KsK+pi-
33mode==10: KsK-pi+
34mode==11: K+K-eta
35mode==12: 2(pi+pi-)
36mode==13: pi+pi-2pi0
37mode==14: K+K-pi+pi-
38mode==15: K+K-2pi0
39mode==16: 2(K+K-)
40mode==17: 2(pi+pi-)pi0
41mode==18: 2(pi+pi-)eta
42mode==19: K+K-pi+pi-pi0
43mode==20: K+K-pi+pi-eta
44mode==21: 3(pi+pi-)
45mode==22: 2(pi+pi-pi0)
46mode==23: phi eta
47mode==24: phi pi0
48mode==25: K+K*-
49mode==26: K-K*+
50mode==27: K_SK*0-bar
51mode==28: K*0(892)K+pi-
52mode==29: K*0(892)K-pi+
53mode==30: K*+K-pi0
54mode==31: K*-K+pi0
55mode==32: K_2*(1430)0 K+pi-
56mode==33: K_2*(1430)0 K-pi+
57mode==34: K+K-rho
58mode==35: phi pi+pi-
59mode==36: phi f0(980)
60mode==37: eta pi+pi-
61mode==38: omega pi+ pi-
62mode==39: omega f0(980)
63mode==40: eta' pi+ pi-
64mode==41: f_1(1285)pi+pi-
65mode==42: omega K+K-
66mode==43: omega pi+pi-pi0
67mode==44: Sigma+ Sigma- (Sigma0 Sigma0-bar SU(3) extention )
68mode==45: K+K-
69mode==46: K_S K_L
70mode==47: omega eta
71mode==48: phi eta'
72mode==49: phi K+K-
73mode==50: ppbar pi0
74mode==51: ppbar eta
75mode==52: omega pi0
76
77mode==70: D0 anti-D0
78mode==71: D+D-
79mode==72: D+D*-
80mode==73: D-D*+
81mode==74: D*+D*-
82mode==68: D0 anti-D*0
83mode==69: anti-D0 D*0
84mode==67: D*0 anti-D*0
85
86mode==59: D0 anti-D0 pi0
87mode==63: D+D-pi0
88mode==75: D0D-pi+
89mode==76: anti-D0 D+pi-
90
91mode==77: D0D*-pi+
92mode==78: anti-D0 D*+pi-
93mode==65: D-D*0pi+
94mode==66: D+ anti-D*0pi-
95mode==60: anti-D0 D*0pi0
96mode==61: D0 anti-D*0pi0
97mode==62: D+D*-pi0
98mode==64: D-D*+pi0
99
100mode==90: pi+pi- J/psi
101mode==92: pi0pi0 J/psi
102mode==91: pi+pi- psi(2S)
103mode==79: pi0pi0 psi(2S)
104mode==80: eta J/psi
105mode==81: pi+pi- h_c
106mode==82: pi0pi0 h_c
107mode==83: K+K- J/psi
108mode==84: KsKs J/psi
109
110mode==93: D_s+ D_s-
111mode==94: D_s^{*+}D_s^-
112mode==95: D_s^{*-}D_s^+
113mode==85: D_s^{*+}D_s^{*-}
114
115mode==96: Lambda_c+ anti-Lamba_c-
116mode==97: K- D_s+ anti-D0
117mode==111: e+e- ->Zc(3900)+ pi-
118mode==112: e+e- ->Zc(3900)- pi+
119//---
120For cross section, see Ref. Rev.Mod. Phys.83,1545
121*************************************/
122//------------------------------------------------------------------------
123#include "EvtConExc.hh"
124#include "../EvtGenBase/EvtDecayBase.hh"
125#include "../EvtGenBase/EvtDecayTable.hh"
126#include "../EvtGenBase/EvtGenKine.hh"
127#include "../EvtGenBase/EvtGlobalSet.hh"
128#include "../EvtGenBase/EvtHelSys.hh"
129#include "../EvtGenBase/EvtPDL.hh"
130#include "../EvtGenBase/EvtParticle.hh"
131#include "../EvtGenBase/EvtParticleFactory.hh"
132#include "../EvtGenBase/EvtPatches.hh"
133#include "../EvtGenBase/EvtRandom.hh"
134#include "../EvtGenBase/EvtReport.hh"
135#include "../EvtGenBase/EvtVector4R.hh"
136#include "EvtPhsp.hh"
137#include "TCanvas.h"
138#include "TGraphErrors.h"
139#include "TMultiGraph.h"
140#include "TPostScript.h"
141#include "TStyle.h"
142#include <fstream>
143#include <iomanip>
144#include <iostream>
145#include <istream>
146#include <math.h>
147#include <stdio.h>
148#include <stdlib.h>
149#include <string>
150#include <strstream>
151using namespace std;
152
153int EvtConExc::nconexcdecays = 0;
154EvtDecayBasePtr* EvtConExc::conexcdecays = 0;
155int EvtConExc::ntable = 0;
156
157int EvtConExc::ncommand = 0;
158int EvtConExc::lcommand = 0;
159std::string* EvtConExc::commands = 0;
160double EvtConExc::AF[600];
161double EvtConExc::AA[600];
162double EvtConExc::MH[600];
163
164double EvtConExc::mlow = 0;
165double EvtConExc::mup = 0;
166
167double EvtConExc::SetMthr = 0;
168
169extern "C" {
170extern double dgauss_( double ( *fun )( double* ), double*, double*, double* );
171}
172
173extern "C" {
174extern double divdif_( float*, float*, int*, float*, int* );
175}
176
177extern "C" {
178extern void polint_( float*, float*, int*, float*, float* );
179}
180
181extern "C" {
182extern void hadr5n12_( float*, float*, float*, float*, float*, float* );
183}
184
185double Rad2difXs( double* m );
186double Rad2difXs_er( double* m );
187
188EvtXsection* EvtConExc::myxsection;
189EvtXsection* EvtConExc::staxsection;
190double EvtConExc::_cms;
191double EvtConExc::XS_max;
192
193double EvtConExc::_xs0 = 0;
194double EvtConExc::_xs1 = 0;
195double EvtConExc::_er0 = 0;
196double EvtConExc::_er1 = 0;
197int EvtConExc::_nevt = 0;
198
199std::ofstream oa;
200int EvtConExc::getMode;
201int EvtConExc::conexcmode = -1;
202int EvtConExc::multimode = 0;
203std::vector<std::vector<double>> EvtConExc::VXS;
204
205EvtConExc::~EvtConExc() {
206 // if(_mode==74110)checkEvtRatio();
207 if ( myFisr.size() ) OutStaISR();
208 if ( mydbg )
209 {
210 // xs->Write();
211 myfile->Write();
212 }
213
214 delete myxsection;
215 delete staxsection;
216 gamH->deleteTree();
217
218 // the deletion of commands is really uggly!
219
220 int i;
221 if ( nconexcdecays == 0 )
222 {
223 delete[] commands;
224 commands = 0;
225 return;
226 }
227
228 for ( i = 0; i < nconexcdecays; i++ )
229 {
230 if ( conexcdecays[i] == this )
231 {
232 conexcdecays[i] = conexcdecays[nconexcdecays - 1];
233 nconexcdecays--;
234 if ( nconexcdecays == 0 )
235 {
236 delete[] commands;
237 commands = 0;
238 }
239 return;
240 }
241 }
242}
243
244void EvtConExc::getName( std::string& model_name ) { model_name = "ConExc"; }
245
246EvtDecayBase* EvtConExc::clone() { return new EvtConExc; }
247
248void EvtConExc::init() {
249 mydbg = false;
250 myFisr.clear();
251 ReadVP();
252 VXS.resize( 120 );
253 for ( int i = 0; i < 120; i++ ) { VXS[i].resize( 600 ); }
254 _mode = getArg( 0 );
255 for ( int i = 0; i < 97; i++ )
256 {
257 if ( 53 <= i && i <= 58 ) continue;
258 if ( 86 <= i && i <= 89 ) continue;
259 if ( _mode == 74110 || _mode == -100 ) vmode.push_back( i );
260 }
261 if ( _mode == 74110 || _mode == -100 )
262 {
263 vmode.push_back( 74100 );
264 vmode.push_back( 74101 );
265 vmode.push_back( 74102 );
266 vmode.push_back( 74103 );
267 vmode.push_back( 74104 );
268 vmode.push_back( 74105 );
269 vmode.push_back( 74110 );
270 }
271
272 std::cout << "==== Parameters set by users=====" << std::endl;
273 for ( int i = 0; i < ncommand; i++ ) { std::cout << commands[i] << std::endl; }
274 std::cout << "===================================" << std::endl;
275 InitPars();
276 std::cout << threshold << " " << beamEnergySpread << std::endl;
277 //----------------
278 // check that there are 0 arguments
279 // checkNArg(1);
280 // Vector ISR process
281 VISR = 0;
282 if ( getNDaug() == 2 )
283 {
284 if ( getDaugs()[0] == EvtPDL::getId( "gamma" ) &&
285 getDaugs()[1] == EvtPDL::getId( "gamma*" ) ||
286 getDaugs()[0] == EvtPDL::getId( "gamma*" ) &&
287 getDaugs()[1] == EvtPDL::getId( "gamma" ) )
288 VISR = 1;
289 }
290 else if ( getNDaug() > 2 )
291 {
292 std::cout << "Bad daughter specified" << std::endl;
293 abort();
294 }
295
296 // cmsspread=0.0015; //CMC energy spread
297 testflag = 0;
298 getResMass();
299 if ( getArg( 0 ) == -1 )
300 {
301 radflag = 0;
302 mydbg = false;
303 _mode = getArg( 0 );
304 pdgcode = getArg( 1 );
305 pid = EvtPDL::evtIdFromStdHep( pdgcode );
306 nson = getNArg() - 2;
307 std::cout << "The decay daughters are " << std::endl;
308 for ( int i = 2; i < getNArg(); i++ )
309 {
310 son[i - 2] = EvtPDL::evtIdFromStdHep( getArg( i ) );
311 std::cout << EvtPDL::name( son[i - 2] ) << " ";
312 }
313 std::cout << std::endl;
314 }
315 else if ( getArg( 0 ) == -2 )
316 {
317 radflag = 0;
318 mydbg = false;
319 _mode = getArg( 0 );
320 nson = getNArg() - 1;
321 for ( int i = 1; i < getNArg(); i++ )
322 { son[i - 1] = EvtPDL::evtIdFromStdHep( getArg( i ) ); }
323 }
324 else if ( getArg( 0 ) == -100 )
325 {
326 _mode = getArg( 0 );
327 multimode = 1;
328 vmd.clear();
329 for ( int i = 1; i < getNArg(); i++ ) { vmd.push_back( getArg( i ) ); }
330 std::cout << " multi-exclusive mode " << std::endl;
331 for ( int i = 1; i < getNArg(); i++ ) { std::cout << getArg( i ) << " "; }
332 std::cout << std::endl;
333 if ( vmd[vmd.size() - 1] == 74110 )
334 {
335 std::cout << "74110 is not allowd for multimode simulation" << std::endl;
336 abort();
337 }
338 }
339 else if ( getNArg() == 1 )
340 {
341 _mode = getArg( 0 );
342 radflag = 0;
343 mydbg = false;
344 }
345 else if ( getNArg() == 2 )
346 {
347 _mode = getArg( 0 );
348 radflag = getArg( 1 );
349 mydbg = false;
350 }
351 else if ( getNArg() == 3 )
352 {
353 _mode = getArg( 0 );
354 radflag = getArg( 1 );
355 mydbg = getArg( 2 );
356 }
357 else
358 {
359 std::cout << "ConExc:number of parameter should be 1,2 or 3, but you set " << getNArg()
360 << std::endl;
361 ::abort();
362 }
363 //--- debugging
364 // std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<"
365 // "<<EvtPDL::name(son[1])<<std::endl;
366
367 gamId = EvtPDL::getId( std::string( "gamma" ) );
368 init_mode( _mode );
369 XS_max = -1;
370 //-----debugging to make root file
371 mydbg = false; // designer need to remove this line
372 if ( mydbg )
373 {
374 myfile = new TFile( "mydbg.root", "RECREATE" );
375 xs = new TTree( "xs", "momentum of rad. gamma and hadrons" );
376 xs->Branch( "imode", &imode, "imode/I" );
377 xs->Branch( "pgam", pgam, "pgam[4]/D" );
378 xs->Branch( "phds", phds, "phds[4]/D" );
379 xs->Branch( "ph1", ph1, "ph1[4]/D" );
380 xs->Branch( "ph2", ph2, "ph2[4]/D" );
381 xs->Branch( "mhds", &mhds, "mhds/D" );
382 xs->Branch( "mass1", &mass1, "mass1/D" );
383 xs->Branch( "mass2", &mass2, "mass2/D" );
384 xs->Branch( "costheta", &costheta, "costheta/D" );
385 xs->Branch( "cosp", &cosp, "cosp/D" );
386 xs->Branch( "cosk", &cosk, "cosk/D" );
387 xs->Branch( "sumxs", &sumxs, "sumxs/D" );
388 xs->Branch( "selectmode", &selectmode, "selectmode/D" );
389 }
390 //--- prepare the output information
391 EvtId parentId = EvtPDL::getId( std::string( "vpho" ) );
392 EvtPDL::reSetWidth( parentId, 0.0 );
393 double parentMass = EvtPDL::getMass( parentId );
394 std::cout << "parent mass = " << parentMass << std::endl;
395
396 EvtVector4R p_init( parentMass, 0.0, 0.0, 0.0 );
397 EvtParticle* p = EvtParticleFactory::particleFactory( parentId, p_init );
398 theparent = p;
399 myxsection = new EvtXsection( _mode );
400 staxsection = new EvtXsection( _mode );
401 if ( _mode == -100 )
402 {
403 myxsection->setModes( vmd );
404 myxsection->ini_data_multimode();
405 staxsection->setModes( vmd );
406 staxsection->ini_data_multimode();
407 }
408 double mth = 0;
409 double mx = EvtPDL::getMeanMass( parentId ); // p->getP4().mass();
410 if ( _mode == -1 )
411 {
412 myxsection->setBW( pdgcode );
413 staxsection->setBW( pdgcode );
414 for ( int i = 0; i < nson; i++ ) { mth += EvtPDL::getMeanMass( son[i] ); }
415 if ( mx < mth )
416 {
417 std::cout << "The CMS energy is lower than the threshold of the final states"
418 << std::endl;
419 abort();
420 }
421 }
422 else if ( _mode == -2 ) { mth = myxsection->getXlw(); }
423 else { mth = myxsection->getXlw(); }
424 _cms = mx;
425 _unit = myxsection->getUnit();
426
427 std::cout << "The specified mode is " << _mode << std::endl;
428 EvtConExc::getMode = _mode;
429 // std::cout<<"getXlw= "<<mth<<std::endl;
430
431 Egamcut = 0.001; // set photon energy cut according to the BESIII detector
432 double Esig = 0.0001; // sigularity energy
433 double x = 2 * Egamcut / _cms;
434 double s = _cms * _cms;
435 double mhscut = sqrt( s * ( 1 - x ) );
436
437 // for vaccum polarization
438 float fe, fst2, fder, ferrder, fdeg, ferrdeg;
439 fe = _cms;
440 // using the vacuum pol. value as given by
441 // http://www-com.physik.hu-berlin.de/~fjeger/software.html
442 vph = getVP( _cms );
443 if ( 3.0943 < _cms && _cms < 3.102 )
444 vph = 1; // For J/psi, the vacuum factor is included in the resonance
445 std::cout << "sum of leptonic, hadronic contributions(u,d,s,c,b) to vacuum polarization: "
446 << vph << " @" << fe << "GeV" << std::endl;
447
448 // caculate the xs for ISR to psiprime, J/psi and phi narrow resonance.
449 double totxsIRS = 0;
450 init_Br_ee();
451 double mthrld = EvtPDL::getMeanMass( daugs[0] ); // threshold mass of hadrons
452 for ( int jj = 1; jj < _ndaugs; jj++ ) { mthrld += EvtPDL::getMeanMass( daugs[jj] ); }
453
454 ISRXS.clear();
455 ISRM.clear();
456 ISRFLAG.clear();
457 for ( int ii = 0; ii < 3; ii++ )
458 {
459 double mR = EvtPDL::getMeanMass( ResId[ii] );
460 if ( mx < mR || _mode != 74110 ) continue;
461 double narRxs = Ros_xs( mx, BR_ee[ii], ResId[ii] );
462 std::cout << "The corss section for gamma " << EvtPDL::name( ResId[ii] ).c_str()
463 << " is: " << narRxs << "*Br (nb)." << std::endl;
464 ISRXS.push_back( narRxs );
465 ISRM.push_back( mR );
466 ISRFLAG.push_back( 1. );
467 ISRID.push_back( ResId[ii] );
468 totxsIRS += narRxs;
469 }
470 std::cout << "==========================================================================="
471 << std::endl;
472
473 //-- calculated with MC integration method
474 double mhdL = myxsection->getXlw();
475 if ( mhdL < SetMthr ) mhdL = SetMthr;
476 std::cout << "SetMthr= " << SetMthr << std::endl;
477 EgamH = ( s - mhdL * mhdL ) / 2 / sqrt( s );
478 int nmc = 1000000;
479 _xs0 = gamHXSection( s, Esig, Egamcut, nmc );
480 _er0 = _er1; // stored when calculate _xs0
481 std::cout << "_er0= " << _er0 << std::endl;
482 _xs1 = gamHXSection( s, Egamcut, EgamH, nmc );
483 std::cout << "_xs1= " << _xs0 << std::endl;
484 double xs1_err = _er1;
485
486 //--- sigularity point x from 0 to 0.0001GeV
487 double xb = 2 * Esig / _cms;
488 double sig_xs = SoftPhoton_xs( s, xb ) * ( myxsection->getXsection( mx ) );
489 _xs0 += sig_xs;
490 std::cout << "_xs0= " << _xs0 << std::endl;
491
492 // prepare the observed cross section with interpotation function divdif in CERNLIB
493 testflag = 1;
494 int Nstp = 598;
495 double stp = ( EgamH - Egamcut ) / Nstp;
496 for ( int i = 0; i < Nstp; i++ )
497 { // points within Egam(max) ~ Egamcut=25MeV
498 double Eg0 = EgamH - i * stp;
499 double Eg1 = EgamH - ( i + 1 ) * stp;
500 int nmc = 20000;
501 int nint = 100;
502 // double xsi= gamHXSection(s,Eg1,Eg0,nmc);
503 double xsi = trapezoid( s, Eg1, Eg0, nint, myxsection );
504 AA[i] = ( Eg1 + Eg0 ) / 2;
505 double mhi = sqrt( _cms * _cms - 2 * _cms * AA[i] );
506 double mh2 = sqrt( _cms * _cms - 2 * _cms * Eg1 );
507 double binwidth = mh2 - mhi;
508 // if(_mode==74110) xsi += addNarrowRXS(mhi,binwidth);
509 if ( i == 0 ) AF[0] = xsi;
510 if ( i >= 1 ) AF[i] = AF[i - 1] + xsi;
511 RadXS[i] = xsi / stp;
512 }
513 // add the interval 0~Esig
514 AA[598] = Egamcut;
515 AA[599] = 0; // AA[i]: E_gamma
516 AF[598] = AF[597];
517 AF[599] = AF[598] + _xs0;
518 RadXS[599] = _xs0;
519
520 // prepare observed cross section for each mode
521 std::cout << "mhdL = " << mhdL << ", SetMthr= " << SetMthr << ", EgamH = " << EgamH
522 << std::endl;
523 for ( int i = 0; i < vmode.size(); i++ )
524 {
525 if ( _mode == 74110 || _mode == -100 ) mk_VXS( Esig, Egamcut, EgamH, i );
526 }
527 if ( _mode == 74110 || _mode == -100 ) writeDecayTabel();
528 // after mk_VXS, restore the myxsection to _mode
529 if ( _mode != -100 )
530 {
531 delete myxsection;
532 myxsection = new EvtXsection( _mode );
533 }
534 // debugging VXS
535 /*
536 double xtmp=0;
537 for(int i=0;i<vmode.size();i++){
538 xtmp+=VXS[i][599];
539 for(int j=0;j<600;j++){
540 std::cout<<VXS[i][j]<<" ";
541 }
542 std::cout<<std::endl;
543 }
544 */
545 // output observed cross section for the gamma Resonance mode
546 // std::cout<<"vmode.size======="<<vmode.size()<<std::endl;
547 for ( int i = 0; i < vmode.size(); i++ )
548 {
549 int md = vmode[i];
550 if ( md < 74100 || md > 74106 ) continue;
551 std::string partname = "";
552 if ( md == 74100 ) { partname = "J/psi"; }
553 else if ( md == 74101 ) { partname = "psi(2S)"; }
554 else if ( md == 74102 ) { partname = "psi(3770)"; }
555 else if ( md == 74103 ) { partname = "phi"; }
556 else if ( md == 74104 ) { partname = "omega"; }
557 else if ( md == 74105 ) { partname = "rho0"; }
558 else if ( md == 74106 ) { partname = "rho(3S)0"; }
559 else { ; }
560 std::cout << "The observed cross section for gamma " << partname << ": " << VXS[i][599]
561 << " nb" << std::endl;
562 }
563 //--
564
565 for ( int i = 0; i < 600; i++ )
566 {
567 MH[i] = sqrt( _cms * _cms - 2 * _cms * AA[i] );
568 // std::cout<<i<<" Egamma "<<AA[i]<<" Mhadons "<<MH[i]<<std::endl;
569 }
570 // for debugging
571 // for(int i=0;i<600;i++){double s=_cms*_cms; double
572 // mhds=sqrt(s-2*_cms*AA[i]);std::cout<<"Mhds="<<mhds<<" Egam="<<AA[i]<<"
573 // "<<AF[i]<<std::endl;} std::cout<<"VXS[79][599]: "<<VXS[79][599]<<" VXS[79][598]=
574 // "<<VXS[79][598]<<std::endl;
575 std::cout << "EgamH = " << EgamH << ", obsvXS = " << _xs0 + _xs1 << ", _xs1 = " << _xs1
576 << ", _xs0 = " << _xs0 << std::endl;
577 std::cout << "EgamH = " << EgamH << ", AF[599] = " << AF[599] << ", AF[598] = " << AF[598]
578 << ", _xs0 = " << _xs0 << std::endl;
579 // double Etst=0.0241838;
580 // std::cout<<"Etst="<<Etst<<" "<<gamHXSection(s,Etst,EgamH,nmc)<<" "<<
581 // lgr(AA,AF,600,Etst)<<std::endl; abort();
582
583 // for information output
584 double bornXS = myxsection->getXsection( mx ); // std::cout<<"BornXS= "<<bornXS<<std::endl;
585 double bornXS_er = myxsection->getErr( mx );
586 double obsvXS = _xs0 + _xs1; // gamHXSection(mth ,mx);
587 double obsvXS_er = _er0 + xs1_err;
588 double corr = obsvXS / bornXS;
589 double corr_er = corr * sqrt( bornXS_er * bornXS_er / bornXS / bornXS +
590 obsvXS_er * obsvXS_er / obsvXS / obsvXS );
591
592 if ( bornXS == 0 )
593 {
594 std::cout << "EvtConExc::init : The Born cross section at this point is zero!"
595 << std::endl;
596 abort();
597 }
598 if ( obsvXS == 0 )
599 {
600 std::cout << "EvtConExc::init : The observed cross section at this point is zero!"
601 << std::endl;
602 abort();
603 }
604
605 //_xs0 += bornXS;// events with very soft photon (Egam<0.025) are taken as the born process
606 //_er0 = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
607
608 std::cout << "" << std::endl;
609 std::cout << "========= Generation using cross section (Egamcut=" << Egamcut
610 << " GeV) ==============" << std::endl;
611 std::cout << " sqrt(s)= " << mx << " GeV " << std::endl;
612 if ( _mode >= 0 || _mode == -2 || _mode == -100 )
613 {
614 std::cout << "The generated born cross section (sigma_b): " << _xs0 << "+/-" << _er0
615 << " in Unit " << _unit.c_str() << std::endl;
616 std::cout << "The generated radiative correction cross section (sigma_isr): " << _xs1
617 << "+/-" << xs1_err << " in Unit " << _unit.c_str() << std::endl;
618 std::cout
619 << "The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"
620 << std::endl;
621 std::cout << "The radiative correction factor f_ISR*f_vacuum= sigma_obs/sigma_born(s) = "
622 << corr << "+/-" << fabs( corr_er ) << "," << std::endl;
623 std::cout << "which is calcualted with these quantities:" << std::endl;
624 std::cout << "the observed cross section is " << obsvXS << "+/-" << obsvXS_er
625 << _unit.c_str() << std::endl;
626 std::cout << "and the Born cross section (s) is " << bornXS << "+/-" << bornXS_er
627 << _unit.c_str() << std::endl;
628 std::cout << "and the vaccum polarziation factor (lep. + hadr.) 1/|1-Pi|^2= " << vph
629 << std::endl;
630 std::cout << "Within the range " << myxsection->getXlw() << " ~ " << myxsection->getXup()
631 << " GeV, " << myxsection->getMsg().c_str() << std::endl;
632 std::cout << "==========================================================================="
633 << std::endl;
634 }
635 else if ( _mode == -1 )
636 {
637 std::cout << "The generated born cross section (sigma_b): " << _xs0 << " *Br_ee"
638 << " in Unit " << _unit.c_str() << std::endl;
639 std::cout << "The generated radiative correction cross section (sigma_isr): " << _xs1
640 << "*Br_ee in Unit " << _unit.c_str() << std::endl;
641 std::cout
642 << "The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"
643 << std::endl;
644 std::cout << "The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "
645 << corr << "+/-" << fabs( corr_er ) << std::endl;
646 std::cout << "==========================================================================="
647 << std::endl;
648 }
649 std::cout << std::endl;
650 std::cout << std::endl;
651
652 // out put born cross section for each mode
653 if ( _mode == 74110 )
654 {
655 std::cout << " Summary of Born cross section for each exclusive mode" << std::endl;
656 std::cout << "Mode index "
657 << " Born cross section at " << mx << " GeV "
658 << " final states" << std::endl;
659 PrintXS( mx );
660 std::cout << "End of summary Born cross section" << std::endl;
661 }
662
663 findMaxXS( p );
664 _mhdL = myxsection->getXlw();
665 //--debugging
666 // std::cout<<"maxXS= "<<maxXS<<std::endl;
667
668 if ( _xs0 == 0 && _xs1 == 0 )
669 {
670 std::cout << "EvtConExc::ini() has zero cross section" << std::endl;
671 abort();
672 }
673
674 std::cout << std::endl;
675 std::cout << std::endl;
676 mlow = myxsection->getXlw();
677 mup = myxsection->getXup();
678 //--- for debugging
679 if ( mydbg && _mode != 74110 )
680 {
681 int nd = myxsection->getYY().size();
682 double xx[10000], yy[10000], xer[10000], yer[10000];
683 for ( int i = 0; i < nd; i++ )
684 {
685 xx[i] = myxsection->getXX()[i];
686 yy[i] = myxsection->getYY()[i];
687 yer[i] = myxsection->getEr()[i];
688 xer[i] = 0.;
689 // std::cout<<"yy "<<yy[i]<<std::endl;
690 }
691 myth = new TH1F( "myth", "Exp.data", 200, xx[0], xx[nd] );
692 for ( int i = 0; i < nd; i++ ) { myth->Fill( xx[i], yy[i] ); }
693 myth->SetError( yer );
694 myth->Write();
695
696 // std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
697 }
698 else if ( mydbg && _mode == 74110 )
699 {
700 int nd = myxsection->getYY().size();
701 double xx[10000], yy[10000], xer[10000], yer[10000], ysum[10000], yysum[10000];
702 for ( int i = 0; i < nd; i++ )
703 {
704 xx[i] = myxsection->getXX()[i];
705 yy[i] = myxsection->getYY()[i];
706 yer[i] = myxsection->getEr()[i];
707 xer[i] = 0.;
708 std::cout << "yy " << yy[i] << std::endl;
709 }
710
711#include "sty.C"
712 double xhigh = 5.0;
713 double xlow = myxsection->getXlw();
714 TGraphErrors* Gdt = new TGraphErrors( nd, xx, yy, xer, yer );
715
716 myth = new TH1F( "myth", "Exp.data", 600, xlow, xhigh );
717 Xsum = new TH1F( "sumxs", "sum of exclusive xs", 600, xlow, xhigh );
718 double mstp = ( xhigh - xlow ) / 600;
719 for ( int i = 0; i < 600; i++ )
720 {
721 double mx = i * mstp + xlow;
722 double xsi = myxsection->getXsection( mx );
723 if ( xsi == 0 ) continue;
724 myth->Fill( mx, xsi );
725 // std::cout<<mx<<" "<<xsi<<std::endl;
726 }
727
728 for ( int i = 0; i < 600; i++ )
729 {
730 double mx = i * mstp + xlow;
731 ysum[i] = sumExc( mx );
732 if ( ysum[i] == 0 ) continue;
733 Xsum->Fill( mx, ysum[i] );
734 // std::cout<<mx<<" "<<ysum[i]<<std::endl;
735 }
736
737 for ( int i = 0; i < nd; i++ ) { yysum[i] = sumExc( xx[i] ); }
738
739 myth->SetError( yer );
740 myth->Write();
741 Xsum->Write();
742
743 TGraphErrors* Gsum = new TGraphErrors( nd, xx, yysum, xer, yer );
744 // draw graph
745 double ex[600] = { 600 * 0 };
746 double ey[600], ta[600];
747 double exz[600] = { 600 * 0 };
748 double eyz[600] = { 600 * 0 };
749 for ( int i = 0; i < 600; i++ ) { ey[i] = AF[i] * 0.001; }
750 TGraphErrors* gr0 = new TGraphErrors( 600, MH, AF, ex, ey );
751 TGraphErrors* gr1 = new TGraphErrors( 600, MH, RadXS, exz, eyz );
752 TPostScript* ps = new TPostScript( "xsobs.ps", 111 );
753 TCanvas* c1 = new TCanvas( "c1", "TGraphs for data", 200, 10, 600, 400 );
754 gPad->SetLogy( 1 );
755 // c1->SetLogy(1);
756 gStyle->SetCanvasColor( 0 );
757 gStyle->SetStatBorderSize( 1 );
758 c1->Divide( 1, 1 );
759
760 c1->Update();
761 ps->NewPage();
762 c1->Draw();
763 c1->cd( 1 );
764 c1->SetLogy( 1 );
765 gr0->SetMarkerStyle( 10 );
766 gr0->Draw( "AP" );
767 gr0->GetXaxis()->SetTitle( "M_{hadrons} (GeV)" );
768 gr0->GetYaxis()->SetTitle( "observed cross section (nb)" );
769 gr0->Write();
770
771 c1->Update();
772 ps->NewPage();
773 c1->Draw();
774 c1->cd( 1 );
775 c1->SetLogy( 1 );
776 gr1->SetMarkerStyle( 10 );
777 gr1->Draw( "AP" );
778 gr1->GetXaxis()->SetTitle( "M_{hadrons} (GeV)" );
779 gr1->GetYaxis()->SetTitle( "Rad*BornXS" );
780 gr1->Write();
781
782 //--------- imposing graphs to a pad
783 TMultiGraph* mg = new TMultiGraph();
784 mg->SetTitle( "Exclusion graphs" );
785 Gdt->SetMarkerStyle( 2 );
786 Gdt->SetMarkerSize( 1 );
787 Gsum->SetLineColor( 2 );
788 Gsum->SetLineWidth( 1 );
789 mg->Add( Gdt );
790 mg->Add( Gsum );
791
792 c1->Update();
793 ps->NewPage();
794 c1->Draw();
795 c1->cd( 1 );
796 gPad->SetLogy( 1 );
797 mg->Draw( "APL" );
798 mg->GetXaxis()->SetTitle( "M_{hadrons} (GeV)" );
799 mg->GetYaxis()->SetTitle( "observed cross section (nb)" );
800 mg->Write();
801 //-------
802
803 c1->Update();
804 ps->NewPage();
805 c1->Draw();
806 c1->cd( 1 );
807 Gdt->SetMarkerStyle( 2 );
808 Gdt->Draw( "AP" );
809 Gdt->GetXaxis()->SetTitle( "M_{hadrons} (GeV)" );
810 Gdt->GetYaxis()->SetTitle( "observed cross section (nb)" );
811 Gdt->Write();
812
813 c1->Update();
814 ps->NewPage();
815 c1->Draw();
816 c1->cd( 1 );
817 Gsum->SetMarkerStyle( 2 );
818 Gsum->SetMarkerSize( 1 );
819 Gsum->Draw( "AP" );
820 Gsum->SetLineWidth( 0 );
821 Gsum->GetXaxis()->SetTitle( "M_{hadrons} (GeV)" );
822 Gsum->GetYaxis()->SetTitle( "observed cross section (nb)" );
823 Gsum->Write();
824
825 c1->Update();
826 ps->NewPage();
827 c1->Draw();
828 c1->cd( 1 );
829 // gPad-> SetLogx(1);
830 Gdt->SetMarkerStyle( 2 );
831 Gdt->SetMarkerSize( 1 );
832 Gdt->SetMaximum( 8000.0 );
833 Gsum->SetLineColor( 2 );
834 Gsum->SetLineWidth( 1.5 );
835 Gsum->Draw( "AL" ); // A: draw axis
836 Gdt->Draw( "P" ); // superpose to the Gsum, without A-option
837 Gsum->GetXaxis()->SetTitle( "M_{hadrons} (GeV)" );
838 Gsum->GetYaxis()->SetTitle( "cross section (nb)" );
839 Gsum->Write();
840
841 ps->Close();
842 } // if(mydbg)_block
843 //-------------------------
844}
845
846//--
847void EvtConExc::init_mode( int mode ) {
848 if ( mode == 0 )
849 {
850 _ndaugs = 2;
851 daugs[0] = EvtPDL::getId( std::string( "p+" ) );
852 daugs[1] = EvtPDL::getId( std::string( "anti-p-" ) );
853 }
854 else if ( mode == 1 )
855 {
856 _ndaugs = 2;
857 daugs[0] = EvtPDL::getId( std::string( "n0" ) );
858 daugs[1] = EvtPDL::getId( std::string( "anti-n0" ) );
859 }
860 else if ( mode == 2 )
861 {
862 _ndaugs = 2;
863 daugs[0] = EvtPDL::getId( std::string( "Lambda0" ) );
864 daugs[1] = EvtPDL::getId( std::string( "anti-Lambda0" ) );
865 }
866 else if ( mode == 3 )
867 {
868 _ndaugs = 2;
869 daugs[0] = EvtPDL::getId( std::string( "Sigma0" ) );
870 daugs[1] = EvtPDL::getId( std::string( "anti-Sigma0" ) );
871 }
872 else if ( mode == 4 )
873 {
874 _ndaugs = 2;
875 daugs[0] = EvtPDL::getId( std::string( "Lambda0" ) );
876 daugs[1] = EvtPDL::getId( std::string( "anti-Sigma0" ) );
877 }
878 else if ( mode == 5 )
879 {
880 _ndaugs = 2;
881 daugs[0] = EvtPDL::getId( std::string( "Sigma0" ) );
882 daugs[1] = EvtPDL::getId( std::string( "anti-Lambda0" ) );
883 }
884 else if ( mode == 6 )
885 {
886 _ndaugs = 2;
887 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
888 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
889 }
890 else if ( mode == 7 )
891 {
892 _ndaugs = 3;
893 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
894 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
895 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
896 }
897 else if ( mode == 8 )
898 {
899 _ndaugs = 3;
900 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
901 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
902 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
903 }
904 else if ( mode == 9 )
905 {
906 _ndaugs = 3;
907 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
908 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
909 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
910 }
911 else if ( mode == 10 )
912 {
913 _ndaugs = 3;
914 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
915 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
916 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
917 }
918 else if ( mode == 11 )
919 {
920 _ndaugs = 3;
921 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
922 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
923 daugs[2] = EvtPDL::getId( std::string( "eta" ) );
924 }
925 else if ( mode == 12 )
926 {
927 _ndaugs = 4;
928 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
929 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
930 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
931 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
932 }
933 else if ( mode == 13 )
934 {
935 _ndaugs = 4;
936 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
937 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
938 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
939 daugs[3] = EvtPDL::getId( std::string( "pi0" ) );
940 }
941 else if ( mode == 14 )
942 {
943 _ndaugs = 4;
944 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
945 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
946 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
947 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
948 }
949 else if ( mode == 15 )
950 {
951 _ndaugs = 4;
952 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
953 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
954 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
955 daugs[3] = EvtPDL::getId( std::string( "pi0" ) );
956 }
957 else if ( mode == 16 )
958 {
959 _ndaugs = 4;
960 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
961 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
962 daugs[2] = EvtPDL::getId( std::string( "K+" ) );
963 daugs[3] = EvtPDL::getId( std::string( "K-" ) );
964 }
965 else if ( mode == 17 )
966 {
967 _ndaugs = 5;
968 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
969 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
970 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
971 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
972 daugs[4] = EvtPDL::getId( std::string( "pi0" ) );
973 }
974 else if ( mode == 18 )
975 {
976 _ndaugs = 5;
977 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
978 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
979 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
980 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
981 daugs[4] = EvtPDL::getId( std::string( "eta" ) );
982 }
983 else if ( mode == 19 )
984 {
985 _ndaugs = 5;
986 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
987 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
988 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
989 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
990 daugs[4] = EvtPDL::getId( std::string( "pi0" ) );
991 }
992 else if ( mode == 20 )
993 {
994 _ndaugs = 5;
995 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
996 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
997 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
998 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
999 daugs[4] = EvtPDL::getId( std::string( "eta" ) );
1000 }
1001 else if ( mode == 21 )
1002 {
1003 _ndaugs = 6;
1004 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1005 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1006 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1007 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1008 daugs[4] = EvtPDL::getId( std::string( "pi+" ) );
1009 daugs[5] = EvtPDL::getId( std::string( "pi-" ) );
1010 }
1011 else if ( mode == 22 )
1012 {
1013 _ndaugs = 6;
1014 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1015 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1016 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1017 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1018 daugs[4] = EvtPDL::getId( std::string( "pi0" ) );
1019 daugs[5] = EvtPDL::getId( std::string( "pi0" ) );
1020 }
1021 else if ( mode == 23 )
1022 {
1023 _ndaugs = 2;
1024 daugs[0] = EvtPDL::getId( std::string( "eta" ) );
1025 daugs[1] = EvtPDL::getId( std::string( "phi" ) );
1026 }
1027 else if ( mode == 24 )
1028 {
1029 _ndaugs = 2;
1030 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1031 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1032 }
1033 else if ( mode == 25 )
1034 {
1035 _ndaugs = 2;
1036 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1037 daugs[1] = EvtPDL::getId( std::string( "K*-" ) );
1038 }
1039 else if ( mode == 26 )
1040 {
1041 _ndaugs = 2;
1042 daugs[0] = EvtPDL::getId( std::string( "K-" ) );
1043 daugs[1] = EvtPDL::getId( std::string( "K*+" ) );
1044 }
1045 else if ( mode == 27 )
1046 {
1047 _ndaugs = 2;
1048 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
1049 daugs[1] = EvtPDL::getId( std::string( "anti-K*0" ) );
1050 }
1051 else if ( mode == 28 )
1052 {
1053 _ndaugs = 3;
1054 daugs[0] = EvtPDL::getId( std::string( "anti-K*0" ) );
1055 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1056 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1057 }
1058 else if ( mode == 29 )
1059 {
1060 _ndaugs = 3;
1061 daugs[0] = EvtPDL::getId( std::string( "K*0" ) );
1062 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1063 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1064 }
1065 else if ( mode == 30 )
1066 {
1067 _ndaugs = 3;
1068 daugs[0] = EvtPDL::getId( std::string( "K*+" ) );
1069 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1070 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1071 }
1072 else if ( mode == 31 )
1073 {
1074 _ndaugs = 3;
1075 daugs[0] = EvtPDL::getId( std::string( "K*-" ) );
1076 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1077 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1078 }
1079 else if ( mode == 32 )
1080 {
1081 _ndaugs = 3;
1082 daugs[0] = EvtPDL::getId( std::string( "anti-K_2*0" ) );
1083 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1084 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1085 }
1086 else if ( mode == 33 )
1087 {
1088 _ndaugs = 3;
1089 daugs[0] = EvtPDL::getId( std::string( "K_2*0" ) );
1090 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1091 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1092 }
1093 else if ( mode == 34 )
1094 {
1095 _ndaugs = 3;
1096 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1097 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1098 daugs[2] = EvtPDL::getId( std::string( "rho0" ) );
1099 }
1100 else if ( mode == 35 )
1101 {
1102 _ndaugs = 3;
1103 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1104 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1105 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1106 }
1107 else if ( mode == 36 )
1108 {
1109 _ndaugs = 2;
1110 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1111 daugs[1] = EvtPDL::getId( std::string( "f_0" ) );
1112 }
1113 else if ( mode == 37 )
1114 {
1115 _ndaugs = 3;
1116 daugs[0] = EvtPDL::getId( std::string( "eta" ) );
1117 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1118 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1119 }
1120 else if ( mode == 38 )
1121 {
1122 _ndaugs = 3;
1123 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1124 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1125 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1126 }
1127 else if ( mode == 39 )
1128 {
1129 _ndaugs = 2;
1130 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1131 daugs[1] = EvtPDL::getId( std::string( "f_0" ) );
1132 }
1133 else if ( mode == 40 )
1134 {
1135 _ndaugs = 3;
1136 daugs[0] = EvtPDL::getId( std::string( "eta'" ) );
1137 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1138 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1139 }
1140 else if ( mode == 41 )
1141 {
1142 _ndaugs = 3;
1143 daugs[0] = EvtPDL::getId( std::string( "f_1" ) );
1144 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1145 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1146 }
1147 else if ( mode == 42 )
1148 {
1149 _ndaugs = 3;
1150 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1151 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1152 daugs[2] = EvtPDL::getId( std::string( "K-" ) );
1153 }
1154 else if ( mode == 43 )
1155 {
1156 _ndaugs = 4;
1157 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1158 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1159 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1160 daugs[3] = EvtPDL::getId( std::string( "pi0" ) );
1161 }
1162 else if ( mode == 44 )
1163 {
1164 _ndaugs = 2;
1165 daugs[0] = EvtPDL::getId( std::string( "Sigma-" ) );
1166 daugs[1] = EvtPDL::getId( std::string( "anti-Sigma+" ) );
1167 }
1168 else if ( mode == 45 )
1169 {
1170 _ndaugs = 2;
1171 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1172 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1173 }
1174 else if ( mode == 46 )
1175 {
1176 _ndaugs = 2;
1177 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
1178 daugs[1] = EvtPDL::getId( std::string( "K_L0" ) );
1179 }
1180 else if ( mode == 47 )
1181 {
1182 _ndaugs = 2;
1183 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1184 daugs[1] = EvtPDL::getId( std::string( "eta" ) );
1185 }
1186 else if ( mode == 48 )
1187 {
1188 _ndaugs = 2;
1189 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1190 daugs[1] = EvtPDL::getId( std::string( "eta'" ) );
1191 }
1192 else if ( mode == 49 )
1193 {
1194 _ndaugs = 3;
1195 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1196 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1197 daugs[2] = EvtPDL::getId( std::string( "K-" ) );
1198 }
1199 else if ( mode == 50 )
1200 {
1201 _ndaugs = 3;
1202 daugs[0] = EvtPDL::getId( std::string( "p+" ) );
1203 daugs[1] = EvtPDL::getId( std::string( "anti-p-" ) );
1204 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1205 }
1206 else if ( mode == 51 )
1207 {
1208 _ndaugs = 3;
1209 daugs[0] = EvtPDL::getId( std::string( "p+" ) );
1210 daugs[1] = EvtPDL::getId( std::string( "anti-p-" ) );
1211 daugs[2] = EvtPDL::getId( std::string( "eta" ) );
1212 }
1213 else if ( mode == 52 )
1214 {
1215 _ndaugs = 2;
1216 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1217 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1218 }
1219 else if ( mode == 59 )
1220 {
1221 _ndaugs = 3;
1222 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1223 daugs[1] = EvtPDL::getId( std::string( "D0" ) );
1224 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1225 }
1226 else if ( mode == 60 )
1227 {
1228 _ndaugs = 3;
1229 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1230 daugs[1] = EvtPDL::getId( std::string( "D*0" ) );
1231 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1232 }
1233 else if ( mode == 61 )
1234 {
1235 _ndaugs = 3;
1236 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
1237 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
1238 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1239 }
1240 else if ( mode == 62 )
1241 {
1242 _ndaugs = 3;
1243 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
1244 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
1245 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1246 }
1247 else if ( mode == 63 )
1248 {
1249 _ndaugs = 3;
1250 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1251 daugs[1] = EvtPDL::getId( std::string( "D+" ) );
1252 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1253 }
1254 else if ( mode == 64 )
1255 {
1256 _ndaugs = 3;
1257 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1258 daugs[1] = EvtPDL::getId( std::string( "D*+" ) );
1259 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1260 }
1261 else if ( mode == 65 )
1262 {
1263 _ndaugs = 3;
1264 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1265 daugs[1] = EvtPDL::getId( std::string( "D*0" ) );
1266 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1267 }
1268 else if ( mode == 66 )
1269 {
1270 _ndaugs = 3;
1271 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
1272 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
1273 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1274 }
1275 else if ( mode == 67 )
1276 {
1277 _ndaugs = 2;
1278 daugs[0] = EvtPDL::getId( std::string( "D*0" ) );
1279 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
1280 }
1281 else if ( mode == 68 )
1282 {
1283 _ndaugs = 2;
1284 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
1285 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
1286 }
1287 else if ( mode == 69 )
1288 {
1289 _ndaugs = 2;
1290 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1291 daugs[1] = EvtPDL::getId( std::string( "D*0" ) );
1292 }
1293 else if ( mode == 70 )
1294 {
1295 _ndaugs = 2;
1296 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
1297 daugs[1] = EvtPDL::getId( std::string( "anti-D0" ) );
1298 }
1299 else if ( mode == 71 )
1300 {
1301 _ndaugs = 2;
1302 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
1303 daugs[1] = EvtPDL::getId( std::string( "D-" ) );
1304 }
1305 else if ( mode == 72 )
1306 {
1307 _ndaugs = 2;
1308 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
1309 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
1310 }
1311 else if ( mode == 73 )
1312 {
1313 _ndaugs = 2;
1314 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1315 daugs[1] = EvtPDL::getId( std::string( "D*+" ) );
1316 }
1317 else if ( mode == 74 )
1318 {
1319 _ndaugs = 2;
1320 daugs[0] = EvtPDL::getId( std::string( "D*+" ) );
1321 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
1322 }
1323 else if ( mode == 75 )
1324 {
1325 _ndaugs = 3;
1326 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
1327 daugs[1] = EvtPDL::getId( std::string( "D-" ) );
1328 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1329 }
1330 else if ( mode == 76 )
1331 {
1332 _ndaugs = 3;
1333 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1334 daugs[1] = EvtPDL::getId( std::string( "D+" ) );
1335 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1336 }
1337 else if ( mode == 77 )
1338 {
1339 _ndaugs = 3;
1340 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
1341 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
1342 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1343 }
1344 else if ( mode == 78 )
1345 {
1346 _ndaugs = 3;
1347 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1348 daugs[1] = EvtPDL::getId( std::string( "D*+" ) );
1349 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1350 }
1351 else if ( mode == 79 )
1352 {
1353 _ndaugs = 3;
1354 daugs[0] = EvtPDL::getId( std::string( "psi(2S)" ) );
1355 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1356 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1357 }
1358 else if ( mode == 80 )
1359 {
1360 _ndaugs = 2;
1361 daugs[0] = EvtPDL::getId( std::string( "eta" ) );
1362 daugs[1] = EvtPDL::getId( std::string( "J/psi" ) );
1363 }
1364 else if ( mode == 81 )
1365 {
1366 _ndaugs = 3;
1367 daugs[0] = EvtPDL::getId( std::string( "h_c" ) );
1368 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1369 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1370 }
1371 else if ( mode == 82 )
1372 {
1373 _ndaugs = 3;
1374 daugs[0] = EvtPDL::getId( std::string( "h_c" ) );
1375 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1376 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1377 }
1378 else if ( mode == 83 )
1379 {
1380 _ndaugs = 3;
1381 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
1382 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1383 daugs[2] = EvtPDL::getId( std::string( "K-" ) );
1384 }
1385 else if ( mode == 84 )
1386 {
1387 _ndaugs = 3;
1388 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
1389 daugs[1] = EvtPDL::getId( std::string( "K_S0" ) );
1390 daugs[2] = EvtPDL::getId( std::string( "K_S0" ) );
1391 }
1392 else if ( mode == 85 )
1393 {
1394 _ndaugs = 2;
1395 daugs[0] = EvtPDL::getId( std::string( "D_s*+" ) );
1396 daugs[1] = EvtPDL::getId( std::string( "D_s*-" ) );
1397 }
1398 else if ( mode == 90 )
1399 {
1400 _ndaugs = 3;
1401 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
1402 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1403 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1404 }
1405 else if ( mode == 91 )
1406 {
1407 _ndaugs = 3;
1408 daugs[0] = EvtPDL::getId( std::string( "psi(2S)" ) );
1409 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1410 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1411 }
1412 else if ( mode == 92 )
1413 {
1414 _ndaugs = 3;
1415 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
1416 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1417 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1418 }
1419 else if ( mode == 93 )
1420 {
1421 _ndaugs = 2;
1422 daugs[0] = EvtPDL::getId( std::string( "D_s+" ) );
1423 daugs[1] = EvtPDL::getId( std::string( "D_s-" ) );
1424 }
1425 else if ( mode == 94 )
1426 {
1427 _ndaugs = 2;
1428 daugs[0] = EvtPDL::getId( std::string( "D_s*+" ) );
1429 daugs[1] = EvtPDL::getId( std::string( "D_s-" ) );
1430 }
1431 else if ( mode == 95 )
1432 {
1433 _ndaugs = 2;
1434 daugs[0] = EvtPDL::getId( std::string( "D_s*-" ) );
1435 daugs[1] = EvtPDL::getId( std::string( "D_s+" ) );
1436 }
1437 else if ( mode == 96 )
1438 {
1439 _ndaugs = 2;
1440 daugs[0] = EvtPDL::getId( std::string( "Lambda_c+" ) );
1441 daugs[1] = EvtPDL::getId( std::string( "anti-Lambda_c-" ) );
1442 }
1443 else if ( mode == 97 )
1444 {
1445 _ndaugs = 3;
1446 daugs[0] = EvtPDL::getId( std::string( "K-" ) );
1447 daugs[1] = EvtPDL::getId( std::string( "D_s+" ) );
1448 daugs[2] = EvtPDL::getId( std::string( "anti-D0" ) );
1449 }
1450 else if ( mode == 111 )
1451 {
1452 _ndaugs = 2;
1453 daugs[0] = EvtPDL::getId( std::string( "Zc(4020)+" ) );
1454 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1455 }
1456 else if ( mode == 112 )
1457 {
1458 _ndaugs = 2;
1459 daugs[0] = EvtPDL::getId( std::string( "Zc(4020)-" ) );
1460 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1461 }
1462 else if ( mode == 74100 )
1463 {
1464 _ndaugs = 1;
1465 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
1466 }
1467 else if ( mode == 74101 )
1468 {
1469 _ndaugs = 1;
1470 daugs[0] = EvtPDL::getId( std::string( "psi(2S)" ) );
1471 }
1472 else if ( mode == 74102 )
1473 {
1474 _ndaugs = 1;
1475 daugs[0] = EvtPDL::getId( std::string( "psi(3770)" ) );
1476 }
1477 else if ( mode == 74103 )
1478 {
1479 _ndaugs = 1;
1480 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1481 }
1482 else if ( mode == 74104 )
1483 {
1484 _ndaugs = 1;
1485 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1486 }
1487 else if ( mode == 74105 )
1488 {
1489 _ndaugs = 1;
1490 daugs[0] = EvtPDL::getId( std::string( "rho0" ) );
1491 }
1492 else if ( mode == 74106 )
1493 {
1494 _ndaugs = 1;
1495 daugs[0] = EvtPDL::getId( std::string( "rho(3S)0" ) );
1496 }
1497 else if ( mode == 74107 )
1498 {
1499 _ndaugs = 1;
1500 daugs[0] = EvtPDL::getId( std::string( "omega(2S)" ) );
1501 }
1502 else if ( mode == 74110 )
1503 {
1504 _ndaugs = 1;
1505 daugs[0] = EvtPDL::getId( std::string( "gamma*" ) );
1506 EvtId myvpho = EvtPDL::getId( std::string( "gamma*" ) );
1507 EvtPDL::reSetMass( myvpho, mhds * 0.9 ); // for calculating maxXS, it will be resent when
1508 // this mode is selected
1509 _modeFlag.clear();
1510 _modeFlag.push_back( 74110 ); // R-value generator tag
1511 _modeFlag.push_back( 0 ); // to contain the mode selected
1512 }
1513 else if ( mode == -1 )
1514 {
1515 _ndaugs = nson;
1516 for ( int i = 0; i < nson; i++ ) { daugs[i] = son[i]; }
1517 std::cout << "The paraent particle: " << EvtPDL::name( pid ) << std::endl;
1518 }
1519 else if ( mode == -2 )
1520 {
1521 _ndaugs = nson;
1522 for ( int i = 0; i < nson; i++ ) { daugs[i] = son[i]; }
1523 }
1524 else if ( mode == -100 )
1525 {
1526 _ndaugs = 1;
1527 daugs[0] = EvtPDL::getId( std::string( "gamma*" ) );
1528 _modeFlag.clear();
1529 _modeFlag.push_back( -100 ); // R-value generator tag
1530 _modeFlag.push_back( 0 ); // to contain the mode selected
1531 return;
1532 }
1533 else if ( mode == 10000 )
1534 { // use for check ISR
1535 _ndaugs = 2;
1536 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1537 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1538 }
1539 else
1540 {
1541 std::cout << "Bad mode_index number " << mode << ", please refer to the manual."
1542 << std::endl;
1543 ::abort();
1544 }
1545 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1546
1547 if ( VISR )
1548 {
1549 _ndaugs = 2;
1550 daugs[0] = getDaugs()[0]; // EvtPDL::getId(std::string("gamma"));
1551 daugs[1] = getDaugs()[1]; // EvtPDL::getId(std::string("gamma*"));
1552 }
1553
1554 double fmth = 0;
1555 for ( int i = 0; i < _ndaugs; i++ )
1556 {
1557 double xmi = EvtPDL::getMinMass( daugs[i] );
1558 fmth += xmi;
1559 }
1560 myxsection = new EvtXsection( mode );
1561 Mthr = myxsection->getXlw();
1562 if ( !( _mode == 74110 || _mode == -100 ) )
1563 {
1564 if ( Mthr < fmth )
1565 {
1566 std::cout << "Low end of cross section: (" << Mthr << " < (mass of final state)" << fmth
1567 << ") GeV." << std::endl;
1568 abort();
1569 }
1570 }
1571}
1572
1573//--
1574std::vector<EvtId> EvtConExc::get_mode( int mode ) {
1575 int _ndaugs;
1576 EvtId daugs[100];
1577 if ( mode == 0 )
1578 {
1579 _ndaugs = 2;
1580 daugs[0] = EvtPDL::getId( std::string( "p+" ) );
1581 daugs[1] = EvtPDL::getId( std::string( "anti-p-" ) );
1582 }
1583 else if ( mode == 1 )
1584 {
1585 _ndaugs = 2;
1586 daugs[0] = EvtPDL::getId( std::string( "n0" ) );
1587 daugs[1] = EvtPDL::getId( std::string( "anti-n0" ) );
1588 }
1589 else if ( mode == 2 )
1590 {
1591 _ndaugs = 2;
1592 daugs[0] = EvtPDL::getId( std::string( "Lambda0" ) );
1593 daugs[1] = EvtPDL::getId( std::string( "anti-Lambda0" ) );
1594 }
1595 else if ( mode == 3 )
1596 {
1597 _ndaugs = 2;
1598 daugs[0] = EvtPDL::getId( std::string( "Sigma0" ) );
1599 daugs[1] = EvtPDL::getId( std::string( "anti-Sigma0" ) );
1600 }
1601 else if ( mode == 4 )
1602 {
1603 _ndaugs = 2;
1604 daugs[0] = EvtPDL::getId( std::string( "Lambda0" ) );
1605 daugs[1] = EvtPDL::getId( std::string( "anti-Sigma0" ) );
1606 }
1607 else if ( mode == 5 )
1608 {
1609 _ndaugs = 2;
1610 daugs[0] = EvtPDL::getId( std::string( "Sigma0" ) );
1611 daugs[1] = EvtPDL::getId( std::string( "anti-Lambda0" ) );
1612 }
1613 else if ( mode == 6 )
1614 {
1615 _ndaugs = 2;
1616 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1617 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1618 }
1619 else if ( mode == 7 )
1620 {
1621 _ndaugs = 3;
1622 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1623 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1624 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1625 }
1626 else if ( mode == 8 )
1627 {
1628 _ndaugs = 3;
1629 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1630 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1631 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1632 }
1633 else if ( mode == 9 )
1634 {
1635 _ndaugs = 3;
1636 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
1637 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1638 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1639 }
1640 else if ( mode == 10 )
1641 {
1642 _ndaugs = 3;
1643 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
1644 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1645 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1646 }
1647 else if ( mode == 11 )
1648 {
1649 _ndaugs = 3;
1650 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1651 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1652 daugs[2] = EvtPDL::getId( std::string( "eta" ) );
1653 }
1654 else if ( mode == 12 )
1655 {
1656 _ndaugs = 4;
1657 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1658 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1659 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1660 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1661 }
1662 else if ( mode == 13 )
1663 {
1664 _ndaugs = 4;
1665 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1666 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1667 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1668 daugs[3] = EvtPDL::getId( std::string( "pi0" ) );
1669 }
1670 else if ( mode == 14 )
1671 {
1672 _ndaugs = 4;
1673 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1674 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1675 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1676 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1677 }
1678 else if ( mode == 15 )
1679 {
1680 _ndaugs = 4;
1681 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1682 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1683 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1684 daugs[3] = EvtPDL::getId( std::string( "pi0" ) );
1685 }
1686 else if ( mode == 16 )
1687 {
1688 _ndaugs = 4;
1689 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1690 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1691 daugs[2] = EvtPDL::getId( std::string( "K+" ) );
1692 daugs[3] = EvtPDL::getId( std::string( "K-" ) );
1693 }
1694 else if ( mode == 17 )
1695 {
1696 _ndaugs = 5;
1697 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1698 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1699 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1700 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1701 daugs[4] = EvtPDL::getId( std::string( "pi0" ) );
1702 }
1703 else if ( mode == 18 )
1704 {
1705 _ndaugs = 5;
1706 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1707 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1708 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1709 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1710 daugs[4] = EvtPDL::getId( std::string( "eta" ) );
1711 }
1712 else if ( mode == 19 )
1713 {
1714 _ndaugs = 5;
1715 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1716 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1717 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1718 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1719 daugs[4] = EvtPDL::getId( std::string( "pi0" ) );
1720 }
1721 else if ( mode == 20 )
1722 {
1723 _ndaugs = 5;
1724 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1725 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1726 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1727 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1728 daugs[4] = EvtPDL::getId( std::string( "eta" ) );
1729 }
1730 else if ( mode == 21 )
1731 {
1732 _ndaugs = 6;
1733 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1734 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1735 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1736 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1737 daugs[4] = EvtPDL::getId( std::string( "pi+" ) );
1738 daugs[5] = EvtPDL::getId( std::string( "pi-" ) );
1739 }
1740 else if ( mode == 22 )
1741 {
1742 _ndaugs = 6;
1743 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
1744 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1745 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1746 daugs[3] = EvtPDL::getId( std::string( "pi-" ) );
1747 daugs[4] = EvtPDL::getId( std::string( "pi0" ) );
1748 daugs[5] = EvtPDL::getId( std::string( "pi0" ) );
1749 }
1750 else if ( mode == 23 )
1751 {
1752 _ndaugs = 2;
1753 daugs[0] = EvtPDL::getId( std::string( "eta" ) );
1754 daugs[1] = EvtPDL::getId( std::string( "phi" ) );
1755 }
1756 else if ( mode == 24 )
1757 {
1758 _ndaugs = 2;
1759 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1760 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1761 }
1762 else if ( mode == 25 )
1763 {
1764 _ndaugs = 2;
1765 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1766 daugs[1] = EvtPDL::getId( std::string( "K*-" ) );
1767 }
1768 else if ( mode == 26 )
1769 {
1770 _ndaugs = 2;
1771 daugs[0] = EvtPDL::getId( std::string( "K-" ) );
1772 daugs[1] = EvtPDL::getId( std::string( "K*+" ) );
1773 }
1774 else if ( mode == 27 )
1775 {
1776 _ndaugs = 2;
1777 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
1778 daugs[1] = EvtPDL::getId( std::string( "anti-K*0" ) );
1779 }
1780 else if ( mode == 28 )
1781 {
1782 _ndaugs = 3;
1783 daugs[0] = EvtPDL::getId( std::string( "anti-K*0" ) );
1784 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1785 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1786 }
1787 else if ( mode == 29 )
1788 {
1789 _ndaugs = 3;
1790 daugs[0] = EvtPDL::getId( std::string( "K*0" ) );
1791 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1792 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1793 }
1794 else if ( mode == 30 )
1795 {
1796 _ndaugs = 3;
1797 daugs[0] = EvtPDL::getId( std::string( "K*+" ) );
1798 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1799 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1800 }
1801 else if ( mode == 31 )
1802 {
1803 _ndaugs = 3;
1804 daugs[0] = EvtPDL::getId( std::string( "K*-" ) );
1805 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1806 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1807 }
1808 else if ( mode == 32 )
1809 {
1810 _ndaugs = 3;
1811 daugs[0] = EvtPDL::getId( std::string( "anti-K_2*0" ) );
1812 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1813 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1814 }
1815 else if ( mode == 33 )
1816 {
1817 _ndaugs = 3;
1818 daugs[0] = EvtPDL::getId( std::string( "K_2*0" ) );
1819 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1820 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1821 }
1822 else if ( mode == 34 )
1823 {
1824 _ndaugs = 3;
1825 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1826 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1827 daugs[2] = EvtPDL::getId( std::string( "rho0" ) );
1828 }
1829 else if ( mode == 35 )
1830 {
1831 _ndaugs = 3;
1832 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1833 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
1834 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1835 }
1836 else if ( mode == 36 )
1837 {
1838 _ndaugs = 2;
1839 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1840 daugs[1] = EvtPDL::getId( std::string( "f_0" ) );
1841 }
1842 else if ( mode == 37 )
1843 {
1844 _ndaugs = 3;
1845 daugs[0] = EvtPDL::getId( std::string( "eta" ) );
1846 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1847 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1848 }
1849 else if ( mode == 38 )
1850 {
1851 _ndaugs = 3;
1852 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1853 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1854 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1855 }
1856 else if ( mode == 39 )
1857 {
1858 _ndaugs = 2;
1859 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1860 daugs[1] = EvtPDL::getId( std::string( "f_0" ) );
1861 }
1862 else if ( mode == 40 )
1863 {
1864 _ndaugs = 3;
1865 daugs[0] = EvtPDL::getId( std::string( "eta'" ) );
1866 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1867 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1868 }
1869 else if ( mode == 41 )
1870 {
1871 _ndaugs = 3;
1872 daugs[0] = EvtPDL::getId( std::string( "f_1" ) );
1873 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1874 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1875 }
1876 else if ( mode == 42 )
1877 {
1878 _ndaugs = 3;
1879 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1880 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1881 daugs[2] = EvtPDL::getId( std::string( "K-" ) );
1882 }
1883 else if ( mode == 43 )
1884 {
1885 _ndaugs = 4;
1886 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1887 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
1888 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
1889 daugs[3] = EvtPDL::getId( std::string( "pi0" ) );
1890 }
1891 else if ( mode == 44 )
1892 {
1893 _ndaugs = 2;
1894 daugs[0] = EvtPDL::getId( std::string( "Sigma-" ) );
1895 daugs[1] = EvtPDL::getId( std::string( "anti-Sigma+" ) );
1896 }
1897 else if ( mode == 45 )
1898 {
1899 _ndaugs = 2;
1900 daugs[0] = EvtPDL::getId( std::string( "K+" ) );
1901 daugs[1] = EvtPDL::getId( std::string( "K-" ) );
1902 }
1903 else if ( mode == 46 )
1904 {
1905 _ndaugs = 2;
1906 daugs[0] = EvtPDL::getId( std::string( "K_S0" ) );
1907 daugs[1] = EvtPDL::getId( std::string( "K_L0" ) );
1908 }
1909 else if ( mode == 47 )
1910 {
1911 _ndaugs = 2;
1912 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1913 daugs[1] = EvtPDL::getId( std::string( "eta" ) );
1914 }
1915 else if ( mode == 48 )
1916 {
1917 _ndaugs = 2;
1918 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1919 daugs[1] = EvtPDL::getId( std::string( "eta'" ) );
1920 }
1921 else if ( mode == 49 )
1922 {
1923 _ndaugs = 3;
1924 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
1925 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
1926 daugs[2] = EvtPDL::getId( std::string( "K-" ) );
1927 }
1928 else if ( mode == 50 )
1929 {
1930 _ndaugs = 3;
1931 daugs[0] = EvtPDL::getId( std::string( "p+" ) );
1932 daugs[1] = EvtPDL::getId( std::string( "anti-p-" ) );
1933 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1934 }
1935 else if ( mode == 51 )
1936 {
1937 _ndaugs = 3;
1938 daugs[0] = EvtPDL::getId( std::string( "p+" ) );
1939 daugs[1] = EvtPDL::getId( std::string( "anti-p-" ) );
1940 daugs[2] = EvtPDL::getId( std::string( "eta" ) );
1941 }
1942 else if ( mode == 52 )
1943 {
1944 _ndaugs = 2;
1945 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
1946 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
1947 }
1948 else if ( mode == 59 )
1949 {
1950 _ndaugs = 3;
1951 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1952 daugs[1] = EvtPDL::getId( std::string( "D0" ) );
1953 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1954 }
1955 else if ( mode == 60 )
1956 {
1957 _ndaugs = 3;
1958 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
1959 daugs[1] = EvtPDL::getId( std::string( "D*0" ) );
1960 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1961 }
1962 else if ( mode == 61 )
1963 {
1964 _ndaugs = 3;
1965 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
1966 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
1967 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1968 }
1969 else if ( mode == 62 )
1970 {
1971 _ndaugs = 3;
1972 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
1973 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
1974 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1975 }
1976 else if ( mode == 63 )
1977 {
1978 _ndaugs = 3;
1979 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1980 daugs[1] = EvtPDL::getId( std::string( "D+" ) );
1981 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1982 }
1983 else if ( mode == 64 )
1984 {
1985 _ndaugs = 3;
1986 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1987 daugs[1] = EvtPDL::getId( std::string( "D*+" ) );
1988 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
1989 }
1990 else if ( mode == 65 )
1991 {
1992 _ndaugs = 3;
1993 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
1994 daugs[1] = EvtPDL::getId( std::string( "D*0" ) );
1995 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
1996 }
1997 else if ( mode == 66 )
1998 {
1999 _ndaugs = 3;
2000 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
2001 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
2002 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
2003 }
2004 else if ( mode == 67 )
2005 {
2006 _ndaugs = 2;
2007 daugs[0] = EvtPDL::getId( std::string( "D*0" ) );
2008 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
2009 }
2010 else if ( mode == 68 )
2011 {
2012 _ndaugs = 2;
2013 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
2014 daugs[1] = EvtPDL::getId( std::string( "anti-D*0" ) );
2015 }
2016 else if ( mode == 69 )
2017 {
2018 _ndaugs = 2;
2019 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
2020 daugs[1] = EvtPDL::getId( std::string( "D*0" ) );
2021 }
2022 else if ( mode == 70 )
2023 {
2024 _ndaugs = 2;
2025 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
2026 daugs[1] = EvtPDL::getId( std::string( "anti-D0" ) );
2027 }
2028 else if ( mode == 71 )
2029 {
2030 _ndaugs = 2;
2031 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
2032 daugs[1] = EvtPDL::getId( std::string( "D-" ) );
2033 }
2034 else if ( mode == 72 )
2035 {
2036 _ndaugs = 2;
2037 daugs[0] = EvtPDL::getId( std::string( "D+" ) );
2038 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
2039 }
2040 else if ( mode == 73 )
2041 {
2042 _ndaugs = 2;
2043 daugs[0] = EvtPDL::getId( std::string( "D-" ) );
2044 daugs[1] = EvtPDL::getId( std::string( "D*+" ) );
2045 }
2046 else if ( mode == 74 )
2047 {
2048 _ndaugs = 2;
2049 daugs[0] = EvtPDL::getId( std::string( "D*+" ) );
2050 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
2051 }
2052 else if ( mode == 75 )
2053 {
2054 _ndaugs = 3;
2055 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
2056 daugs[1] = EvtPDL::getId( std::string( "D-" ) );
2057 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
2058 }
2059 else if ( mode == 76 )
2060 {
2061 _ndaugs = 3;
2062 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
2063 daugs[1] = EvtPDL::getId( std::string( "D+" ) );
2064 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
2065 }
2066 else if ( mode == 77 )
2067 {
2068 _ndaugs = 3;
2069 daugs[0] = EvtPDL::getId( std::string( "D0" ) );
2070 daugs[1] = EvtPDL::getId( std::string( "D*-" ) );
2071 daugs[2] = EvtPDL::getId( std::string( "pi+" ) );
2072 }
2073 else if ( mode == 78 )
2074 {
2075 _ndaugs = 3;
2076 daugs[0] = EvtPDL::getId( std::string( "anti-D0" ) );
2077 daugs[1] = EvtPDL::getId( std::string( "D*+" ) );
2078 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
2079 }
2080 else if ( mode == 79 )
2081 {
2082 _ndaugs = 3;
2083 daugs[0] = EvtPDL::getId( std::string( "psi(2S)" ) );
2084 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
2085 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
2086 }
2087 else if ( mode == 80 )
2088 {
2089 _ndaugs = 2;
2090 daugs[0] = EvtPDL::getId( std::string( "eta" ) );
2091 daugs[1] = EvtPDL::getId( std::string( "J/psi" ) );
2092 }
2093 else if ( mode == 81 )
2094 {
2095 _ndaugs = 3;
2096 daugs[0] = EvtPDL::getId( std::string( "h_c" ) );
2097 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
2098 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
2099 }
2100 else if ( mode == 82 )
2101 {
2102 _ndaugs = 3;
2103 daugs[0] = EvtPDL::getId( std::string( "h_c" ) );
2104 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
2105 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
2106 }
2107 else if ( mode == 83 )
2108 {
2109 _ndaugs = 3;
2110 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
2111 daugs[1] = EvtPDL::getId( std::string( "K+" ) );
2112 daugs[2] = EvtPDL::getId( std::string( "K-" ) );
2113 }
2114 else if ( mode == 84 )
2115 {
2116 _ndaugs = 3;
2117 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
2118 daugs[1] = EvtPDL::getId( std::string( "K_S0" ) );
2119 daugs[2] = EvtPDL::getId( std::string( "K_S0" ) );
2120 }
2121 else if ( mode == 85 )
2122 {
2123 _ndaugs = 2;
2124 daugs[0] = EvtPDL::getId( std::string( "D_s*+" ) );
2125 daugs[1] = EvtPDL::getId( std::string( "D_s*-" ) );
2126 }
2127 else if ( mode == 90 )
2128 {
2129 _ndaugs = 3;
2130 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
2131 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
2132 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
2133 }
2134 else if ( mode == 91 )
2135 {
2136 _ndaugs = 3;
2137 daugs[0] = EvtPDL::getId( std::string( "psi(2S)" ) );
2138 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
2139 daugs[2] = EvtPDL::getId( std::string( "pi-" ) );
2140 }
2141 else if ( mode == 92 )
2142 {
2143 _ndaugs = 3;
2144 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
2145 daugs[1] = EvtPDL::getId( std::string( "pi0" ) );
2146 daugs[2] = EvtPDL::getId( std::string( "pi0" ) );
2147 }
2148 else if ( mode == 93 )
2149 {
2150 _ndaugs = 2;
2151 daugs[0] = EvtPDL::getId( std::string( "D_s+" ) );
2152 daugs[1] = EvtPDL::getId( std::string( "D_s-" ) );
2153 }
2154 else if ( mode == 94 )
2155 {
2156 _ndaugs = 2;
2157 daugs[0] = EvtPDL::getId( std::string( "D_s*+" ) );
2158 daugs[1] = EvtPDL::getId( std::string( "D_s-" ) );
2159 }
2160 else if ( mode == 95 )
2161 {
2162 _ndaugs = 2;
2163 daugs[0] = EvtPDL::getId( std::string( "D_s*-" ) );
2164 daugs[1] = EvtPDL::getId( std::string( "D_s+" ) );
2165 }
2166 else if ( mode == 96 )
2167 {
2168 _ndaugs = 2;
2169 daugs[0] = EvtPDL::getId( std::string( "Lambda_c+" ) );
2170 daugs[1] = EvtPDL::getId( std::string( "anti-Lambda_c-" ) );
2171 }
2172 else if ( mode == 97 )
2173 {
2174 _ndaugs = 3;
2175 daugs[0] = EvtPDL::getId( std::string( "K-" ) );
2176 daugs[1] = EvtPDL::getId( std::string( "D_s+" ) );
2177 daugs[2] = EvtPDL::getId( std::string( "anti-D0" ) );
2178 }
2179 else if ( mode == 111 )
2180 {
2181 _ndaugs = 2;
2182 daugs[0] = EvtPDL::getId( std::string( "Zc(4020)+" ) );
2183 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
2184 }
2185 else if ( mode == 112 )
2186 {
2187 _ndaugs = 2;
2188 daugs[0] = EvtPDL::getId( std::string( "Zc(4020)-" ) );
2189 daugs[1] = EvtPDL::getId( std::string( "pi+" ) );
2190 }
2191 else if ( mode == 74100 )
2192 {
2193 _ndaugs = 2;
2194 daugs[0] = EvtPDL::getId( std::string( "J/psi" ) );
2195 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2196 }
2197 else if ( mode == 74101 )
2198 {
2199 _ndaugs = 2;
2200 daugs[0] = EvtPDL::getId( std::string( "psi(2S)" ) );
2201 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2202 }
2203 else if ( mode == 74102 )
2204 {
2205 _ndaugs = 2;
2206 daugs[0] = EvtPDL::getId( std::string( "psi(3770)" ) );
2207 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2208 }
2209 else if ( mode == 74103 )
2210 {
2211 _ndaugs = 2;
2212 daugs[0] = EvtPDL::getId( std::string( "phi" ) );
2213 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2214 }
2215 else if ( mode == 74104 )
2216 {
2217 _ndaugs = 2;
2218 daugs[0] = EvtPDL::getId( std::string( "omega" ) );
2219 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2220 }
2221 else if ( mode == 74105 )
2222 {
2223 _ndaugs = 2;
2224 daugs[0] = EvtPDL::getId( std::string( "rho0" ) );
2225 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2226 }
2227 else if ( mode == 74106 )
2228 {
2229 _ndaugs = 2;
2230 daugs[0] = EvtPDL::getId( std::string( "rho(3S)0" ) );
2231 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2232 }
2233 else if ( mode == 74107 )
2234 {
2235 _ndaugs = 2;
2236 daugs[0] = EvtPDL::getId( std::string( "omega(2S)" ) );
2237 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2238 }
2239 else if ( mode == 74110 )
2240 {
2241 _ndaugs = 2;
2242 daugs[0] = EvtPDL::getId( std::string( "gamma*" ) );
2243 daugs[1] = EvtPDL::getId( std::string( "gamma" ) );
2244 EvtId myvpho = EvtPDL::getId( std::string( "gamma*" ) );
2245 EvtPDL::reSetMass( myvpho, mhds * 0.9 ); // for calculating maxXS, it will be resent when
2246 // this mode is selected
2247 _modeFlag.clear();
2248 _modeFlag.push_back( 74110 ); // R-value generator tag
2249 _modeFlag.push_back( 0 ); // to contain the mode selected
2250 }
2251 else if ( mode == -1 )
2252 {
2253 _ndaugs = nson;
2254 for ( int i = 0; i < nson; i++ ) { daugs[i] = son[i]; }
2255 std::cout << "The paraent particle: " << EvtPDL::name( pid ) << std::endl;
2256 }
2257 else if ( mode == -2 )
2258 {
2259 _ndaugs = nson;
2260 for ( int i = 0; i < nson; i++ ) { daugs[i] = son[i]; }
2261 }
2262 else if ( mode == -100 )
2263 {
2264 _ndaugs = 1;
2265 daugs[0] = EvtPDL::getId( std::string( "gamma*" ) );
2266 _modeFlag.clear();
2267 _modeFlag.push_back( -100 ); // R-value generator tag
2268 _modeFlag.push_back( 0 ); // to contain the mode selected
2269 }
2270 else if ( mode == 10000 )
2271 { // use for check ISR
2272 _ndaugs = 2;
2273 daugs[0] = EvtPDL::getId( std::string( "pi+" ) );
2274 daugs[1] = EvtPDL::getId( std::string( "pi-" ) );
2275 }
2276 else
2277 {
2278 std::cout << "Bad mode_index number " << mode << ", please refer to the manual."
2279 << std::endl;
2280 ::abort();
2281 }
2282 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
2283
2284 if ( VISR )
2285 {
2286 _ndaugs = 2;
2287 daugs[0] = getDaugs()[0]; // EvtPDL::getId(std::string("gamma"));
2288 daugs[1] = getDaugs()[1]; // EvtPDL::getId(std::string("gamma*"));
2289 }
2290
2291 std::vector<EvtId> theDaugs;
2292 for ( int i = 0; i < _ndaugs; i++ ) { theDaugs.push_back( daugs[i] ); }
2293 if ( theDaugs.size() ) { return theDaugs; }
2294 else
2295 {
2296 std::cout << "EvtConExc::get_mode: zero size" << std::endl;
2297 abort();
2298 }
2299}
2300
2301void EvtConExc::initProbMax() { noProbMax(); }
2302
2303void EvtConExc::decay( EvtParticle* p ) {
2304 // std::cout<<"_cms= "<<_cms<<" mode= "<<_mode<<std::endl;
2305 EvtId myvpho = EvtPDL::getId( std::string( "vpho" ) );
2306 if ( myvpho != p->getId() )
2307 {
2308 std::cout << "Parent needs to be vpho, but found " << EvtPDL::name( p->getId() )
2309 << std::endl;
2310 abort();
2311 }
2312 //-- to read _cms from database or to consider beam energy spread
2313 double mvpho = EvtPDL::getMeanMass( myvpho );
2314 bool beamflag = 0;
2315 if ( mvpho != EvtConExc::_cms )
2316 { // if read cms energy from database
2317 EvtConExc::_cms = mvpho;
2318 beamflag = 1;
2319 }
2320
2321 /*
2322 if(beamEnergySpread!=0){
2323 cmsspread=sqrt(2.)*beamEnergySpread;
2324 _cms = energySpread(_cms,cmsspread);
2325 beamflag=1;
2326 std::cout<<"test_cms "<<_cms<<std::endl;
2327 }
2328 */
2329
2330 if ( beamflag ) calAF( _cms );
2331
2332 // debugging
2333 // double mvpho=EvtPDL::getMeanMass(myvpho);
2334 // std::cout<<"myvpho= "<<mvpho<<std::endl;
2335
2336 // for fill root tree
2337 EvtVector4R vgam, hd1, hd2, vhds;
2338
2339 // first for Rscan generator with _mode=74110
2340 if ( _mode == 74110 )
2341 {
2342 // preparation of mode sampling
2343 std::vector<int> vmod;
2344 vmod.clear();
2345 int mn[83] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 16, 17, 18, 19, 20, 21, 22,
2346 23, 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 40, 41, 44, 45,
2347 46, // 12,14, 30,31,39,42 and 43 is
2348 // removed
2349 50, 51, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
2350 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100,
2351 74101, 74102, 74103, 74104, 74105, 74110 };
2352 int mn2[84] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 20, 21,
2353 22, 23, 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 40, 41, 44, 45,
2354 46, // mode 14,30,31,42, 43 are
2355 // removed
2356 50, 51, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
2357 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100,
2358 74101, 74102, 74103, 74104, 74105, 74110 };
2359 int mn3[77] = {
2360 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
2361 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
2362 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 44, 45, 46, 50, 51,
2363 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83,
2364 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100, 74101, 74102, 74110 }; // remove 43,
2365 // 74103,74104,74105,
2366 // this is included
2367 // in PHOKHARA
2368 double mx = p->getP4().mass();
2369
2370 if ( _cms > 2.5 )
2371 {
2372 for ( int i = 0; i < 83; i++ ) { vmod.push_back( mn[i] ); }
2373 }
2374 else
2375 {
2376 for ( int i = 0; i < 84; i++ ) { vmod.push_back( mn2[i] ); }
2377 }
2378
2379 // for(int i=0;i<76;i++){vmod.push_back(mn3[i]);}
2380
2381 int mymode;
2382 double pm = EvtRandom::Flat( 0., 1 );
2383
2384 // std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
2385 //--
2386 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross
2387 // section within Egam=(0.025-Egamma^max GeV)
2388 if ( pm < _xs0 / ( _xs0 + _xs1 ) )
2389 { // without ISR photon for mode=74110
2390
2391 mhds = _cms;
2392 mymode = selectMode( vmod, mhds );
2393 _selectedMode = mymode;
2394 std::cout << "A selected mode is " << mymode << " with Mhds= " << mhds
2395 << std::endl; // debugging
2396 delete myxsection;
2397 myxsection = new EvtXsection( mymode );
2398 init_mode( mymode );
2399 resetResMass();
2400
2401 if ( _ndaugs > 1 )
2402 { // for e+e- -> exclusive decays
2403 checkA:
2404 p->makeDaughters( _ndaugs, daugs );
2405 double totMass = 0;
2406 for ( int i = 0; i < _ndaugs; i++ )
2407 {
2408 EvtParticle* di = p->getDaug( i );
2409 double xmi = EvtPDL::getMass( di->getId() );
2410 di->setMass( xmi );
2411 totMass += xmi;
2412 }
2413 if ( totMass > p->mass() ) goto checkA;
2414 p->initializePhaseSpace( _ndaugs, daugs );
2415 if ( !checkdecay( p ) ) goto checkA;
2416 vhds = p->getP4();
2417 if ( _cms > 2.5 && !angularSampling( p ) ) goto checkA;
2418 if ( _cms < 2.5 && !photonSampling( p ) ) goto checkA;
2419 }
2420 else
2421 { // for e+e- -> vector resonace, add a very soft photon
2422 EvtId mydaugs[2];
2423 mydaugs[0] = daugs[0];
2424 EvtPDL::reSetMass( mydaugs[0], mhds - 0.001 );
2425 // EvtPDL::reSetWidth(mydaugs[0],0);
2426 mydaugs[1] = gamId; // ISR photon
2427 p->makeDaughters( 2, mydaugs );
2428 checkB:
2429 double totMass = 0;
2430 for ( int i = 0; i < 2; i++ )
2431 {
2432 EvtParticle* di = p->getDaug( i );
2433 double xmi = EvtPDL::getMass( di->getId() );
2434 di->setMass( xmi );
2435 totMass += xmi;
2436 }
2437 // std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<"
2438 // "<<checkdecay(p)<<std::endl;
2439 if ( totMass > p->mass() ) goto checkB;
2440 p->initializePhaseSpace( 2, mydaugs );
2441
2442 if ( !checkdecay( p ) ) goto checkB;
2443 vhds = p->getDaug( 0 )->getP4();
2444 ;
2445 vgam = p->getDaug( 1 )->getP4();
2446 }
2447 //-----
2448 }
2449 else
2450 { // with ISR photon for mode=74110
2451 Sampling_mhds:
2452 mhds = Mhad_sampling( MH, AF );
2453 // std::cout<<"SetMthr= "<<SetMthr<<std::endl;
2454 if ( mhds < SetMthr ) goto Sampling_mhds;
2455 double xeng = 1 - mhds * mhds / ( _cms * _cms ); // photon energy ratio
2456 double theta = ISR_ang_sampling( xeng ); // ISR photon polar angle in rad unit
2457
2458 if ( mydbg ) mass2 = mhds;
2459
2460 // generating events
2461 ModeSelection:
2462 EvtGlobalSet::ConExcPythia = 1;
2463 mymode = selectMode( vmod, mhds );
2464 if ( mymode == -10 ) goto Sampling_mhds;
2465 conexcmode = mymode; // save for model REXC ( see EvtRexc.cc) decay
2466 if ( mhds < 2.3 && mymode == 74110 )
2467 { goto ModeSelection; } // E<2.3 GeV, fully exclusive production
2468 std::cout << "A selected mode is " << mymode << " with Mhds= " << mhds
2469 << std::endl; // debugging
2470 _selectedMode = mymode;
2471 delete myxsection;
2472 myxsection = new EvtXsection( mymode );
2473 init_mode( mymode );
2474 EvtId myvpho = EvtPDL::getId( std::string( "vgam" ) );
2475 EvtPDL::reSetMass( myvpho, mhds ); // set to continum cms energy
2476 EvtPDL::reSetWidth( myvpho, 0 );
2477
2478 // debugging
2479 // for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<"
2480 // "<<EvtPDL::name(daugs[i])<<std::endl;}
2481
2482 // decay e+e- ->gamma_ISR gamma*
2483 EvtId mydaugs[2];
2484 if ( _ndaugs > 1 )
2485 { // gamma* -> exclusive decays
2486 resetResMass();
2487 mydaugs[0] = myvpho;
2488 }
2489 else
2490 { // vector resonance decays
2491 resetResMass();
2492 mydaugs[0] = daugs[0];
2493 EvtPDL::reSetMass( mydaugs[0], mhds );
2494 // EvtPDL::reSetWidth(mydaugs[0],0);
2495 }
2496 mydaugs[1] = gamId; // ISR photon
2497 int maxflag = 0;
2498 int trycount = 0;
2499 ISRphotonAng_sampling:
2500 double totMass = 0;
2501 p->makeDaughters( 2, mydaugs );
2502 for ( int i = 0; i < 2; i++ )
2503 {
2504 EvtParticle* di = p->getDaug( i );
2505 double xmi = EvtPDL::getMass( di->getId() );
2506 di->setMass( xmi );
2507 totMass += xmi;
2508 }
2509 if ( totMass > p->mass() ) goto ISRphotonAng_sampling;
2510 // std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2511 double weight1 = p->initializePhaseSpace( 2, mydaugs );
2512 if ( !checkdecay( p ) ) goto ISRphotonAng_sampling;
2513 // set the ISR photon angular distribution
2514 SetP4Rvalue( p, mhds, xeng, theta ); // end of decay e+e- -> vpho gamma_ISR
2515
2516 if ( maxflag == 0 )
2517 {
2518 findMaxXS( mhds );
2519 maxflag = 1;
2520 }
2521 vhds = p->getDaug( 0 )->getP4();
2522 vgam = p->getDaug( 1 )->getP4();
2523 double gx = vgam.get( 1 );
2524 double gy = vgam.get( 2 );
2525 double sintheta = sqrt( gx * gx + gy * gy ) / vgam.d3mag();
2526 bool gam_ang = gam_sampling( mhds, sintheta );
2527 trycount++;
2528
2529 } // with and without ISR sampling block
2530 // finish event generation
2531 // for debugging
2532 if ( mydbg )
2533 {
2534 pgam[0] = vgam.get( 0 );
2535 pgam[1] = vgam.get( 1 );
2536 pgam[2] = vgam.get( 2 );
2537 pgam[3] = vgam.get( 3 );
2538
2539 phds[0] = vhds.get( 0 );
2540 phds[1] = vhds.get( 1 );
2541 phds[2] = vhds.get( 2 );
2542 phds[3] = vhds.get( 3 );
2543 costheta = vgam.get( 3 ) / vgam.d3mag();
2544 selectmode = mymode;
2545 xs->Fill();
2546 // std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
2547 }
2548 _modeFlag[1] = _selectedMode;
2549 p->setIntFlag( _modeFlag );
2550 return;
2551 } // end block if(_mode==74110)
2552
2553 // ======================================================
2554 // multi-exclusive mode
2555 //=======================================================
2556 if ( _mode == -100 )
2557 {
2558 int mymode;
2559 double pm = EvtRandom::Flat( 0., 1 );
2560
2561 // std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
2562 //--
2563 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross
2564 // section within Egam=(0.025-Egamma^max GeV)
2565 if ( pm < _xs0 / ( _xs0 + _xs1 ) )
2566 { // without ISR photon for mode=74110
2567
2568 mhds = _cms;
2569 mymode = selectMode( vmd, mhds );
2570 _selectedMode = mymode;
2571 std::cout << "A selected mode is " << mymode << " with Mhds= " << mhds
2572 << std::endl; // debugging
2573 delete myxsection;
2574 myxsection = new EvtXsection( mymode );
2575 init_mode( mymode );
2576 resetResMass();
2577
2578 if ( _ndaugs > 1 )
2579 { // for e+e- -> exclusive decays
2580 checkAA:
2581 p->makeDaughters( _ndaugs, daugs );
2582 double totMass = 0;
2583 for ( int i = 0; i < _ndaugs; i++ )
2584 {
2585 EvtParticle* di = p->getDaug( i );
2586 double xmi = EvtPDL::getMass( di->getId() );
2587 di->setMass( xmi );
2588 totMass += xmi;
2589 }
2590 if ( totMass > p->mass() ) goto checkAA;
2591 p->initializePhaseSpace( _ndaugs, daugs );
2592 if ( !checkdecay( p ) ) goto checkAA;
2593 vhds = p->getP4();
2594 bool byn_ang;
2595 EvtVector4R pd1 = p->getDaug( 0 )->getP4();
2596 EvtVector4R pcm( _cms, 0, 0, 0 );
2597 // p->printTree();
2598 if ( _ndaugs == 2 )
2599 { // angular distribution for the hadrons only for two-body decays
2600 byn_ang = hadron_angle_sampling( pd1, pcm );
2601 if ( !byn_ang ) goto checkAA;
2602 }
2603 }
2604 else
2605 { // for e+e- -> vector resonace, add a very soft photon
2606 EvtId mydaugs[2];
2607 mydaugs[0] = daugs[0];
2608 EvtPDL::reSetMass( mydaugs[0], mhds - 0.001 );
2609 // EvtPDL::reSetWidth(mydaugs[0],0);
2610 mydaugs[1] = gamId; // ISR photon
2611 p->makeDaughters( 2, mydaugs );
2612 checkBA:
2613 double totMass = 0;
2614 for ( int i = 0; i < 2; i++ )
2615 {
2616 EvtParticle* di = p->getDaug( i );
2617 double xmi = EvtPDL::getMass( di->getId() );
2618 di->setMass( xmi );
2619 totMass += xmi;
2620 }
2621 // std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<"
2622 // "<<checkdecay(p)<<std::endl;
2623 if ( totMass > p->mass() ) goto checkBA;
2624 p->initializePhaseSpace( 2, mydaugs );
2625
2626 if ( !checkdecay( p ) ) goto checkBA;
2627 vhds = p->getDaug( 0 )->getP4();
2628 ;
2629 vgam = p->getDaug( 1 )->getP4();
2630 }
2631 //-----
2632 }
2633 else
2634 { // with ISR photon for mode=-100
2635 Sampling_mhds_A:
2636 mhds = Mhad_sampling( MH, AF );
2637 // std::cout<<"SetMthr= "<<SetMthr<<std::endl;
2638 if ( mhds < SetMthr ) goto Sampling_mhds_A;
2639 double xeng = 1 - mhds * mhds / ( _cms * _cms ); // photon energy ratio
2640 double theta = ISR_ang_sampling( xeng ); // ISR photon polar angle in rad unit
2641
2642 if ( mydbg ) mass2 = mhds;
2643
2644 // generating events
2645 ModeSelection_A:
2646 mymode = selectMode( vmd, mhds );
2647 if ( mymode == -10 ) goto Sampling_mhds_A;
2648 conexcmode = mymode; // save for model REXC ( see EvtRexc.cc) decay
2649 std::cout << "A selected mode is " << mymode << " with Mhds= " << mhds
2650 << std::endl; // debugging
2651 _selectedMode = mymode;
2652 delete myxsection;
2653 myxsection = new EvtXsection( mymode );
2654 init_mode( mymode );
2655 EvtId myvpho = EvtPDL::getId( std::string( "vgam" ) );
2656 EvtPDL::reSetMass( myvpho, mhds ); // set to continum cms energy
2657 EvtPDL::reSetWidth( myvpho, 0 );
2658
2659 // debugging
2660 // for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<"
2661 // "<<EvtPDL::name(daugs[i])<<std::endl;}
2662
2663 // decay e+e- ->gamma_ISR gamma*
2664 EvtId mydaugs[2];
2665 if ( _ndaugs > 1 )
2666 { // gamma* -> exclusive decays
2667 resetResMass();
2668 mydaugs[0] = myvpho;
2669 }
2670 else
2671 { // vector resonance decays
2672 resetResMass();
2673 mydaugs[0] = daugs[0];
2674 EvtPDL::reSetMass( mydaugs[0], mhds );
2675 // EvtPDL::reSetWidth(mydaugs[0],0);
2676 }
2677 mydaugs[1] = gamId; // ISR photon
2678 int maxflag = 0;
2679 int trycount = 0;
2680 ISRphotonAng_sampling_A:
2681 double totMass = 0;
2682 p->makeDaughters( 2, mydaugs );
2683 for ( int i = 0; i < 2; i++ )
2684 {
2685 EvtParticle* di = p->getDaug( i );
2686 double xmi = EvtPDL::getMass( di->getId() );
2687 di->setMass( xmi );
2688 totMass += xmi;
2689 }
2690 if ( totMass > p->mass() ) goto ISRphotonAng_sampling_A;
2691 // std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2692 double weight1 = p->initializePhaseSpace( 2, mydaugs );
2693 if ( !checkdecay( p ) ) goto ISRphotonAng_sampling_A;
2694 // set the ISR photon angular distribution
2695 SetP4Rvalue( p, mhds, xeng, theta ); // end of decay e+e- -> vpho gamma_ISR
2696
2697 //--debugging
2698 // p->printTree();
2699
2700 if ( maxflag == 0 )
2701 {
2702 findMaxXS( mhds );
2703 maxflag = 1;
2704 }
2705 vhds = p->getDaug( 0 )->getP4();
2706 vgam = p->getDaug( 1 )->getP4();
2707 double gx = vgam.get( 1 );
2708 double gy = vgam.get( 2 );
2709 double sintheta = sqrt( gx * gx + gy * gy ) / vgam.d3mag();
2710 bool gam_ang = gam_sampling( mhds, sintheta );
2711 trycount++;
2712
2713 } // with and without ISR sampling block
2714 _modeFlag[0] = -100;
2715 _modeFlag[1] = _selectedMode;
2716 p->setIntFlag( _modeFlag );
2717 // finish event generation
2718 // for debugging
2719 if ( mydbg )
2720 {
2721 pgam[0] = vgam.get( 0 );
2722 pgam[1] = vgam.get( 1 );
2723 pgam[2] = vgam.get( 2 );
2724 pgam[3] = vgam.get( 3 );
2725
2726 phds[0] = vhds.get( 0 );
2727 phds[1] = vhds.get( 1 );
2728 phds[2] = vhds.get( 2 );
2729 phds[3] = vhds.get( 3 );
2730 costheta = vgam.get( 3 ) / vgam.d3mag();
2731 selectmode = mymode;
2732 xs->Fill();
2733 // std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
2734 }
2735 return;
2736 } // end block if(_mode==-100)
2737
2738 //=======================================================
2739 // e+e- -> gamma_ISR gamma*
2740 //=======================================================
2741 if ( VISR )
2742 {
2743 EvtId gid = EvtPDL::getId( "gamma*" );
2744 double pm = EvtRandom::Flat( 0., 1 );
2745
2746 if ( pm < _xs0 / ( _xs0 + _xs1 ) )
2747 { // with a very soft ISR photon
2748 EvtPDL::reSetMass( gid, ( _cms - 0.0001 ) );
2749 mhds = _cms - 0.0001;
2750 }
2751 else
2752 {
2753 mhds = Mhad_sampling( MH, AF );
2754 EvtPDL::reSetMass( gid, mhds );
2755 }
2756 // debugging
2757 std::cout << "generatedMass " << EvtPDL::getMeanMass( gid ) << std::endl;
2758 double xeng = 1 - mhds * mhds / ( _cms * _cms ); // photon energy ratio
2759 double theta = ISR_ang_sampling( xeng ); // ISR photon polar angle in rad unit
2760 p->makeDaughters( 2, daugs );
2761 for ( int i = 0; i < 2; i++ )
2762 {
2763 EvtParticle* di = p->getDaug( i );
2764 // if(EvtPDL::name(di->getId())=="gamma*") di->setVectorSpinDensity();
2765 double xmi = EvtPDL::getMeanMass( di->getId() );
2766 di->setMass( xmi );
2767 }
2768 p->initializePhaseSpace( 2, daugs );
2769 SetP4( p, mhds, xeng, theta );
2770 return;
2771 }
2772
2773 //========================================================
2774 //=== for other mode generation with _mode != 74110
2775 //========================================================
2776
2777 if ( ( _xs0 + _xs1 ) == 0 )
2778 {
2779 std::cout << "EvtConExc::zero total xsection" << std::endl;
2780 ::abort();
2781 }
2782 double pm = EvtRandom::Flat( 0., 1 );
2783 double xsratio = _xs0 / ( _xs0 + _xs1 );
2784 int iflag = 2; // flag = 0 for ee->hadrons, 1 for ee->gamma_ISR hadrons, 2: mix 0 and 1 case
2785 if ( getArg( 0 ) != -2 )
2786 { // exclude DIY case
2787 if ( getNArg() == 3 && radflag == 1 )
2788 {
2789 iflag = 1;
2790 xsratio = 0;
2791 } // only generating gamma hadrons mode
2792 else if ( getNArg() == 3 && radflag == 0 )
2793 {
2794 iflag = 0;
2795 xsratio = 1;
2796 }
2797 }
2798
2799 // std::cout<<"xsratio= "<<xsratio<<std::endl;
2800
2801 if ( pm < xsratio )
2802 { // no ISR photon
2803 masscheck:
2804 double summassx = 0;
2805 p->makeDaughters( _ndaugs, daugs );
2806 for ( int i = 0; i < _ndaugs; i++ )
2807 {
2808 EvtParticle* di = p->getDaug( i );
2809 double xmi = EvtPDL::getMass( di->getId() );
2810 di->setMass( xmi );
2811 summassx += xmi;
2812 // std::cout<<"PID and mass "<<di->getId()<<" "<<xmi<<std::endl;
2813 }
2814 if ( summassx > p->mass() ) goto masscheck;
2815 angular_hadron:
2816 p->initializePhaseSpace( _ndaugs, daugs );
2817 bool byn_ang;
2818 EvtVector4R pd1 = p->getDaug( 0 )->getP4();
2819 EvtVector4R pcm( _cms, 0, 0, 0 );
2820 if ( _ndaugs == 2 )
2821 { // angular distribution for the hadrons only for two-body decays
2822 byn_ang = hadron_angle_sampling( pd1, pcm );
2823 if ( !byn_ang ) goto angular_hadron;
2824 }
2825
2826 // for histogram
2827 cosp = pd1.get( 3 ) / pd1.d3mag();
2828 mhds = _cms;
2829 // std::cout<<"cosp, mhds "<<cosp<<" "<<mhds<<std::endl;
2830 // p->printTree();
2831 }
2832 else
2833 { // sampling m_hadrons and decay e+e- ->gamma gamma*
2834 double mhdr = Mhad_sampling( MH, AF );
2835 double xeng = 1 - mhdr * mhdr / ( _cms * _cms ); // photon energy ratio
2836 double theta = ISR_ang_sampling( xeng ); // ISR photon polar angle in rad unit
2837 EvtId gid = EvtPDL::getId( "vhdr" );
2838 EvtPDL::reSetMass( gid, mhdr );
2839 int ndaugs = 2;
2840 EvtId mydaugs[2];
2841 mydaugs[0] = EvtPDL::getId( std::string( "gamma" ) );
2842 mydaugs[1] = EvtPDL::getId( std::string( "vhdr" ) );
2843
2844 // for histogram
2845 mhds = mhdr;
2846 costheta = cos( theta );
2847 //--
2848
2849 p->makeDaughters( 2, mydaugs );
2850 for ( int i = 0; i < 2; i++ )
2851 {
2852 EvtParticle* di = p->getDaug( i );
2853 double xmi = EvtPDL::getMass( di->getId() );
2854 // if(EvtPDL::name(di->getId())=="vhdr") di->setVectorSpinDensity();
2855 di->setMass( xmi );
2856 }
2857 p->initializePhaseSpace( 2, mydaugs );
2858 SetP4( p, mhdr, xeng, theta ); // end of decay e+e- -> gamma_ISR gamma*
2859 // p->printTree();
2860
2861 //----
2862 } // end of gamma gamma*, gamma*->hadrons generation
2863 // p->printTree();
2864
2865 //-----------------------------------------
2866 // End of decays for all topology
2867 //----------------------------------------
2868 //================================= event analysis
2869 if ( _nevt == 0 )
2870 {
2871 std::cout << "The decay chain: " << std::endl;
2872 p->printTree();
2873 }
2874 _nevt++;
2875 //--- for debuggint to make root file
2876 if ( mydbg ) { xs->Fill(); }
2877
2878 //----
2879 return;
2880}
2881
2882bool EvtConExc::hadron_angle_sampling( EvtVector4R ppi, EvtVector4R pcm ) {
2883 EvtVector4R pbst = -1 * pcm;
2884 pbst.set( 0, pcm.get( 0 ) );
2885 EvtVector4R p4i = boostTo( ppi, pbst );
2886 if ( ( _mode >= 0 && _mode <= 5 ) || _mode == 44 || _mode == 96 )
2887 { // ee->two baryon mode, VP, SP, mode=-2 is excluded
2888 bool byn_ang = baryon_sampling( pcm, p4i );
2889 return byn_ang;
2890 }
2891 else if ( _mode == 6 || _mode == 45 || _mode == 46 || _mode == 70 || _mode == 71 ||
2892 _mode == 93 )
2893 { // ee->PP (pi+pi-,..) mode
2894 bool msn_ang = meson_sampling( pcm, p4i );
2895 return msn_ang;
2896 }
2897 else if ( _mode == 23 || _mode == 24 || _mode == 25 || _mode == 26 || _mode == 27 ||
2898 _mode == 36 || _mode == 39 || _mode == 47 || _mode == 48 || _mode == 52 ||
2899 _mode == 68 || _mode == 69 || _mode == 72 || _mode == 73 || _mode == 80 ||
2900 _mode == 94 || _mode == 95 )
2901 {
2902 bool msn_ang = VP_sampling( pcm, p4i );
2903 return msn_ang;
2904 }
2905 else if ( _mode == -2 )
2906 {
2907 EvtSpinType::spintype type0 = EvtPDL::getSpinType( daugs[0] );
2908 EvtSpinType::spintype type1 = EvtPDL::getSpinType( daugs[1] );
2909 if ( type0 == EvtSpinType::SCALAR && type1 == EvtSpinType::SCALAR )
2910 {
2911 bool msn_ang = meson_sampling( pcm, p4i );
2912 return msn_ang;
2913 }
2914 else if ( type0 == EvtSpinType::VECTOR && type1 == EvtSpinType::SCALAR ||
2915 type1 == EvtSpinType::VECTOR && type0 == EvtSpinType::SCALAR )
2916 {
2917 bool msn_ang = VP_sampling( pcm, p4i );
2918 return msn_ang;
2919 }
2920 }
2921 return true;
2922}
2923
2924double EvtConExc::gamHXSection( EvtParticle* p, double El, double Eh,
2925 int nmc ) { // El, Eh : the lower and higher energy of
2926 // radiative photons
2927 // std::cout<<"nmc= "<<nmc<<std::endl;
2928 gamId = EvtPDL::getId( std::string( "gamma" ) );
2929 EvtId xvp = EvtPDL::getId( std::string( "xvpho" ) );
2930 EvtVector4R p_init( p->getP4().mass(), 0., 0., 0. );
2931 double totxs = 0;
2932 maxXS = -1;
2933 _er1 = 0;
2934 Rad2Xs = 0;
2935 for ( int ii = 0; ii < nmc; ii++ )
2936 { // integarted the gamma hadrons xsection
2937 gamH = EvtParticleFactory::particleFactory( xvp, p_init );
2938 int gamdaugs = _ndaugs + 1;
2939
2940 EvtId theDaugs[10];
2941 for ( int i = 0; i < _ndaugs; i++ ) { theDaugs[i] = daugs[i]; }
2942 theDaugs[_ndaugs] = gamId;
2943
2944 gamH->makeDaughters( gamdaugs, theDaugs );
2945
2946 for ( int i = 0; i < gamdaugs; i++ )
2947 {
2948 EvtParticle* di = gamH->getDaug( i );
2949 double xmi = EvtPDL::getMass( di->getId() );
2950 di->setMass( xmi );
2951 }
2952 loop:
2953 double weight = gamH->initializePhaseSpace( gamdaugs, theDaugs );
2954 //-- slection:theta_gamma > 1 degree
2955 EvtVector4R pgam = gamH->getDaug( _ndaugs )->getP4();
2956 double pmag = pgam.d3mag();
2957 double pz = pgam.get( 3 );
2958 double sintheta = sqrt( pmag * pmag - pz * pz ) / pmag;
2959 double onedegree = 1. / 180. * 3.1415926;
2960 // if(sintheta < onedegree) {goto loop;}
2961 if ( pmag < El || pmag > Eh ) { goto loop; }
2962 //--------
2963 // std::cout<<"pmag= "<<pmag<<std::endl;
2964
2965 double thexs = difgamXs( gamH ); // prepare the photon angular distribution
2966 Rad2Xs += Rad2difXs( gamH );
2967 if ( thexs > maxXS ) { maxXS = thexs; }
2968 double s = p_init.mass2();
2969 double x = 2 * pgam.get( 0 ) / sqrt( s );
2970 double rad1xs = Rad1difXs( gamH ); // first order xs by Eq(4)hep-ph/9910523
2971 totxs += rad1xs;
2972 _er1 += differ;
2973 gamH->deleteDaughters();
2974 } // for int_i block
2975 _er1 /= nmc;
2976
2977 Rad2Xs /= nmc; // the leading second order correction
2978 totxs /= nmc; // first order correction XS
2979
2980 // return totxs; // first order correction XS
2981 return Rad2Xs; // the leading second order correction
2982}
2983
2984double EvtConExc::gamHXSection( double s, double El, double Eh,
2985 int nmc ) { // El, Eh : the lower and higher energy of
2986 // radiative photons
2987 //--for narrow resonance psi(2S),J/psi, phi, which can not calculated with the MC integration
2988 // double mxL = sqrt( s-2*Eh*sqrt(s) ); //low mass
2989 // double mxH = sqrt( s-2*El*sqrt(s) ); //high mass
2990 // double sigv = narrowRXS(mxL,mxH);
2991 //--
2992
2993 double totxs = 0;
2994 maxXS = -1;
2995 _er1 = 0;
2996 Rad2Xs = 0;
2997 double xL = 2 * El / sqrt( s );
2998 double xH = 2 * Eh / sqrt( s );
2999 for ( int i = 0; i < nmc; i++ )
3000 { // It is found that the narrow resonance can not included in this MC integartion
3001 double rdm = EvtRandom::Flat( 0., 1. ); // set angular cut 1^o
3002 double x = xL + ( xH - xL ) * rdm;
3003 Rad2Xs += Rad2difXs( s, x );
3004 _er1 += differ2; // stored when calculate Rad2Xs
3005 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
3006 }
3007 _er1 /= nmc;
3008 _er1 *= ( xH - xL );
3009 // std::cout<<"_er1= "<<_er1<<std::endl;
3010 Rad2Xs /= nmc; // the leading second order correction
3011 double thexs = Rad2Xs * ( xH - xL ); // xH-xL is the Jacobi factor
3012 // std::cout<<"thexs= "<<thexs<<std::endl;
3013 // if(sigv != -1) thexs += sigv; //add narrow resonance ISR cross section
3014 return thexs;
3015}
3016
3017double EvtConExc::gamHXSection( double El, double Eh ) { // using Gaussian integration
3018 // subroutine from Cern lib
3019 std::cout << " " << std::endl;
3020 extern double Rad2difXs( double* x );
3021 extern double Rad2difXs2( double* x );
3022 double eps = 0.01;
3023 double Rad2Xs;
3024 if ( _mode != -2 ) { Rad2Xs = dgauss_( &Rad2difXs, &El, &Eh, &eps ); }
3025 else { Rad2Xs = dgauss_( &Rad2difXs2, &El, &Eh, &eps ); }
3026 if ( Rad2Xs == 0 )
3027 {
3028 for ( int iter = 0; iter < 10; iter++ )
3029 {
3030 eps += 0.01;
3031 if ( _mode != -2 ) { Rad2Xs = dgauss_( &Rad2difXs, &El, &Eh, &eps ); }
3032 else { Rad2Xs = dgauss_( &Rad2difXs2, &El, &Eh, &eps ); }
3033 if ( !Rad2Xs ) return Rad2Xs;
3034 }
3035 }
3036 return Rad2Xs; // the leading second order correction
3037}
3038
3039double EvtConExc::gamHXSection( double El, double Eh,
3040 int mode ) { // using Gaussian integration subroutine from Cern
3041 // lib
3042 std::cout << " " << std::endl;
3043 extern double Rad2difXs( double* x );
3044 extern double Rad2difXs2( double* x );
3045 double eps = 0.01;
3046 double Rad2Xs;
3047 if ( mode != -2 ) { Rad2Xs = dgauss_( &Rad2difXs, &El, &Eh, &eps ); }
3048 else { Rad2Xs = dgauss_( &Rad2difXs2, &El, &Eh, &eps ); }
3049 if ( Rad2Xs == 0 )
3050 {
3051 for ( int iter = 0; iter < 10; iter++ )
3052 {
3053 eps += 0.01;
3054 if ( mode != -2 ) { Rad2Xs = dgauss_( &Rad2difXs, &El, &Eh, &eps ); }
3055 else { Rad2Xs = dgauss_( &Rad2difXs2, &El, &Eh, &eps ); }
3056 if ( !Rad2Xs ) return Rad2Xs;
3057 }
3058 }
3059 return Rad2Xs; // the leading second order correction
3060}
3061
3062double EvtConExc::gamHXSection_er( double El, double Eh ) { // using Gaussian integration
3063 // subroutine from Cern lib
3064 std::cout << " " << std::endl;
3065 extern double Rad2difXs_er( double* x );
3066 extern double Rad2difXs_er2( double* x );
3067 double eps = 0.01;
3068 double Rad2Xs;
3069 if ( _mode != -2 ) { Rad2Xs = dgauss_( &Rad2difXs_er, &El, &Eh, &eps ); }
3070 else { Rad2Xs = dgauss_( &Rad2difXs_er2, &El, &Eh, &eps ); }
3071 if ( Rad2Xs == 0 )
3072 {
3073 for ( int iter = 0; iter < 10; iter++ )
3074 {
3075 eps += 0.01;
3076 if ( _mode != -2 ) { Rad2Xs = dgauss_( &Rad2difXs_er, &El, &Eh, &eps ); }
3077 else { Rad2Xs = dgauss_( &Rad2difXs_er2, &El, &Eh, &eps ); }
3078 if ( !Rad2Xs ) return Rad2Xs;
3079 ;
3080 }
3081 }
3082 return Rad2Xs; // the leading second order correction
3083}
3084
3085void EvtConExc::findMaxXS( EvtParticle* p ) {
3086 // std::cout<<"nmc= "<<nmc<<std::endl;
3087 gamId = EvtPDL::getId( std::string( "gamma" ) );
3088 EvtId xvp = EvtPDL::getId( std::string( "xvpho" ) );
3089 EvtVector4R p_init( p->getP4().mass(), 0., 0., 0. );
3090 double totxs = 0;
3091 maxXS = -1;
3092 _er1 = 0;
3093 Rad2Xs = 0;
3094 int nmc = 50000;
3095 for ( int ii = 0; ii < nmc; ii++ )
3096 { // integarted the gamma hadrons xsection
3097 gamH = EvtParticleFactory::particleFactory( xvp, p_init );
3098 int gamdaugs = _ndaugs + 1;
3099
3100 EvtId theDaugs[10];
3101 for ( int i = 0; i < _ndaugs; i++ ) { theDaugs[i] = daugs[i]; }
3102 theDaugs[_ndaugs] = gamId;
3103
3104 gamH->makeDaughters( gamdaugs, theDaugs );
3105 loop:
3106 double totm = 0;
3107 for ( int i = 0; i < gamdaugs; i++ )
3108 {
3109 EvtParticle* di = gamH->getDaug( i );
3110 double xmi = EvtPDL::getMass( di->getId() );
3111 di->setMass( xmi );
3112 totm += xmi;
3113 }
3114
3115 // std::cout<<totm<<" "<<p_init.get(0)<<std::endl;
3116 if ( totm >= p_init.get( 0 ) ) goto loop;
3117
3118 double weight = gamH->initializePhaseSpace( gamdaugs, theDaugs );
3119 double thexs = difgamXs( gamH ); // prepare the photon angular distribution
3120 EvtVector4R p4gam = gamH->getDaug( _ndaugs )->getP4();
3121 double costheta = p4gam.get( 3 ) / p4gam.d3mag();
3122 double sintheta = sqrt( 1 - costheta * costheta );
3123 bool acut = ( sintheta > 0.04 ) && ( sintheta < 0.96 ); // set photon auglar cut 2^o
3124 if ( thexs > maxXS && acut ) { maxXS = thexs; }
3125 // gamH->deleteTree();
3126 }
3127}
3128
3129double EvtConExc::difgamXs( EvtParticle* p ) { // dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3130 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
3131 EvtVector4R p0 = p->getDaug( 0 )->getP4();
3132 for ( int i = 1; i < _ndaugs; i++ ) { p0 += p->getDaug( i )->getP4(); }
3133 EvtVector4R pgam = p->getDaug( _ndaugs )->getP4();
3134 double mhs = p0.mass();
3135 double egam = pgam.get( 0 );
3136 double sin2theta = pgam.get( 3 ) / pgam.d3mag();
3137 sin2theta = 1 - sin2theta * sin2theta;
3138
3139 double cms = p->getP4().mass();
3140 double alpha = 1. / 137.;
3141 double pi = 3.1415926;
3142 double x = 2 * egam / cms;
3143 double wsx = alpha / pi / x * ( ( 2 - 2 * x + x * x ) / sin2theta - x * x / 2. );
3144 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<",
3145 // "<<sin2theta<<", "<<cms<<",
3146 // "<<x<<", "<<wsx<<std::endl;
3147
3148 double xs = myxsection->getXsection( mhs );
3149 double difxs = 2 * mhs / cms / cms * wsx * xs;
3150 differ = 2 * mhs / cms / cms * wsx * ( myxsection->getErr( mhs ) );
3151 return difxs;
3152}
3153
3154bool EvtConExc::gam_sampling( EvtParticle* p ) { // photon angular distribution sampling
3155 double pm = EvtRandom::Flat( 0., 1 );
3156 double xs = difgamXs( p );
3157 double xsratio = xs / maxXS;
3158 if ( pm < xsratio ) { return true; }
3159 else { return false; }
3160}
3161
3162bool EvtConExc::gam_sampling( double mhds, double sintheta ) {
3163 double pm = EvtRandom::Flat( 0., 1 );
3164 double xs = difgamXs( mhds, sintheta );
3165 double xsratio = xs / maxXS;
3166 if ( pm < xsratio ) { return true; }
3167 else { return false; }
3168}
3169
3170bool EvtConExc::xs_sampling( double xs ) {
3171 double pm = EvtRandom::Flat( 0., 1 );
3172 // std::cout<<"Rdm= "<<pm<<std::endl;
3173 if ( pm < xs / XS_max ) { return true; }
3174 else { return false; }
3175}
3176
3177bool EvtConExc::xs_sampling( double xs, double xs0 ) {
3178 double pm = EvtRandom::Flat( 0., 1 );
3179 // std::cout<<"Rdm= "<<pm<<std::endl;
3180 if ( pm < xs / ( xs0 * 1.1 ) ) { return true; }
3181 else { return false; }
3182}
3183
3184bool EvtConExc::baryon_sampling( EvtVector4R pcm, EvtVector4R pi ) {
3185 EvtHelSys angles( pcm, pi ); // using helicity sys.angles
3186 double theta = angles.getHelAng( 1 );
3187 double phi = angles.getHelAng( 2 );
3188 double costheta = cos( theta ); // using helicity angles in parent system
3189
3190 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
3191 // std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
3192 double alpha = baryonAng( _cms );
3193 if ( _mode == 96 ) { alpha = -0.34; }
3194 double pm = EvtRandom::Flat( 0., 1 );
3195 double ang = 1 + alpha * costheta * costheta;
3196 double myang;
3197 if ( alpha >= 0 ) { myang = 1 + alpha; }
3198 else { myang = 1; }
3199 if ( pm < ang / myang ) { return true; }
3200 else { return false; }
3201}
3202
3203bool EvtConExc::meson_sampling( EvtVector4R pcm, EvtVector4R pi ) {
3204 EvtHelSys angles( pcm, pi ); // using helicity sys.angles
3205 double theta = angles.getHelAng( 1 );
3206 double phi = angles.getHelAng( 2 );
3207 double costheta = cos( theta ); // using helicity angles in parent system
3208
3209 double pm = EvtRandom::Flat( 0., 1 );
3210 double ang = 1 - costheta * costheta;
3211 if ( pm < ang / 1. ) { return true; }
3212 else { return false; }
3213}
3214
3215bool EvtConExc::VP_sampling( EvtVector4R pcm, EvtVector4R pi ) {
3216 EvtHelSys angles( pcm, pi ); // using helicity sys.angles
3217 double theta = angles.getHelAng( 1 );
3218 double phi = angles.getHelAng( 2 );
3219 double costheta = cos( theta ); // using helicity angles in parent system
3220
3221 double pm = EvtRandom::Flat( 0., 1 );
3222 double ang = 1 + costheta * costheta;
3223 if ( pm < ang / 2. ) { return true; }
3224 else { return false; }
3225}
3226
3227double EvtConExc::Rad1( double s, double x ) { // first order correction
3228 // radiator correction upto second order, see Int. J. of Moder.Phys. A
3229 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
3230 double alpha = 1. / 137.;
3231 double pi = 3.1415926;
3232 double me = 0.5 * 0.001; // e mass
3233 double sme = sqrt( s ) / me;
3234 double L = 2 * log( sme );
3235 double wsx = 2 * alpha / pi / x * ( L - 1 ) * ( 1 - x + x * x / 2 );
3236 return wsx;
3237}
3238
3239double EvtConExc::Rad2( double s, double x ) {
3240 // radiator correction upto second order, see Int. J. of Moder.Phys. A, hep-ph/9910523, Eq(8)
3241 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
3242 double alpha = 1. / 137.;
3243 double pi = 3.1415926;
3244 double me = 0.5 * 0.001; // e mass
3245 double xi2 = 1.64493407;
3246 double xi3 = 1.2020569;
3247 double sme = sqrt( s ) / me;
3248 double L = 2 * log( sme );
3249 double beta = 2 * alpha / pi * ( L - 1 );
3250 double delta2 = ( 9. / 8. - 2 * xi2 ) * L * L - ( 45. / 16. - 5.5 * xi2 - 3 * xi3 ) * L -
3251 6. / 5. * xi2 * xi2 - 4.5 * xi3 - 6 * xi2 * log( 2. ) + 3. / 8. * xi2 +
3252 57. / 12.;
3253 double ap = alpha / pi;
3254 double Delta = 1 + ap * ( 1.5 * L + 1. / 3. * pi * pi - 2 ) + ap * ap * delta2;
3255 double wsx = Delta * beta * pow( x, beta - 1 ) - 0.5 * beta * ( 2 - x );
3256 double wsx2 = ( 2 - x ) * ( 3 * log( 1 - x ) - 4 * log( x ) ) - 4 * log( 1 - x ) / x - 6 + x;
3257 wsx = wsx + beta * beta / 8. * wsx2;
3258 double mymx = sqrt( s * ( 1 - x ) );
3259 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
3260 return wsx * getVP( mymx ); // vph is vaccum polarization factor
3261}
3262
3263double EvtConExc::Rad2difXs( EvtParticle* p ) { // leading second order xsection
3264 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
3265 double summass = p->getDaug( 0 )->getP4().mass();
3266 EvtVector4R v4p = p->getDaug( 0 )->getP4();
3267 for ( int i = 1; i < _ndaugs; i++ )
3268 {
3269 summass += p->getDaug( i )->getP4().mass();
3270 v4p += p->getDaug( i )->getP4();
3271 }
3272
3273 double Egam = p->getDaug( _ndaugs )->getP4().d3mag();
3274 double cms = p->getP4().mass();
3275 double mrg = cms - summass;
3276 double pm = EvtRandom::Flat( 0., 1 );
3277 double mhs = pm * mrg + summass;
3278
3279 double s = cms * cms;
3280 double x = 2 * Egam / cms;
3281 // double mhs = sqrt( s*(1-x) );
3282 double wsx = Rad2( s, x );
3283
3284 // std::cout<<"random : "<<pm<<std::endl;
3285
3286 double xs = myxsection->getXsection( mhs );
3287 if ( xs > XS_max ) { XS_max = xs; }
3288 double difxs = 2 * mhs / cms / cms * wsx * xs;
3289 differ2 = 2 * mhs / cms / cms * wsx * ( myxsection->getErr( mhs ) );
3290 differ *= mrg; // Jacobi factor for variable m_hds
3291 difxs *= mrg;
3292 return difxs;
3293}
3294
3295double EvtConExc::Rad2difXs( double s, double x ) { // leading second order xsection
3296 double wsx = Rad2( s, x );
3297 double mhs = sqrt( s * ( 1 - x ) );
3298 double xs = myxsection->getXsection( mhs );
3299 // if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
3300 if ( xs > XS_max ) { XS_max = xs; }
3301 double difxs = wsx * xs;
3302 differ2 = wsx * ( myxsection->getErr( mhs ) );
3303 return difxs;
3304}
3305
3306double EvtConExc::Rad2difXs( double s, double x,
3307 EvtXsection* myxsection ) { // leading second order xsection
3308 double wsx = Rad2( s, x );
3309 double mhs = sqrt( s * ( 1 - x ) );
3310 double xs = myxsection->getXsection( mhs );
3311 // if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
3312 if ( xs > XS_max ) { XS_max = xs; }
3313 double difxs = wsx * xs;
3314 differ2 = wsx * ( myxsection->getErr( mhs ) );
3315 return difxs;
3316}
3317
3318double EvtConExc::Rad1difXs( EvtParticle* p ) { // // first order xsection
3319 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
3320 double summass = p->getDaug( 0 )->getP4().mass();
3321 for ( int i = 1; i < _ndaugs; i++ ) { summass += p->getDaug( i )->getP4().mass(); }
3322
3323 double cms = p->getP4().mass();
3324 double mrg = cms - summass;
3325 double pm = EvtRandom::Flat( 0., 1 );
3326 double mhs = pm * mrg + summass;
3327
3328 double s = cms * cms;
3329 double x = 1 - mhs * mhs / s;
3330 double wsx = Rad1( s, x );
3331
3332 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<",
3333 // "<<sin2theta<<", "<<cms<<",
3334 // "<<x<<", "<<wsx<<std::endl;
3335
3336 double xs = myxsection->getXsection( mhs );
3337 if ( xs > XS_max ) { XS_max = xs; }
3338 double difxs = 2 * mhs / cms / cms * wsx * xs;
3339
3340 differ = 2 * mhs / cms / cms * wsx * ( myxsection->getErr( mhs ) );
3341 differ *= mrg; // Jacobi factor for variable m_hds
3342 difxs *= mrg;
3343 return difxs;
3344}
3345
3346void EvtConExc::init_Br_ee() {
3347 // double psipp_ee=9.6E-06;
3348 double psip_ee = 7.73E-03;
3349 double jsi_ee = 5.94 * 0.01;
3350 double phi_ee = 2.954E-04;
3351 // double omega_ee=7.28E-05;
3352 // double rho_ee = 4.72E-05;
3353 EvtId psppId = EvtPDL::getId( std::string( "psi(3770)" ) );
3354 EvtId pspId = EvtPDL::getId( std::string( "psi(2S)" ) );
3355 EvtId jsiId = EvtPDL::getId( std::string( "J/psi" ) );
3356 EvtId phiId = EvtPDL::getId( std::string( "phi" ) );
3357 EvtId omegaId = EvtPDL::getId( std::string( "omega" ) );
3358 EvtId rhoId = EvtPDL::getId( std::string( "rho0" ) );
3359 BR_ee.clear();
3360 ResId.clear();
3361
3362 // BR_ee.push_back(rho_ee);
3363 // BR_ee.push_back(omega_ee);
3364 BR_ee.push_back( phi_ee );
3365 BR_ee.push_back( jsi_ee );
3366 BR_ee.push_back( psip_ee );
3367 // BR_ee.push_back(psipp_ee);
3368
3369 // ResId.push_back(rhoId);
3370 // ResId.push_back(omegaId);
3371 ResId.push_back( phiId );
3372 ResId.push_back( jsiId );
3373 ResId.push_back( pspId );
3374 // ResId.push_back(psppId);
3375}
3376
3377double EvtConExc::Ros_xs( double mx, double bree, EvtId pid ) { // in unit nb
3378 double pi = 3.1415926;
3379 double s = mx * mx;
3380 double width = EvtPDL::getWidth( pid );
3381 double mass = EvtPDL::getMeanMass( pid );
3382 double xv = 1 - mass * mass / s;
3383 double sigma = 12 * pi * pi * bree * width / mass / s;
3384 sigma *= Rad2( s, xv );
3385 double nbar = 3.89379304 * 100000; // ! GeV-2 = 3.89379304*10^5 nbar
3386 return sigma * nbar;
3387}
3388
3389double Rad2difXs( double* mhs ) { // leading second order xsection, mhs: the mass of final
3390 // states
3391 double cms = EvtConExc::_cms;
3392 double s = cms * cms;
3393 double x = 1 - ( *mhs ) * ( *mhs ) / s;
3394 EvtConExc myconexc;
3395 double wsx;
3396 double dhs = ( *mhs );
3397 double xs = EvtConExc::myxsection->getXsection( dhs );
3398 wsx = myconexc.Rad2( s, x );
3399 if ( xs > EvtConExc::XS_max ) { EvtConExc::XS_max = xs; }
3400 double difxs = 2 * dhs / cms / cms * wsx * xs;
3401 std::ofstream oa;
3402 oa << x << std::endl; // without this line, the running will breaking, I don't know why
3403 return difxs;
3404}
3405double Rad2difXs_er( double* mhs ) { // leading second order xsection, mhs: the mass of final
3406 // states
3407 double cms = EvtConExc::_cms;
3408 double s = cms * cms;
3409 double x = 1 - ( *mhs ) * ( *mhs ) / s;
3410 EvtConExc myconexc;
3411 double wsx;
3412 double xs = EvtConExc::myxsection->getErr( *mhs );
3413 wsx = myconexc.Rad2( s, x );
3414 double differ2 = 2 * ( *mhs ) / cms / cms * wsx * xs;
3415 std::ofstream oa;
3416 oa << x << std::endl;
3417 return differ2;
3418}
3419
3420double Rad2difXs2( double* mhs ) { // leading second order xsection, mhs: the mass of final
3421 // states
3422 double cms = EvtConExc::_cms;
3423 double s = cms * cms;
3424 double x = 1 - ( *mhs ) * ( *mhs ) / s;
3425 EvtConExc myconexc;
3426 double wsx;
3427 double dhs = ( *mhs );
3428 double xs = EvtConExc::myxsection->getXsection( dhs );
3429 wsx = myconexc.Rad2( s, x );
3430 if ( xs > EvtConExc::XS_max ) { EvtConExc::XS_max = xs; }
3431 double difxs = 2 * dhs / cms / cms * wsx * xs;
3432 oa << x << std::endl; // without this line, the running will breaking, I don't know why
3433 return difxs;
3434}
3435
3436double Rad2difXs_er2( double* mhs ) { // leading second order xsection, mhs: the mass of final
3437 // states
3438 double cms = EvtConExc::_cms;
3439 double s = cms * cms;
3440 double x = 1 - ( *mhs ) * ( *mhs ) / s;
3441 EvtConExc myconexc;
3442 double wsx;
3443 double xs = EvtConExc::myxsection->getErr( *mhs );
3444 wsx = myconexc.Rad2( s, x );
3445 double differ2 = 2 * ( *mhs ) / cms / cms * wsx * xs;
3446 oa << x << std::endl;
3447 return differ2;
3448}
3449
3450double EvtConExc::SoftPhoton_xs( double s, double b ) {
3451 double alpha = 1. / 137.;
3452 double pi = 3.1415926;
3453 double me = 0.5 * 0.001; // e mass
3454 double xi2 = 1.64493407;
3455 double xi3 = 1.2020569;
3456 double sme = sqrt( s ) / me;
3457 double L = 2 * log( sme );
3458 double beta = 2 * alpha / pi * ( L - 1 );
3459 double delta2 = ( 9. / 8. - 2 * xi2 ) * L * L - ( 45. / 16. - 5.5 * xi2 - 3 * xi3 ) * L -
3460 6. / 5. * xi2 * xi2 - 4.5 * xi3 - 6 * xi2 * log( 2. ) + 3. / 8. * xi2 +
3461 57. / 12.;
3462 double ap = alpha / pi;
3463 double Delta = 1 + ap * ( 1.5 * L + 1. / 3. * pi * pi - 2 ) + ap * ap * delta2;
3464
3465 double beta2 = beta * beta;
3466 double b2 = b * b;
3467
3468 double xs =
3469 ( -32 * b * beta + 8 * pow( b, 2 ) * beta - 10 * b * pow( beta, 2 ) +
3470 pow( b, 2 ) * pow( beta, 2 ) + 32 * pow( b, beta ) * Delta -
3471 6 * ( 3 - 4 * b + pow( b, 2 ) ) * pow( beta, 2 ) * log( 1 - b ) -
3472 32 * b * pow( beta, 2 ) * log( b ) + 8 * pow( b, 2 ) * pow( beta, 2 ) * log( b ) +
3473 16 * pow( beta, 2 ) * Li2( b ) ) /
3474 32.;
3475 double mymx = sqrt( s * ( 1 - b ) );
3476 // return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
3477 return xs * getVP( mymx ); // vph :the vaccum polarzation factor
3478}
3479
3480double EvtConExc::Li2( double x ) {
3481 // double li2= -x +x*x/4. - x*x*x/9; // wangwp: may be not correct!
3482 double li2 = +x + x * x / 4. + x * x * x / 9; // corrected by wangwp
3483 return li2;
3484}
3485
3486double EvtConExc::lgr( double* x, double* y, int n, double t ) {
3487 int n0 = -1;
3488 double z;
3489 for ( int i = 0; i < n - 1; i++ )
3490 {
3491 if ( x[i] >= t && t > x[i + 1] )
3492 {
3493 n0 = i;
3494 break;
3495 }
3496 }
3497 if ( n0 == -1 ) { return 0.0; }
3498 else
3499 {
3500 double k = ( y[n0] - y[n0 + 1] ) / ( x[n0] - x[n0 + 1] );
3501 z = y[n0 + 1] + k * ( t - x[n0 + 1] );
3502 }
3503 // std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
3504 // std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
3505 return z;
3506}
3507
3508bool EvtConExc::islgr( double* x, double* y, int n, double t ) {
3509 int n0 = -1;
3510 double z;
3511 for ( int i = 0; i < n - 1; i++ )
3512 {
3513 if ( x[i] >= t && t > x[i + 1] )
3514 {
3515 n0 = i;
3516 break;
3517 }
3518 }
3519 double xstotal = y[599];
3520 if ( n0 == -1 ) { return 0; }
3521 else
3522 {
3523 double p1 = y[n0] / xstotal;
3524 double p2 = y[n0 + 1] / xstotal;
3525 double pm = EvtRandom::Flat( 0., 1 );
3526 if ( p1 < pm && pm <= p2 ) { return 1; }
3527 else { return 0; }
3528 }
3529}
3530
3531double EvtConExc::LLr( double* x, double* y, int n, double t ) {
3532 int n0 = -1;
3533 double z;
3534 if ( t == x[n - 1] ) return y[n - 1];
3535 for ( int i = 0; i < n - 1; i++ )
3536 {
3537 if ( x[i] <= t && t < x[i + 1] )
3538 {
3539 n0 = i;
3540 break;
3541 }
3542 }
3543 if ( n0 == -1 ) { return 0.0; }
3544 else
3545 {
3546 double k = ( y[n0] - y[n0 + 1] ) / ( x[n0] - x[n0 + 1] );
3547 z = y[n0 + 1] + k * ( t - x[n0 + 1] );
3548 }
3549 // std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
3550 // std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
3551 return z;
3552}
3553
3554double EvtConExc::Mhad_sampling( double* x, double* y ) { // sample ISR hadrons, return Mhadron
3555 // x=MH: array contains the Mhad
3556 // y=AF: array contains the Mhad*Xsection
3557 // n: specify how many bins in the hadron sampling
3558 int n = 598; // AF[599] is the total xs including Ecms point
3559 double pm = EvtRandom::Flat( 0., 1 );
3560 double xrt = 1;
3561 int mybin = 1;
3562 double xst = y[n];
3563 for ( int i = 0; i < n; i++ )
3564 {
3565 if ( ( y[i] / xst ) < pm && pm <= ( y[i + 1] / xst ) )
3566 {
3567 mybin = i + 1;
3568 break;
3569 }
3570 }
3571 pm = EvtRandom::Flat( 0., 1 );
3572 if ( pm > 1 ) { std::cout << "random larger than 1: " << pm << std::endl; }
3573 double mhd = x[mybin - 1] + ( x[mybin] - x[mybin - 1] ) * pm;
3574
3575 if ( ( mhd - _cms ) > 0 )
3576 {
3577 std::cout << "selected mhd larger than Ecms: " << mhd << " > " << x[mybin] << std::endl;
3578 abort();
3579 }
3580 if ( mhd < _mhdL )
3581 {
3582 std::cout << "the sample mhassrons is less than the low end of XS" << mhd << " <" << _mhdL
3583 << std::endl;
3584 for ( int i = 0; i < 598; i++ )
3585 { std::cout << i << " " << x[i] << " " << y[i] << std::endl; }
3586 abort();
3587 }
3588 return mhd;
3589}
3590
3591double EvtConExc::ISR_ang_integrate( double x, double theta ) {
3592 // see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
3593 double costheta = cos( theta );
3594 double eb = _cms / 2;
3595 double cos2 = costheta * costheta;
3596 double sin2 = 1 - cos2;
3597 double me = 0.51 * 0.001;
3598 double L = 2 * log( _cms / me );
3599 double meE2 = me * me / eb / eb;
3600 double hpi = L - 1;
3601 double hq1 = meE2 / 2 * costheta / ( sin2 + meE2 * cos2 );
3602 double hq2 =
3603 -0.5 * log( ( 1 + sqrt( 1 - meE2 ) * costheta ) / ( 1 - sqrt( 1 - meE2 ) * costheta ) );
3604 double hq3 =
3605 x * x * costheta / 2 / ( x * x - 2 * x + 2 ) * ( 1 - meE2 / ( sin2 + meE2 * cos2 ) );
3606 double hq = ( L - 1 ) / 2. + hq1 + hq2 + hq3;
3607 hq /= hpi; // Eq (11) in arXif:hep-ph/9910523
3608 return hq;
3609}
3610
3611double EvtConExc::ISR_ang_sampling( double x ) {
3612 double xx[180], yy[180];
3613 double dgr2rad = 1. / 180. * 3.1415926;
3614 xx[0] = 0;
3615 xx[1] = 5 * dgr2rad; // first bin is 5 degree
3616 int nc = 2;
3617 for ( int i = 6; i <= 175; i++ )
3618 { // last bin is 5 degree
3619 xx[nc] = i * dgr2rad;
3620 nc++;
3621 }
3622 xx[nc] = 180 * dgr2rad;
3623 for ( int j = 0; j <= nc; j++ ) { yy[j] = ISR_ang_integrate( x, xx[j] ); }
3624 double pm = EvtRandom::Flat( 0., 1 );
3625 int mypt = 1;
3626 for ( int j = 1; j <= nc; j++ )
3627 {
3628 if ( yy[j - 1] < pm && pm <= yy[j] )
3629 {
3630 mypt = j;
3631 break;
3632 }
3633 }
3634 pm = EvtRandom::Flat( 0., 1 );
3635 double mytheta = xx[mypt - 1] + ( xx[mypt] - xx[mypt - 1] ) * pm;
3636 return mytheta; // theta in rad unit
3637}
3638
3639void EvtConExc::SetP4( EvtParticle* p, double mhdr, double xeng,
3640 double theta ) { // set the gamma and gamma* momentum according sampled
3641 // results
3642 double pm = EvtRandom::Flat( 0., 1 );
3643 double phi = 2 * 3.1415926 * pm;
3644 double gamE = _cms / 2 * xeng;
3645 double hdrE = sqrt( mhdr * mhdr + gamE * gamE );
3646 double px = gamE * sin( theta ) * cos( phi );
3647 double py = gamE * sin( theta ) * sin( phi );
3648 double pz = gamE * cos( theta );
3649 EvtVector4R p4[2];
3650 p4[0].set( gamE, px, py, pz ); // ISR photon
3651 p4[1].set( hdrE, -px, -py, -pz ); // mhdraons
3652
3653 for ( int i = 0; i < 2; i++ )
3654 {
3655 EvtId myid = p->getDaug( i )->getId();
3656 p->getDaug( i )->init( myid, p4[i] );
3657 if ( EvtPDL::name( myid ) == "gamma*" )
3658 {
3659 if ( ( _cms - mhdr ) < 0.002 )
3660 {
3661 EvtVector4R PG( _cms, 0, 0, 0 );
3662 p->getDaug( i )->init( myid, PG );
3663 }
3664 }
3665 }
3666}
3667
3668void EvtConExc::SetP4Rvalue( EvtParticle* p, double mhdr, double xeng,
3669 double theta ) { // set the gamma and gamma* momentum according
3670 // sampled results
3671 double pm = EvtRandom::Flat( 0., 1 );
3672 double phi = 2 * 3.1415926 * pm;
3673 double gamE = _cms / 2 * xeng;
3674 double hdrE = sqrt( mhdr * mhdr + gamE * gamE );
3675 double px = gamE * sin( theta ) * cos( phi );
3676 double py = gamE * sin( theta ) * sin( phi );
3677 double pz = gamE * cos( theta );
3678 EvtVector4R p4[2];
3679 p4[0].set( hdrE, px, py, pz ); // mhdraons
3680 p4[1].set( gamE, -px, -py, -pz ); // ISR photon
3681 for ( int i = 0; i < 2; i++ )
3682 {
3683 EvtId myid = p->getDaug( i )->getId();
3684 p->getDaug( i )->init( myid, p4[i] );
3685 }
3686}
3687
3688void EvtConExc::findMaxXS( double mhds ) { // dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3689 maxXS = -1;
3690 for ( int i = 0; i < 90000; i++ )
3691 {
3692 double x = 1 - ( mhds / _cms ) * ( mhds / _cms );
3693 double cos = EvtRandom::Flat( 0.0006, 0.9994 ); // set cut on ISR gamma 2^o
3694 double sin2theta = sqrt( 1 - cos * cos );
3695 double alpha = 1. / 137.;
3696 double pi = 3.1415926;
3697 double wsx = alpha / pi / x * ( ( 2 - 2 * x + x * x ) / sin2theta - x * x / 2. );
3698 double xs = myxsection->getXsection( mhds );
3699 double difxs = 2 * mhds / _cms / _cms * wsx * xs;
3700 if ( difxs > maxXS ) maxXS = difxs;
3701 }
3702}
3703
3704double EvtConExc::difgamXs( double mhds,
3705 double sintheta ) { // dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3706 double x = 1 - ( mhds / _cms ) * ( mhds / _cms );
3707 double sin2theta = sintheta * sintheta;
3708 double alpha = 1. / 137.;
3709 double pi = 3.1415926;
3710 double wsx = alpha / pi / x * ( ( 2 - 2 * x + x * x ) / sin2theta - x * x / 2. );
3711 double xs = myxsection->getXsection( mhds );
3712 double difxs = 2 * mhds / _cms / _cms * wsx * xs;
3713 return difxs;
3714}
3715
3716int EvtConExc::selectMode( std::vector<int> vmod, double mhds ) {
3717 // first get xsection for each mode in vmod array
3718 // std::cout<<"vmod.size, mhds "<<vmod.size()<<" "<<mhds<<std::endl;
3719 std::vector<int> excmod;
3720 excmod.push_back( 0 );
3721 excmod.push_back( 1 );
3722 excmod.push_back( 2 );
3723 excmod.push_back( 6 );
3724 excmod.push_back( 7 );
3725 excmod.push_back( 12 );
3726 excmod.push_back( 13 );
3727 excmod.push_back( 45 );
3728 excmod.push_back( 46 );
3729 double uscale = 1;
3730 std::vector<double> vxs;
3731 vxs.clear();
3732 for ( int i = 0; i < vmod.size(); i++ )
3733 {
3734 int imod = vmod[i];
3735 delete myxsection;
3736 myxsection = new EvtXsection( imod );
3737 if ( myxsection->getUnit() == "pb" ) { uscale = 0.001; }
3738 else { uscale = 1; }
3739 // if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
3740
3741 double myxs = uscale * myxsection->getXsection( mhds );
3742
3743 if ( i == 0 ) { vxs.push_back( myxs ); }
3744 else if ( imod == 74110 )
3745 { // for continuum process
3746 double xcon = myxs - vxs[i - 1];
3747 if ( xcon < 0 ) { xcon = vxs[i - 1]; }
3748 else { xcon = myxs; }
3749 if ( mhds < 2.0 )
3750 xcon = vxs[i - 1]; // continuum cut: _cms energy less than 1.1, sampling phi,rho0 and
3751 // omega resonance, see parp(2)=1.08 at Pythia.F
3752 vxs.push_back( xcon );
3753 }
3754 else { vxs.push_back( myxs + vxs[i - 1] ); }
3755 } // for_i_block
3756 // degugging
3757 // for(int i=0;i<vxs.size();i++){std::cout<<"Mhds="<<mhds<<" Mode="<<vmod[i]<<" xs_i=
3758 // "<<vxs[i]<<" totalxs "<< vxs[vxs.size()-1]<<std::endl;}
3759
3760 double totxs = vxs[vxs.size() - 1];
3761 // if(totxs==0){std::cout<<"totalXs=0, vxs.size()= "<<vxs.size()<<std::endl;}
3762 int icount = 0;
3763 int themode;
3764mode_iter:
3765 double pm = EvtRandom::Flat( 0., 1 ); // std::cout<<"pm= "<<pm<<" mhds= "<<mhds<<std::endl;
3766 if ( pm <= vxs[0] / totxs )
3767 {
3768 themode = vmod[0];
3769 std::vector<EvtId> theDaug = get_mode( themode );
3770 EvtId daugs[100];
3771 for ( int i = 0; i < theDaug.size(); i++ ) { daugs[i] = theDaug[i]; }
3772 int exmode = EvtDecayTable::inChannelList( theparent->getId(), theDaug.size(), daugs );
3773 if ( _mode != 74110 ) { return themode; }
3774 else if ( exmode == -1 ) { return themode; }
3775 else { goto mycount; }
3776 }
3777
3778 for ( int j = 1; j < vxs.size(); j++ )
3779 {
3780 double x0 = vxs[j - 1] / totxs;
3781 double x1 = vxs[j] / totxs; // std::cout<<"j,x0,x1 "<<j<<" "<<x0<<" "<<x1<<std::endl;
3782 if ( x0 < pm && pm <= x1 )
3783 {
3784 themode = vmod[j];
3785 std::vector<EvtId> theDaug = get_mode( themode );
3786 EvtId daugs[100];
3787 for ( int i = 0; i < theDaug.size(); i++ ) { daugs[i] = theDaug[i]; }
3788 int exmode = EvtDecayTable::inChannelList( theparent->getId(), theDaug.size(), daugs );
3789 if ( _mode != 74110 ) { return themode; }
3790 else if ( exmode == -1 ) { return themode; }
3791 else { break; }
3792 }
3793 }
3794
3795mycount:
3796 icount++;
3797 if ( icount < 20000 ) goto mode_iter;
3798 // debugging
3799 // for(int i=0;i<vxs.size();i++){std::cout<<"Random="<<pm<<" Mode,Mhad "<<vmod[i]<<",
3800 // "<<mhds<<" xs_i "<<vxs[i]<<" xsi/totalxs "<<vxs[i]/totxs<<std::endl;}
3801 std::cout << "EvtConExc::selectMode can not find a mode with 20000 iteration for Mhadrons="
3802 << mhds << std::endl;
3803 return -10;
3804 // abort();
3805}
3806
3807void EvtConExc::resetResMass() {
3808 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
3809 EvtId myres = EvtPDL::getId( std::string( "J/psi" ) );
3810 EvtPDL::reSetMass( myres, mjsi );
3811 // EvtPDL::reSetWidth(myres,wjsi);
3812
3813 myres = EvtPDL::getId( std::string( "psi(2S)" ) );
3814 EvtPDL::reSetMass( myres, mpsip );
3815 // EvtPDL::reSetWidth(myres,wpsip);
3816
3817 myres = EvtPDL::getId( std::string( "psi(3770)" ) );
3818 EvtPDL::reSetMass( myres, mpsipp );
3819 // EvtPDL::reSetWidth(myres,wpsipp);
3820
3821 myres = EvtPDL::getId( std::string( "phi" ) );
3822 EvtPDL::reSetMass( myres, mphi );
3823 // EvtPDL::reSetWidth(myres,wphi);
3824
3825 myres = EvtPDL::getId( std::string( "omega" ) );
3826 EvtPDL::reSetMass( myres, momega );
3827 // EvtPDL::reSetWidth(myres,womega);
3828
3829 myres = EvtPDL::getId( std::string( "rho0" ) );
3830 EvtPDL::reSetMass( myres, mrho0 );
3831 // EvtPDL::reSetWidth(myres,wrho0);
3832
3833 myres = EvtPDL::getId( std::string( "rho(3S)0" ) );
3834 EvtPDL::reSetMass( myres, mrho3s );
3835 // EvtPDL::reSetWidth(myres,wrho3s);
3836
3837 myres = EvtPDL::getId( std::string( "omega(2S)" ) );
3838 EvtPDL::reSetMass( myres, momega2s );
3839 // EvtPDL::reSetWidth(myres,womega2s);
3840}
3841
3842void EvtConExc::getResMass() {
3843 EvtId myres = EvtPDL::getId( std::string( "J/psi" ) );
3844 mjsi = EvtPDL::getMeanMass( myres );
3845 wjsi = EvtPDL::getWidth( myres );
3846
3847 myres = EvtPDL::getId( std::string( "psi(2S)" ) );
3848 mpsip = EvtPDL::getMeanMass( myres );
3849 wpsip = EvtPDL::getWidth( myres );
3850
3851 myres = EvtPDL::getId( std::string( "psi(3770)" ) );
3852 mpsipp = EvtPDL::getMeanMass( myres );
3853 wpsipp = EvtPDL::getWidth( myres );
3854
3855 myres = EvtPDL::getId( std::string( "phi" ) );
3856 mphi = EvtPDL::getMeanMass( myres );
3857 wphi = EvtPDL::getWidth( myres );
3858
3859 myres = EvtPDL::getId( std::string( "omega" ) );
3860 momega = EvtPDL::getMeanMass( myres );
3861 womega = EvtPDL::getWidth( myres );
3862
3863 myres = EvtPDL::getId( std::string( "rho0" ) );
3864 mrho0 = EvtPDL::getMeanMass( myres );
3865 wrho0 = EvtPDL::getWidth( myres );
3866
3867 myres = EvtPDL::getId( std::string( "rho(3S)0" ) );
3868 mrho3s = EvtPDL::getMeanMass( myres );
3869 wrho3s = EvtPDL::getWidth( myres );
3870
3871 myres = EvtPDL::getId( std::string( "omega(2S)" ) );
3872 momega2s = EvtPDL::getMeanMass( myres );
3873 womega2s = EvtPDL::getWidth( myres );
3874}
3875
3876void EvtConExc::showResMass() {
3877 std::cout << "J/psi: " << mjsi << " " << wjsi << std::endl;
3878 std::cout << "psipr: " << mpsip << " " << wpsip << std::endl;
3879 std::cout << "psipp: " << mpsipp << " " << wpsipp << std::endl;
3880 std::cout << "phi : " << mphi << " " << wphi << std::endl;
3881 std::cout << "omega: " << momega << " " << womega << std::endl;
3882 std::cout << "rho0 : " << mrho0 << " " << wrho0 << std::endl;
3883 std::cout << "rho(3S)0: " << mrho3s << " " << wrho3s << std::endl;
3884 std::cout << "omega(2S): " << momega2s << " " << womega2s << std::endl;
3885}
3886
3887bool EvtConExc::checkdecay( EvtParticle* p ) {
3888 int nson = p->getNDaug();
3889 double massflag = 1;
3890 for ( int i = 0; i < nson; i++ )
3891 {
3892 if ( EvtPDL::getStdHep( p->getDaug( i )->getId() ) == 22 ) continue;
3893 massflag *= p->getDaug( i )->mass();
3894 }
3895 if ( massflag == 0 )
3896 {
3897 std::cout << "Zero mass!" << std::endl;
3898 return 0;
3899 }
3900 else { return 1; }
3901}
3902
3903double EvtConExc::sumExc( double mx ) { // summation all cross section of exclusive decays
3904 std::vector<int> vmod;
3905 vmod.clear();
3906 int mn[83] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 16, 17, 18, 19, 20, 21, 22, 23,
3907 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 40, 41, 44, 45,
3908 46, // 12,14, 30,31,39,42 and 43 is
3909 // removed
3910 50, 51, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
3911 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100,
3912 74101, 74102, 74103, 74104, 74105, 74110 };
3913 int mn2[84] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 20, 21, 22,
3914 23, 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 40, 41, 44, 45,
3915 46, // mode 14,30,31,42, 43 are
3916 // removed
3917 50, 51, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
3918 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100,
3919 74101, 74102, 74103, 74104, 74105, 74110 };
3920
3921 if ( _cms > 2.5 )
3922 {
3923 for ( int i = 0; i < 83; i++ ) { vmod.push_back( mn[i] ); }
3924 }
3925 else
3926 {
3927 for ( int i = 0; i < 84; i++ ) { vmod.push_back( mn2[i] ); }
3928 }
3929
3930 // for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
3931
3932 double xsum = 0;
3933 double uscale = 1;
3934 for ( int i = 0; i < vmod.size(); i++ )
3935 {
3936 int mymode = vmod[i];
3937 if ( mymode == 74110 ) continue;
3938 delete myxsection;
3939 myxsection = new EvtXsection( mymode );
3940 init_mode( mymode );
3941 if ( myxsection->getUnit() == "pb" ) { uscale = 0.001; }
3942 else { uscale = 1; }
3943 xsum += uscale * myxsection->getXsection( mx );
3944 // if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
3945 }
3946 return xsum;
3947}
3948
3949void EvtConExc::PrintXS( double mx ) { // print xsection mode by mode at mx energy point
3950 std::vector<int> vmod;
3951 vmod.clear();
3952 int mn[83] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 15, 16, 17, 18, 19, 20, 21, 22, 23,
3953 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 40, 41, 44, 45,
3954 46, // 12,14, 30,31,39,42 and 43 is
3955 // removed
3956 50, 51, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
3957 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100,
3958 74101, 74102, 74103, 74104, 74105, 74110 };
3959 int mn2[84] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18, 19, 20, 21, 22,
3960 23, 24, 25, 26, 27, 28, 29, 32, 33, 34, 35, 36, 37, 38, 40, 41, 44, 45,
3961 46, // mode 14,30,31,42, 43 are
3962 // removed
3963 50, 51, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75,
3964 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 90, 91, 92, 93, 94, 95, 96, 74100,
3965 74101, 74102, 74103, 74104, 74105, 74110 };
3966
3967 if ( _cms > 2.5 )
3968 {
3969 for ( int i = 0; i < 83; i++ ) { vmod.push_back( mn[i] ); }
3970 }
3971 else
3972 {
3973 for ( int i = 0; i < 84; i++ ) { vmod.push_back( mn2[i] ); }
3974 }
3975
3976 // for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
3977
3978 double xsum = 0;
3979 double uscale = 1;
3980 double hxs = 0;
3981 for ( int i = 0; i < vmod.size(); i++ )
3982 {
3983 int mymode = vmod[i];
3984 delete myxsection;
3985 myxsection = new EvtXsection( mymode );
3986 init_mode( mymode );
3987 if ( myxsection->getUnit() == "pb" ) { uscale = 0.001; }
3988 else { uscale = 1; }
3989 double mixs = uscale * myxsection->getXsection( mx );
3990 if ( mymode == 74110 )
3991 {
3992 hxs = mixs;
3993 continue;
3994 }
3995 xsum += mixs;
3996 // if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
3997 std::cout << setw( 5 ) << mymode << " : ";
3998 std::cout << setw( 10 ) << mixs << " nb ";
3999 for ( int ii = 0; ii < _ndaugs; ii++ )
4000 { std::cout << setw( 10 ) << EvtPDL::name( daugs[ii] ) << " "; }
4001 std::cout << std::endl;
4002 }
4003 std::cout << setw( 6 ) << "LUARLW " << setw( 10 ) << hxs - xsum << " nb " << setw( 10 )
4004 << " anything " << std::endl;
4005 std::cout << "Total hadron cross section here is " << hxs << " nb" << std::endl;
4006}
4007
4008bool EvtConExc::angularSampling( EvtParticle* par ) {
4009 bool tagp, tagk;
4010 tagk = 0;
4011 tagp = 0;
4012 int nds = par->getNDaug();
4013 for ( int i = 0; i < par->getNDaug(); i++ )
4014 {
4015 EvtId di = par->getDaug( i )->getId();
4016 EvtVector4R p4i = par->getDaug( i )->getP4Lab();
4017 int pdgcode = EvtPDL::getStdHep( di );
4018 double alpha = 1;
4019 /*
4020 if(fabs(pdgcode)==2212){//proton and anti-proton
4021 alpha = baryonAng(p4i.mass());
4022 cosp = cos(p4i.theta());
4023 tagp=1;
4024 }else if(fabs(pdgcode)==321||fabs(pdgcode)==211){
4025 alpha=1;
4026 cosk = cos(p4i.theta());
4027 tagk=1;
4028 }else continue;
4029 */
4030 double angmax = 1 + alpha;
4031 double costheta = cos( p4i.theta() );
4032 double ang = 1 + alpha * costheta * costheta;
4033 double xratio = ang / angmax;
4034 double xi = EvtRandom::Flat( 0., 1 );
4035 // std::cout<<"pdgcode "<<pdgcode<<std::endl;
4036 // std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
4037 if ( xi > xratio ) return false;
4038 } // loop over duaghters
4039 // if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
4040 // if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
4041
4042 return true;
4043}
4044
4045double EvtConExc::baryonAng( double mx ) {
4046 double mp = 0.938;
4047 double u = 0.938 / mx;
4048 u = u * u;
4049 double u2 = ( 1 + u ) * ( 1 + u );
4050 double uu = u * ( 1 + 6 * u );
4051 double alpha = ( u2 - uu ) / ( u2 + uu );
4052 return alpha;
4053}
4054
4055bool EvtConExc::photonSampling( EvtParticle* par ) {
4056 bool tagp, tagk;
4057 tagk = 0;
4058 tagp = 0;
4059 int nds = par->getNDaug();
4060 for ( int i = 0; i < par->getNDaug(); i++ )
4061 {
4062 EvtId di = par->getDaug( i )->getId();
4063 EvtVector4R p4i = par->getDaug( i )->getP4Lab();
4064 int pdgcode = EvtPDL::getStdHep( di );
4065 double alpha = 0;
4066
4067 if ( pdgcode == 111 || pdgcode == 221 || pdgcode == 331 )
4068 { // for photon
4069 alpha = 1;
4070 }
4071 else continue;
4072
4073 double angmax = 1 + alpha;
4074 double costheta = cos( p4i.theta() );
4075 double ang = 1 + alpha * costheta * costheta;
4076 double xratio = ang / angmax;
4077 double xi = EvtRandom::Flat( 0., 1 );
4078 // std::cout<<"pdgcode "<<pdgcode<<std::endl;
4079 // std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
4080 if ( xi > xratio ) return false;
4081 } // loop over duaghters
4082 // if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
4083 // if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
4084
4085 return true;
4086}
4087
4088double EvtConExc::narrowRXS( double mxL, double mxH ) {
4089 // br for ee
4090 double psippee, psipee, jsiee, phiee, omegaee, rhoee;
4091 double kev2Gev = 0.000001;
4092 psippee = 0.262 * kev2Gev;
4093 psipee = 2.36 * kev2Gev;
4094 jsiee = 5.55 * kev2Gev;
4095 phiee = 4.266 * 0.001 * 2.954 * 0.0001;
4096 omegaee = 0.6 * kev2Gev;
4097 rhoee = 7.04 * kev2Gev;
4098 double s = _cms * _cms;
4099 double sigv = 0;
4100 double mv = 0;
4101 double widee = 0;
4102 double xpi = 12 * 3.1415926 * 3.1415926;
4103 if ( mxL <= 3.0957 && 3.0957 <= mxH || mxL <= 3.0983 && 3.0983 <= mxH )
4104 {
4105 widee = jsiee;
4106 mv = 3.096916;
4107 }
4108 else if ( mxL <= 3.6847 && 3.6847 <= mxH || mxL <= 3.6873 && 3.6873 <= mxH )
4109 {
4110 widee = psipee;
4111 mv = 3.686109;
4112 }
4113 else if ( mxL <= 1.01906 && 1.01906 <= mxH || mxL <= 1.01986 && 1.01986 <= mxH )
4114 {
4115 widee = phiee;
4116 mv = 1.01946;
4117 }
4118 else { return -1.0; }
4119
4120 if ( _cms <= mv ) return -1.;
4121 double xv = 1 - mv * mv / s;
4122 sigv += xpi * widee / mv / s * Rad2( s, xv );
4123 double unic = 3.89379304 * 100000; // in unit nb
4124 return sigv * unic; // vaccum factor has included in the Rad2
4125}
4126
4127double EvtConExc::addNarrowRXS( double mhi, double binwidth ) {
4128 double myxs = 0;
4129 for ( int i = 0; i < ISRXS.size(); i++ )
4130 { // std::cout<<"ISRM "<<ISRM[i]<<std::endl;
4131 if ( ( ISRM[i] - mhi ) < binwidth && ISRFLAG[i] )
4132 {
4133 ISRFLAG[i] = 0;
4134 myxs = ISRXS[i];
4135 }
4136 }
4137 // std::cout<<mhi<<" "<<binwidth<<" "<<myxs<<std::endl;
4138 return myxs;
4139}
4140
4141double EvtConExc::selectMass() {
4142 double pm, mhdL, mhds;
4143 pm = EvtRandom::Flat( 0., 1 );
4144 mhdL = _mhdL;
4145 mhds = pm * ( _cms - mhdL ) + mhdL; // non narrow resonance
4146 std::vector<double> sxs;
4147 for ( int i = 0; i < ISRID.size(); i++ )
4148 {
4149 double ri = ISRXS[i] / AF[598];
4150 sxs.push_back( ri );
4151 }
4152 for ( int j = 0; j < sxs.size(); j++ )
4153 {
4154 if ( j > 0 ) sxs[j] += sxs[j - 1];
4155 }
4156 pm = EvtRandom::Flat( 0., 1 );
4157 if ( pm < sxs[0] )
4158 {
4159 mhds = EvtPDL::getMass( ISRID[0] ); // narrow resonance
4160 }
4161 for ( int j = 1; j < sxs.size(); j++ )
4162 {
4163 double x0 = sxs[j - 1];
4164 double x1 = sxs[j];
4165 if ( x0 < pm && pm <= x1 ) mhds = EvtPDL::getMass( ISRID[j] ); // narrow resonance
4166 }
4167 return mhds;
4168}
4169
4170double EvtConExc::getVP( double mx ) {
4171 double xx, r1, i1;
4172 double x1, y1, y2;
4173 xx = 0;
4174 int mytag = 1;
4175 for ( int i = 0; i < 4001; i++ )
4176 {
4177 x1 = vpx[i];
4178 y1 = vpr[i];
4179 y2 = vpi[i];
4180 if ( xx <= mx && mx < x1 )
4181 {
4182 xx = x1;
4183 r1 = y1;
4184 i1 = y2;
4185 mytag = 0;
4186 break;
4187 }
4188 xx = x1;
4189 r1 = y1;
4190 i1 = y2;
4191 }
4192 if ( mytag == 1 )
4193 {
4194 std::cout << "No vacuum polarization value found" << std::endl;
4195 abort();
4196 }
4197 EvtComplex cvp( r1, i1 );
4198 double thevp = abs( 1. / ( 1 - cvp ) );
4199 if ( 3.0933 < mx && mx < 3.1035 ) return 1.; // J/psi peak excluded
4200 if ( 3.6810 < mx && mx < 3.6913 ) return 1.; // psi(2S) peak excluded
4201 return thevp * thevp;
4202}
4203
4204void EvtConExc::ReadVP() {
4205 vpx.clear();
4206 vpr.clear();
4207 vpi.clear();
4208 int vpsize = 4001;
4209 vpx.resize( vpsize );
4210 vpr.resize( vpsize );
4211 vpi.resize( vpsize );
4212 std::string locvp = getenv( "BESEVTGENROOT" );
4213 locvp += "/share/vp.dat";
4214 ifstream m_inputFile;
4215 m_inputFile.open( locvp.c_str() );
4216 double xx, r1, i1;
4217 double x1, y1, y2;
4218 xx = 0;
4219 int mytag = 1;
4220 for ( int i = 0; i < 4001; i++ ) { m_inputFile >> vpx[i] >> vpr[i] >> vpi[i]; }
4221}
4222
4223void EvtConExc::mk_VXS( double Esig, double Egamcut, double EgamH, int midx ) {
4224 // the mode index is put into vmode
4225 // initial
4226 // midx: mode index in vmode[midx] and VXS[midx][bin]
4227 int mode = vmode[midx];
4228 double uscale;
4229
4230 EvtXsection* myxsection = new EvtXsection( mode );
4231 if ( myxsection->getUnit() == "pb" ) { uscale = 0.001; }
4232 else { uscale = 1; }
4233 double s = _cms * _cms;
4234 double mlow = myxsection->getXlw();
4235 if ( _cms <= mlow )
4236 {
4237 for ( int i = 0; i < 600; i++ ) { VXS[midx][i] = 0; }
4238 return;
4239 }
4240
4241 int nint = 50;
4242 int nmc = 800;
4243 // double myxs0 =uscale*gamHXSection(s,Esig,Egamcut,nmc);
4244 double myxs0 = uscale * trapezoid( s, Esig, Egamcut, nint, myxsection );
4245 double xb = 2 * Esig / _cms;
4246 double sig_xs = uscale * SoftPhoton_xs( s, xb ) * ( myxsection->getXsection( _cms ) );
4247 myxs0 += sig_xs;
4248
4249 int Nstp = 598;
4250 double stp = ( EgamH - Egamcut ) / Nstp;
4251 for ( int i = 0; i < Nstp; i++ )
4252 { // points within Egam(max) ~ Egamcut=25MeV
4253 double Eg0 = EgamH - i * stp;
4254 double Eg1 = EgamH - ( i + 1 ) * stp;
4255 // double xsi= uscale*gamHXSection(s,Eg1,Eg0,nmc);
4256 double xsi = uscale * trapezoid( s, Eg1, Eg0, nint, myxsection );
4257 if ( i == 0 ) VXS[midx][0] = xsi;
4258 if ( i >= 1 ) VXS[midx][i] = VXS[midx][i - 1] + xsi;
4259 }
4260 VXS[midx][598] = VXS[midx][597];
4261 VXS[midx][599] = VXS[midx][598] + myxs0;
4262 // std::cout<<"mode "<<mode<<" "<<uscale<<" "<<VXS[midx][599]<<std::endl;
4263}
4264
4265int EvtConExc::get_mode_index( int mode ) {
4266 int i = 0;
4267 for ( int j = 0; i < vmode.size(); j++ )
4268 {
4269 if ( mode == vmode[j] ) return j;
4270 }
4271 std::cout << " EvtConExc::get_mode_index: no index is found in vmode" << std::endl;
4272 abort();
4273}
4274
4275double EvtConExc::getObsXsection( double mhds, int mode ) {
4276 double hwid = ( AA[0] - AA[1] ) / 2.;
4277 double s = _cms * _cms;
4278 int idx = get_mode_index( mode );
4279 double Egam = ( s - mhds * mhds ) / 2 / _cms;
4280 for ( int i = 0; i < 600; i++ )
4281 {
4282 if ( ( AA[i] - hwid ) <= Egam && Egam < ( AA[i] + hwid ) ) return VXS[idx][i];
4283 }
4284
4285std:
4286 cout << "EvtConExc::getObsXsection : no observed xs is found " << endl;
4287 abort();
4288}
4289
4290double EvtConExc::Egam2Mhds( double Egam ) {
4291 double s = _cms * _cms;
4292 double mhds = sqrt( s - 2 * Egam * _cms );
4293 return mhds;
4294}
4295
4296void EvtConExc::writeDecayTabel() {
4297 ofstream oa, ob;
4298 oa.open( "_pkhr.dec" );
4299 ob.open( "obsxs.dat" );
4300 //
4301 int im = getModeIndex( 74110 );
4302
4303 double xscon = VXS[im][599];
4304 // std::cout<<"tsize "<<tsize<<" "<<xscon<<" "<<VXS[70][599]<<std::endl;
4305 std::vector<int> Vmode;
4306 Vmode.push_back( 6 ); // 1:pi+pi-
4307 Vmode.push_back( 13 ); // 2:pi+pi-2pi0
4308 Vmode.push_back( 12 ); // 3:2(pi+pi-)
4309 Vmode.push_back( 0 ); // 4:ppbar
4310 Vmode.push_back( 1 ); // 5:nnbar
4311 Vmode.push_back( 45 ); // 6:K+K-
4312 Vmode.push_back( 46 ); // 7:K0K0bar
4313 Vmode.push_back( 7 ); // 8:pi+pi-pi0
4314 Vmode.push_back( 2 ); // 9:Lambda anti-Lambda
4315 Vmode.push_back( 37 ); // 10: pi+pi- eta
4316 std::vector<std::string> vmdl;
4317 vmdl.push_back( "PHOKHARA_pipi" );
4318 vmdl.push_back( "PHOKHARA_pi0pi0pipi" );
4319 vmdl.push_back( "PHOKHARA_4pi" );
4320 vmdl.push_back( "PHOKHARA_ppbar" );
4321 vmdl.push_back( "PHOKHARA_nnbar" );
4322 vmdl.push_back( "PHOKHARA_KK" );
4323 vmdl.push_back( "PHOKHARA_K0K0" );
4324 vmdl.push_back( "PHOKHARA_pipipi0" );
4325 vmdl.push_back( "PHOKHARA_LLB" );
4326 vmdl.push_back( "PHOKHARA_pipieta" );
4327
4328 ob << "mode_index observed cross /nb" << std::endl;
4329 for ( int i = 0; i < vmode.size(); i++ )
4330 {
4331 ob << vmode[i] << " " << VXS[getModeIndex( vmode[i] )][599] << std::endl;
4332 // cout<<"vmode["<<i<<"] = "<<vmode[i]<<", VXS["<<getModeIndex(vmode[i])<<"][599] =
4333 // "<<XS[getModeIndex(vmode[i])][599]<<std::endl; // wangwp
4334 }
4335
4336 for ( int i = 0; i < Vmode.size(); i++ )
4337 {
4338 // std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
4339 xscon -= VXS[getModeIndex( Vmode[i] )][599];
4340 }
4341 // debugging
4342 for ( int i = 0; i < Vmode.size(); i++ )
4343 {
4344 // std::cout<<Vmode[i]<<"-th mode: "<<VXS[getModeIndex(Vmode[i])][599]<<" "<<std::endl;
4345 }
4346
4347 oa << "LundaPar PARJ(11)=0.611798" << std::endl;
4348 oa << "LundaPar PARJ(12)=7.92527E-12" << std::endl;
4349 oa << "LundaPar PARJ(14)=0.244495" << std::endl;
4350 oa << "LundaPar PARJ(15)=1.16573E-15" << std::endl;
4351 oa << "LundaPar PARJ(16)=0.436516" << std::endl;
4352 oa << "LundaPar PARJ(17)=0.530517" << std::endl;
4353 oa << "LundaPar PARJ(1)=0.0651577" << std::endl;
4354 oa << "LundaPar PARJ(2)=0.260378" << std::endl;
4355 oa << "LundaPar PARJ(21)=0.0664835" << std::endl;
4356 oa << "LundaPar RALPA(15)=0.576687" << std::endl;
4357 oa << "LundaPar RALPA(16)=0.364796" << std::endl;
4358 oa << "LundaPar RALPA(17)=3.19486E-06" << std::endl;
4359 oa << "noPhotos" << std::endl;
4360 oa << "Particle vpho " << _cms << " 0" << std::endl;
4361 oa << "Decay vpho" << std::endl;
4362 oa << "0 pi+ pi- PHSP ;" << std::endl;
4363 oa << "0 pi0 pi0 pi- pi+ PHSP ;" << std::endl;
4364 oa << "0 pi+ pi- pi- pi+ PHSP ;" << std::endl;
4365 oa << "0 anti-p- p+ PHSP ;" << std::endl;
4366 oa << "0 anti-n0 n0 PHSP ;" << std::endl;
4367 oa << "0 K+ K- PHSP ;" << std::endl;
4368 oa << "0 K_S0 K_L0 PHSP ;" << std::endl;
4369 oa << "0 pi+ pi- pi0 PHSP ;" << std::endl;
4370 oa << "0 Lambda0 anti-Lambda0 PHSP ;" << std::endl;
4371 oa << "0 pi+ pi- eta PHSP ;" << std::endl;
4372 oa << "0 gamma phi PHSP;" << std::endl;
4373 oa << "0 gamma rho0 PHSP;" << std::endl;
4374 oa << "0 gamma omega PHSP;" << std::endl;
4375 oa << xscon << " ConExc 74110;" << std::endl;
4376 for ( int i = 0; i < Vmode.size(); i++ )
4377 { oa << VXS[getModeIndex( Vmode[i] )][599] << " " << vmdl[i] << " ;" << std::endl; }
4378 oa << "Enddecay" << std::endl;
4379 oa << "End" << std::endl;
4380}
4381
4382void EvtConExc::checkEvtRatio() {
4383 ofstream oa;
4384 oa.open( "_evtR.dat" );
4385 //
4386 int im = getModeIndex( 74110 );
4387 double xscon = VXS[im][599];
4388 double xscon0 = xscon;
4389 oa << "Ecms= " << _cms << " GeV" << std::endl;
4390 oa << "The total observed xs: " << xscon << " nb" << std::endl;
4391 oa << "=== mode =========== ratio ===========" << std::endl;
4392 for ( int i = 0; i < vmode.size(); i++ )
4393 {
4394 // std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
4395 int themode = getModeIndex( vmode[i] );
4396 if ( vmode[i] == 74110 ) continue;
4397 xscon -= VXS[themode][599];
4398 }
4399 if ( xscon < 0 ) xscon = 0;
4400 // debugging
4401 for ( int i = 0; i < vmode.size(); i++ )
4402 {
4403 int themode = getModeIndex( vmode[i] );
4404 if ( vmode[i] == 74110 ) continue;
4405 oa << vmode[i] << "-th mode: " << 100 * VXS[themode][599] / xscon0 << " % " << std::endl;
4406 }
4407 oa << "74110-th mode: " << 100 * xscon / xscon0 << " % " << std::endl;
4408}
4409
4410int EvtConExc::getModeIndex( int m ) {
4411 for ( int i = 0; i < vmode.size(); i++ )
4412 {
4413 if ( m == vmode[i] ) return i;
4414 }
4415 std::cout << "EvtConExc::getModeIndex: no index in vmode found " << std::endl;
4416 abort();
4417}
4418
4419std::string EvtConExc::commandName() { return std::string( "ConExcPar" ); }
4420
4421void EvtConExc::command( std::string cmd ) {
4422 if ( ncommand == lcommand )
4423 {
4424 lcommand = 10 + 2 * lcommand;
4425 std::string* newcommands = new std::string[lcommand];
4426 int i;
4427 for ( i = 0; i < ncommand; i++ ) { newcommands[i] = commands[i]; }
4428 delete[] commands;
4429 commands = newcommands;
4430 }
4431 commands[ncommand] = cmd;
4432 ncommand++;
4433}
4434
4435void EvtConExc::store( EvtDecayBase* jsdecay ) {
4436
4437 if ( nconexcdecays == ntable )
4438 {
4439
4440 EvtDecayBasePtr* newconexcdecays = new EvtDecayBasePtr[2 * ntable + 10];
4441 int i;
4442 for ( i = 0; i < ntable; i++ ) { newconexcdecays[i] = conexcdecays[i]; }
4443 ntable = 2 * ntable + 10;
4444 delete[] conexcdecays;
4445 conexcdecays = newconexcdecays;
4446 }
4447
4448 conexcdecays[nconexcdecays++] = jsdecay;
4449}
4450
4451std::vector<std::string> EvtConExc::split( std::string str, std::string pattern ) {
4452 std::string::size_type pos;
4453 std::vector<std::string> result;
4454 str += pattern;
4455 int size = str.size();
4456
4457 for ( int i = 0; i < size; i++ )
4458 {
4459 pos = str.find( pattern, i );
4460 if ( pos < size )
4461 {
4462 std::string s = str.substr( i, pos - i );
4463 result.push_back( s );
4464 i = pos + pattern.size() - 1;
4465 }
4466 }
4467 return result;
4468}
4469
4470void EvtConExc::InitPars() {
4471 threshold = 0;
4472 beamEnergySpread = 0;
4473 std::string pattern = "=";
4474 for ( int i = 0; i < ncommand; i++ )
4475 {
4476 std::vector<std::string> result = split( commands[i], pattern );
4477 if ( result[0] == "threshold" ) { threshold = atof( result[1].c_str() ); }
4478 else if ( result[0] == "beamEnergySpread" )
4479 { beamEnergySpread = atof( result[1].c_str() ); }
4480 else
4481 {
4482 std::cout << "Your parameter in decay card \"" << result[0] << "\" incorect"
4483 << std::endl;
4484 std::cout << "Possible words: threshold, beamEnergySpread" << std::endl;
4485 abort();
4486 }
4487 }
4488}
4489
4490double EvtConExc::energySpread( double mu, double sigma ) {
4491 // mu: mean value in Gaussian
4492 // sigma: variance in Gaussian
4493rloop:
4494 int n = 12;
4495 double ri = 0;
4496 for ( int i = 0; i < n; i++ )
4497 {
4498 double pm = EvtRandom::Flat( 0., 1 );
4499 ri += pm;
4500 }
4501 double eta = sqrt( n * 12.0 ) * ( ri / 12 - 0.5 );
4502 double xsig = sigma * eta + mu; // smapling Gaussion value
4503 double mx0 = staxsection->getXlw();
4504 double mx1 = staxsection->getXup();
4505 if ( xsig < mx0 || xsig > mx1 ) goto rloop;
4506 return xsig;
4507}
4508
4509void EvtConExc::calAF( double myecms ) {
4510
4511 double _cms = myecms;
4512 double Esig = 0.0001; // sigularity engy
4513 double x = 2 * Egamcut / _cms;
4514 double s = _cms * _cms;
4515
4516 double mhdL = staxsection->getXlw();
4517 double EgamH = ( s - mhdL * mhdL ) / 2 / sqrt( s );
4518
4519 int nint = 50;
4520 int nmc = 1000;
4521 _xs0 = trapezoid( s, Esig, Egamcut, nint ); // std::cout<<"Int: _xs0= "<<_xs0<<std::endl;
4522 //--- sigularity point x from 0 to 0.0001GeV
4523 double xb = 2 * Esig / _cms;
4524 double sig_xs = SoftPhoton_xs( s, xb ) * ( staxsection->getXsection( _cms ) );
4525 _xs0 += sig_xs;
4526
4527 int Nstp = 598;
4528 double stp = ( EgamH - Egamcut ) / Nstp;
4529 for ( int i = 0; i < Nstp; i++ )
4530 { // points within Egam(max) ~ Egamcut=25MeV
4531 double Eg0 = EgamH - i * stp;
4532 double Eg1 = EgamH - ( i + 1 ) * stp;
4533 double xsi = trapezoid( s, Eg1, Eg0, nint ); // std::cout<<"Int xsi= " <<xsi<<std::endl;
4534
4535 AA[i] = ( Eg1 + Eg0 ) / 2;
4536 double mhi = sqrt( _cms * _cms - 2 * _cms * AA[i] );
4537 double mh2 = sqrt( _cms * _cms - 2 * _cms * Eg1 );
4538 double binwidth = mh2 - mhi;
4539 if ( i == 0 ) AF[0] = xsi;
4540 if ( i >= 1 ) AF[i] = AF[i - 1] + xsi;
4541 RadXS[i] = xsi / stp;
4542 }
4543 // add the interval 0~Esig
4544 AA[598] = Egamcut;
4545 AA[599] = 0; // AA[i]: E_gamma
4546 AF[598] = AF[597];
4547 AF[599] = AF[598] + _xs0;
4548 //--
4549 for ( int i = 0; i < 600; i++ ) { MH[i] = sqrt( _cms * _cms - 2 * _cms * AA[i] ); }
4550 //--debugging
4551 double bornXS = staxsection->getXsection( _cms );
4552 if ( bornXS == 0 )
4553 {
4554 std::cout << "EvtConExc::calAF: bad BornXS at " << _cms << " GeV" << std::endl;
4555 abort();
4556 }
4557 double fisr = AF[599] / bornXS;
4558 myFisr.push_back( fisr );
4559 // std::cout<<"fisr= "<<fisr<<std::endl;
4560}
4561
4562void EvtConExc::OutStaISR() {
4563 int ntot = myFisr.size();
4564 double mymu = 0;
4565 for ( int i = 0; i < ntot; i++ ) { mymu += myFisr[i]; }
4566 mymu /= ntot;
4567 double mysig = 0;
4568 for ( int i = 0; i < ntot; i++ ) { mysig += ( myFisr[i] - mymu ) * ( myFisr[i] - mymu ); }
4569 mysig /= ntot;
4570 mysig = sqrt( mysig );
4571 std::cout << "========= Iteration times over ISR factor: " << ntot << std::endl;
4572 std::cout << " ISR factor * VP factor= observedXS/BornXS: " << mymu << " + " << mysig
4573 << std::endl;
4574}
4575
4576double EvtConExc::trapezoid( double s, double El, double Eh,
4577 int n ) // integration with trapezoid method
4578{
4579 double xL = 2 * El / sqrt( s );
4580 double xH = 2 * Eh / sqrt( s );
4581 double sum = 0.0;
4582 double gaps = ( xH - xL ) / double( n ); // interval
4583 for ( int i = 0; i < n; i++ )
4584 {
4585 sum += ( gaps / 2.0 ) * ( Rad2difXs( s, xL + i * gaps, staxsection ) +
4586 Rad2difXs( s, xL + ( i + 1 ) * gaps, staxsection ) );
4587 }
4588 return sum;
4589}
4590
4591double EvtConExc::trapezoid( double s, double El, double Eh, int n,
4592 EvtXsection* myxs ) // integration with trapezoid method
4593{
4594 double xL = 2 * El / sqrt( s );
4595 double xH = 2 * Eh / sqrt( s );
4596 double sum = 0.0;
4597 double gaps = ( xH - xL ) / double( n ); // interval
4598 for ( int i = 0; i < n; i++ )
4599 {
4600 sum += ( gaps / 2.0 ) * ( Rad2difXs( s, xL + i * gaps, myxs ) +
4601 Rad2difXs( s, xL + ( i + 1 ) * gaps, myxs ) );
4602 }
4603 return sum;
4604}
p2
double p2[4]
Definition BabayagaNLO.cxx:41
p1
double p1[4]
Definition BabayagaNLO.cxx:41
mass
double mass
Definition CosmicGenerator.cxx:128
n
const Int_t n
Definition DataBase/tau_mode.c:57
Rad2difXs2
double Rad2difXs2(double *mhs)
Definition EvtConExc.cc:3420
hadr5n12_
void hadr5n12_(float *, float *, float *, float *, float *, float *)
divdif_
double divdif_(float *, float *, int *, float *, int *)
oa
std::ofstream oa
Definition EvtConExc.cc:199
Rad2difXs
double Rad2difXs(double *m)
Definition EvtConExc.cc:3389
Rad2difXs_er2
double Rad2difXs_er2(double *mhs)
Definition EvtConExc.cc:3436
dgauss_
double dgauss_(double(*fun)(double *), double *, double *, double *)
Rad2difXs_er
double Rad2difXs_er(double *m)
Definition EvtConExc.cc:3405
polint_
void polint_(float *, float *, int *, float *, float *)
EvtDecayBasePtr
EvtDecayBase * EvtDecayBasePtr
Definition EvtDecayTable.hh:32
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:67
alpha
double alpha
Definition EvtPhokharaDef.hh:17
pi
double pi
Definition EvtPhokharaDef.hh:17
me
double me
Definition EvtPhokharaDef.hh:17
xlow
DOUBLE_PRECISION xlow[20]
Definition EvtPhokharaDef.hh:56
mp
double mp
Definition EvtPhokharaDef.hh:17
mstp
int mstp[200]
Definition EvtPycont.cc:52
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
eps
EvtTensor3C eps(const EvtVector3R &v)
Definition EvtTensor3C.cc:269
s
XmlRpcServer s
Definition HelloServer.cpp:11
sin
double sin(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:185
cos
double cos(const BesAngle a)
Definition InstallArea/x86_64-el9-gcc13-dbg/include/MdcGeom/BesAngle.h:187
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition KarFin.h:34
phiee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rh2ee DOUBLE PRECISION omepi DOUBLE PRECISION om3ee DOUBLE PRECISION phiee
Definition RRes.h:37
rhoee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rhoee
Definition RRes.h:34
EvtConExc::checkdecay
bool checkdecay(EvtParticle *p)
Definition EvtConExc.cc:3887
EvtConExc::Li2
double Li2(double x)
Definition EvtConExc.cc:3480
EvtConExc::findMaxXS
void findMaxXS(EvtParticle *p)
Definition EvtConExc.cc:3085
EvtConExc::getMode
static int getMode
Definition EvtConExc.hh:196
EvtConExc::conexcmode
static int conexcmode
Definition EvtConExc.hh:214
EvtConExc::addNarrowRXS
double addNarrowRXS(double mhi, double binwidth)
Definition EvtConExc.cc:4127
EvtConExc::init
void init()
Definition EvtConExc.cc:248
EvtConExc::narrowRXS
double narrowRXS(double mxL, double mxH)
Definition EvtConExc.cc:4088
EvtConExc::VP_sampling
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:3215
EvtConExc::command
void command(std::string cmd)
Definition EvtConExc.cc:4421
EvtConExc::init_Br_ee
void init_Br_ee()
Definition EvtConExc.cc:3346
EvtConExc::gamHXSection_er
double gamHXSection_er(double El, double Eh)
Definition EvtConExc.cc:3062
EvtConExc::islgr
bool islgr(double *x, double *y, int n, double t)
Definition EvtConExc.cc:3508
EvtConExc::showResMass
void showResMass()
Definition EvtConExc.cc:3876
EvtConExc::lgr
double lgr(double *x, double *y, int n, double t)
Definition EvtConExc.cc:3486
EvtConExc::ReadVP
void ReadVP()
Definition EvtConExc.cc:4204
EvtConExc::baryonAng
double baryonAng(double mx)
Definition EvtConExc.cc:4045
EvtConExc::InitPars
void InitPars()
Definition EvtConExc.cc:4470
EvtConExc::Ros_xs
double Ros_xs(double mx, double bree, EvtId pid)
Definition EvtConExc.cc:3377
EvtConExc::OutStaISR
void OutStaISR()
Definition EvtConExc.cc:4562
EvtConExc::Rad1
double Rad1(double s, double x)
Definition EvtConExc.cc:3227
EvtConExc::staxsection
static EvtXsection * staxsection
Definition EvtConExc.hh:191
EvtConExc::myxsection
static EvtXsection * myxsection
Definition EvtConExc.hh:191
EvtConExc::SetP4Rvalue
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition EvtConExc.cc:3668
EvtConExc::getModeIndex
int getModeIndex(int m)
Definition EvtConExc.cc:4410
EvtConExc::meson_sampling
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:3203
EvtConExc::multimode
static int multimode
Definition EvtConExc.hh:214
EvtConExc::photonSampling
bool photonSampling(EvtParticle *part)
Definition EvtConExc.cc:4055
EvtConExc::ISR_ang_integrate
double ISR_ang_integrate(double x, double theta)
Definition EvtConExc.cc:3591
EvtConExc::SetMthr
static double SetMthr
Definition EvtConExc.hh:194
EvtConExc::xs_sampling
bool xs_sampling(double xs)
Definition EvtConExc.cc:3170
EvtConExc::selectMode
int selectMode(std::vector< int > vmod, double mhds)
Definition EvtConExc.cc:3716
EvtConExc::gamHXSection
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition EvtConExc.cc:2924
EvtConExc::Rad1difXs
double Rad1difXs(EvtParticle *p)
Definition EvtConExc.cc:3318
EvtConExc::sumExc
double sumExc(double mx)
Definition EvtConExc.cc:3903
EvtConExc::init_mode
void init_mode(int mode)
Definition EvtConExc.cc:847
EvtConExc::baryon_sampling
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:3184
EvtConExc::get_mode_index
int get_mode_index(int mode)
Definition EvtConExc.cc:4265
EvtConExc::LLr
double LLr(double *x, double *y, int n, double t)
Definition EvtConExc.cc:3531
EvtConExc::Rad2difXs
double Rad2difXs(EvtParticle *p)
Definition EvtConExc.cc:3263
EvtConExc::SoftPhoton_xs
double SoftPhoton_xs(double s, double b)
Definition EvtConExc.cc:3450
EvtConExc::angularSampling
bool angularSampling(EvtParticle *part)
Definition EvtConExc.cc:4008
EvtConExc::difgamXs
double difgamXs(EvtParticle *p)
Definition EvtConExc.cc:3129
EvtConExc::ISR_ang_sampling
double ISR_ang_sampling(double x)
Definition EvtConExc.cc:3611
EvtConExc::checkEvtRatio
void checkEvtRatio()
Definition EvtConExc.cc:4382
EvtConExc::initProbMax
void initProbMax()
Definition EvtConExc.cc:2301
EvtConExc::Mhad_sampling
double Mhad_sampling(double *x, double *y)
Definition EvtConExc.cc:3554
EvtConExc::_cms
static double _cms
Definition EvtConExc.hh:192
EvtConExc::writeDecayTabel
void writeDecayTabel()
Definition EvtConExc.cc:4296
EvtConExc::energySpread
double energySpread(double mu, double sigma)
Definition EvtConExc.cc:4490
EvtConExc::get_mode
std::vector< EvtId > get_mode(int mode)
Definition EvtConExc.cc:1574
EvtConExc::mlow
static double mlow
Definition EvtConExc.hh:239
EvtConExc::resetResMass
void resetResMass()
Definition EvtConExc.cc:3807
EvtConExc::getObsXsection
double getObsXsection(double mhds, int mode)
Definition EvtConExc.cc:4275
EvtConExc::mk_VXS
void mk_VXS(double Esig, double Egamcut, double EgamH, int midx)
Definition EvtConExc.cc:4223
EvtConExc::calAF
void calAF(double myecms)
Definition EvtConExc.cc:4509
EvtConExc::hadron_angle_sampling
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition EvtConExc.cc:2882
EvtConExc::~EvtConExc
virtual ~EvtConExc()
Definition EvtConExc.cc:205
EvtConExc::getResMass
void getResMass()
Definition EvtConExc.cc:3842
EvtConExc::gam_sampling
bool gam_sampling(EvtParticle *p)
Definition EvtConExc.cc:3154
EvtConExc::Rad2
double Rad2(double s, double x)
Definition EvtConExc.cc:3239
EvtConExc::selectMass
double selectMass()
Definition EvtConExc.cc:4141
EvtConExc::getVP
double getVP(double cms)
Definition EvtConExc.cc:4170
EvtConExc::PrintXS
void PrintXS(double mx)
Definition EvtConExc.cc:3949
EvtConExc::trapezoid
double trapezoid(double s, double a, double b, int n)
Definition EvtConExc.cc:4576
EvtConExc::getName
void getName(std::string &name)
Definition EvtConExc.cc:244
EvtConExc::decay
void decay(EvtParticle *p)
Definition EvtConExc.cc:2303
EvtConExc::split
std::vector< std::string > split(std::string str, std::string pattern)
Definition EvtConExc.cc:4451
EvtConExc::commandName
std::string commandName()
Definition EvtConExc.cc:4419
EvtConExc::XS_max
static double XS_max
Definition EvtConExc.hh:193
EvtConExc::Egam2Mhds
double Egam2Mhds(double Egam)
Definition EvtConExc.cc:4290
EvtConExc::mup
static double mup
Definition EvtConExc.hh:239
EvtConExc::clone
EvtDecayBase * clone()
Definition EvtConExc.cc:246
EvtConExc::SetP4
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition EvtConExc.cc:3639
EvtDecayBase::getArg
double getArg(int j)
Definition EvtDecayBase.cc:520
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition EvtDecayBase.hh:68
EvtDecayTable::inChannelList
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition EvtDecayTable.cc:831
EvtGlobalSet::ConExcPythia
static bool ConExcPythia
Definition EvtGlobalSet.hh:19
EvtHelSys::getHelAng
double getHelAng(int i)
Definition EvtHelSys.cc:52
EvtId
Definition EvtId.hh:27
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition EvtPDL.hh:59
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition EvtPDL.hh:61
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition EvtPDL.hh:43
EvtPDL::reSetWidth
static void reSetWidth(EvtId i, double width)
Definition EvtPDL.hh:75
EvtPDL::reSetMass
static void reSetMass(EvtId i, double mass)
Definition EvtPDL.hh:74
EvtPDL::evtIdFromStdHep
static EvtId evtIdFromStdHep(int stdhep)
Definition EvtPDL.cc:232
EvtPDL::name
static std::string name(EvtId i)
Definition EvtPDL.hh:70
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition EvtPDL.hh:56
EvtPDL::getSpinType
static EvtSpinType::spintype getSpinType(EvtId i)
Definition EvtPDL.hh:66
EvtPDL::getMass
static double getMass(EvtId i)
Definition EvtPDL.hh:44
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:272
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition EvtParticleFactory.cc:46
EvtParticle::setMass
void setMass(double m)
Definition EvtParticle.hh:368
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition EvtParticle.cc:1150
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getP4Lab
EvtVector4R getP4Lab()
Definition EvtParticle.cc:694
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:100
EvtParticle::setIntFlag
void setIntFlag(std::vector< int > vi)
Definition EvtParticle.hh:150
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:108
EvtParticle::printTree
void printTree() const
Definition EvtParticle.cc:883
EvtParticle::getNDaug
int getNDaug() const
Definition EvtParticle.cc:112
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:80
EvtParticle::mass
double mass() const
Definition EvtParticle.cc:114
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cc:69
EvtVector4R::mass
double mass() const
Definition EvtVector4R.cc:40
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:171
EvtVector4R::d3mag
double d3mag() const
Definition EvtVector4R.cc:168
EvtVector4R::mass2
double mass2() const
Definition EvtVector4R.hh:110
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:173
EvtVector4R::theta
double theta()
Definition EvtVector4R.cc:224
t
int t()
Definition t.c:1