Loading...
Searching...
No Matches
EvtConExc Class Reference
#include <EvtConExc.hh>
Inheritance diagram for EvtConExc:
Public Member Functions | |
| EvtConExc () | |
| virtual | ~EvtConExc () |
| void | getName (std::string &name) |
| EvtDecayBase * | clone () |
| void | initProbMax () |
| void | init () |
| void | init_Br_ee () |
| void | decay (EvtParticle *p) |
| void | init_mode (int mode) |
| double | gamHXSection (EvtParticle *p, double El, double Eh, int nmc=100000) |
| double | gamHXSection (double s, double El, double Eh, int nmc=100000) |
| double | gamHXSection (double El, double Eh) |
| double | gamHXSection (double El, double Eh, int mode) |
| double | gamHXSection_er (double El, double Eh) |
| void | findMaxXS (EvtParticle *p) |
| double | difgamXs (EvtParticle *p) |
| double | difgamXs (double mhds, double sintheta) |
| double | Mhad_sampling (double *x, double *y) |
| double | ISR_ang_integrate (double x, double theta) |
| double | ISR_ang_sampling (double x) |
| bool | hadron_angle_sampling (EvtVector4R ppi, EvtVector4R pcm) |
| void | SetP4 (EvtParticle *part, double mhdr, double xeng, double theta) |
| void | SetP4Rvalue (EvtParticle *part, double mhdr, double xeng, double theta) |
| bool | gam_sampling (EvtParticle *p) |
| bool | xs_sampling (double xs) |
| bool | xs_sampling (double xs, double xs1) |
| bool | baryon_sampling (EvtVector4R pcm, EvtVector4R pi) |
| bool | meson_sampling (EvtVector4R pcm, EvtVector4R pi) |
| bool | VP_sampling (EvtVector4R pcm, EvtVector4R pi) |
| bool | angularSampling (EvtParticle *part) |
| bool | photonSampling (EvtParticle *part) |
| double | baryonAng (double mx) |
| double | Rad1 (double s, double x) |
| double | Rad2 (double s, double x) |
| double | Rad2difXs (EvtParticle *p) |
| double | Rad2difXs (double s, double x) |
| double | Rad2difXs (double s, double x, EvtXsection *myxsection) |
| double | Rad1difXs (EvtParticle *p) |
| double | Ros_xs (double mx, double bree, EvtId pid) |
| double | Li2 (double x) |
| double | SoftPhoton_xs (double s, double b) |
| double | lgr (double *x, double *y, int n, double t) |
| bool | islgr (double *x, double *y, int n, double t) |
| double | LLr (double *x, double *y, int n, double t) |
| int | selectMode (std::vector< int > vmod, double mhds) |
| void | findMaxXS (double mhds) |
| bool | gam_sampling (double mhds, double sintheta) |
| void | resetResMass () |
| void | getResMass () |
| bool | checkdecay (EvtParticle *p) |
| double | sumExc (double mx) |
| void | PrintXS (double mx) |
| void | showResMass () |
| int | getNdaugs () |
| EvtId * | getDaugId () |
| int | getSelectedMode () |
| double | narrowRXS (double mxL, double mxH) |
| double | selectMass () |
| double | addNarrowRXS (double mhi, double binwidth) |
| void | ReadVP () |
| double | getVP (double cms) |
| void | mk_VXS (double Esig, double Egamcut, double EgamH, int midx) |
| int | get_mode_index (int mode) |
| double | getObsXsection (double mhds, int mode) |
| double | Egam2Mhds (double Egam) |
| std::vector< EvtId > | get_mode (int mode) |
| void | writeDecayTabel () |
| void | checkEvtRatio () |
| int | getModeIndex (int m) |
| std::string | commandName () |
| void | command (std::string cmd) |
| std::vector< std::string > | split (std::string str, std::string pattern) |
| void | InitPars () |
| double | energySpread (double mu, double sigma) |
| void | calAF (double myecms) |
| void | OutStaISR () |
| double | trapezoid (double s, double a, double b, int n) |
| double | trapezoid (double s, double El, double Eh, int n, EvtXsection *myxs) |
| Public Member Functions inherited from EvtDecayIncoherent | |
| void | makeDecay (EvtParticle *p) |
| virtual | ~EvtDecayIncoherent () |
| void | setDaughterSpinDensity (int daughter) |
| int | isDaughterSpinDensitySet (int daughter) |
| Public Member Functions inherited from EvtDecayBase | |
| double | getProbMax (double prob) |
| double | resetProbMax (double prob) |
| EvtDecayBase () | |
| virtual | ~EvtDecayBase () |
| virtual bool | matchingDecay (const EvtDecayBase &other) const |
| EvtId | getParentId () |
| double | getBranchingFraction () |
| void | disableCheckQ () |
| void | checkQ () |
| int | getNDaug () |
| EvtId * | getDaugs () |
| EvtId | getDaug (int i) |
| int | getNArg () |
| int | getPHOTOS () |
| void | setPHOTOS () |
| void | setVerbose () |
| void | setSummary () |
| double * | getArgs () |
| std::string * | getArgsStr () |
| double | getArg (int j) |
| std::string | getArgStr (int j) |
| std::string | getModelName () |
| int | getDSum () |
| int | summary () |
| int | verbose () |
| void | saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr) |
| void | printSummary () |
| void | setProbMax (double prbmx) |
| void | noProbMax () |
| void | checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1) |
| void | checkNDaug (int d1, int d2=-1) |
| void | checkSpinParent (EvtSpinType::spintype sp) |
| void | checkSpinDaughter (int d1, EvtSpinType::spintype sp) |
| virtual int | nRealDaughters () |
Static Public Attributes | |
| static EvtXsection * | myxsection |
| static EvtXsection * | staxsection |
| static double | _cms |
| static double | XS_max |
| static double | SetMthr = 0 |
| static int | getMode |
| static int | conexcmode = -1 |
| static int | multimode = 0 |
| static double | mlow = 0 |
| static double | mup = 0 |
Additional Inherited Members | |
| Static Public Member Functions inherited from EvtDecayBase | |
| static void | findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10]) |
| static void | findMass (EvtParticle *p) |
| static double | findMaxMass (EvtParticle *p) |
| Protected Member Functions inherited from EvtDecayBase | |
| bool | daugsDecayedByParentModel () |
| Protected Attributes inherited from EvtDecayBase | |
| bool | _daugsDecayedByParentModel |
Detailed Description
Definition at line 138 of file EvtConExc.hh.
Constructor & Destructor Documentation
◆ EvtConExc()
|
inline |
◆ ~EvtConExc()
|
virtual |
Definition at line 205 of file EvtConExc.cc.
205 {
206 // if(_mode==74110)checkEvtRatio();
208 if ( mydbg )
209 {
210 // xs->Write();
211 myfile->Write();
212 }
213
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}
Member Function Documentation
◆ addNarrowRXS()
| double EvtConExc::addNarrowRXS | ( | double | mhi, |
| double | binwidth ) |
Definition at line 4127 of file EvtConExc.cc.
4127 {
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}
◆ angularSampling()
| bool EvtConExc::angularSampling | ( | EvtParticle * | part | ) |
Definition at line 4008 of file EvtConExc.cc.
4008 {
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();
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 */
4033 double xratio = ang / angmax;
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}
double cos(const BesAngle a)
Referenced by decay().
◆ baryon_sampling()
| bool EvtConExc::baryon_sampling | ( | EvtVector4R | pcm, |
| EvtVector4R | pi ) |
Definition at line 3184 of file EvtConExc.cc.
3184 {
3186 double theta = angles.getHelAng( 1 );
3187 double phi = angles.getHelAng( 2 );
3189
3190 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
3191 // std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
3196 double myang;
3198 else { myang = 1; }
3199 if ( pm < ang / myang ) { return true; }
3200 else { return false; }
3201}
Referenced by hadron_angle_sampling().
◆ baryonAng()
| double EvtConExc::baryonAng | ( | double | mx | ) |
Definition at line 4045 of file EvtConExc.cc.
4045 {
4047 double u = 0.938 / mx;
4048 u = u * u;
4049 double u2 = ( 1 + u ) * ( 1 + u );
4050 double uu = u * ( 1 + 6 * u );
4053}
Referenced by baryon_sampling().
◆ calAF()
| void EvtConExc::calAF | ( | double | myecms | ) |
Definition at line 4509 of file EvtConExc.cc.
4509 {
4510
4512 double Esig = 0.0001; // sigularity engy
4515
4518
4519 int nint = 50;
4520 int nmc = 1000;
4522 //--- sigularity point x from 0 to 0.0001GeV
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;
4534
4535 AA[i] = ( Eg1 + Eg0 ) / 2;
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 //--
4550 //--debugging
4552 if ( bornXS == 0 )
4553 {
4555 abort();
4556 }
4557 double fisr = AF[599] / bornXS;
4558 myFisr.push_back( fisr );
4559 // std::cout<<"fisr= "<<fisr<<std::endl;
4560}
double trapezoid(double s, double a, double b, int n)
Definition EvtConExc.cc:4576
Referenced by decay().
◆ checkdecay()
| bool EvtConExc::checkdecay | ( | EvtParticle * | p | ) |
Definition at line 3887 of file EvtConExc.cc.
3887 {
3889 double massflag = 1;
3890 for ( int i = 0; i < nson; i++ )
3891 {
3894 }
3895 if ( massflag == 0 )
3896 {
3897 std::cout << "Zero mass!" << std::endl;
3898 return 0;
3899 }
3900 else { return 1; }
3901}
Referenced by decay().
◆ checkEvtRatio()
| void EvtConExc::checkEvtRatio | ( | ) |
Definition at line 4382 of file EvtConExc.cc.
4382 {
4385 //
4387 double xscon = VXS[im][599];
4388 double xscon0 = xscon;
4392 for ( int i = 0; i < vmode.size(); i++ )
4393 {
4394 // std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
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 {
4404 if ( vmode[i] == 74110 ) continue;
4406 }
4408}
std::ofstream ofstream
◆ clone()
|
virtual |
◆ command()
|
virtual |
Reimplemented from EvtDecayBase.
Definition at line 4421 of file EvtConExc.cc.
4421 {
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}
◆ commandName()
|
virtual |
Reimplemented from EvtDecayBase.
Definition at line 4419 of file EvtConExc.cc.
4419{ return std::string( "ConExcPar" ); }
◆ decay()
|
virtual |
Implements EvtDecayBase.
Definition at line 2303 of file EvtConExc.cc.
2303 {
2304 // std::cout<<"_cms= "<<_cms<<" mode= "<<_mode<<std::endl;
2307 {
2309 << std::endl;
2310 abort();
2311 }
2312 //-- to read _cms from database or to consider beam energy spread
2314 bool beamflag = 0;
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
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
2369
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;
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
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 );
2410 di->setMass( xmi );
2411 totMass += xmi;
2412 }
2414 p->initializePhaseSpace( _ndaugs, daugs );
2416 vhds = p->getP4();
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 );
2434 di->setMass( xmi );
2435 totMass += xmi;
2436 }
2437 // std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<"
2438 // "<<checkdecay(p)<<std::endl;
2440 p->initializePhaseSpace( 2, mydaugs );
2441
2444 ;
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;
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;
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;
2473 init_mode( mymode );
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 );
2506 di->setMass( xmi );
2507 totMass += xmi;
2508 }
2510 // std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2513 // set the ISR photon angular distribution
2515
2516 if ( maxflag == 0 )
2517 {
2518 findMaxXS( mhds );
2519 maxflag = 1;
2520 }
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 );
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;
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
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 );
2587 di->setMass( xmi );
2588 totMass += xmi;
2589 }
2591 p->initializePhaseSpace( _ndaugs, daugs );
2593 vhds = p->getP4();
2594 bool byn_ang;
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 );
2618 di->setMass( xmi );
2619 totMass += xmi;
2620 }
2621 // std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<"
2622 // "<<checkdecay(p)<<std::endl;
2624 p->initializePhaseSpace( 2, mydaugs );
2625
2628 ;
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;
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;
2649 std::cout << "A selected mode is " << mymode << " with Mhds= " << mhds
2650 << std::endl; // debugging
2651 _selectedMode = mymode;
2654 init_mode( mymode );
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 );
2687 di->setMass( xmi );
2688 totMass += xmi;
2689 }
2691 // std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2694 // set the ISR photon angular distribution
2696
2697 //--debugging
2698 // p->printTree();
2699
2700 if ( maxflag == 0 )
2701 {
2702 findMaxXS( mhds );
2703 maxflag = 1;
2704 }
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 );
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 {
2745
2746 if ( pm < _xs0 / ( _xs0 + _xs1 ) )
2747 { // with a very soft ISR photon
2749 mhds = _cms - 0.0001;
2750 }
2751 else
2752 {
2753 mhds = Mhad_sampling( MH, AF );
2754 EvtPDL::reSetMass( gid, mhds );
2755 }
2756 // debugging
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();
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 }
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
2786 { // exclude DIY case
2788 {
2789 iflag = 1;
2790 xsratio = 0;
2791 } // only generating gamma hadrons mode
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 );
2810 di->setMass( xmi );
2811 summassx += xmi;
2812 // std::cout<<"PID and mass "<<di->getId()<<" "<<xmi<<std::endl;
2813 }
2815 angular_hadron:
2816 p->initializePhaseSpace( _ndaugs, daugs );
2817 bool byn_ang;
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
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*
2838 EvtPDL::reSetMass( gid, mhdr );
2839 int ndaugs = 2;
2840 EvtId mydaugs[2];
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 );
2854 // if(EvtPDL::name(di->getId())=="vhdr") di->setVectorSpinDensity();
2855 di->setMass( xmi );
2856 }
2857 p->initializePhaseSpace( 2, mydaugs );
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}
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition EvtConExc.cc:3668
int selectMode(std::vector< int > vmod, double mhds)
Definition EvtConExc.cc:3716
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition EvtConExc.cc:2882
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition EvtConExc.cc:3639
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1042
◆ difgamXs() [1/2]
| double EvtConExc::difgamXs | ( | double | mhds, |
| double | sintheta ) |
Definition at line 3704 of file EvtConExc.cc.
3705 { // dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3707 double sin2theta = sintheta * sintheta;
3713 return difxs;
3714}
◆ difgamXs() [2/2]
| double EvtConExc::difgamXs | ( | EvtParticle * | p | ) |
Definition at line 3129 of file EvtConExc.cc.
3129 { // dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3130 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
3135 double egam = pgam.get( 0 );
3136 double sin2theta = pgam.get( 3 ) / pgam.d3mag();
3137 sin2theta = 1 - sin2theta * sin2theta;
3138
3144 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<",
3145 // "<<sin2theta<<", "<<cms<<",
3146 // "<<x<<", "<<wsx<<std::endl;
3147
3149 double difxs = 2 * mhs / cms / cms * wsx * xs;
3150 differ = 2 * mhs / cms / cms * wsx * ( myxsection->getErr( mhs ) );
3151 return difxs;
3152}
Referenced by findMaxXS(), gam_sampling(), gam_sampling(), and gamHXSection().
◆ Egam2Mhds()
| double EvtConExc::Egam2Mhds | ( | double | Egam | ) |
◆ energySpread()
| double EvtConExc::energySpread | ( | double | mu, |
| double | sigma ) |
Definition at line 4490 of file EvtConExc.cc.
4490 {
4491 // mu: mean value in Gaussian
4492 // sigma: variance in Gaussian
4493rloop:
4495 double ri = 0;
4497 {
4499 ri += pm;
4500 }
4502 double xsig = sigma * eta + mu; // smapling Gaussion value
4505 if ( xsig < mx0 || xsig > mx1 ) goto rloop;
4506 return xsig;
4507}
◆ findMaxXS() [1/2]
| void EvtConExc::findMaxXS | ( | double | mhds | ) |
Definition at line 3688 of file EvtConExc.cc.
3688 { // dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3689 maxXS = -1;
3690 for ( int i = 0; i < 90000; i++ )
3691 {
3700 if ( difxs > maxXS ) maxXS = difxs;
3701 }
3702}
◆ findMaxXS() [2/2]
| void EvtConExc::findMaxXS | ( | EvtParticle * | p | ) |
Definition at line 3085 of file EvtConExc.cc.
3085 {
3086 // std::cout<<"nmc= "<<nmc<<std::endl;
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 );
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
3120 EvtVector4R p4gam = gamH->getDaug( _ndaugs )->getP4();
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}
*********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
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition EvtParticleFactory.cc:46
◆ gam_sampling() [1/2]
| bool EvtConExc::gam_sampling | ( | double | mhds, |
| double | sintheta ) |
Definition at line 3162 of file EvtConExc.cc.
3162 {
3165 double xsratio = xs / maxXS;
3166 if ( pm < xsratio ) { return true; }
3167 else { return false; }
3168}
◆ gam_sampling() [2/2]
| bool EvtConExc::gam_sampling | ( | EvtParticle * | p | ) |
Definition at line 3154 of file EvtConExc.cc.
3154 { // photon angular distribution sampling
3157 double xsratio = xs / maxXS;
3158 if ( pm < xsratio ) { return true; }
3159 else { return false; }
3160}
Referenced by decay().
◆ gamHXSection() [1/4]
| double EvtConExc::gamHXSection | ( | double | El, |
| double | Eh ) |
Definition at line 3017 of file EvtConExc.cc.
3017 { // using Gaussian integration
3018 // subroutine from Cern lib
3019 std::cout << " " << std::endl;
3023 double Rad2Xs;
3026 if ( Rad2Xs == 0 )
3027 {
3029 {
3030 eps += 0.01;
3033 if ( !Rad2Xs ) return Rad2Xs;
3034 }
3035 }
3036 return Rad2Xs; // the leading second order correction
3037}
double dgauss_(double(*fun)(double *), double *, double *, double *)
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:83
◆ gamHXSection() [2/4]
| double EvtConExc::gamHXSection | ( | double | El, |
| double | Eh, | ||
| int | mode ) |
Definition at line 3039 of file EvtConExc.cc.
3040 { // using Gaussian integration subroutine from Cern
3041 // lib
3042 std::cout << " " << std::endl;
3046 double Rad2Xs;
3049 if ( Rad2Xs == 0 )
3050 {
3052 {
3053 eps += 0.01;
3056 if ( !Rad2Xs ) return Rad2Xs;
3057 }
3058 }
3059 return Rad2Xs; // the leading second order correction
3060}
◆ gamHXSection() [3/4]
| double EvtConExc::gamHXSection | ( | double | s, |
| double | El, | ||
| double | Eh, | ||
| int | nmc = 100000 ) |
Definition at line 2984 of file EvtConExc.cc.
2985 { // 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;
2999 for ( int i = 0; i < nmc; i++ )
3000 { // It is found that the narrow resonance can not included in this MC integartion
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}
◆ gamHXSection() [4/4]
| double EvtConExc::gamHXSection | ( | EvtParticle * | p, |
| double | El, | ||
| double | Eh, | ||
| int | nmc = 100000 ) |
Definition at line 2924 of file EvtConExc.cc.
2925 { // El, Eh : the lower and higher energy of
2926 // radiative photons
2927 // std::cout<<"nmc= "<<nmc<<std::endl;
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 );
2950 di->setMass( xmi );
2951 }
2952 loop:
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
2966 Rad2Xs += Rad2difXs( gamH );
2967 if ( thexs > maxXS ) { maxXS = thexs; }
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}
Referenced by init().
◆ gamHXSection_er()
| double EvtConExc::gamHXSection_er | ( | double | El, |
| double | Eh ) |
Definition at line 3062 of file EvtConExc.cc.
3062 { // using Gaussian integration
3063 // subroutine from Cern lib
3064 std::cout << " " << std::endl;
3068 double Rad2Xs;
3071 if ( Rad2Xs == 0 )
3072 {
3074 {
3075 eps += 0.01;
3078 if ( !Rad2Xs ) return Rad2Xs;
3079 ;
3080 }
3081 }
3082 return Rad2Xs; // the leading second order correction
3083}
◆ get_mode()
| std::vector< EvtId > EvtConExc::get_mode | ( | int | mode | ) |
Definition at line 1574 of file EvtConExc.cc.
1574 {
1575 int _ndaugs;
1576 EvtId daugs[100];
1577 if ( mode == 0 )
1578 {
1579 _ndaugs = 2;
1582 }
1583 else if ( mode == 1 )
1584 {
1585 _ndaugs = 2;
1588 }
1589 else if ( mode == 2 )
1590 {
1591 _ndaugs = 2;
1594 }
1595 else if ( mode == 3 )
1596 {
1597 _ndaugs = 2;
1600 }
1601 else if ( mode == 4 )
1602 {
1603 _ndaugs = 2;
1606 }
1607 else if ( mode == 5 )
1608 {
1609 _ndaugs = 2;
1612 }
1613 else if ( mode == 6 )
1614 {
1615 _ndaugs = 2;
1618 }
1619 else if ( mode == 7 )
1620 {
1621 _ndaugs = 3;
1625 }
1626 else if ( mode == 8 )
1627 {
1628 _ndaugs = 3;
1632 }
1633 else if ( mode == 9 )
1634 {
1635 _ndaugs = 3;
1639 }
1640 else if ( mode == 10 )
1641 {
1642 _ndaugs = 3;
1646 }
1647 else if ( mode == 11 )
1648 {
1649 _ndaugs = 3;
1653 }
1654 else if ( mode == 12 )
1655 {
1656 _ndaugs = 4;
1661 }
1662 else if ( mode == 13 )
1663 {
1664 _ndaugs = 4;
1669 }
1670 else if ( mode == 14 )
1671 {
1672 _ndaugs = 4;
1677 }
1678 else if ( mode == 15 )
1679 {
1680 _ndaugs = 4;
1685 }
1686 else if ( mode == 16 )
1687 {
1688 _ndaugs = 4;
1693 }
1694 else if ( mode == 17 )
1695 {
1696 _ndaugs = 5;
1702 }
1703 else if ( mode == 18 )
1704 {
1705 _ndaugs = 5;
1711 }
1712 else if ( mode == 19 )
1713 {
1714 _ndaugs = 5;
1720 }
1721 else if ( mode == 20 )
1722 {
1723 _ndaugs = 5;
1729 }
1730 else if ( mode == 21 )
1731 {
1732 _ndaugs = 6;
1739 }
1740 else if ( mode == 22 )
1741 {
1742 _ndaugs = 6;
1749 }
1750 else if ( mode == 23 )
1751 {
1752 _ndaugs = 2;
1755 }
1756 else if ( mode == 24 )
1757 {
1758 _ndaugs = 2;
1761 }
1762 else if ( mode == 25 )
1763 {
1764 _ndaugs = 2;
1767 }
1768 else if ( mode == 26 )
1769 {
1770 _ndaugs = 2;
1773 }
1774 else if ( mode == 27 )
1775 {
1776 _ndaugs = 2;
1779 }
1780 else if ( mode == 28 )
1781 {
1782 _ndaugs = 3;
1786 }
1787 else if ( mode == 29 )
1788 {
1789 _ndaugs = 3;
1793 }
1794 else if ( mode == 30 )
1795 {
1796 _ndaugs = 3;
1800 }
1801 else if ( mode == 31 )
1802 {
1803 _ndaugs = 3;
1807 }
1808 else if ( mode == 32 )
1809 {
1810 _ndaugs = 3;
1814 }
1815 else if ( mode == 33 )
1816 {
1817 _ndaugs = 3;
1821 }
1822 else if ( mode == 34 )
1823 {
1824 _ndaugs = 3;
1828 }
1829 else if ( mode == 35 )
1830 {
1831 _ndaugs = 3;
1835 }
1836 else if ( mode == 36 )
1837 {
1838 _ndaugs = 2;
1841 }
1842 else if ( mode == 37 )
1843 {
1844 _ndaugs = 3;
1848 }
1849 else if ( mode == 38 )
1850 {
1851 _ndaugs = 3;
1855 }
1856 else if ( mode == 39 )
1857 {
1858 _ndaugs = 2;
1861 }
1862 else if ( mode == 40 )
1863 {
1864 _ndaugs = 3;
1868 }
1869 else if ( mode == 41 )
1870 {
1871 _ndaugs = 3;
1875 }
1876 else if ( mode == 42 )
1877 {
1878 _ndaugs = 3;
1882 }
1883 else if ( mode == 43 )
1884 {
1885 _ndaugs = 4;
1890 }
1891 else if ( mode == 44 )
1892 {
1893 _ndaugs = 2;
1896 }
1897 else if ( mode == 45 )
1898 {
1899 _ndaugs = 2;
1902 }
1903 else if ( mode == 46 )
1904 {
1905 _ndaugs = 2;
1908 }
1909 else if ( mode == 47 )
1910 {
1911 _ndaugs = 2;
1914 }
1915 else if ( mode == 48 )
1916 {
1917 _ndaugs = 2;
1920 }
1921 else if ( mode == 49 )
1922 {
1923 _ndaugs = 3;
1927 }
1928 else if ( mode == 50 )
1929 {
1930 _ndaugs = 3;
1934 }
1935 else if ( mode == 51 )
1936 {
1937 _ndaugs = 3;
1941 }
1942 else if ( mode == 52 )
1943 {
1944 _ndaugs = 2;
1947 }
1948 else if ( mode == 59 )
1949 {
1950 _ndaugs = 3;
1954 }
1955 else if ( mode == 60 )
1956 {
1957 _ndaugs = 3;
1961 }
1962 else if ( mode == 61 )
1963 {
1964 _ndaugs = 3;
1968 }
1969 else if ( mode == 62 )
1970 {
1971 _ndaugs = 3;
1975 }
1976 else if ( mode == 63 )
1977 {
1978 _ndaugs = 3;
1982 }
1983 else if ( mode == 64 )
1984 {
1985 _ndaugs = 3;
1989 }
1990 else if ( mode == 65 )
1991 {
1992 _ndaugs = 3;
1996 }
1997 else if ( mode == 66 )
1998 {
1999 _ndaugs = 3;
2003 }
2004 else if ( mode == 67 )
2005 {
2006 _ndaugs = 2;
2009 }
2010 else if ( mode == 68 )
2011 {
2012 _ndaugs = 2;
2015 }
2016 else if ( mode == 69 )
2017 {
2018 _ndaugs = 2;
2021 }
2022 else if ( mode == 70 )
2023 {
2024 _ndaugs = 2;
2027 }
2028 else if ( mode == 71 )
2029 {
2030 _ndaugs = 2;
2033 }
2034 else if ( mode == 72 )
2035 {
2036 _ndaugs = 2;
2039 }
2040 else if ( mode == 73 )
2041 {
2042 _ndaugs = 2;
2045 }
2046 else if ( mode == 74 )
2047 {
2048 _ndaugs = 2;
2051 }
2052 else if ( mode == 75 )
2053 {
2054 _ndaugs = 3;
2058 }
2059 else if ( mode == 76 )
2060 {
2061 _ndaugs = 3;
2065 }
2066 else if ( mode == 77 )
2067 {
2068 _ndaugs = 3;
2072 }
2073 else if ( mode == 78 )
2074 {
2075 _ndaugs = 3;
2079 }
2080 else if ( mode == 79 )
2081 {
2082 _ndaugs = 3;
2086 }
2087 else if ( mode == 80 )
2088 {
2089 _ndaugs = 2;
2092 }
2093 else if ( mode == 81 )
2094 {
2095 _ndaugs = 3;
2099 }
2100 else if ( mode == 82 )
2101 {
2102 _ndaugs = 3;
2106 }
2107 else if ( mode == 83 )
2108 {
2109 _ndaugs = 3;
2113 }
2114 else if ( mode == 84 )
2115 {
2116 _ndaugs = 3;
2120 }
2121 else if ( mode == 85 )
2122 {
2123 _ndaugs = 2;
2126 }
2127 else if ( mode == 90 )
2128 {
2129 _ndaugs = 3;
2133 }
2134 else if ( mode == 91 )
2135 {
2136 _ndaugs = 3;
2140 }
2141 else if ( mode == 92 )
2142 {
2143 _ndaugs = 3;
2147 }
2148 else if ( mode == 93 )
2149 {
2150 _ndaugs = 2;
2153 }
2154 else if ( mode == 94 )
2155 {
2156 _ndaugs = 2;
2159 }
2160 else if ( mode == 95 )
2161 {
2162 _ndaugs = 2;
2165 }
2166 else if ( mode == 96 )
2167 {
2168 _ndaugs = 2;
2171 }
2172 else if ( mode == 97 )
2173 {
2174 _ndaugs = 3;
2178 }
2179 else if ( mode == 111 )
2180 {
2181 _ndaugs = 2;
2184 }
2185 else if ( mode == 112 )
2186 {
2187 _ndaugs = 2;
2190 }
2191 else if ( mode == 74100 )
2192 {
2193 _ndaugs = 2;
2196 }
2197 else if ( mode == 74101 )
2198 {
2199 _ndaugs = 2;
2202 }
2203 else if ( mode == 74102 )
2204 {
2205 _ndaugs = 2;
2208 }
2209 else if ( mode == 74103 )
2210 {
2211 _ndaugs = 2;
2214 }
2215 else if ( mode == 74104 )
2216 {
2217 _ndaugs = 2;
2220 }
2221 else if ( mode == 74105 )
2222 {
2223 _ndaugs = 2;
2226 }
2227 else if ( mode == 74106 )
2228 {
2229 _ndaugs = 2;
2232 }
2233 else if ( mode == 74107 )
2234 {
2235 _ndaugs = 2;
2238 }
2239 else if ( mode == 74110 )
2240 {
2241 _ndaugs = 2;
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]; }
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;
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;
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;
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}
Referenced by selectMode().
◆ get_mode_index()
| int EvtConExc::get_mode_index | ( | int | mode | ) |
Definition at line 4265 of file EvtConExc.cc.
4265 {
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}
Referenced by getObsXsection().
◆ getDaugId()
|
inline |
Definition at line 212 of file EvtConExc.hh.
212{ return daugs; }
◆ getModeIndex()
| int EvtConExc::getModeIndex | ( | int | m | ) |
Definition at line 4410 of file EvtConExc.cc.
4410 {
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}
Referenced by checkEvtRatio(), and writeDecayTabel().
◆ getName()
|
virtual |
◆ getNdaugs()
|
inline |
Definition at line 211 of file EvtConExc.hh.
211{ return _ndaugs; }
◆ getObsXsection()
| double EvtConExc::getObsXsection | ( | double | mhds, |
| int | mode ) |
Definition at line 4275 of file EvtConExc.cc.
4275 {
4276 double hwid = ( AA[0] - AA[1] ) / 2.;
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}
◆ getResMass()
| void EvtConExc::getResMass | ( | ) |
Definition at line 3842 of file EvtConExc.cc.
3842 {
3844 mjsi = EvtPDL::getMeanMass( myres );
3845 wjsi = EvtPDL::getWidth( myres );
3846
3848 mpsip = EvtPDL::getMeanMass( myres );
3849 wpsip = EvtPDL::getWidth( myres );
3850
3852 mpsipp = EvtPDL::getMeanMass( myres );
3853 wpsipp = EvtPDL::getWidth( myres );
3854
3856 mphi = EvtPDL::getMeanMass( myres );
3857 wphi = EvtPDL::getWidth( myres );
3858
3860 momega = EvtPDL::getMeanMass( myres );
3861 womega = EvtPDL::getWidth( myres );
3862
3864 mrho0 = EvtPDL::getMeanMass( myres );
3865 wrho0 = EvtPDL::getWidth( myres );
3866
3868 mrho3s = EvtPDL::getMeanMass( myres );
3869 wrho3s = EvtPDL::getWidth( myres );
3870
3872 momega2s = EvtPDL::getMeanMass( myres );
3873 womega2s = EvtPDL::getWidth( myres );
3874}
Referenced by init().
◆ getSelectedMode()
|
inline |
Definition at line 213 of file EvtConExc.hh.
213{ return _selectedMode; }
◆ getVP()
| double EvtConExc::getVP | ( | double | cms | ) |
Definition at line 4170 of file EvtConExc.cc.
4170 {
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 );
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}
Referenced by init(), Rad2(), and SoftPhoton_xs().
◆ hadron_angle_sampling()
| bool EvtConExc::hadron_angle_sampling | ( | EvtVector4R | ppi, |
| EvtVector4R | pcm ) |
Definition at line 2882 of file EvtConExc.cc.
2882 {
2883 EvtVector4R pbst = -1 * pcm;
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
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
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 {
2903 return msn_ang;
2904 }
2905 else if ( _mode == -2 )
2906 {
2910 {
2912 return msn_ang;
2913 }
2916 {
2918 return msn_ang;
2919 }
2920 }
2921 return true;
2922}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:67
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:3203
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:3184
Referenced by decay().
◆ init()
|
virtual |
Reimplemented from EvtDecayBase.
Definition at line 248 of file EvtConExc.cc.
248 {
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;
283 {
288 VISR = 1;
289 }
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();
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;
309 {
312 }
313 std::cout << std::endl;
314 }
316 {
317 radflag = 0;
318 mydbg = false;
319 _mode = getArg( 0 );
320 nson = getNArg() - 1;
323 }
325 {
326 _mode = getArg( 0 );
327 multimode = 1;
328 vmd.clear();
330 std::cout << " multi-exclusive mode " << std::endl;
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 }
340 {
341 _mode = getArg( 0 );
342 radflag = 0;
343 mydbg = false;
344 }
346 {
347 _mode = getArg( 0 );
348 radflag = getArg( 1 );
349 mydbg = false;
350 }
352 {
353 _mode = getArg( 0 );
354 radflag = getArg( 1 );
355 mydbg = getArg( 2 );
356 }
357 else
358 {
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
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
392 EvtPDL::reSetWidth( parentId, 0.0 );
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;
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;
410 if ( _mode == -1 )
411 {
412 myxsection->setBW( pdgcode );
413 staxsection->setBW( pdgcode );
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 }
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
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
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();
453
454 ISRXS.clear();
455 ISRM.clear();
456 ISRFLAG.clear();
457 for ( int ii = 0; ii < 3; ii++ )
458 {
460 if ( mx < mR || _mode != 74110 ) continue;
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
478 int nmc = 1000000;
480 _er0 = _er1; // stored when calculate _xs0
481 std::cout << "_er0= " << _er0 << std::endl;
483 std::cout << "_xs1= " << _xs0 << std::endl;
484 double xs1_err = _er1;
485
486 //--- sigularity point x from 0 to 0.0001GeV
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);
504 AA[i] = ( Eg1 + Eg0 ) / 2;
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
522 << std::endl;
523 for ( int i = 0; i < vmode.size(); i++ )
524 {
526 }
528 // after mk_VXS, restore the myxsection to _mode
529 if ( _mode != -100 )
530 {
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 {
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
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;
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;
678 //--- for debugging
679 if ( mydbg && _mode != 74110 )
680 {
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 {
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;
714 TGraphErrors* Gdt = new TGraphErrors( nd, xx, yy, xer, yer );
715
719 for ( int i = 0; i < 600; i++ )
720 {
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 {
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
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}
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition EvtConExc.cc:2924
void mk_VXS(double Esig, double Egamcut, double EgamH, int midx)
Definition EvtConExc.cc:4223
◆ init_Br_ee()
| void EvtConExc::init_Br_ee | ( | ) |
Definition at line 3346 of file EvtConExc.cc.
3346 {
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;
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}
Referenced by init().
◆ init_mode()
| void EvtConExc::init_mode | ( | int | mode | ) |
Definition at line 847 of file EvtConExc.cc.
847 {
848 if ( mode == 0 )
849 {
850 _ndaugs = 2;
853 }
854 else if ( mode == 1 )
855 {
856 _ndaugs = 2;
859 }
860 else if ( mode == 2 )
861 {
862 _ndaugs = 2;
865 }
866 else if ( mode == 3 )
867 {
868 _ndaugs = 2;
871 }
872 else if ( mode == 4 )
873 {
874 _ndaugs = 2;
877 }
878 else if ( mode == 5 )
879 {
880 _ndaugs = 2;
883 }
884 else if ( mode == 6 )
885 {
886 _ndaugs = 2;
889 }
890 else if ( mode == 7 )
891 {
892 _ndaugs = 3;
896 }
897 else if ( mode == 8 )
898 {
899 _ndaugs = 3;
903 }
904 else if ( mode == 9 )
905 {
906 _ndaugs = 3;
910 }
911 else if ( mode == 10 )
912 {
913 _ndaugs = 3;
917 }
918 else if ( mode == 11 )
919 {
920 _ndaugs = 3;
924 }
925 else if ( mode == 12 )
926 {
927 _ndaugs = 4;
932 }
933 else if ( mode == 13 )
934 {
935 _ndaugs = 4;
940 }
941 else if ( mode == 14 )
942 {
943 _ndaugs = 4;
948 }
949 else if ( mode == 15 )
950 {
951 _ndaugs = 4;
956 }
957 else if ( mode == 16 )
958 {
959 _ndaugs = 4;
964 }
965 else if ( mode == 17 )
966 {
967 _ndaugs = 5;
973 }
974 else if ( mode == 18 )
975 {
976 _ndaugs = 5;
982 }
983 else if ( mode == 19 )
984 {
985 _ndaugs = 5;
991 }
992 else if ( mode == 20 )
993 {
994 _ndaugs = 5;
1000 }
1001 else if ( mode == 21 )
1002 {
1003 _ndaugs = 6;
1010 }
1011 else if ( mode == 22 )
1012 {
1013 _ndaugs = 6;
1020 }
1021 else if ( mode == 23 )
1022 {
1023 _ndaugs = 2;
1026 }
1027 else if ( mode == 24 )
1028 {
1029 _ndaugs = 2;
1032 }
1033 else if ( mode == 25 )
1034 {
1035 _ndaugs = 2;
1038 }
1039 else if ( mode == 26 )
1040 {
1041 _ndaugs = 2;
1044 }
1045 else if ( mode == 27 )
1046 {
1047 _ndaugs = 2;
1050 }
1051 else if ( mode == 28 )
1052 {
1053 _ndaugs = 3;
1057 }
1058 else if ( mode == 29 )
1059 {
1060 _ndaugs = 3;
1064 }
1065 else if ( mode == 30 )
1066 {
1067 _ndaugs = 3;
1071 }
1072 else if ( mode == 31 )
1073 {
1074 _ndaugs = 3;
1078 }
1079 else if ( mode == 32 )
1080 {
1081 _ndaugs = 3;
1085 }
1086 else if ( mode == 33 )
1087 {
1088 _ndaugs = 3;
1092 }
1093 else if ( mode == 34 )
1094 {
1095 _ndaugs = 3;
1099 }
1100 else if ( mode == 35 )
1101 {
1102 _ndaugs = 3;
1106 }
1107 else if ( mode == 36 )
1108 {
1109 _ndaugs = 2;
1112 }
1113 else if ( mode == 37 )
1114 {
1115 _ndaugs = 3;
1119 }
1120 else if ( mode == 38 )
1121 {
1122 _ndaugs = 3;
1126 }
1127 else if ( mode == 39 )
1128 {
1129 _ndaugs = 2;
1132 }
1133 else if ( mode == 40 )
1134 {
1135 _ndaugs = 3;
1139 }
1140 else if ( mode == 41 )
1141 {
1142 _ndaugs = 3;
1146 }
1147 else if ( mode == 42 )
1148 {
1149 _ndaugs = 3;
1153 }
1154 else if ( mode == 43 )
1155 {
1156 _ndaugs = 4;
1161 }
1162 else if ( mode == 44 )
1163 {
1164 _ndaugs = 2;
1167 }
1168 else if ( mode == 45 )
1169 {
1170 _ndaugs = 2;
1173 }
1174 else if ( mode == 46 )
1175 {
1176 _ndaugs = 2;
1179 }
1180 else if ( mode == 47 )
1181 {
1182 _ndaugs = 2;
1185 }
1186 else if ( mode == 48 )
1187 {
1188 _ndaugs = 2;
1191 }
1192 else if ( mode == 49 )
1193 {
1194 _ndaugs = 3;
1198 }
1199 else if ( mode == 50 )
1200 {
1201 _ndaugs = 3;
1205 }
1206 else if ( mode == 51 )
1207 {
1208 _ndaugs = 3;
1212 }
1213 else if ( mode == 52 )
1214 {
1215 _ndaugs = 2;
1218 }
1219 else if ( mode == 59 )
1220 {
1221 _ndaugs = 3;
1225 }
1226 else if ( mode == 60 )
1227 {
1228 _ndaugs = 3;
1232 }
1233 else if ( mode == 61 )
1234 {
1235 _ndaugs = 3;
1239 }
1240 else if ( mode == 62 )
1241 {
1242 _ndaugs = 3;
1246 }
1247 else if ( mode == 63 )
1248 {
1249 _ndaugs = 3;
1253 }
1254 else if ( mode == 64 )
1255 {
1256 _ndaugs = 3;
1260 }
1261 else if ( mode == 65 )
1262 {
1263 _ndaugs = 3;
1267 }
1268 else if ( mode == 66 )
1269 {
1270 _ndaugs = 3;
1274 }
1275 else if ( mode == 67 )
1276 {
1277 _ndaugs = 2;
1280 }
1281 else if ( mode == 68 )
1282 {
1283 _ndaugs = 2;
1286 }
1287 else if ( mode == 69 )
1288 {
1289 _ndaugs = 2;
1292 }
1293 else if ( mode == 70 )
1294 {
1295 _ndaugs = 2;
1298 }
1299 else if ( mode == 71 )
1300 {
1301 _ndaugs = 2;
1304 }
1305 else if ( mode == 72 )
1306 {
1307 _ndaugs = 2;
1310 }
1311 else if ( mode == 73 )
1312 {
1313 _ndaugs = 2;
1316 }
1317 else if ( mode == 74 )
1318 {
1319 _ndaugs = 2;
1322 }
1323 else if ( mode == 75 )
1324 {
1325 _ndaugs = 3;
1329 }
1330 else if ( mode == 76 )
1331 {
1332 _ndaugs = 3;
1336 }
1337 else if ( mode == 77 )
1338 {
1339 _ndaugs = 3;
1343 }
1344 else if ( mode == 78 )
1345 {
1346 _ndaugs = 3;
1350 }
1351 else if ( mode == 79 )
1352 {
1353 _ndaugs = 3;
1357 }
1358 else if ( mode == 80 )
1359 {
1360 _ndaugs = 2;
1363 }
1364 else if ( mode == 81 )
1365 {
1366 _ndaugs = 3;
1370 }
1371 else if ( mode == 82 )
1372 {
1373 _ndaugs = 3;
1377 }
1378 else if ( mode == 83 )
1379 {
1380 _ndaugs = 3;
1384 }
1385 else if ( mode == 84 )
1386 {
1387 _ndaugs = 3;
1391 }
1392 else if ( mode == 85 )
1393 {
1394 _ndaugs = 2;
1397 }
1398 else if ( mode == 90 )
1399 {
1400 _ndaugs = 3;
1404 }
1405 else if ( mode == 91 )
1406 {
1407 _ndaugs = 3;
1411 }
1412 else if ( mode == 92 )
1413 {
1414 _ndaugs = 3;
1418 }
1419 else if ( mode == 93 )
1420 {
1421 _ndaugs = 2;
1424 }
1425 else if ( mode == 94 )
1426 {
1427 _ndaugs = 2;
1430 }
1431 else if ( mode == 95 )
1432 {
1433 _ndaugs = 2;
1436 }
1437 else if ( mode == 96 )
1438 {
1439 _ndaugs = 2;
1442 }
1443 else if ( mode == 97 )
1444 {
1445 _ndaugs = 3;
1449 }
1450 else if ( mode == 111 )
1451 {
1452 _ndaugs = 2;
1455 }
1456 else if ( mode == 112 )
1457 {
1458 _ndaugs = 2;
1461 }
1462 else if ( mode == 74100 )
1463 {
1464 _ndaugs = 1;
1466 }
1467 else if ( mode == 74101 )
1468 {
1469 _ndaugs = 1;
1471 }
1472 else if ( mode == 74102 )
1473 {
1474 _ndaugs = 1;
1476 }
1477 else if ( mode == 74103 )
1478 {
1479 _ndaugs = 1;
1481 }
1482 else if ( mode == 74104 )
1483 {
1484 _ndaugs = 1;
1486 }
1487 else if ( mode == 74105 )
1488 {
1489 _ndaugs = 1;
1491 }
1492 else if ( mode == 74106 )
1493 {
1494 _ndaugs = 1;
1496 }
1497 else if ( mode == 74107 )
1498 {
1499 _ndaugs = 1;
1501 }
1502 else if ( mode == 74110 )
1503 {
1504 _ndaugs = 1;
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]; }
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;
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;
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;
1552 }
1553
1554 double fmth = 0;
1555 for ( int i = 0; i < _ndaugs; i++ )
1556 {
1558 fmth += xmi;
1559 }
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}
◆ InitPars()
| void EvtConExc::InitPars | ( | ) |
Definition at line 4470 of file EvtConExc.cc.
4470 {
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 );
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}
std::vector< std::string > split(std::string str, std::string pattern)
Definition EvtConExc.cc:4451
Referenced by init().
◆ initProbMax()
|
virtual |
◆ islgr()
| bool EvtConExc::islgr | ( | double * | x, |
| double * | y, | ||
| int | n, | ||
| double | t ) |
Definition at line 3508 of file EvtConExc.cc.
3508 {
3509 int n0 = -1;
3510 double z;
3512 {
3514 {
3515 n0 = i;
3516 break;
3517 }
3518 }
3519 double xstotal = y[599];
3520 if ( n0 == -1 ) { return 0; }
3521 else
3522 {
3527 else { return 0; }
3528 }
3529}
◆ ISR_ang_integrate()
| double EvtConExc::ISR_ang_integrate | ( | double | x, |
| double | theta ) |
Definition at line 3591 of file EvtConExc.cc.
3591 {
3592 // see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
3595 double cos2 = costheta * costheta;
3596 double sin2 = 1 - cos2;
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 =
3606 double hq = ( L - 1 ) / 2. + hq1 + hq2 + hq3;
3607 hq /= hpi; // Eq (11) in arXif:hep-ph/9910523
3608 return hq;
3609}
Referenced by ISR_ang_sampling().
◆ ISR_ang_sampling()
| double EvtConExc::ISR_ang_sampling | ( | double | x | ) |
Definition at line 3611 of file EvtConExc.cc.
3611 {
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;
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}
double ISR_ang_integrate(double x, double theta)
Definition EvtConExc.cc:3591
Referenced by decay().
◆ lgr()
| double EvtConExc::lgr | ( | double * | x, |
| double * | y, | ||
| int | n, | ||
| double | t ) |
Definition at line 3486 of file EvtConExc.cc.
3486 {
3487 int n0 = -1;
3488 double z;
3490 {
3492 {
3493 n0 = i;
3494 break;
3495 }
3496 }
3497 if ( n0 == -1 ) { return 0.0; }
3498 else
3499 {
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}
◆ Li2()
| double EvtConExc::Li2 | ( | double | x | ) |
Definition at line 3480 of file EvtConExc.cc.
3480 {
3481 // double li2= -x +x*x/4. - x*x*x/9; // wangwp: may be not correct!
3483 return li2;
3484}
Referenced by SoftPhoton_xs().
◆ LLr()
| double EvtConExc::LLr | ( | double * | x, |
| double * | y, | ||
| int | n, | ||
| double | t ) |
Definition at line 3531 of file EvtConExc.cc.
3531 {
3532 int n0 = -1;
3533 double z;
3536 {
3538 {
3539 n0 = i;
3540 break;
3541 }
3542 }
3543 if ( n0 == -1 ) { return 0.0; }
3544 else
3545 {
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}
◆ meson_sampling()
| bool EvtConExc::meson_sampling | ( | EvtVector4R | pcm, |
| EvtVector4R | pi ) |
Definition at line 3203 of file EvtConExc.cc.
3203 {
3205 double theta = angles.getHelAng( 1 );
3206 double phi = angles.getHelAng( 2 );
3208
3210 double ang = 1 - costheta * costheta;
3211 if ( pm < ang / 1. ) { return true; }
3212 else { return false; }
3213}
Referenced by hadron_angle_sampling().
◆ Mhad_sampling()
| double EvtConExc::Mhad_sampling | ( | double * | x, |
| double * | y ) |
Definition at line 3554 of file EvtConExc.cc.
3554 { // 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
3560 double xrt = 1;
3561 int mybin = 1;
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; }
3574
3576 {
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++ )
3586 abort();
3587 }
3588 return mhd;
3589}
Referenced by decay().
◆ mk_VXS()
| void EvtConExc::mk_VXS | ( | double | Esig, |
| double | Egamcut, | ||
| double | EgamH, | ||
| int | midx ) |
Definition at line 4223 of file EvtConExc.cc.
4223 {
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
4232 else { uscale = 1; }
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);
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);
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}
Referenced by init().
◆ narrowRXS()
| double EvtConExc::narrowRXS | ( | double | mxL, |
| double | mxH ) |
Definition at line 4088 of file EvtConExc.cc.
4088 {
4089 // br for ee
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;
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
4123 double unic = 3.89379304 * 100000; // in unit nb
4124 return sigv * unic; // vaccum factor has included in the Rad2
4125}
***************************************************************************************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
***************************************************************************************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
◆ OutStaISR()
| void EvtConExc::OutStaISR | ( | ) |
Definition at line 4562 of file EvtConExc.cc.
4562 {
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}
Referenced by ~EvtConExc().
◆ photonSampling()
| bool EvtConExc::photonSampling | ( | EvtParticle * | part | ) |
Definition at line 4055 of file EvtConExc.cc.
4055 {
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();
4066
4067 if ( pdgcode == 111 || pdgcode == 221 || pdgcode == 331 )
4068 { // for photon
4069 alpha = 1;
4070 }
4071 else continue;
4072
4076 double xratio = ang / angmax;
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}
Referenced by decay().
◆ PrintXS()
| void EvtConExc::PrintXS | ( | double | mx | ) |
Definition at line 3949 of file EvtConExc.cc.
3949 { // 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
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];
3986 init_mode( mymode );
3988 else { uscale = 1; }
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++ )
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}
Referenced by init().
◆ Rad1()
| double EvtConExc::Rad1 | ( | double | s, |
| double | x ) |
Definition at line 3227 of file EvtConExc.cc.
3227 { // 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
3234 double L = 2 * log( sme );
3236 return wsx;
3237}
Referenced by Rad1difXs().
◆ Rad1difXs()
| double EvtConExc::Rad1difXs | ( | EvtParticle * | p | ) |
Definition at line 3318 of file EvtConExc.cc.
3318 { // // first order xsection
3319 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
3322
3324 double mrg = cms - summass;
3326 double mhs = pm * mrg + summass;
3327
3331
3332 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<",
3333 // "<<sin2theta<<", "<<cms<<",
3334 // "<<x<<", "<<wsx<<std::endl;
3335
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}
Referenced by gamHXSection().
◆ Rad2()
| double EvtConExc::Rad2 | ( | double | s, |
| double | x ) |
Definition at line 3239 of file EvtConExc.cc.
3239 {
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
3245 double xi2 = 1.64493407;
3246 double xi3 = 1.2020569;
3248 double L = 2 * log( sme );
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.;
3257 wsx = wsx + beta * beta / 8. * wsx2;
3259 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
3261}
Referenced by narrowRXS(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), and Ros_xs().
◆ Rad2difXs() [1/3]
| double EvtConExc::Rad2difXs | ( | double | s, |
| double | x ) |
Definition at line 3295 of file EvtConExc.cc.
3295 { // leading second order xsection
3299 // if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
3301 double difxs = wsx * xs;
3302 differ2 = wsx * ( myxsection->getErr( mhs ) );
3303 return difxs;
3304}
◆ Rad2difXs() [2/3]
| double EvtConExc::Rad2difXs | ( | double | s, |
| double | x, | ||
| EvtXsection * | myxsection ) |
Definition at line 3306 of file EvtConExc.cc.
3307 { // leading second order xsection
3311 // if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
3313 double difxs = wsx * xs;
3314 differ2 = wsx * ( myxsection->getErr( mhs ) );
3315 return difxs;
3316}
◆ Rad2difXs() [3/3]
| double EvtConExc::Rad2difXs | ( | EvtParticle * | p | ) |
Definition at line 3263 of file EvtConExc.cc.
3263 { // leading second order xsection
3264 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
3267 for ( int i = 1; i < _ndaugs; i++ )
3268 {
3271 }
3272
3275 double mrg = cms - summass;
3277 double mhs = pm * mrg + summass;
3278
3281 // double mhs = sqrt( s*(1-x) );
3283
3284 // std::cout<<"random : "<<pm<<std::endl;
3285
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}
Referenced by gamHXSection(), gamHXSection(), gamHXSection(), gamHXSection(), trapezoid(), and trapezoid().
◆ ReadVP()
| void EvtConExc::ReadVP | ( | ) |
Definition at line 4204 of file EvtConExc.cc.
4204 {
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}
std::ifstream ifstream
Referenced by init().
◆ resetResMass()
| void EvtConExc::resetResMass | ( | ) |
Definition at line 3807 of file EvtConExc.cc.
3807 {
3808 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
3810 EvtPDL::reSetMass( myres, mjsi );
3811 // EvtPDL::reSetWidth(myres,wjsi);
3812
3814 EvtPDL::reSetMass( myres, mpsip );
3815 // EvtPDL::reSetWidth(myres,wpsip);
3816
3818 EvtPDL::reSetMass( myres, mpsipp );
3819 // EvtPDL::reSetWidth(myres,wpsipp);
3820
3822 EvtPDL::reSetMass( myres, mphi );
3823 // EvtPDL::reSetWidth(myres,wphi);
3824
3826 EvtPDL::reSetMass( myres, momega );
3827 // EvtPDL::reSetWidth(myres,womega);
3828
3830 EvtPDL::reSetMass( myres, mrho0 );
3831 // EvtPDL::reSetWidth(myres,wrho0);
3832
3834 EvtPDL::reSetMass( myres, mrho3s );
3835 // EvtPDL::reSetWidth(myres,wrho3s);
3836
3838 EvtPDL::reSetMass( myres, momega2s );
3839 // EvtPDL::reSetWidth(myres,womega2s);
3840}
Referenced by decay().
◆ Ros_xs()
| double EvtConExc::Ros_xs | ( | double | mx, |
| double | bree, | ||
| EvtId | pid ) |
Definition at line 3377 of file EvtConExc.cc.
3377 { // in unit nb
3385 double nbar = 3.89379304 * 100000; // ! GeV-2 = 3.89379304*10^5 nbar
3386 return sigma * nbar;
3387}
Referenced by init().
◆ selectMass()
| double EvtConExc::selectMass | ( | ) |
Definition at line 4141 of file EvtConExc.cc.
4141 {
4142 double pm, mhdL, mhds;
4143 pm = EvtRandom::Flat( 0., 1 );
4144 mhdL = _mhdL;
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 {
4160 }
4161 for ( int j = 1; j < sxs.size(); j++ )
4162 {
4163 double x0 = sxs[j - 1];
4164 double x1 = sxs[j];
4166 }
4167 return mhds;
4168}
◆ selectMode()
| int EvtConExc::selectMode | ( | std::vector< int > | vmod, |
| double | mhds ) |
Definition at line 3716 of file EvtConExc.cc.
3716 {
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];
3738 else { uscale = 1; }
3739 // if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
3740
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:
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]; }
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]; }
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}
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition EvtDecayTable.cc:831
Referenced by decay().
◆ SetP4()
| void EvtConExc::SetP4 | ( | EvtParticle * | part, |
| double | mhdr, | ||
| double | xeng, | ||
| double | theta ) |
Definition at line 3639 of file EvtConExc.cc.
3640 { // set the gamma and gamma* momentum according sampled
3641 // results
3643 double phi = 2 * 3.1415926 * pm;
3645 double hdrE = sqrt( mhdr * mhdr + gamE * gamE );
3649 EvtVector4R p4[2];
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] );
3658 {
3660 {
3661 EvtVector4R PG( _cms, 0, 0, 0 );
3662 p->getDaug( i )->init( myid, PG );
3663 }
3664 }
3665 }
3666}
double sin(const BesAngle a)
Referenced by decay().
◆ SetP4Rvalue()
| void EvtConExc::SetP4Rvalue | ( | EvtParticle * | part, |
| double | mhdr, | ||
| double | xeng, | ||
| double | theta ) |
Definition at line 3668 of file EvtConExc.cc.
3669 { // set the gamma and gamma* momentum according
3670 // sampled results
3672 double phi = 2 * 3.1415926 * pm;
3674 double hdrE = sqrt( mhdr * mhdr + gamE * gamE );
3678 EvtVector4R p4[2];
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}
Referenced by decay().
◆ showResMass()
| void EvtConExc::showResMass | ( | ) |
Definition at line 3876 of file EvtConExc.cc.
3876 {
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}
◆ SoftPhoton_xs()
| double EvtConExc::SoftPhoton_xs | ( | double | s, |
| double | b ) |
Definition at line 3450 of file EvtConExc.cc.
3450 {
3454 double xi2 = 1.64493407;
3455 double xi3 = 1.2020569;
3457 double L = 2 * log( sme );
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.;
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.;
3476 // return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
3478}
◆ split()
| std::vector< std::string > EvtConExc::split | ( | std::string | str, |
| std::string | pattern ) |
Definition at line 4451 of file EvtConExc.cc.
4451 {
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}
Referenced by InitPars().
◆ sumExc()
| double EvtConExc::sumExc | ( | double | mx | ) |
Definition at line 3903 of file EvtConExc.cc.
3903 { // 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
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;
3940 init_mode( mymode );
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}
Referenced by init().
◆ trapezoid() [1/2]
| double EvtConExc::trapezoid | ( | double | s, |
| double | a, | ||
| double | b, | ||
| int | n ) |
Definition at line 4576 of file EvtConExc.cc.
4578{
4581 double sum = 0.0;
4584 {
4587 }
4588 return sum;
4589}
◆ trapezoid() [2/2]
| double EvtConExc::trapezoid | ( | double | s, |
| double | El, | ||
| double | Eh, | ||
| int | n, | ||
| EvtXsection * | myxs ) |
Definition at line 4591 of file EvtConExc.cc.
4593{
4596 double sum = 0.0;
4599 {
4602 }
4603 return sum;
4604}
◆ VP_sampling()
| bool EvtConExc::VP_sampling | ( | EvtVector4R | pcm, |
| EvtVector4R | pi ) |
Definition at line 3215 of file EvtConExc.cc.
3215 {
3217 double theta = angles.getHelAng( 1 );
3218 double phi = angles.getHelAng( 2 );
3220
3222 double ang = 1 + costheta * costheta;
3223 if ( pm < ang / 2. ) { return true; }
3224 else { return false; }
3225}
Referenced by hadron_angle_sampling().
◆ writeDecayTabel()
| void EvtConExc::writeDecayTabel | ( | ) |
Definition at line 4296 of file EvtConExc.cc.
4296 {
4299 ob.open( "obsxs.dat" );
4300 //
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 {
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
4376 for ( int i = 0; i < Vmode.size(); i++ )
4380}
Referenced by init().
◆ xs_sampling() [1/2]
| bool EvtConExc::xs_sampling | ( | double | xs | ) |
Definition at line 3170 of file EvtConExc.cc.
3170 {
3172 // std::cout<<"Rdm= "<<pm<<std::endl;
3174 else { return false; }
3175}
◆ xs_sampling() [2/2]
| bool EvtConExc::xs_sampling | ( | double | xs, |
| double | xs1 ) |
Definition at line 3177 of file EvtConExc.cc.
3177 {
3179 // std::cout<<"Rdm= "<<pm<<std::endl;
3180 if ( pm < xs / ( xs0 * 1.1 ) ) { return true; }
3181 else { return false; }
3182}
Member Data Documentation
◆ _cms
|
static |
Definition at line 192 of file EvtConExc.hh.
Referenced by baryon_sampling(), calAF(), checkEvtRatio(), decay(), difgamXs(), Egam2Mhds(), findMaxXS(), getObsXsection(), init(), ISR_ang_integrate(), Mhad_sampling(), mk_VXS(), narrowRXS(), PrintXS(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMass(), SetP4(), SetP4Rvalue(), sumExc(), and writeDecayTabel().
◆ conexcmode
|
static |
Definition at line 214 of file EvtConExc.hh.
Referenced by decay(), and EvtRexc::decay().
◆ getMode
|
static |
Definition at line 196 of file EvtConExc.hh.
Referenced by init().
◆ mlow
|
static |
Definition at line 239 of file EvtConExc.hh.
Referenced by EvtDecay::energySpread(), init(), and mk_VXS().
◆ multimode
|
static |
Definition at line 214 of file EvtConExc.hh.
Referenced by EvtRexc::decay(), and init().
◆ mup
|
static |
Definition at line 239 of file EvtConExc.hh.
Referenced by EvtDecay::energySpread(), and init().
◆ myxsection
|
static |
Definition at line 191 of file EvtConExc.hh.
Referenced by decay(), difgamXs(), difgamXs(), findMaxXS(), init(), init_mode(), mk_VXS(), PrintXS(), Rad1difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMode(), sumExc(), and ~EvtConExc().
◆ SetMthr
|
static |
Definition at line 194 of file EvtConExc.hh.
Referenced by decay(), init(), and EvtDecay::initialize().
◆ staxsection
|
static |
Definition at line 191 of file EvtConExc.hh.
Referenced by calAF(), energySpread(), init(), trapezoid(), and ~EvtConExc().
◆ XS_max
|
static |
Definition at line 193 of file EvtConExc.hh.
Referenced by init(), Rad1difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), and xs_sampling().
The documentation for this class was generated from the following files:
