BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Static Public Attributes | List of all members
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 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
virtual void getName (std::string &name)=0
 
virtual void decay (EvtParticle *p)=0
 
virtual void makeDecay (EvtParticle *p)=0
 
virtual EvtDecayBase * clone ()=0
 
virtual void init ()
 
virtual void initProbMax ()
 
virtual std::string commandName ()
 
virtual void command (std::string cmd)
 
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 88 of file EvtConExc.hh.

Constructor & Destructor Documentation

◆ EvtConExc()

EvtConExc::EvtConExc ( )
inline

Definition at line 92 of file EvtConExc.hh.

92 {
93 } //constructor

Referenced by clone().

◆ ~EvtConExc()

EvtConExc::~EvtConExc ( )
virtual

Definition at line 183 of file EvtConExc.cc.

183 {
184 if(_mode==74110)checkEvtRatio();
185 if(myFisr.size())OutStaISR();
186 if(mydbg){
187 // xs->Write();
188 myfile->Write();
189 }
190 delete myxsection;
191 delete staxsection;
192 gamH->deleteTree();
193
194 //the deletion of commands is really uggly!
195 int i;
196 if (nconexcdecays==0) {
197 delete [] commands;
198 commands=0;
199 return;
200 }
201
202 for(i=0;i<nconexcdecays;i++){
203 if (conexcdecays[i]==this){
204 conexcdecays[i]=conexcdecays[nconexcdecays-1];
205 nconexcdecays--;
206 if (nconexcdecays==0) {
207 delete [] commands;
208 commands=0;
209 }
210 return;
211 }
212 }
213
214
215}
EvtConExc::OutStaISR
void OutStaISR()
Definition: EvtConExc.cc:3698
EvtConExc::staxsection
static EvtXsection * staxsection
Definition: EvtConExc.hh:139
EvtConExc::myxsection
static EvtXsection * myxsection
Definition: EvtConExc.hh:139
EvtConExc::checkEvtRatio
void checkEvtRatio()
Definition: EvtConExc.cc:3513
EvtParticle::deleteTree
void deleteTree()
Definition: EvtParticle.cc:557

Member Function Documentation

◆ addNarrowRXS()

double EvtConExc::addNarrowRXS ( double  mhi,
double  binwidth 
)

Definition at line 3284 of file EvtConExc.cc.

3284 {
3285 double myxs = 0;
3286 for(int i=0;i<ISRXS.size();i++){ // std::cout<<"ISRM "<<ISRM[i]<<std::endl;
3287 if( (ISRM[i]-mhi)<binwidth && ISRFLAG[i]){ISRFLAG[i]=0; myxs = ISRXS[i];}
3288 }
3289 //std::cout<<mhi<<" "<<binwidth<<" "<<myxs<<std::endl;
3290 return myxs;
3291}

◆ angularSampling()

bool EvtConExc::angularSampling ( EvtParticle *  part)

Definition at line 3173 of file EvtConExc.cc.

3173 {
3174 bool tagp,tagk;
3175 tagk=0;
3176 tagp=0;
3177 int nds = par->getNDaug();
3178 for(int i=0;i<par->getNDaug();i++){
3179 EvtId di=par->getDaug(i)->getId();
3180 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
3181 int pdgcode =EvtPDL::getStdHep( di );
3182 double alpha=1;
3183 /*
3184 if(fabs(pdgcode)==2212){//proton and anti-proton
3185 alpha = baryonAng(p4i.mass());
3186 cosp = cos(p4i.theta());
3187 tagp=1;
3188 }else if(fabs(pdgcode)==321||fabs(pdgcode)==211 ){
3189 alpha=1;
3190 cosk = cos(p4i.theta());
3191 tagk=1;
3192 }else continue;
3193 */
3194 double angmax = 1+alpha;
3195 double costheta = cos(p4i.theta());
3196 double ang=1+alpha*costheta*costheta;
3197 double xratio = ang/angmax;
3198 double xi=EvtRandom::Flat(0.,1);
3199 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
3200 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
3201 if(xi>xratio) return false;
3202 }//loop over duaghters
3203 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
3204 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
3205
3206 return true;
3207}
alpha
const double alpha
Definition: FastVertexFit.cxx:4
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
EvtId
Definition: EvtId.hh:27
EvtId::getId
int getId() const
Definition: EvtId.hh:41
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:73
EvtVector4R::theta
double theta()
Definition: EvtVector4R.cc:249

Referenced by decay().

◆ baryon_sampling()

bool EvtConExc::baryon_sampling ( EvtVector4R  pcm,
EvtVector4R  pi 
)

Definition at line 2488 of file EvtConExc.cc.

2488 {
2489 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2490 double theta=angles.getHelAng(1);
2491 double phi =angles.getHelAng(2);
2492 double costheta=cos(theta); //using helicity angles in parent system
2493
2494 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
2495 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
2496 double alpha = baryonAng(_cms);
2497 if(_mode ==96){alpha=-0.34;}
2498 double pm= EvtRandom::Flat(0.,1);
2499 double ang = 1 + alpha*costheta*costheta;
2500 double myang;
2501 if(alpha>=0){myang=1+alpha;}else{myang=1;}
2502 if(pm< ang/myang) {return true;}else{return false;}
2503}
EvtConExc::baryonAng
double baryonAng(double mx)
Definition: EvtConExc.cc:3209
EvtConExc::_cms
static double _cms
Definition: EvtConExc.hh:140

Referenced by hadron_angle_sampling().

◆ baryonAng()

double EvtConExc::baryonAng ( double  mx)

Definition at line 3209 of file EvtConExc.cc.

3209 {
3210 double mp=0.938;
3211 double u = 0.938/mx;
3212 u = u*u;
3213 double u2 = (1+u)*(1+u);
3214 double uu = u*(1+6*u);
3215 double alpha = (u2-uu)/(u2+uu);
3216 return alpha;
3217}
mp
const double mp
Definition: incllambda.cxx:45

Referenced by baryon_sampling().

◆ calAF()

void EvtConExc::calAF ( double  myecms)

Definition at line 3647 of file EvtConExc.cc.

3647 {
3648
3649
3650 double _cms=myecms;
3651 double Esig = 0.0001; //sigularity engy
3652 double x = 2*Egamcut/_cms;
3653 double s = _cms*_cms;
3654
3655 double mhdL=staxsection->getXlw();
3656 double EgamH = (s-mhdL*mhdL)/2/sqrt(s);
3657
3658 int nint=50;
3659 int nmc= 1000;
3660 _xs0 = trapezoid(s,Esig,Egamcut,nint); // std::cout<<"Int: _xs0= "<<_xs0<<std::endl;
3661 //--- sigularity point x from 0 to 0.0001GeV
3662 double xb= 2*Esig/_cms;
3663 double sig_xs = SoftPhoton_xs(s,xb)*(staxsection->getXsection(_cms));
3664 _xs0 += sig_xs;
3665
3666 int Nstp=598;
3667 double stp=(EgamH - Egamcut)/Nstp;
3668 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
3669 double Eg0=EgamH - i*stp;
3670 double Eg1=EgamH - (i+1)*stp;
3671 double xsi= trapezoid(s,Eg1,Eg0,nint); // std::cout<<"Int xsi= " <<xsi<<std::endl;
3672
3673 AA[i]=(Eg1+Eg0)/2;
3674 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
3675 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
3676 double binwidth = mh2-mhi;
3677 if(i==0) AF[0]=xsi;
3678 if(i>=1) AF[i]=AF[i-1]+xsi;
3679 RadXS[i]=xsi/stp;
3680 }
3681 //add the interval 0~Esig
3682 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
3683 AF[598]=AF[597];
3684 AF[599]=AF[598]+ _xs0;
3685 //--
3686 for(int i=0;i<600;i++){
3687 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
3688 }
3689 //--debugging
3690 double bornXS = staxsection->getXsection(_cms);
3691 if(bornXS==0){std::cout<<"EvtConExc::calAF: bad BornXS at "<<_cms<<" GeV"<<std::endl;abort();}
3692 double fisr=AF[599]/bornXS;
3693 myFisr.push_back(fisr);
3694 //std::cout<<"fisr= "<<fisr<<std::endl;
3695}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
s
XmlRpcServer s
Definition: HelloServer.cpp:11
EvtConExc::SoftPhoton_xs
double SoftPhoton_xs(double s, double b)
Definition: EvtConExc.cc:2746
EvtConExc::trapezoid
double trapezoid(double s, double a, double b, int n)
Definition: EvtConExc.cc:3712
EvtXsection::getXlw
double getXlw()
Definition: EvtXsection.hh:156
EvtXsection::getXsection
double getXsection(double mx)
Definition: EvtXsection.cc:12008

Referenced by decay().

◆ checkdecay()

bool EvtConExc::checkdecay ( EvtParticle *  p)

Definition at line 3123 of file EvtConExc.cc.

3123 {
3124 int nson=p->getNDaug();
3125 double massflag=1;
3126 for(int i=0;i<nson;i++){
3127 if( EvtPDL::getStdHep(p->getDaug(i)->getId())==22 ) continue;
3128 massflag *= p->getDaug(i)->mass();
3129 }
3130 if(massflag==0){
3131 std::cout<<"Zero mass!"<<std::endl;
3132 return 0;
3133 }else{return 1;}
3134}
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle::getNDaug
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:127

Referenced by decay().

◆ checkEvtRatio()

void EvtConExc::checkEvtRatio ( )

Definition at line 3513 of file EvtConExc.cc.

3513 {
3514 ofstream oa;
3515 oa.open("_evtR.dat");
3516 //
3517 int im=getModeIndex(74110);
3518 double xscon=VXS[im][599];
3519 double xscon0=xscon;
3520 oa<<"Ecms= "<<_cms<<" GeV"<<std::endl;
3521 oa<<"The total observed xs: "<<xscon<<" nb"<<std::endl;
3522 oa<<"=== mode =========== ratio ==========="<<std::endl;
3523 for(int i=0;i<vmode.size();i++){
3524 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
3525 int themode=getModeIndex(vmode[i]);
3526 if(vmode[i]==74110) continue;
3527 xscon -= VXS[themode ][599];
3528 }
3529 if(xscon<0) xscon=0;
3530 //debugging
3531 for(int i=0;i<vmode.size();i++){
3532 int themode=getModeIndex(vmode[i]);
3533 if(vmode[i]==74110) continue;
3534 oa<<vmode[i]<<"-th mode: "<<100*VXS[themode][599]/xscon0<<" % "<<std::endl;
3535 }
3536 oa<<"74110-th mode: "<<100*xscon/xscon0<<" % "<<std::endl;
3537
3538
3539}
oa
std::ofstream oa
Definition: EvtConExc.cc:177
EvtConExc::getModeIndex
int getModeIndex(int m)
Definition: EvtConExc.cc:3541

Referenced by ~EvtConExc().

◆ clone()

EvtDecayBase * EvtConExc::clone ( )
virtual

Implements EvtDecayBase.

Definition at line 223 of file EvtConExc.cc.

223 {
224
225 return new EvtConExc;
226
227}

◆ command()

void EvtConExc::command ( std::string  cmd)
virtual

Reimplemented from EvtDecayBase.

Definition at line 3555 of file EvtConExc.cc.

3555 {
3556 if (ncommand==lcommand){
3557 lcommand=10+2*lcommand;
3558 std::string* newcommands=new std::string[lcommand];
3559 int i;
3560 for(i=0;i<ncommand;i++){
3561 newcommands[i]=commands[i];
3562 }
3563 delete [] commands;
3564 commands=newcommands;
3565 }
3566 commands[ncommand]=cmd;
3567 ncommand++;
3568}

◆ commandName()

std::string EvtConExc::commandName ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 3549 of file EvtConExc.cc.

3549 {
3550
3551 return std::string("ConExcPar");
3552
3553}

◆ decay()

void EvtConExc::decay ( EvtParticle *  p)
virtual

Implements EvtDecayBase.

Definition at line 1722 of file EvtConExc.cc.

1722 {
1723 //std::cout<<"_cms= "<<_cms<<" mode= "<<_mode<<std::endl;
1724 EvtId myvpho = EvtPDL::getId(std::string("vpho"));
1725 if(myvpho != p->getId()){
1726 std::cout<<"Parent needs to be vpho, but found "<<EvtPDL::name(p->getId())<<std::endl;
1727 abort();
1728 }
1729 //-- to read _cms from database or to consider beam energy spread
1730 double mvpho=EvtPDL::getMeanMass(myvpho);
1731 bool beamflag=0;
1732 if(mvpho!=EvtConExc::_cms){ //if read cms energy from database
1733 EvtConExc::_cms=mvpho;
1734 beamflag=1;
1735 }
1736
1737 /*
1738 if(beamEnergySpread!=0){
1739 cmsspread=sqrt(2.)*beamEnergySpread;
1740 _cms = energySpread(_cms,cmsspread);
1741 beamflag=1;
1742 std::cout<<"test_cms "<<_cms<<std::endl;
1743 }
1744 */
1745
1746 if(beamflag) calAF(_cms);
1747
1748
1749 //debugging
1750 //double mvpho=EvtPDL::getMeanMass(myvpho);
1751 //std::cout<<"myvpho= "<<mvpho<<std::endl;
1752
1753 //for fill root tree
1754 EvtVector4R vgam,hd1,hd2,vhds;
1755
1756 //first for Rscan generator with _mode=74110
1757 if(_mode==74110){
1758 //preparation of mode sampling
1759 std::vector<int> vmod; vmod.clear();
1760 int mn[82]={0,1,2,3,4,5,6,7,8,9,10,11,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46, // 12,14, 30,31,39,42 and 43 is removed
1761 50,51,
1762 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1763 90,91,92,93,94,95,96,
1764 74100,74101,74102,74103,74104,74105,74110};
1765 int mn2[83]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46,//mode 14,30,31,42, 43 are removed
1766 50,51,
1767 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1768 90,91,92,93,94,95,96,
1769 74100,74101,74102,74103,74104,74105,74110};
1770 int mn3[76]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,44,45,46,
1771 50,51,
1772 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1773 90,91,92,93,94,95,96,
1774 74100,74101,74102,74110}; //remove 43, 74103,74104,74105, this is included in PHOKHARA
1775 double mx = p->getP4().mass();
1776
1777 if(_cms>2.5){
1778 for(int i=0;i<82;i++){vmod.push_back(mn[i]);}
1779 }else{
1780 for(int i=0;i<83;i++){vmod.push_back(mn2[i]);}
1781 }
1782
1783 //for(int i=0;i<76;i++){vmod.push_back(mn3[i]);}
1784
1785 int mymode;
1786 double pm= EvtRandom::Flat(0.,1);
1787
1788 //std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
1789 //--
1790 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross section within Egam=(0.025-Egamma^max GeV)
1791 if(pm <_xs0/(_xs0 + _xs1) ){//without ISR photon for mode=74110
1792
1793 mhds=_cms;
1794 mymode=selectMode(vmod,mhds);
1795 _selectedMode = mymode;
1796 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1797 delete myxsection;
1798 myxsection =new EvtXsection (mymode);
1799 init_mode(mymode);
1800 resetResMass();
1801
1802 if(_ndaugs>1){//for e+e- -> exclusive decays
1803 checkA:
1804 p->makeDaughters(_ndaugs,daugs);
1805 double totMass=0;
1806 for(int i=0;i< _ndaugs;i++){
1807 EvtParticle* di=p->getDaug(i);
1808 double xmi=EvtPDL::getMass(di->getId());
1809 di->setMass(xmi);
1810 totMass+=xmi;
1811 }
1812 if(totMass>p->mass()) goto checkA;
1813 p->initializePhaseSpace( _ndaugs , daugs);
1814 if(!checkdecay(p)) goto checkA;
1815 vhds = p->getP4();
1816 if(_cms>2.5 && !angularSampling(p)) goto checkA;
1817 if(_cms<2.5 && !photonSampling(p)) goto checkA;
1818 }else{// for e+e- -> vector resonace, add a very soft photon
1819 EvtId mydaugs[2];
1820 mydaugs[0]=daugs[0];
1821 EvtPDL::reSetMass(mydaugs[0],mhds-0.001);
1822 //EvtPDL::reSetWidth(mydaugs[0],0);
1823 mydaugs[1]=gamId; //ISR photon
1824 p->makeDaughters(2,mydaugs);
1825 checkB:
1826 double totMass=0;
1827 for(int i=0;i< 2;i++){
1828 EvtParticle* di=p->getDaug(i);
1829 double xmi=EvtPDL::getMass(di->getId());
1830 di->setMass(xmi);
1831 totMass+=xmi;
1832 }
1833 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
1834 if(totMass>p->mass()) goto checkB;
1835 p->initializePhaseSpace(2,mydaugs);
1836
1837 if(!checkdecay(p)) goto checkB;
1838 vhds = p->getDaug(0)->getP4();;
1839 vgam = p->getDaug(1)->getP4();
1840 }
1841 //-----
1842 }else{// with ISR photon for mode=74110
1843Sampling_mhds:
1844 mhds=Mhad_sampling(MH,AF);
1845 //std::cout<<"SetMthr= "<<SetMthr<<std::endl;
1846 if(mhds<SetMthr) goto Sampling_mhds;
1847 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
1848 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1849
1850 if(mydbg) mass2=mhds;
1851
1852 //generating events
1853 ModeSelection:
1854 EvtGlobalSet::ConExcPythia=1;
1855 mymode=selectMode(vmod,mhds);
1856 if(mymode==-10) goto Sampling_mhds;
1857 conexcmode = mymode; //save for model REXC ( see EvtRexc.cc) decay
1858 if(mhds<2.3 && mymode==74110) {goto ModeSelection;} // E<2.3 GeV, fully exclusive production
1859 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1860 _selectedMode = mymode;
1861 delete myxsection;
1862 myxsection =new EvtXsection (mymode);
1863 init_mode(mymode);
1864 EvtId myvpho = EvtPDL::getId(std::string("vgam"));
1865 EvtPDL::reSetMass(myvpho,mhds); //set to continum cms energy
1866 EvtPDL::reSetWidth(myvpho,0);
1867
1868 //debugging
1869 //for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<" "<<EvtPDL::name(daugs[i])<<std::endl;}
1870
1871 //decay e+e- ->gamma_ISR gamma*
1872 EvtId mydaugs[2];
1873 if(_ndaugs>1) { //gamma* -> exclusive decays
1874 resetResMass();
1875 mydaugs[0]=myvpho;
1876 }else{// vector resonance decays
1877 resetResMass();
1878 mydaugs[0]=daugs[0];
1879 EvtPDL::reSetMass(mydaugs[0],mhds);
1880 //EvtPDL::reSetWidth(mydaugs[0],0);
1881 }
1882 mydaugs[1]=gamId; //ISR photon
1883 int maxflag=0;
1884 int trycount=0;
1885 ISRphotonAng_sampling:
1886 double totMass=0;
1887 p->makeDaughters(2,mydaugs);
1888 for(int i=0;i< 2;i++){
1889 EvtParticle* di=p->getDaug(i);
1890 double xmi=EvtPDL::getMass(di->getId());
1891 di->setMass(xmi);
1892 totMass+=xmi;
1893 }
1894 if(totMass>p->mass()) goto ISRphotonAng_sampling;
1895 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
1896 double weight1 = p->initializePhaseSpace(2,mydaugs);
1897 if(!checkdecay(p)) goto ISRphotonAng_sampling;
1898 //set the ISR photon angular distribution
1899 SetP4Rvalue(p,mhds,xeng,theta); //end of decay e+e- -> vpho gamma_ISR
1900
1901 if(maxflag==0) {findMaxXS(mhds); maxflag=1;}
1902 vhds = p->getDaug(0)->getP4();
1903 vgam=p->getDaug(1)->getP4();
1904 double gx=vgam.get(1);
1905 double gy=vgam.get(2);
1906 double sintheta= sqrt(gx*gx+gy*gy)/vgam.d3mag();
1907 bool gam_ang=gam_sampling(mhds,sintheta);
1908 trycount++;
1909
1910 } // with and without ISR sampling block
1911 // finish event generation
1912 // for debugging
1913 if(mydbg){
1914 pgam[0]=vgam.get(0);
1915 pgam[1]=vgam.get(1);
1916 pgam[2]=vgam.get(2);
1917 pgam[3]=vgam.get(3);
1918
1919 phds[0]=vhds.get(0);
1920 phds[1]=vhds.get(1);
1921 phds[2]=vhds.get(2);
1922 phds[3]=vhds.get(3);
1923 costheta = vgam.get(3)/vgam.d3mag();
1924 selectmode = mymode;
1925 xs->Fill();
1926 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
1927 }
1928 _modeFlag[1]= _selectedMode;
1929 p->setIntFlag(_modeFlag);
1930 return;
1931 }//end block if(_mode==74110)
1932
1933 // ======================================================
1934 // multi-exclusive mode
1935 //=======================================================
1936 if(_mode==-100){
1937 int mymode;
1938 double pm= EvtRandom::Flat(0.,1);
1939
1940 //std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
1941 //--
1942 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross section within Egam=(0.025-Egamma^max GeV)
1943 if(pm <_xs0/(_xs0 + _xs1) ){//without ISR photon for mode=74110
1944
1945 mhds=_cms;
1946 mymode=selectMode(vmd,mhds);
1947 _selectedMode = mymode;
1948 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1949 delete myxsection;
1950 myxsection =new EvtXsection (mymode);
1951 init_mode(mymode);
1952 resetResMass();
1953
1954 if(_ndaugs>1){//for e+e- -> exclusive decays
1955 checkAA:
1956 p->makeDaughters(_ndaugs,daugs);
1957 double totMass=0;
1958 for(int i=0;i< _ndaugs;i++){
1959 EvtParticle* di=p->getDaug(i);
1960 double xmi=EvtPDL::getMass(di->getId());
1961 di->setMass(xmi);
1962 totMass+=xmi;
1963 }
1964 if(totMass>p->mass()) goto checkAA;
1965 p->initializePhaseSpace( _ndaugs , daugs);
1966 if(!checkdecay(p)) goto checkAA;
1967 vhds = p->getP4();
1968 bool byn_ang;
1969 EvtVector4R pd1 = p->getDaug(0)->getP4();
1970 EvtVector4R pcm(_cms,0,0,0);
1971 //p->printTree();
1972 if(_ndaugs==2){//angular distribution for the hadrons only for two-body decays
1973 byn_ang=hadron_angle_sampling(pd1, pcm);
1974 if(!byn_ang) goto checkAA;
1975 }
1976
1977 }else{// for e+e- -> vector resonace, add a very soft photon
1978 EvtId mydaugs[2];
1979 mydaugs[0]=daugs[0];
1980 EvtPDL::reSetMass(mydaugs[0],mhds-0.001);
1981 //EvtPDL::reSetWidth(mydaugs[0],0);
1982 mydaugs[1]=gamId; //ISR photon
1983 p->makeDaughters(2,mydaugs);
1984 checkBA:
1985 double totMass=0;
1986 for(int i=0;i< 2;i++){
1987 EvtParticle* di=p->getDaug(i);
1988 double xmi=EvtPDL::getMass(di->getId());
1989 di->setMass(xmi);
1990 totMass+=xmi;
1991 }
1992 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
1993 if(totMass>p->mass()) goto checkBA;
1994 p->initializePhaseSpace(2,mydaugs);
1995
1996 if(!checkdecay(p)) goto checkBA;
1997 vhds = p->getDaug(0)->getP4();;
1998 vgam = p->getDaug(1)->getP4();
1999 }
2000 //-----
2001 }else{// with ISR photon for mode=-100
2002Sampling_mhds_A:
2003 mhds=Mhad_sampling(MH,AF);
2004 //std::cout<<"SetMthr= "<<SetMthr<<std::endl;
2005 if(mhds<SetMthr) goto Sampling_mhds_A;
2006 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
2007 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
2008
2009 if(mydbg) mass2=mhds;
2010
2011 //generating events
2012 ModeSelection_A:
2013 mymode=selectMode(vmd,mhds);
2014 if(mymode==-10) goto Sampling_mhds_A;
2015 conexcmode = mymode; //save for model REXC ( see EvtRexc.cc) decay
2016 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
2017 _selectedMode = mymode;
2018 delete myxsection;
2019 myxsection =new EvtXsection (mymode);
2020 init_mode(mymode);
2021 EvtId myvpho = EvtPDL::getId(std::string("vgam"));
2022 EvtPDL::reSetMass(myvpho,mhds); //set to continum cms energy
2023 EvtPDL::reSetWidth(myvpho,0);
2024
2025 //debugging
2026 // for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<" "<<EvtPDL::name(daugs[i])<<std::endl;}
2027
2028 //decay e+e- ->gamma_ISR gamma*
2029 EvtId mydaugs[2];
2030 if(_ndaugs>1) { //gamma* -> exclusive decays
2031 resetResMass();
2032 mydaugs[0]=myvpho;
2033 }else{// vector resonance decays
2034 resetResMass();
2035 mydaugs[0]=daugs[0];
2036 EvtPDL::reSetMass(mydaugs[0],mhds);
2037 //EvtPDL::reSetWidth(mydaugs[0],0);
2038 }
2039 mydaugs[1]=gamId; //ISR photon
2040 int maxflag=0;
2041 int trycount=0;
2042 ISRphotonAng_sampling_A:
2043 double totMass=0;
2044 p->makeDaughters(2,mydaugs);
2045 for(int i=0;i< 2;i++){
2046 EvtParticle* di=p->getDaug(i);
2047 double xmi=EvtPDL::getMass(di->getId());
2048 di->setMass(xmi);
2049 totMass+=xmi;
2050 }
2051 if(totMass>p->mass()) goto ISRphotonAng_sampling_A;
2052 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2053 double weight1 = p->initializePhaseSpace(2,mydaugs);
2054 if(!checkdecay(p)) goto ISRphotonAng_sampling_A;
2055 //set the ISR photon angular distribution
2056 SetP4Rvalue(p,mhds,xeng,theta); //end of decay e+e- -> vpho gamma_ISR
2057
2058 //--debugging
2059 //p->printTree();
2060
2061 if(maxflag==0) {findMaxXS(mhds); maxflag=1;}
2062 vhds = p->getDaug(0)->getP4();
2063 vgam=p->getDaug(1)->getP4();
2064 double gx=vgam.get(1);
2065 double gy=vgam.get(2);
2066 double sintheta= sqrt(gx*gx+gy*gy)/vgam.d3mag();
2067 bool gam_ang=gam_sampling(mhds,sintheta);
2068 trycount++;
2069
2070 } // with and without ISR sampling block
2071 _modeFlag[0]=-100;
2072 _modeFlag[1]= _selectedMode;
2073 p->setIntFlag(_modeFlag);
2074 // finish event generation
2075 // for debugging
2076 if(mydbg){
2077 pgam[0]=vgam.get(0);
2078 pgam[1]=vgam.get(1);
2079 pgam[2]=vgam.get(2);
2080 pgam[3]=vgam.get(3);
2081
2082 phds[0]=vhds.get(0);
2083 phds[1]=vhds.get(1);
2084 phds[2]=vhds.get(2);
2085 phds[3]=vhds.get(3);
2086 costheta = vgam.get(3)/vgam.d3mag();
2087 selectmode = mymode;
2088 xs->Fill();
2089 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
2090 }
2091 return;
2092 }//end block if(_mode==-100)
2093
2094 //=======================================================
2095 // e+e- -> gamma_ISR gamma*
2096 //=======================================================
2097 if(VISR){
2098 EvtId gid=EvtPDL::getId("gamma*");
2099 double pm= EvtRandom::Flat(0.,1);
2100
2101 if(pm <_xs0/(_xs0 + _xs1) ){//with a very soft ISR photon
2102 EvtPDL::reSetMass(gid,(_cms-0.0001) );
2103 mhds = _cms-0.0001;
2104 }else{
2105 mhds=Mhad_sampling(MH,AF);
2106 EvtPDL::reSetMass(gid,mhds);
2107 }
2108 //debugging
2109 std::cout<<"generatedMass "<<EvtPDL::getMeanMass(gid)<<std::endl;
2110 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
2111 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
2112 p->makeDaughters(2,daugs);
2113 for(int i=0;i< 2;i++){
2114 EvtParticle* di=p->getDaug(i);
2115 //if(EvtPDL::name(di->getId())=="gamma*") di->setVectorSpinDensity();
2116 double xmi=EvtPDL::getMeanMass(di->getId());
2117 di->setMass(xmi);
2118 }
2119 p->initializePhaseSpace(2,daugs);
2120 SetP4(p,mhds,xeng,theta);
2121 return;
2122 }
2123
2124
2125 //========================================================
2126 //=== for other mode generation with _mode != 74110
2127 //========================================================
2128
2129 if((_xs0 + _xs1)==0) {std::cout<<"EvtConExc::zero total xsection"<<std::endl;::abort();}
2130 double pm= EvtRandom::Flat(0.,1);
2131 double xsratio = _xs0/(_xs0 + _xs1);
2132 int iflag=2; //flag = 0 for ee->hadrons, 1 for ee->gamma_ISR hadrons, 2: mix 0 and 1 case
2133 if(getArg(0)!= -2){// exclude DIY case
2134 if(getNArg()==3 && radflag==1){iflag=1;xsratio=0;} // only generating gamma hadrons mode
2135 else if(getNArg()==3 && radflag==0) {iflag=0;xsratio=1;}
2136 }
2137
2138 // std::cout<<"xsratio= "<<xsratio<<std::endl;
2139
2140 if(pm<xsratio ){// no ISR photon
2141 masscheck:
2142 double summassx=0;
2143 p->makeDaughters(_ndaugs,daugs);
2144 for(int i=0;i< _ndaugs;i++){
2145 EvtParticle* di=p->getDaug(i);
2146 double xmi=EvtPDL::getMass(di->getId());
2147 di->setMass(xmi);
2148 summassx += xmi;
2149 //std::cout<<"PID and mass "<<di->getId()<<" "<<xmi<<std::endl;
2150 }
2151 if(summassx>p->mass()) goto masscheck;
2152 angular_hadron:
2153 p->initializePhaseSpace(_ndaugs,daugs);
2154 bool byn_ang;
2155 EvtVector4R pd1 = p->getDaug(0)->getP4();
2156 EvtVector4R pcm(_cms,0,0,0);
2157 if(_ndaugs==2){//angular distribution for the hadrons only for two-body decays
2158 byn_ang=hadron_angle_sampling(pd1, pcm);
2159 if(!byn_ang) goto angular_hadron;
2160 }
2161
2162 //for histogram
2163 cosp = pd1.get(3)/pd1.d3mag();
2164 mhds = _cms;
2165 //std::cout<<"cosp, mhds "<<cosp<<" "<<mhds<<std::endl;
2166 //p->printTree();
2167
2168 }else{// sampling m_hadrons and decay e+e- ->gamma gamma*
2169 double mhdr=Mhad_sampling(MH,AF);
2170 double xeng=1-mhdr*mhdr/(_cms*_cms); //photon energy ratio
2171 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
2172 EvtId gid =EvtPDL::getId("vhdr");
2173 EvtPDL::reSetMass(gid,mhdr);
2174 int ndaugs =2;
2175 EvtId mydaugs[2];
2176 mydaugs[0] =EvtPDL::getId(std::string("gamma"));
2177 mydaugs[1]=EvtPDL::getId(std::string("vhdr"));
2178
2179 //for histogram
2180 mhds = mhdr;
2181 costheta = cos(theta);
2182 //--
2183
2184 p->makeDaughters(2,mydaugs);
2185 for(int i=0;i< 2;i++){
2186 EvtParticle* di=p->getDaug(i);
2187 double xmi=EvtPDL::getMass(di->getId());
2188 //if(EvtPDL::name(di->getId())=="vhdr") di->setVectorSpinDensity();
2189 di->setMass(xmi);
2190 }
2191 p->initializePhaseSpace(2,mydaugs);
2192 SetP4(p,mhdr,xeng,theta); //end of decay e+e- -> gamma_ISR gamma*
2193 //p->printTree();
2194
2195 //----
2196 }//end of gamma gamma*, gamma*->hadrons generation
2197 // p->printTree();
2198
2199 //-----------------------------------------
2200 // End of decays for all topology
2201 //----------------------------------------
2202 //================================= event analysis
2203 if(_nevt ==0 ){
2204 std::cout<<"The decay chain: "<<std::endl;
2205 p->printTree();
2206 }
2207 _nevt++;
2208 //--- for debuggint to make root file
2209 if(mydbg){
2210 xs->Fill();
2211 }
2212
2213 //----
2214 return ;
2215}
EvtConExc::checkdecay
bool checkdecay(EvtParticle *p)
Definition: EvtConExc.cc:3123
EvtConExc::findMaxXS
void findMaxXS(EvtParticle *p)
Definition: EvtConExc.cc:2390
EvtConExc::conexcmode
static int conexcmode
Definition: EvtConExc.hh:161
EvtConExc::SetP4Rvalue
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:2917
EvtConExc::photonSampling
bool photonSampling(EvtParticle *part)
Definition: EvtConExc.cc:3219
EvtConExc::SetMthr
static double SetMthr
Definition: EvtConExc.hh:142
EvtConExc::selectMode
int selectMode(std::vector< int > vmod, double mhds)
Definition: EvtConExc.cc:2961
EvtConExc::init_mode
void init_mode(int mode)
Definition: EvtConExc.cc:710
EvtConExc::angularSampling
bool angularSampling(EvtParticle *part)
Definition: EvtConExc.cc:3173
EvtConExc::ISR_ang_sampling
double ISR_ang_sampling(double x)
Definition: EvtConExc.cc:2869
EvtConExc::Mhad_sampling
double Mhad_sampling(double *x, double *y)
Definition: EvtConExc.cc:2823
EvtConExc::resetResMass
void resetResMass()
Definition: EvtConExc.cc:3041
EvtConExc::calAF
void calAF(double myecms)
Definition: EvtConExc.cc:3647
EvtConExc::hadron_angle_sampling
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition: EvtConExc.cc:2217
EvtConExc::gam_sampling
bool gam_sampling(EvtParticle *p)
Definition: EvtConExc.cc:2462
EvtConExc::SetP4
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:2893
EvtDecayBase::getArg
double getArg(int j)
Definition: EvtDecayBase.cc:564
EvtGlobalSet::ConExcPythia
static bool ConExcPythia
Definition: EvtGlobalSet.hh:20
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45
EvtPDL::reSetWidth
static void reSetWidth(EvtId i, double width)
Definition: EvtPDL.hh:69
EvtPDL::reSetMass
static void reSetMass(EvtId i, double mass)
Definition: EvtPDL.hh:68
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
EvtPDL::getMass
static double getMass(EvtId i)
Definition: EvtPDL.hh:46
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
EvtParticle::setMass
void setMass(double m)
Definition: EvtParticle.hh:372
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1177
EvtParticle::setIntFlag
void setIntFlag(std::vector< int > vi)
Definition: EvtParticle.hh:151
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
EvtParticle::printTree
void printTree() const
Definition: EvtParticle.cc:897
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1071
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cc:39
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:179
EvtVector4R::d3mag
double d3mag() const
Definition: EvtVector4R.cc:186

◆ difgamXs() [1/2]

double EvtConExc::difgamXs ( double  mhds,
double  sintheta 
)

Definition at line 2950 of file EvtConExc.cc.

2950 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2951 double x = 1-(mhds/_cms)*(mhds/_cms);
2952 double sin2theta = sintheta*sintheta;
2953 double alpha = 1./137.;
2954 double pi = 3.1415926;
2955 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2956 double xs = myxsection->getXsection(mhds);
2957 double difxs = 2*mhds/_cms/_cms * wsx *xs;
2958 return difxs;
2959}

◆ difgamXs() [2/2]

double EvtConExc::difgamXs ( EvtParticle *  p)

Definition at line 2435 of file EvtConExc.cc.

2435 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2436 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2437 EvtVector4R p0 = p->getDaug(0)->getP4();
2438 for(int i=1;i<_ndaugs;i++){
2439 p0 += p->getDaug(i)->getP4();
2440 }
2441 EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
2442 double mhs = p0.mass();
2443 double egam = pgam.get(0);
2444 double sin2theta = pgam.get(3)/pgam.d3mag();
2445 sin2theta = 1-sin2theta*sin2theta;
2446
2447 double cms = p->getP4().mass();
2448 double alpha = 1./137.;
2449 double pi = 3.1415926;
2450 double x = 2*egam/cms;
2451 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2452 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2453
2454 double xs = myxsection->getXsection(mhs);
2455 double difxs = 2*mhs/cms/cms * wsx *xs;
2456 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2457 return difxs;
2458
2459}
EvtXsection::getErr
double getErr(double mx)
Definition: EvtXsection.cc:12014

Referenced by findMaxXS(), gam_sampling(), and gamHXSection().

◆ Egam2Mhds()

double EvtConExc::Egam2Mhds ( double  Egam)

Definition at line 3423 of file EvtConExc.cc.

3423 {
3424 double s=_cms*_cms;
3425 double mhds = sqrt( s - 2*Egam*_cms);
3426 return mhds;
3427}

◆ energySpread()

double EvtConExc::energySpread ( double  mu,
double  sigma 
)

Definition at line 3629 of file EvtConExc.cc.

3629 {
3630 //mu: mean value in Gaussian
3631 //sigma: variance in Gaussian
3632 rloop:
3633 int n=12;
3634 double ri=0;
3635 for(int i=0;i<n;i++){
3636 double pm= EvtRandom::Flat(0.,1);
3637 ri += pm;
3638 }
3639 double eta=sqrt(n*12.0)*(ri/12-0.5);
3640 double xsig= sigma*eta+mu;//smapling Gaussion value
3641 double mx0=staxsection->getXlw();
3642 double mx1=staxsection->getXup();
3643 if(xsig<mx0 || xsig>mx1) goto rloop;
3644 return xsig;
3645}
n
const Int_t n
Definition: DataBase/tau_mode.c:65
EvtXsection::getXup
double getXup()
Definition: EvtXsection.hh:155

◆ findMaxXS() [1/2]

void EvtConExc::findMaxXS ( double  mhds)

Definition at line 2935 of file EvtConExc.cc.

2935 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2936 maxXS=-1;
2937 for(int i=0;i<90000;i++){
2938 double x = 1-(mhds/_cms)*(mhds/_cms);
2939 double cos = EvtRandom::Flat(0.0006,0.9994); //set cut on ISR gamma 2^o
2940 double sin2theta =sqrt(1-cos*cos);
2941 double alpha = 1./137.;
2942 double pi = 3.1415926;
2943 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2944 double xs = myxsection->getXsection(mhds);
2945 double difxs = 2*mhds/_cms/_cms * wsx *xs;
2946 if(difxs>maxXS)maxXS=difxs;
2947 }
2948}

◆ findMaxXS() [2/2]

void EvtConExc::findMaxXS ( EvtParticle *  p)

Definition at line 2390 of file EvtConExc.cc.

2390 {
2391 //std::cout<<"nmc= "<<nmc<<std::endl;
2392 gamId = EvtPDL::getId(std::string("gamma"));
2393 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
2394 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
2395 double totxs = 0;
2396 maxXS=-1;
2397 _er1=0;
2398 Rad2Xs =0;
2399 int nmc = 50000;
2400 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
2401 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
2402 int gamdaugs = _ndaugs+1;
2403
2404 EvtId theDaugs[10];
2405 for(int i=0;i<_ndaugs;i++){
2406 theDaugs[i] = daugs[i];
2407 }
2408 theDaugs[_ndaugs]=gamId;
2409
2410 gamH->makeDaughters(gamdaugs,theDaugs);
2411 loop:
2412 double totm=0;
2413 for(int i=0;i<gamdaugs;i++){
2414 EvtParticle* di=gamH->getDaug(i);
2415 double xmi=EvtPDL::getMass(di->getId());
2416 di->setMass(xmi);
2417 totm += xmi;
2418 }
2419
2420 //std::cout<<totm<<" "<<p_init.get(0)<<std::endl;
2421 if(totm >= p_init.get(0) ) goto loop;
2422
2423 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
2424 double thexs = difgamXs(gamH); //prepare the photon angular distribution
2425 EvtVector4R p4gam = gamH->getDaug(_ndaugs)->getP4();
2426 double costheta = p4gam.get(3)/p4gam.d3mag();
2427 double sintheta = sqrt(1-costheta*costheta);
2428 bool acut=(sintheta>0.04) && (sintheta<0.96); //set photon auglar cut 2^o
2429 if(thexs>maxXS && acut ) {maxXS=thexs;}
2430 //gamH->deleteTree();
2431 }
2432
2433}
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition: KarFin.h:34
EvtConExc::difgamXs
double difgamXs(EvtParticle *p)
Definition: EvtConExc.cc:2435
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition: EvtParticleFactory.cc:46

Referenced by decay(), and init().

◆ gam_sampling() [1/2]

bool EvtConExc::gam_sampling ( double  mhds,
double  sintheta 
)

Definition at line 2469 of file EvtConExc.cc.

2469 {
2470 double pm= EvtRandom::Flat(0.,1);
2471 double xs = difgamXs( mhds,sintheta );
2472 double xsratio = xs/maxXS;
2473 if(pm<xsratio){return true;}else{return false;}
2474}

◆ gam_sampling() [2/2]

bool EvtConExc::gam_sampling ( EvtParticle *  p)

Definition at line 2462 of file EvtConExc.cc.

2462 {//photon angular distribution sampling
2463 double pm= EvtRandom::Flat(0.,1);
2464 double xs = difgamXs( p );
2465 double xsratio = xs/maxXS;
2466 if(pm<xsratio){return true;}else{return false;}
2467}

Referenced by decay().

◆ gamHXSection() [1/4]

double EvtConExc::gamHXSection ( double  El,
double  Eh 
)

Definition at line 2336 of file EvtConExc.cc.

2336 {// using Gaussian integration subroutine from Cern lib
2337 std::cout<<" "<<std::endl;
2338 extern double Rad2difXs(double *x);
2339 extern double Rad2difXs2(double *x);
2340 double eps = 0.01;
2341 double Rad2Xs;
2342 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2343 if(Rad2Xs==0){
2344 for(int iter=0;iter<10;iter++){
2345 eps += 0.01;
2346 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2347 if(!Rad2Xs) return Rad2Xs;
2348 }
2349 }
2350 return Rad2Xs; // the leading second order correction
2351}
Rad2difXs2
double Rad2difXs2(double *mhs)
Definition: EvtConExc.cc:2717
dgauss_
double dgauss_(double(*fun)(double *), double *, double *, double *)
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
eps
EvtTensor3C eps(const EvtVector3R &v)
Definition: EvtTensor3C.cc:307
EvtConExc::Rad2difXs
double Rad2difXs(EvtParticle *p)
Definition: EvtConExc.cc:2563

◆ gamHXSection() [2/4]

double EvtConExc::gamHXSection ( double  El,
double  Eh,
int  mode 
)

Definition at line 2354 of file EvtConExc.cc.

2354 {// using Gaussian integration subroutine from Cern lib
2355 std::cout<<" "<<std::endl;
2356 extern double Rad2difXs(double *x);
2357 extern double Rad2difXs2(double *x);
2358 double eps = 0.01;
2359 double Rad2Xs;
2360 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2361 if(Rad2Xs==0){
2362 for(int iter=0;iter<10;iter++){
2363 eps += 0.01;
2364 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2365 if(!Rad2Xs) return Rad2Xs;
2366 }
2367 }
2368 return Rad2Xs; // the leading second order correction
2369}

◆ gamHXSection() [3/4]

double EvtConExc::gamHXSection ( double  s,
double  El,
double  Eh,
int  nmc = 100000 
)

Definition at line 2304 of file EvtConExc.cc.

2304 {//El, Eh : the lower and higher energy of radiative photons
2305 //--for narrow resonance psi(2S),J/psi, phi, which can not calculated with the MC integration
2306 //double mxL = sqrt( s-2*Eh*sqrt(s) ); //low mass
2307 //double mxH = sqrt( s-2*El*sqrt(s) ); //high mass
2308 //double sigv = narrowRXS(mxL,mxH);
2309 //--
2310
2311 double totxs = 0;
2312 maxXS=-1;
2313 _er1=0;
2314 Rad2Xs =0;
2315 double xL=2*El/sqrt(s);
2316 double xH=2*Eh/sqrt(s);
2317 for(int i=0;i<nmc;i++){//It is found that the narrow resonance can not included in this MC integartion
2318 double rdm = EvtRandom::Flat(0.,1.);// set angular cut 1^o
2319 double x= xL+ (xH-xL)*rdm;
2320 Rad2Xs += Rad2difXs(s,x);
2321 _er1 += differ2; //stored when calculate Rad2Xs
2322 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
2323 }
2324 _er1/=nmc;
2325 _er1*=(xH-xL);
2326 //std::cout<<"_er1= "<<_er1<<std::endl;
2327 Rad2Xs/=nmc; // the leading second order correction
2328 double thexs= Rad2Xs*(xH-xL); //xH-xL is the Jacobi factor
2329 //std::cout<<"thexs= "<<thexs<<std::endl;
2330 //if(sigv != -1) thexs += sigv; //add narrow resonance ISR cross section
2331 return thexs;
2332}

◆ gamHXSection() [4/4]

double EvtConExc::gamHXSection ( EvtParticle *  p,
double  El,
double  Eh,
int  nmc = 100000 
)

Definition at line 2245 of file EvtConExc.cc.

2245 {//El, Eh : the lower and higher energy of radiative photons
2246 //std::cout<<"nmc= "<<nmc<<std::endl;
2247 gamId = EvtPDL::getId(std::string("gamma"));
2248 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
2249 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
2250 double totxs = 0;
2251 maxXS=-1;
2252 _er1=0;
2253 Rad2Xs =0;
2254 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
2255 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
2256 int gamdaugs = _ndaugs+1;
2257
2258 EvtId theDaugs[10];
2259 for(int i=0;i<_ndaugs;i++){
2260 theDaugs[i] = daugs[i];
2261 }
2262 theDaugs[_ndaugs]=gamId;
2263
2264 gamH->makeDaughters(gamdaugs,theDaugs);
2265
2266 for(int i=0;i<gamdaugs;i++){
2267 EvtParticle* di=gamH->getDaug(i);
2268 double xmi=EvtPDL::getMass(di->getId());
2269 di->setMass(xmi);
2270 }
2271 loop:
2272 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
2273 //-- slection:theta_gamma > 1 degree
2274 EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
2275 double pmag = pgam.d3mag();
2276 double pz = pgam.get(3);
2277 double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
2278 double onedegree = 1./180.*3.1415926;
2279 //if(sintheta < onedegree) {goto loop;}
2280 if(pmag <El || pmag >Eh) {goto loop;}
2281 //--------
2282 // std::cout<<"pmag= "<<pmag<<std::endl;
2283
2284 double thexs = difgamXs(gamH); //prepare the photon angular distribution
2285 Rad2Xs += Rad2difXs( gamH );
2286 if(thexs>maxXS) {maxXS=thexs;}
2287 double s = p_init.mass2();
2288 double x = 2*pgam.get(0)/sqrt(s);
2289 double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
2290 totxs += rad1xs;
2291 _er1 += differ;
2292 gamH->deleteDaughters();
2293 } //for int_i block
2294 _er1/=nmc;
2295
2296 Rad2Xs/=nmc; // the leading second order correction
2297 totxs/=nmc; // first order correction XS
2298
2299 // return totxs; // first order correction XS
2300 return Rad2Xs; // the leading second order correction
2301}
EvtConExc::Rad1difXs
double Rad1difXs(EvtParticle *p)
Definition: EvtConExc.cc:2617
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
EvtParticle::deleteDaughters
void deleteDaughters(bool keepChannel=false)
Definition: EvtParticle.cc:540

Referenced by init().

◆ gamHXSection_er()

double EvtConExc::gamHXSection_er ( double  El,
double  Eh 
)

Definition at line 2372 of file EvtConExc.cc.

2372 {// using Gaussian integration subroutine from Cern lib
2373 std::cout<<" "<<std::endl;
2374 extern double Rad2difXs_er(double *x);
2375 extern double Rad2difXs_er2(double *x);
2376 double eps = 0.01;
2377 double Rad2Xs;
2378 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
2379 if(Rad2Xs==0){
2380 for(int iter=0;iter<10;iter++){
2381 eps += 0.01;
2382 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
2383 if(!Rad2Xs) return Rad2Xs;;
2384 }
2385 }
2386 return Rad2Xs; // the leading second order correction
2387}
Rad2difXs_er2
double Rad2difXs_er2(double *mhs)
Definition: EvtConExc.cc:2732
Rad2difXs_er
double Rad2difXs_er(double *m)
Definition: EvtConExc.cc:2704

◆ get_mode()

std::vector< EvtId > EvtConExc::get_mode ( int  mode)

Definition at line 1210 of file EvtConExc.cc.

1210 {
1211 int _ndaugs;
1212 EvtId daugs[100];
1213 if(mode ==0 ){
1214 _ndaugs =2;
1215 daugs[0]=EvtPDL::getId(std::string("p+"));
1216 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1217 }else if(mode ==1 ){
1218 _ndaugs =2;
1219 daugs[0]=EvtPDL::getId(std::string("n0"));
1220 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
1221 }else if(mode ==2 ){
1222 _ndaugs =2;
1223 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
1224 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
1225 }else if(mode ==3 ){
1226 _ndaugs =2;
1227 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
1228 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
1229 }else if(mode ==4 ){
1230 _ndaugs =2;
1231 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
1232 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
1233 }else if(mode ==5 ){
1234 _ndaugs =2;
1235 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
1236 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
1237 } else if(mode ==6 ){
1238 _ndaugs =2;
1239 daugs[0]=EvtPDL::getId(std::string("pi+"));
1240 daugs[1]=EvtPDL::getId(std::string("pi-"));
1241 }else if(mode ==7 ){
1242 _ndaugs =3;
1243 daugs[0]=EvtPDL::getId(std::string("pi+"));
1244 daugs[1]=EvtPDL::getId(std::string("pi-"));
1245 daugs[2]=EvtPDL::getId(std::string("pi0"));
1246 }else if(mode ==8 ){
1247 _ndaugs =3;
1248 daugs[0]=EvtPDL::getId(std::string("K+"));
1249 daugs[1]=EvtPDL::getId(std::string("K-"));
1250 daugs[2]=EvtPDL::getId(std::string("pi0"));
1251 }else if(mode ==9 ){
1252 _ndaugs =3;
1253 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1254 daugs[1]=EvtPDL::getId(std::string("K+"));
1255 daugs[2]=EvtPDL::getId(std::string("pi-"));
1256 }else if(mode ==10 ){
1257 _ndaugs =3;
1258 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1259 daugs[1]=EvtPDL::getId(std::string("K-"));
1260 daugs[2]=EvtPDL::getId(std::string("pi+"));
1261 }else if(mode ==11 ){
1262 _ndaugs =3;
1263 daugs[0]=EvtPDL::getId(std::string("K+"));
1264 daugs[1]=EvtPDL::getId(std::string("K-"));
1265 daugs[2]=EvtPDL::getId(std::string("eta"));
1266 }else if(mode ==12 ){
1267 _ndaugs =4;
1268 daugs[0]=EvtPDL::getId(std::string("pi+"));
1269 daugs[1]=EvtPDL::getId(std::string("pi-"));
1270 daugs[2]=EvtPDL::getId(std::string("pi+"));
1271 daugs[3]=EvtPDL::getId(std::string("pi-"));
1272 }else if(mode ==13 ){
1273 _ndaugs =4;
1274 daugs[0]=EvtPDL::getId(std::string("pi+"));
1275 daugs[1]=EvtPDL::getId(std::string("pi-"));
1276 daugs[2]=EvtPDL::getId(std::string("pi0"));
1277 daugs[3]=EvtPDL::getId(std::string("pi0"));
1278 }else if(mode ==14 ){
1279 _ndaugs =4;
1280 daugs[0]=EvtPDL::getId(std::string("K+"));
1281 daugs[1]=EvtPDL::getId(std::string("K-"));
1282 daugs[2]=EvtPDL::getId(std::string("pi+"));
1283 daugs[3]=EvtPDL::getId(std::string("pi-"));
1284 }else if(mode ==15 ){
1285 _ndaugs =4;
1286 daugs[0]=EvtPDL::getId(std::string("K+"));
1287 daugs[1]=EvtPDL::getId(std::string("K-"));
1288 daugs[2]=EvtPDL::getId(std::string("pi0"));
1289 daugs[3]=EvtPDL::getId(std::string("pi0"));
1290 }else if(mode ==16 ){
1291 _ndaugs =4;
1292 daugs[0]=EvtPDL::getId(std::string("K+"));
1293 daugs[1]=EvtPDL::getId(std::string("K-"));
1294 daugs[2]=EvtPDL::getId(std::string("K+"));
1295 daugs[3]=EvtPDL::getId(std::string("K-"));
1296 }else if(mode ==17 ){
1297 _ndaugs =5;
1298 daugs[0]=EvtPDL::getId(std::string("pi+"));
1299 daugs[1]=EvtPDL::getId(std::string("pi-"));
1300 daugs[2]=EvtPDL::getId(std::string("pi+"));
1301 daugs[3]=EvtPDL::getId(std::string("pi-"));
1302 daugs[4]=EvtPDL::getId(std::string("pi0"));
1303 }else if(mode ==18 ){
1304 _ndaugs =5;
1305 daugs[0]=EvtPDL::getId(std::string("pi+"));
1306 daugs[1]=EvtPDL::getId(std::string("pi-"));
1307 daugs[2]=EvtPDL::getId(std::string("pi+"));
1308 daugs[3]=EvtPDL::getId(std::string("pi-"));
1309 daugs[4]=EvtPDL::getId(std::string("eta"));
1310 }else if(mode ==19 ){
1311 _ndaugs =5;
1312 daugs[0]=EvtPDL::getId(std::string("K+"));
1313 daugs[1]=EvtPDL::getId(std::string("K-"));
1314 daugs[2]=EvtPDL::getId(std::string("pi+"));
1315 daugs[3]=EvtPDL::getId(std::string("pi-"));
1316 daugs[4]=EvtPDL::getId(std::string("pi0"));
1317 }else if(mode ==20 ){
1318 _ndaugs =5;
1319 daugs[0]=EvtPDL::getId(std::string("K+"));
1320 daugs[1]=EvtPDL::getId(std::string("K-"));
1321 daugs[2]=EvtPDL::getId(std::string("pi+"));
1322 daugs[3]=EvtPDL::getId(std::string("pi-"));
1323 daugs[4]=EvtPDL::getId(std::string("eta"));
1324 }else if(mode ==21 ){
1325 _ndaugs =6;
1326 daugs[0]=EvtPDL::getId(std::string("pi+"));
1327 daugs[1]=EvtPDL::getId(std::string("pi-"));
1328 daugs[2]=EvtPDL::getId(std::string("pi+"));
1329 daugs[3]=EvtPDL::getId(std::string("pi-"));
1330 daugs[4]=EvtPDL::getId(std::string("pi+"));
1331 daugs[5]=EvtPDL::getId(std::string("pi-"));
1332 }else if(mode ==22 ){
1333 _ndaugs =6;
1334 daugs[0]=EvtPDL::getId(std::string("pi+"));
1335 daugs[1]=EvtPDL::getId(std::string("pi-"));
1336 daugs[2]=EvtPDL::getId(std::string("pi+"));
1337 daugs[3]=EvtPDL::getId(std::string("pi-"));
1338 daugs[4]=EvtPDL::getId(std::string("pi0"));
1339 daugs[5]=EvtPDL::getId(std::string("pi0"));
1340 }else if(mode == 23){
1341 _ndaugs =2;
1342 daugs[0]=EvtPDL::getId(std::string("eta"));
1343 daugs[1]=EvtPDL::getId(std::string("phi"));
1344 }else if(mode == 24){
1345 _ndaugs =2;
1346 daugs[0]=EvtPDL::getId(std::string("phi"));
1347 daugs[1]=EvtPDL::getId(std::string("pi0"));
1348 }else if(mode == 25){
1349 _ndaugs =2;
1350 daugs[0]=EvtPDL::getId(std::string("K+"));
1351 daugs[1]=EvtPDL::getId(std::string("K*-"));
1352 }else if(mode == 26){
1353 _ndaugs =2;
1354 daugs[0]=EvtPDL::getId(std::string("K-"));
1355 daugs[1]=EvtPDL::getId(std::string("K*+"));
1356 }else if(mode == 27){
1357 _ndaugs =2;
1358 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1359 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
1360 }else if(mode == 28){
1361 _ndaugs =3;
1362 daugs[0]=EvtPDL::getId(std::string("K*0"));
1363 daugs[1]=EvtPDL::getId(std::string("K+"));
1364 daugs[2]=EvtPDL::getId(std::string("pi-"));
1365 }else if(mode == 29){
1366 _ndaugs =3;
1367 daugs[0]=EvtPDL::getId(std::string("K*0"));
1368 daugs[1]=EvtPDL::getId(std::string("K-"));
1369 daugs[2]=EvtPDL::getId(std::string("pi+"));
1370 }else if(mode == 30){
1371 _ndaugs =3;
1372 daugs[0]=EvtPDL::getId(std::string("K*+"));
1373 daugs[1]=EvtPDL::getId(std::string("K-"));
1374 daugs[2]=EvtPDL::getId(std::string("pi0"));
1375 }else if(mode == 31){
1376 _ndaugs =3;
1377 daugs[0]=EvtPDL::getId(std::string("K*-"));
1378 daugs[1]=EvtPDL::getId(std::string("K+"));
1379 daugs[2]=EvtPDL::getId(std::string("pi0"));
1380 }else if(mode == 32){
1381 _ndaugs =3;
1382 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
1383 daugs[1]=EvtPDL::getId(std::string("K+"));
1384 daugs[2]=EvtPDL::getId(std::string("pi-"));
1385 }else if(mode == 33){
1386 _ndaugs =3;
1387 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
1388 daugs[1]=EvtPDL::getId(std::string("K-"));
1389 daugs[2]=EvtPDL::getId(std::string("pi+"));
1390 }else if(mode == 34){
1391 _ndaugs =3;
1392 daugs[0]=EvtPDL::getId(std::string("K+"));
1393 daugs[1]=EvtPDL::getId(std::string("K-"));
1394 daugs[2]=EvtPDL::getId(std::string("rho0"));
1395 }else if(mode == 35){
1396 _ndaugs =3;
1397 daugs[0]=EvtPDL::getId(std::string("phi"));
1398 daugs[1]=EvtPDL::getId(std::string("pi-"));
1399 daugs[2]=EvtPDL::getId(std::string("pi+"));
1400 }else if(mode == 36){
1401 _ndaugs =2;
1402 daugs[0]=EvtPDL::getId(std::string("phi"));
1403 daugs[1]=EvtPDL::getId(std::string("f_0"));
1404 }else if(mode == 37){
1405 _ndaugs =3;
1406 daugs[0]=EvtPDL::getId(std::string("eta"));
1407 daugs[1]=EvtPDL::getId(std::string("pi+"));
1408 daugs[2]=EvtPDL::getId(std::string("pi-"));
1409 }else if(mode == 38){
1410 _ndaugs =3;
1411 daugs[0]=EvtPDL::getId(std::string("omega"));
1412 daugs[1]=EvtPDL::getId(std::string("pi+"));
1413 daugs[2]=EvtPDL::getId(std::string("pi-"));
1414 }else if(mode == 39){
1415 _ndaugs =2;
1416 daugs[0]=EvtPDL::getId(std::string("omega"));
1417 daugs[1]=EvtPDL::getId(std::string("f_0"));
1418 }else if(mode == 40){
1419 _ndaugs =3;
1420 daugs[0]=EvtPDL::getId(std::string("eta'"));
1421 daugs[1]=EvtPDL::getId(std::string("pi+"));
1422 daugs[2]=EvtPDL::getId(std::string("pi-"));
1423 }else if(mode == 41){
1424 _ndaugs =3;
1425 daugs[0]=EvtPDL::getId(std::string("f_1"));
1426 daugs[1]=EvtPDL::getId(std::string("pi+"));
1427 daugs[2]=EvtPDL::getId(std::string("pi-"));
1428 }else if(mode == 42){
1429 _ndaugs =3;
1430 daugs[0]=EvtPDL::getId(std::string("omega"));
1431 daugs[1]=EvtPDL::getId(std::string("K+"));
1432 daugs[2]=EvtPDL::getId(std::string("K-"));
1433 }else if(mode == 43){
1434 _ndaugs =4;
1435 daugs[0]=EvtPDL::getId(std::string("omega"));
1436 daugs[1]=EvtPDL::getId(std::string("pi+"));
1437 daugs[2]=EvtPDL::getId(std::string("pi-"));
1438 daugs[3]=EvtPDL::getId(std::string("pi0"));
1439 }else if(mode == 44){
1440 _ndaugs =2;
1441 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
1442 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
1443 }else if(mode == 45){
1444 _ndaugs =2;
1445 daugs[0]=EvtPDL::getId(std::string("K+"));
1446 daugs[1]=EvtPDL::getId(std::string("K-"));
1447 }else if(mode == 46){
1448 _ndaugs =2;
1449 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1450 daugs[1]=EvtPDL::getId(std::string("K_L0"));
1451 }else if(mode == 47){
1452 _ndaugs =2;
1453 daugs[0]=EvtPDL::getId(std::string("omega"));
1454 daugs[1]=EvtPDL::getId(std::string("eta"));
1455 }else if(mode == 48){
1456 _ndaugs =2;
1457 daugs[0]=EvtPDL::getId(std::string("phi"));
1458 daugs[1]=EvtPDL::getId(std::string("eta'"));
1459 }else if(mode == 49){
1460 _ndaugs =3;
1461 daugs[0]=EvtPDL::getId(std::string("phi"));
1462 daugs[1]=EvtPDL::getId(std::string("K+"));
1463 daugs[2]=EvtPDL::getId(std::string("K-"));
1464 }else if(mode == 50){
1465 _ndaugs =3;
1466 daugs[0]=EvtPDL::getId(std::string("p+"));
1467 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1468 daugs[2]=EvtPDL::getId(std::string("pi0"));
1469 }else if(mode == 51){
1470 _ndaugs =3;
1471 daugs[0]=EvtPDL::getId(std::string("p+"));
1472 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1473 daugs[2]=EvtPDL::getId(std::string("eta"));
1474 }else if(mode == 52){
1475 _ndaugs =2;
1476 daugs[0]=EvtPDL::getId(std::string("omega"));
1477 daugs[1]=EvtPDL::getId(std::string("pi0"));
1478
1479 }else if(mode == 59){
1480 _ndaugs =3;
1481 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1482 daugs[1]=EvtPDL::getId(std::string("D0"));
1483 daugs[2]=EvtPDL::getId(std::string("pi0"));
1484 }else if(mode == 60){
1485 _ndaugs =3;
1486 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1487 daugs[1]=EvtPDL::getId(std::string("D*0"));
1488 daugs[2]=EvtPDL::getId(std::string("pi0"));
1489 }else if(mode == 61){
1490 _ndaugs =3;
1491 daugs[0]=EvtPDL::getId(std::string("D0"));
1492 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1493 daugs[2]=EvtPDL::getId(std::string("pi0"));
1494 }else if(mode == 62){
1495 _ndaugs =3;
1496 daugs[0]=EvtPDL::getId(std::string("D+"));
1497 daugs[1]=EvtPDL::getId(std::string("D*-"));
1498 daugs[2]=EvtPDL::getId(std::string("pi0"));
1499 }else if(mode == 63){
1500 _ndaugs =3;
1501 daugs[0]=EvtPDL::getId(std::string("D-"));
1502 daugs[1]=EvtPDL::getId(std::string("D+"));
1503 daugs[2]=EvtPDL::getId(std::string("pi0"));
1504 }else if(mode == 64){
1505 _ndaugs =3;
1506 daugs[0]=EvtPDL::getId(std::string("D-"));
1507 daugs[1]=EvtPDL::getId(std::string("D*+"));
1508 daugs[2]=EvtPDL::getId(std::string("pi0"));
1509 }else if(mode == 65){
1510 _ndaugs =3;
1511 daugs[0]=EvtPDL::getId(std::string("D-"));
1512 daugs[1]=EvtPDL::getId(std::string("D*0"));
1513 daugs[2]=EvtPDL::getId(std::string("pi+"));
1514 }else if(mode == 66){
1515 _ndaugs =3;
1516 daugs[0]=EvtPDL::getId(std::string("D+"));
1517 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1518 daugs[2]=EvtPDL::getId(std::string("pi-"));
1519 }else if(mode == 67){
1520 _ndaugs =2;
1521 daugs[0]=EvtPDL::getId(std::string("D*0"));
1522 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1523 }else if(mode == 68){
1524 _ndaugs =2;
1525 daugs[0]=EvtPDL::getId(std::string("D0"));
1526 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1527
1528 }else if(mode == 69){
1529 _ndaugs =2;
1530 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1531 daugs[1]=EvtPDL::getId(std::string("D*0"));
1532
1533 }else if(mode == 70){
1534 _ndaugs =2;
1535 daugs[0]=EvtPDL::getId(std::string("D0"));
1536 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
1537
1538}else if(mode == 71){
1539 _ndaugs =2;
1540 daugs[0]=EvtPDL::getId(std::string("D+"));
1541 daugs[1]=EvtPDL::getId(std::string("D-"));
1542 }else if(mode == 72){
1543 _ndaugs =2;
1544 daugs[0]=EvtPDL::getId(std::string("D+"));
1545 daugs[1]=EvtPDL::getId(std::string("D*-"));
1546
1547 }else if(mode == 73){
1548 _ndaugs =2;
1549 daugs[0]=EvtPDL::getId(std::string("D-"));
1550 daugs[1]=EvtPDL::getId(std::string("D*+"));
1551
1552 }else if(mode == 74){
1553 _ndaugs =2;
1554 daugs[0]=EvtPDL::getId(std::string("D*+"));
1555 daugs[1]=EvtPDL::getId(std::string("D*-"));
1556
1557 }else if(mode == 75){
1558 _ndaugs =3;
1559 daugs[0]=EvtPDL::getId(std::string("D0"));
1560 daugs[1]=EvtPDL::getId(std::string("D-"));
1561 daugs[2]=EvtPDL::getId(std::string("pi+"));
1562 }else if(mode == 76){
1563 _ndaugs =3;
1564 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1565 daugs[1]=EvtPDL::getId(std::string("D+"));
1566 daugs[2]=EvtPDL::getId(std::string("pi-"));
1567 }else if(mode == 77){
1568 _ndaugs =3;
1569 daugs[0]=EvtPDL::getId(std::string("D0"));
1570 daugs[1]=EvtPDL::getId(std::string("D*-"));
1571 daugs[2]=EvtPDL::getId(std::string("pi+"));
1572 }else if(mode == 78){
1573 _ndaugs =3;
1574 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1575 daugs[1]=EvtPDL::getId(std::string("D*+"));
1576 daugs[2]=EvtPDL::getId(std::string("pi-"));
1577 }else if(mode == 79){
1578 _ndaugs =3;
1579 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1580 daugs[1]=EvtPDL::getId(std::string("pi0"));
1581 daugs[2]=EvtPDL::getId(std::string("pi0"));
1582 }else if(mode == 80){
1583 _ndaugs =2;
1584 daugs[0]=EvtPDL::getId(std::string("eta"));
1585 daugs[1]=EvtPDL::getId(std::string("J/psi"));
1586 }else if(mode == 81){
1587 _ndaugs =3;
1588 daugs[0]=EvtPDL::getId(std::string("h_c"));
1589 daugs[1]=EvtPDL::getId(std::string("pi+"));
1590 daugs[2]=EvtPDL::getId(std::string("pi-"));
1591 }else if(mode == 82){
1592 _ndaugs =3;
1593 daugs[0]=EvtPDL::getId(std::string("h_c"));
1594 daugs[1]=EvtPDL::getId(std::string("pi0"));
1595 daugs[2]=EvtPDL::getId(std::string("pi0"));
1596 }else if(mode == 83){
1597 _ndaugs =3;
1598 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1599 daugs[1]=EvtPDL::getId(std::string("K+"));
1600 daugs[2]=EvtPDL::getId(std::string("K-"));
1601 }else if(mode == 84){
1602 _ndaugs =3;
1603 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1604 daugs[1]=EvtPDL::getId(std::string("K_S0"));
1605 daugs[2]=EvtPDL::getId(std::string("K_S0"));
1606 }else if(mode == 90){
1607 _ndaugs =3;
1608 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1609 daugs[1]=EvtPDL::getId(std::string("pi+"));
1610 daugs[2]=EvtPDL::getId(std::string("pi-"));
1611 }else if(mode == 91){
1612 _ndaugs =3;
1613 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1614 daugs[1]=EvtPDL::getId(std::string("pi+"));
1615 daugs[2]=EvtPDL::getId(std::string("pi-"));
1616 }else if(mode == 92){
1617 _ndaugs =3;
1618 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1619 daugs[1]=EvtPDL::getId(std::string("K+"));
1620 daugs[2]=EvtPDL::getId(std::string("K-"));
1621 }else if(mode == 93){
1622 _ndaugs =2;
1623 daugs[0]=EvtPDL::getId(std::string("D_s+"));
1624 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1625 }else if(mode == 94){
1626 _ndaugs =2;
1627 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
1628 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1629 }else if(mode == 95){
1630 _ndaugs =2;
1631 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
1632 daugs[1]=EvtPDL::getId(std::string("D_s+"));
1633 }else if(mode == 96){
1634 _ndaugs =2;
1635 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
1636 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
1637 }else if(mode == 74100){
1638 _ndaugs =2;
1639 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1640 daugs[1]=EvtPDL::getId(std::string("gamma"));
1641 }else if(mode == 74101){
1642 _ndaugs =2;
1643 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1644 daugs[1]=EvtPDL::getId(std::string("gamma"));
1645 }else if(mode == 74102){
1646 _ndaugs =2;
1647 daugs[0]=EvtPDL::getId(std::string("psi(3770)"));
1648 daugs[1]=EvtPDL::getId(std::string("gamma"));
1649 }else if(mode == 74103){
1650 _ndaugs =2;
1651 daugs[0]=EvtPDL::getId(std::string("phi"));
1652 daugs[1]=EvtPDL::getId(std::string("gamma"));
1653 }else if(mode == 74104){
1654 _ndaugs =2;
1655 daugs[0]=EvtPDL::getId(std::string("omega"));
1656 daugs[1]=EvtPDL::getId(std::string("gamma"));
1657 }else if(mode == 74105){
1658 _ndaugs =2;
1659 daugs[0]=EvtPDL::getId(std::string("rho0"));
1660 daugs[1]=EvtPDL::getId(std::string("gamma"));
1661 }else if(mode == 74106){
1662 _ndaugs =2;
1663 daugs[0]=EvtPDL::getId(std::string("rho(3S)0"));
1664 daugs[1]=EvtPDL::getId(std::string("gamma"));
1665 }else if(mode == 74107){
1666 _ndaugs =2;
1667 daugs[0]=EvtPDL::getId(std::string("omega(2S)"));
1668 daugs[1]=EvtPDL::getId(std::string("gamma"));
1669 }else if(mode == 74110){
1670 _ndaugs =2;
1671 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1672 daugs[1]=EvtPDL::getId(std::string("gamma"));
1673 EvtId myvpho = EvtPDL::getId(std::string("gamma*"));
1674 EvtPDL::reSetMass(myvpho,mhds*0.9); //for calculating maxXS, it will be resent when this mode is selected
1675 _modeFlag.clear();
1676 _modeFlag.push_back(74110); //R-value generator tag
1677 _modeFlag.push_back(0); //to contain the mode selected
1678 }else if(mode == -1){
1679 _ndaugs = nson;
1680 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1681 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
1682 }else if(mode == -2){
1683 _ndaugs = nson;
1684 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1685 }else if(mode == -100){
1686 _ndaugs =1;
1687 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1688 _modeFlag.clear();
1689 _modeFlag.push_back(-100); //R-value generator tag
1690 _modeFlag.push_back(0); //to contain the mode selected
1691
1692 }else if(mode == 10000){//use for check ISR
1693 _ndaugs =2;
1694 daugs[0]=EvtPDL::getId(std::string("pi+"));
1695 daugs[1]=EvtPDL::getId(std::string("pi-"));
1696 }else {
1697 std::cout<<"Bad mode_index number " <<mode<<", please refer to the manual."<<std::endl;
1698 ::abort();
1699 }
1700 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1701
1702 if(VISR){
1703 _ndaugs=2;
1704 daugs[0]=getDaugs()[0];//EvtPDL::getId(std::string("gamma"));
1705 daugs[1]=getDaugs()[1];//EvtPDL::getId(std::string("gamma*"));
1706 }
1707
1708 std::vector<EvtId> theDaugs;
1709 for(int i=0;i<_ndaugs;i++){
1710 theDaugs.push_back(daugs[i]);
1711 }
1712 if(theDaugs.size()) {return theDaugs;} else {std::cout<<"EvtConExc::get_mode: zero size"<<std::endl;abort();}
1713}
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65

Referenced by selectMode().

◆ get_mode_index()

int EvtConExc::get_mode_index ( int  mode)

Definition at line 3401 of file EvtConExc.cc.

3401 {
3402 int i=0;
3403 for(int j=0;i<vmode.size();j++){
3404 if(mode==vmode[j]) return j;
3405 }
3406 std::cout<<" EvtConExc::get_mode_index: no index is found in vmode"<<std::endl;
3407 abort();
3408}

Referenced by getObsXsection().

◆ getDaugId()

EvtId * EvtConExc::getDaugId ( )
inline

Definition at line 159 of file EvtConExc.hh.

159{return daugs;}

Referenced by EvtRexc::decay().

◆ getModeIndex()

int EvtConExc::getModeIndex ( int  m)

Definition at line 3541 of file EvtConExc.cc.

3541 {
3542 for (int i=0;i<vmode.size();i++){
3543 if(m==vmode[i]) return i;
3544 }
3545 std::cout<<"EvtConExc::getModeIndex: no index in vmode found "<<std::endl;
3546 abort();
3547}

Referenced by checkEvtRatio(), and writeDecayTabel().

◆ getName()

void EvtConExc::getName ( std::string &  name)
virtual

Implements EvtDecayBase.

Definition at line 217 of file EvtConExc.cc.

217 {
218
219 model_name="ConExc";
220
221}

◆ getNdaugs()

int EvtConExc::getNdaugs ( )
inline

Definition at line 158 of file EvtConExc.hh.

158{return _ndaugs;}

Referenced by EvtRexc::decay().

◆ getObsXsection()

double EvtConExc::getObsXsection ( double  mhds,
int  mode 
)

Definition at line 3410 of file EvtConExc.cc.

3410 {
3411 double hwid=(AA[0]-AA[1])/2.;
3412 double s=_cms*_cms;
3413 int idx=get_mode_index(mode);
3414 double Egam=(s-mhds*mhds)/2/_cms;
3415 for(int i=0;i<600;i++){
3416 if( (AA[i]-hwid)<=Egam && Egam<(AA[i]+hwid) ) return VXS[idx][i];
3417 }
3418
3419 std:cout<<"EvtConExc::getObsXsection : no observed xs is found "<<endl;
3420 abort();
3421}
EvtConExc::get_mode_index
int get_mode_index(int mode)
Definition: EvtConExc.cc:3401

◆ getResMass()

void EvtConExc::getResMass ( )

Definition at line 3076 of file EvtConExc.cc.

3076 {
3077 EvtId myres = EvtPDL::getId(std::string("J/psi"));
3078 mjsi = EvtPDL::getMeanMass(myres);
3079 wjsi = EvtPDL::getWidth(myres);
3080
3081 myres = EvtPDL::getId(std::string("psi(2S)"));
3082 mpsip= EvtPDL::getMeanMass(myres);
3083 wpsip= EvtPDL::getWidth(myres);
3084
3085 myres = EvtPDL::getId(std::string("psi(3770)"));
3086 mpsipp= EvtPDL::getMeanMass(myres);
3087 wpsipp= EvtPDL::getWidth(myres);
3088
3089 myres = EvtPDL::getId(std::string("phi"));
3090 mphi = EvtPDL::getMeanMass(myres);
3091 wphi = EvtPDL::getWidth(myres);
3092
3093 myres = EvtPDL::getId(std::string("omega"));
3094 momega= EvtPDL::getMeanMass(myres);
3095 womega= EvtPDL::getWidth(myres);
3096
3097 myres = EvtPDL::getId(std::string("rho0"));
3098 mrho0 = EvtPDL::getMeanMass(myres);
3099 wrho0 = EvtPDL::getWidth(myres);
3100
3101 myres = EvtPDL::getId(std::string("rho(3S)0"));
3102 mrho3s= EvtPDL::getMeanMass(myres);
3103 wrho3s= EvtPDL::getWidth(myres);
3104
3105
3106 myres = EvtPDL::getId(std::string("omega(2S)"));
3107 momega2s= EvtPDL::getMeanMass(myres);
3108 womega2s= EvtPDL::getWidth(myres);
3109
3110}
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition: EvtPDL.hh:54

Referenced by init().

◆ getSelectedMode()

int EvtConExc::getSelectedMode ( )
inline

Definition at line 160 of file EvtConExc.hh.

160{return _selectedMode;}

◆ getVP()

double EvtConExc::getVP ( double  cms)

Definition at line 3318 of file EvtConExc.cc.

3318 {
3319 double xx,r1,i1;
3320 double x1,y1,y2;
3321 xx=0;
3322 int mytag=1;
3323 for(int i=0;i<4001;i++){
3324 x1=vpx[i];
3325 y1=vpr[i];
3326 y2=vpi[i];
3327 if(xx<=mx && mx <x1){xx=x1; r1=y1; i1=y2; mytag=0; break;}
3328 xx=x1; r1=y1; i1=y2;
3329 }
3330 if(mytag==1){std::cout<<"No vacuum polarization value found"<<std::endl;abort();}
3331 EvtComplex cvp(r1,i1);
3332 double thevp=abs(1./(1-cvp));
3333 if(3.0933<mx && mx<3.1035) return 1.; //J/psi peak excluded
3334 if(3.6810<mx && mx<3.6913) return 1.; //psi(2S) peak excluded
3335 return thevp*thevp;
3336}

Referenced by init(), Rad2(), and SoftPhoton_xs().

◆ hadron_angle_sampling()

bool EvtConExc::hadron_angle_sampling ( EvtVector4R  ppi,
EvtVector4R  pcm 
)

Definition at line 2217 of file EvtConExc.cc.

2217 {
2218 EvtVector4R pbst=-1*pcm;
2219 pbst.set(0,pcm.get(0));
2220 EvtVector4R p4i = boostTo(ppi,pbst);
2221 if( (_mode>=0 && _mode<=5) || _mode ==44|| _mode ==47|| _mode ==48 ||_mode ==72 ||_mode == 73||_mode==94||_mode==95 ||_mode ==96 || _mode>=23 &&_mode<=27 || _mode==36|| _mode ==39 || _mode == 80 ){//ee->two baryon mode, VP, SP, mode=-2 is excluded
2222 bool byn_ang = baryon_sampling(pcm, p4i);
2223 return byn_ang;
2224 }else if(_mode ==6||_mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
2225 bool msn_ang = meson_sampling(pcm,p4i);
2226 return msn_ang;
2227 }else if(_mode==23 || _mode==24 || _mode==25 || _mode==26 || _mode==27 || _mode==36||_mode==47||_mode==48|| _mode==52||_mode==68||_mode==69||_mode==72||_mode==73||_mode==80||_mode==94||_mode==95){
2228 bool msn_ang = VP_sampling(pcm,p4i);
2229 return msn_ang;
2230 }else if(_mode==-2 ){
2231 EvtSpinType::spintype type0 = EvtPDL::getSpinType(daugs[0]);
2232 EvtSpinType::spintype type1 = EvtPDL::getSpinType(daugs[1]);
2233 if(type0==EvtSpinType::SCALAR &&type1==EvtSpinType::SCALAR ){
2234 bool msn_ang = meson_sampling(pcm,p4i);
2235 return msn_ang;
2236 }else if(type0==EvtSpinType::VECTOR &&type1==EvtSpinType::SCALAR || type1==EvtSpinType::VECTOR &&type0==EvtSpinType::SCALAR ){
2237 bool msn_ang = VP_sampling(pcm,p4i);
2238 return msn_ang;
2239 }
2240 }
2241 return true;
2242}
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition: EvtDiracSpinor.cc:75
EvtConExc::VP_sampling
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2516
EvtConExc::meson_sampling
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2505
EvtConExc::baryon_sampling
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2488
EvtPDL::getSpinType
static EvtSpinType::spintype getSpinType(EvtId i)
Definition: EvtPDL.hh:61
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:183

Referenced by decay().

◆ init()

void EvtConExc::init ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 230 of file EvtConExc.cc.

230 {
231 myFisr.clear();
232 ReadVP();
233 VXS.resize(120);
234 for(int i=0;i<120;i++){
235 VXS[i].resize(600);
236 }
237 _mode = getArg(0);
238 for(int i=0;i<97;i++){
239 if(53<=i && i<=58) continue;
240 if(85<=i && i<=89) continue;
241 if(_mode==74110 ||_mode==-100) vmode.push_back(i);
242 }
243 if(_mode==74110||_mode==-100){
244 vmode.push_back(74100);
245 vmode.push_back(74101);
246 vmode.push_back(74102);
247 vmode.push_back(74103);
248 vmode.push_back(74104);
249 vmode.push_back(74105);
250 vmode.push_back(74110);
251 }
252
253 std::cout<<"==== Parameters set by users====="<<std::endl;
254 for(int i=0;i<ncommand;i++){
255 std::cout<<commands[i]<<std::endl;
256 }
257 std::cout<<"==================================="<<std::endl;
258 InitPars();
259 std::cout<<threshold<<" "<<beamEnergySpread<<std::endl;
260 //----------------
261 // check that there are 0 arguments
262 // checkNArg(1);
263 //Vector ISR process
264 VISR=0;
265 if(getNDaug()==2){
266 if(getDaugs()[0]==EvtPDL::getId("gamma") &&getDaugs()[1]==EvtPDL::getId("gamma*") || getDaugs()[0]==EvtPDL::getId("gamma*") &&getDaugs()[1]==EvtPDL::getId("gamma") ) VISR=1;
267 }else if(getNDaug()>2){std::cout<<"Bad daughter specified"<<std::endl;abort();}
268
269 // cmsspread=0.0015; //CMC energy spread
270 testflag=0;
271 getResMass();
272 if(getArg(0) == -1){
273 radflag=0;mydbg=false;
274 _mode = getArg(0);
275 pdgcode = getArg(1);
276 pid = EvtPDL::evtIdFromStdHep(pdgcode );
277 nson = getNArg()-2;
278 std::cout<<"The decay daughters are "<<std::endl;
279 for(int i=2;i<getNArg();i++ ){son[i-2]=EvtPDL::evtIdFromStdHep(getArg(i));std::cout<<EvtPDL::name(son[i-2])<<" ";}
280 std::cout<<std::endl;
281 }else if(getArg(0)==-2){
282 radflag=0;mydbg=false;
283 _mode = getArg(0);
284 nson = getNArg()-1;
285 for(int i=1;i<getNArg();i++ ){son[i-1]=EvtPDL::evtIdFromStdHep(getArg(i));}
286 }else if(getArg(0)==-100){
287 _mode = getArg(0);
288 multimode=1;
289 vmd.clear();
290 for(int i=1;i<getNArg();i++){vmd.push_back(getArg(i));}
291 std::cout<<" multi-exclusive mode "<<std::endl;
292 for(int i=1;i<getNArg();i++){std::cout<<getArg(i)<<" ";}
293 std::cout<<std::endl;
294 if(vmd[vmd.size()-1]==74110){std::cout<<"74110 is not allowd for multimode simulation"<<std::endl;abort();}
295 }else if(getNArg()==1){ _mode = getArg(0); radflag=0;mydbg=false;}
296 else if(getNArg()==2){ _mode = getArg(0); radflag=getArg(1);mydbg=false;}
297 else if(getNArg()==3){ _mode = getArg(0); radflag=getArg(1);mydbg=getArg(2);}
298 else{std::cout<<"ConExc:number of parameter should be 1,2 or 3, but you set "<<getNArg()<<std::endl;::abort(); }
299 //--- debugging
300 //std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<" "<<EvtPDL::name(son[1])<<std::endl;
301
302 gamId = EvtPDL::getId(std::string("gamma"));
303 init_mode(_mode);
304 XS_max=-1;
305 //-----debugging to make root file
306 if(mydbg){
307 myfile = new TFile("mydbg.root","RECREATE");
308 xs=new TTree ("xs","momentum of rad. gamma and hadrons");
309 xs->Branch("imode",&imode,"imode/I");
310 xs->Branch("pgam",pgam,"pgam[4]/D");
311 xs->Branch("phds",phds,"phds[4]/D");
312 xs->Branch("ph1",ph1,"ph1[4]/D");
313 xs->Branch("ph2",ph2,"ph2[4]/D");
314 xs->Branch("mhds",&mhds,"mhds/D");
315 xs->Branch("mass1",&mass1,"mass1/D");
316 xs->Branch("mass2",&mass2,"mass2/D");
317 xs->Branch("costheta",&costheta,"costheta/D");
318 xs->Branch("cosp",&cosp,"cosp/D");
319 xs->Branch("cosk",&cosk,"cosk/D");
320 xs->Branch("sumxs",&sumxs,"sumxs/D");
321 xs->Branch("selectmode",&selectmode,"selectmode/D");
322 }
323 //--- prepare the output information
324 EvtId parentId =EvtPDL::getId(std::string("vpho"));
325 EvtPDL::reSetWidth(parentId,0.0);
326 double parentMass = EvtPDL::getMass(parentId);
327 std::cout<<"parent mass = "<<parentMass<<std::endl;
328
329
330 EvtVector4R p_init(parentMass,0.0,0.0,0.0);
331 EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
332 theparent = p;
333 myxsection =new EvtXsection (_mode);
334 staxsection =new EvtXsection (_mode);
335 if(_mode==-100) {
336 myxsection->setModes(vmd);
337 myxsection->ini_data_multimode();
338 staxsection->setModes(vmd);
339 staxsection->ini_data_multimode();
340 }
341 double mth=0;
342 double mx = EvtPDL::getMeanMass(parentId); //p->getP4().mass();
343 if(_mode ==-1){
344 myxsection->setBW(pdgcode);
345 staxsection->setBW(pdgcode);
346 for(int i=0;i<nson;i++){mth +=EvtPDL::getMeanMass(son[i]);}
347 if(mx<mth){std::cout<<"The CMS energy is lower than the threshold of the final states"<<std::endl;abort();}
348 }else if(_mode == -2){
349 mth=myxsection->getXlw();
350 }else{ mth=myxsection->getXlw();}
351 _cms = mx;
352 _unit = myxsection->getUnit();
353
354 std::cout<<"The specified mode is "<<_mode<<std::endl;
355 EvtConExc::getMode = _mode;
356 //std::cout<<"getXlw= "<<mth<<std::endl;
357
358 Egamcut = 0.001; //set photon energy cut according to the BESIII detector
359 double Esig = 0.0001; //sigularity engy
360 double x = 2*Egamcut/_cms;
361 double s = _cms*_cms;
362 double mhscut = sqrt(s*(1-x));
363
364 //for vaccum polarization
365 float fe,fst2,fder,ferrder,fdeg,ferrdeg;
366 fe=_cms;
367 //using the vacuum pol. value as given by http://www-com.physik.hu-berlin.de/~fjeger/software.html
368 vph=getVP(_cms);
369 if(3.0943<_cms && _cms<3.102) vph=1;// For J/psi, the vacuum factor is included in the resonance
370 std::cout<<"sum of leptonic, hadronic contributions(u,d,s,c,b) to vacuum polarization: "<<vph<<" @"<<fe<<"GeV"<<std::endl;
371
372 //caculate the xs for ISR to psiprime, J/psi and phi narrow resonance.
373 double totxsIRS=0;
374 init_Br_ee();
375 double mthrld = EvtPDL::getMeanMass(daugs[0]); //threshold mass of hadrons
376 for(int jj=1;jj<_ndaugs;jj++){
377 mthrld += EvtPDL::getMeanMass(daugs[jj]);
378 }
379
380 ISRXS.clear();ISRM.clear();ISRFLAG.clear();
381 for(int ii=0;ii<3;ii++){
382 double mR = EvtPDL::getMeanMass(ResId[ii]);
383 if(mx<mR || _mode !=74110) continue;
384 double narRxs=Ros_xs(mx,BR_ee[ii],ResId[ii]);
385 std::cout<<"The corss section for gamma "<<EvtPDL::name(ResId[ii]).c_str()<<" is: "<< narRxs<<"*Br (nb)."<<std::endl;
386 ISRXS.push_back(narRxs);
387 ISRM.push_back(mR);
388 ISRFLAG.push_back(1.);
389 ISRID.push_back(ResId[ii]);
390 totxsIRS += narRxs;
391 }
392 std::cout<<"==========================================================================="<<std::endl;
393
394 //-- calculated with MC integration method
395 double mhdL=myxsection->getXlw();
396 if(mhdL<SetMthr) mhdL=SetMthr;
397 EgamH = (s-mhdL*mhdL)/2/sqrt(s);
398 int nmc=1000000;
399 _xs0 = gamHXSection(s,Esig,Egamcut,nmc);
400 _er0 = _er1; // stored when calculate _xs0
401 std::cout<<"_er0= "<<_er0<<std::endl;
402 _xs1 = gamHXSection(s,Egamcut,EgamH,nmc);
403 double xs1_err = _er1;
404
405
406 //--- sigularity point x from 0 to 0.0001GeV
407 double xb= 2*Esig/_cms;
408 double sig_xs = SoftPhoton_xs(s,xb)*(myxsection->getXsection(mx));
409 _xs0 += sig_xs;
410
411 //prepare the observed cross section with interpotation function divdif in CERNLIB
412 testflag=1;
413 int Nstp=598;
414 double stp=(EgamH - Egamcut)/Nstp;
415 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
416 double Eg0=EgamH - i*stp;
417 double Eg1=EgamH - (i+1)*stp;
418 int nmc=20000;
419 int nint=100;
420 //double xsi= gamHXSection(s,Eg1,Eg0,nmc);
421 double xsi= trapezoid(s,Eg1,Eg0,nint,myxsection);
422 AA[i]=(Eg1+Eg0)/2;
423 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
424 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
425 double binwidth = mh2-mhi;
426 //if(_mode==74110) xsi += addNarrowRXS(mhi,binwidth);
427 if(i==0) AF[0]=xsi;
428 if(i>=1) AF[i]=AF[i-1]+xsi;
429 RadXS[i]=xsi/stp;
430 }
431 //add the interval 0~Esig
432 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
433 AF[598]=AF[597];
434 AF[599]=AF[598]+ _xs0;
435 RadXS[599]=_xs0;
436
437
438 //prepare observed cross section for each mode
439 for(int i=0;i<vmode.size();i++){
440 //std::cout<<"vmode index = "<<i<<std::endl;
441 if(_mode==74110||_mode==-100) mk_VXS(Esig,Egamcut,EgamH,i);
442 }
443 if(_mode==74110||_mode==-100) writeDecayTabel();
444 //after mk_VXS, restore the myxsection to _mode
445 if(_mode !=-100){
446 delete myxsection;
447 myxsection = new EvtXsection (_mode);
448 }
449 //debugging VXS
450 /*
451 double xtmp=0;
452 for(int i=0;i<vmode.size();i++){
453 xtmp+=VXS[i][599];
454 for(int j=0;j<600;j++){
455 std::cout<<VXS[i][j]<<" ";
456 }
457 std::cout<<std::endl;
458 }
459 */
460 //output observed cross section for the gamma Resonance mode
461 //std::cout<<"vmode.size======="<<vmode.size()<<std::endl;
462 for(int i=0;i<vmode.size();i++){
463 int md=vmode[i];
464 if(md<74100 || md>74106) continue;
465 std::string partname="";
466 if(md==74100) {partname="J/psi";}
467 else if(md==74101) {partname="psi(2S)";}
468 else if(md==74102) {partname="psi(3770)";}
469 else if(md==74103) {partname="phi";}
470 else if(md==74104) {partname="omega";}
471 else if(md==74105) {partname="rho0";}
472 else if(md==74106) {partname="rho(3S)0";}
473 else{;}
474 std::cout<<"The observed cross section for gamma "<<partname<<": "<<VXS[i][599]<<" nb"<<std::endl;
475 }
476 //--
477
478 for(int i=0;i<600;i++){
479 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
480 //std::cout<<i<<" Egamma "<<AA[i]<<" Mhadons "<<MH[i]<<std::endl;
481 }
482 //for debugging
483 //for(int i=0;i<600;i++){double s=_cms*_cms; double mhds=sqrt(s-2*_cms*AA[i]);std::cout<<"Mhds="<<mhds<<" Egam="<<AA[i]<<" "<<AF[i]<<std::endl;}
484 std::cout<<"VXS[79][599]: "<<VXS[79][599]<<" VXS[79][598]= "<<VXS[79][598]<<std::endl;
485 std::cout<<"EgamH="<<EgamH<<" "<<_xs0+_xs1<<" AF[599]="<<AF[599]<<" AF[598] "<<AF[598]<<std::endl;
486 //double Etst=0.0241838;
487 //std::cout<<"Etst="<<Etst<<" "<<gamHXSection(s,Etst,EgamH,nmc)<<" "<< lgr(AA,AF,600,Etst)<<std::endl; abort();
488
489 //for information output
490 double bornXS = myxsection->getXsection(mx); // std::cout<<"BornXS= "<<bornXS<<std::endl;
491 double bornXS_er=myxsection->getErr(mx);
492 double obsvXS = _xs0 + _xs1; //gamHXSection(mth ,mx);
493 double obsvXS_er= _er0 + xs1_err;
494 double corr = obsvXS/bornXS;
495 double corr_er =corr*sqrt(bornXS_er*bornXS_er/bornXS/bornXS + obsvXS_er*obsvXS_er/obsvXS/obsvXS);
496
497
498 if(bornXS==0){std::cout<<"EvtConExc::init : The Born cross section at this point is zero!"<<std::endl;abort();}
499 if(obsvXS==0){std::cout<<"EvtConExc::init : The observed cross section at this point is zero!"<<std::endl;abort();}
500
501 //_xs0 += bornXS;// events with very soft photon (Egam<0.025) are taken as the born process
502 //_er0 = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
503
504 std::cout<<""<<std::endl;
505 std::cout<<"========= Generation using cross section (Egamcut="<<Egamcut<<" GeV) =============="<<std::endl;
506 std::cout<<" sqrt(s)= "<<mx<< " GeV "<<std::endl;
507 if(_mode>=0 || _mode ==-2 || _mode ==-100 ){
508 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<"+/-"<<_er0<<" in Unit "<<_unit.c_str()<<std::endl;
509 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"+/-"<<xs1_err<<" in Unit "<<_unit.c_str()<<std::endl;
510 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
511 std::cout<<"The radiative correction factor f_ISR*f_vacuum= sigma_obs/sigma_born(s) = "<<corr<<"+/-"<< fabs(corr_er)<<","<<std::endl;
512 std::cout<<"which is calcualted with these quantities:"<<std::endl;
513 std::cout<<"the observed cross section is "<<obsvXS<<"+/-"<< obsvXS_er<<_unit.c_str()<<std::endl;
514 std::cout<<"and the Born cross section (s) is "<<bornXS<<"+/-"<<bornXS_er<<_unit.c_str()<<std::endl;
515 std::cout<<"and the vaccum polarziation factor (lep. + hadr.) 1/|1-Pi|^2= "<<vph<<std::endl;
516 std::cout<<"Within the range "<<myxsection->getXlw()<<" ~"<<myxsection->getXup()<<" GeV, "<<myxsection->getMsg().c_str()<<std::endl;
517 std::cout<<"==========================================================================="<<std::endl;
518 }else if(_mode==-1){
519 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<" *Br_ee"<<" in Unit "<<_unit.c_str()<<std::endl;
520 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"*Br_ee in Unit "<<_unit.c_str()<<std::endl;
521 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
522 std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<std::endl;
523 std::cout<<"==========================================================================="<<std::endl;
524 }
525 std::cout<<std::endl;
526 std::cout<<std::endl;
527
528 findMaxXS(p);
529 _mhdL=myxsection->getXlw();
530 //--debugging
531 //std::cout<<"maxXS= "<<maxXS<<std::endl;
532
533 if(_xs0 == 0 && _xs1==0){std::cout<<"EvtConExc::ini() has zero cross section"<<std::endl;abort();}
534
535 std::cout<<std::endl;
536 std::cout<<std::endl;
537 mlow=myxsection->getXlw();
538 mup =myxsection->getXup();
539 //--- for debugging
540 if(mydbg && _mode!=74110){
541 int nd = myxsection->getYY().size();
542 double xx[10000],yy[10000],xer[10000],yer[10000];
543 for(int i=0;i<nd;i++){
544 xx[i]=myxsection->getXX()[i];
545 yy[i]=myxsection->getYY()[i];
546 yer[i]=myxsection->getEr()[i];
547 xer[i]=0.;
548 //std::cout<<"yy "<<yy[i]<<std::endl;
549 }
550 myth=new TH1F("myth","Exp.data",200,xx[0],xx[nd]);
551 for(int i=0;i<nd;i++){
552 myth->Fill(xx[i],yy[i]);
553 }
554 myth->SetError(yer);
555 myth->Write();
556
557 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
558 }else
559
560 if(mydbg && _mode==74110 ){
561 int nd = myxsection->getYY().size();
562 double xx[10000],yy[10000],xer[10000],yer[10000],ysum[10000],yysum[10000];
563 for(int i=0;i<nd;i++){
564 xx[i]=myxsection->getXX()[i];
565 yy[i]=myxsection->getYY()[i];
566 yer[i]=myxsection->getEr()[i];
567 xer[i]=0.;
568 //std::cout<<"yy "<<yy[i]<<std::endl;
569 }
570#include "sty.C"
571 double xhigh=5.0;
572 double xlow = myxsection->getXlw();
573 TGraphErrors *Gdt = new TGraphErrors(nd,xx,yy,xer,yer);
574
575 myth=new TH1F("myth","Exp.data",600,xlow,xhigh);
576 Xsum=new TH1F("sumxs","sum of exclusive xs",600,xlow,xhigh);
577 double mstp=(xhigh-xlow)/600;
578 for(int i=0;i<600;i++){
579 double mx=i*mstp+xlow;
580 double xsi = myxsection->getXsection(mx);
581 if(xsi==0) continue;
582 myth->Fill(mx,xsi);
583 //std::cout<<mx<<" "<<xsi<<std::endl;
584 }
585
586 for(int i=0;i<600;i++){
587 double mx=i*mstp+xlow;
588 ysum[i]=sumExc(mx);
589 if(ysum[i]==0) continue;
590 Xsum->Fill(mx,ysum[i]);
591 //std::cout<<mx<<" "<<ysum[i]<<std::endl;
592 }
593
594 for(int i=0;i<nd;i++){
595 yysum[i]=sumExc(xx[i]);
596 }
597
598 myth->SetError(yer);
599 myth->Write();
600 Xsum->Write();
601
602 TGraphErrors *Gsum = new TGraphErrors(nd,xx,yysum,xer,yer);
603 //draw graph
604 double ex[600]={600*0};
605 double ey[600],ta[600];
606 double exz[600]={600*0};
607 double eyz[600]={600*0};
608 for(int i=0;i<600;i++){
609 ey[i]=AF[i]*0.001;
610 }
611 TGraphErrors *gr0 = new TGraphErrors(600,MH,AF,ex,ey);
612 TGraphErrors *gr1 = new TGraphErrors(600,MH,RadXS,exz,eyz);
613 TPostScript *ps = new TPostScript("xsobs.ps",111);
614 TCanvas *c1 = new TCanvas("c1","TGraphs for data",200,10,600,400);
615 gPad-> SetLogy(1);
616 // c1->SetLogy(1);
617 gStyle->SetCanvasColor(0);
618 gStyle->SetStatBorderSize(1);
619 c1->Divide(1,1);
620
621 c1->Update();
622 ps->NewPage();
623 c1->Draw();
624 c1->cd(1);
625 c1->SetLogy(1);
626 gr0->SetMarkerStyle(10);
627 gr0->Draw("AP");
628 gr0->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
629 gr0->GetYaxis()->SetTitle("observed cross section (nb)");
630 gr0->Write();
631
632 c1->Update();
633 ps->NewPage();
634 c1->Draw();
635 c1->cd(1);
636 c1->SetLogy(1);
637 gr1->SetMarkerStyle(10);
638 gr1->Draw("AP");
639 gr1->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
640 gr1->GetYaxis()->SetTitle("Rad*BornXS");
641 gr1->Write();
642
643 //--------- imposing graphs to a pad
644 TMultiGraph *mg = new TMultiGraph();
645 mg->SetTitle("Exclusion graphs");
646 Gdt->SetMarkerStyle(2);
647 Gdt->SetMarkerSize(1);
648 Gsum->SetLineColor(2);
649 Gsum->SetLineWidth(1);
650 mg->Add(Gdt);
651 mg->Add(Gsum);
652
653 c1->Update();
654 ps->NewPage();
655 c1->Draw();
656 c1->cd(1);
657 gPad-> SetLogy(1);
658 mg->Draw("APL");
659 mg->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
660 mg->GetYaxis()->SetTitle("observed cross section (nb)");
661 mg->Write();
662 //-------
663
664 c1->Update();
665 ps->NewPage();
666 c1->Draw();
667 c1->cd(1);
668 Gdt->SetMarkerStyle(2);
669 Gdt->Draw("AP");
670 Gdt->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
671 Gdt->GetYaxis()->SetTitle("observed cross section (nb)");
672 Gdt->Write();
673
674 c1->Update();
675 ps->NewPage();
676 c1->Draw();
677 c1->cd(1);
678 Gsum->SetMarkerStyle(2);
679 Gsum->SetMarkerSize(1);
680 Gsum->Draw("AP");
681 Gsum->SetLineWidth(0);
682 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
683 Gsum->GetYaxis()->SetTitle("observed cross section (nb)");
684 Gsum->Write();
685
686 c1->Update();
687 ps->NewPage();
688 c1->Draw();
689 c1->cd(1);
690 // gPad-> SetLogx(1);
691 Gdt->SetMarkerStyle(2);
692 Gdt->SetMarkerSize(1);
693 Gdt->SetMaximum(8000.0);
694 Gsum->SetLineColor(2);
695 Gsum->SetLineWidth(1.5);
696 Gsum->Draw("AL"); //A: draw axis
697 Gdt->Draw("P"); // superpose to the Gsum, without A-option
698 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
699 Gsum->GetYaxis()->SetTitle("cross section (nb)");
700 Gsum->Write();
701
702 ps->Close();
703 } //if(mydbg)_block
704 //-------------------------
705
706}
xlow
DOUBLE_PRECISION xlow[20]
Definition: EvtPhokharaDef.hh:56
mstp
int mstp[200]
Definition: EvtPycont.cc:54
EvtConExc::getMode
static int getMode
Definition: EvtConExc.hh:144
EvtConExc::init_Br_ee
void init_Br_ee()
Definition: EvtConExc.cc:2645
EvtConExc::ReadVP
void ReadVP()
Definition: EvtConExc.cc:3339
EvtConExc::InitPars
void InitPars()
Definition: EvtConExc.cc:3611
EvtConExc::Ros_xs
double Ros_xs(double mx, double bree, EvtId pid)
Definition: EvtConExc.cc:2677
EvtConExc::multimode
static int multimode
Definition: EvtConExc.hh:161
EvtConExc::gamHXSection
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition: EvtConExc.cc:2245
EvtConExc::sumExc
double sumExc(double mx)
Definition: EvtConExc.cc:3137
EvtConExc::writeDecayTabel
void writeDecayTabel()
Definition: EvtConExc.cc:3429
EvtConExc::mlow
static double mlow
Definition: EvtConExc.hh:184
EvtConExc::mk_VXS
void mk_VXS(double Esig, double Egamcut, double EgamH, int midx)
Definition: EvtConExc.cc:3362
EvtConExc::getResMass
void getResMass()
Definition: EvtConExc.cc:3076
EvtConExc::getVP
double getVP(double cms)
Definition: EvtConExc.cc:3318
EvtConExc::XS_max
static double XS_max
Definition: EvtConExc.hh:141
EvtConExc::mup
static double mup
Definition: EvtConExc.hh:184
EvtPDL::evtIdFromStdHep
static EvtId evtIdFromStdHep(int stdhep)
Definition: EvtPDL.cc:244
EvtXsection::getMsg
std::string getMsg()
Definition: EvtXsection.hh:157
EvtXsection::setModes
void setModes(std::vector< int > vmd)
Definition: EvtXsection.cc:12136
EvtXsection::getEr
std::vector< double > getEr()
Definition: EvtXsection.hh:154
EvtXsection::ini_data_multimode
void ini_data_multimode()
Definition: EvtXsection.cc:11922
EvtXsection::getUnit
std::string getUnit()
Definition: EvtXsection.hh:150
EvtXsection::getYY
std::vector< double > getYY()
Definition: EvtXsection.hh:153
EvtXsection::setBW
void setBW(int pdg)
Definition: EvtXsection.cc:12116
EvtXsection::getXX
std::vector< double > getXX()
Definition: EvtXsection.hh:152

◆ init_Br_ee()

void EvtConExc::init_Br_ee ( )

Definition at line 2645 of file EvtConExc.cc.

2645 {
2646 // double psipp_ee=9.6E-06;
2647 double psip_ee =7.73E-03;
2648 double jsi_ee =5.94*0.01;
2649 double phi_ee =2.954E-04;
2650 // double omega_ee=7.28E-05;
2651 // double rho_ee = 4.72E-05;
2652 EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
2653 EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
2654 EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
2655 EvtId phiId=EvtPDL::getId(std::string("phi"));
2656 EvtId omegaId=EvtPDL::getId(std::string("omega"));
2657 EvtId rhoId=EvtPDL::getId(std::string("rho0"));
2658 BR_ee.clear();
2659 ResId.clear();
2660
2661 //BR_ee.push_back(rho_ee);
2662 //BR_ee.push_back(omega_ee);
2663 BR_ee.push_back(phi_ee);
2664 BR_ee.push_back(jsi_ee);
2665 BR_ee.push_back(psip_ee);
2666 //BR_ee.push_back(psipp_ee);
2667
2668 //ResId.push_back(rhoId);
2669 //ResId.push_back(omegaId);
2670 ResId.push_back(phiId);
2671 ResId.push_back(jsiId);
2672 ResId.push_back(pspId);
2673 //ResId.push_back(psppId);
2674
2675}

Referenced by init().

◆ init_mode()

void EvtConExc::init_mode ( int  mode)

Definition at line 710 of file EvtConExc.cc.

710 {
711 if(mode ==0 ){
712 _ndaugs =2;
713 daugs[0]=EvtPDL::getId(std::string("p+"));
714 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
715 }else if(mode ==1 ){
716 _ndaugs =2;
717 daugs[0]=EvtPDL::getId(std::string("n0"));
718 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
719 }else if(mode ==2 ){
720 _ndaugs =2;
721 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
722 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
723 }else if(mode ==3 ){
724 _ndaugs =2;
725 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
726 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
727 }else if(mode ==4 ){
728 _ndaugs =2;
729 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
730 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
731 }else if(mode ==5 ){
732 _ndaugs =2;
733 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
734 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
735 } else if(mode ==6 ){
736 _ndaugs =2;
737 daugs[0]=EvtPDL::getId(std::string("pi+"));
738 daugs[1]=EvtPDL::getId(std::string("pi-"));
739 }else if(mode ==7 ){
740 _ndaugs =3;
741 daugs[0]=EvtPDL::getId(std::string("pi+"));
742 daugs[1]=EvtPDL::getId(std::string("pi-"));
743 daugs[2]=EvtPDL::getId(std::string("pi0"));
744 }else if(mode ==8 ){
745 _ndaugs =3;
746 daugs[0]=EvtPDL::getId(std::string("K+"));
747 daugs[1]=EvtPDL::getId(std::string("K-"));
748 daugs[2]=EvtPDL::getId(std::string("pi0"));
749 }else if(mode ==9 ){
750 _ndaugs =3;
751 daugs[0]=EvtPDL::getId(std::string("K_S0"));
752 daugs[1]=EvtPDL::getId(std::string("K+"));
753 daugs[2]=EvtPDL::getId(std::string("pi-"));
754 }else if(mode ==10 ){
755 _ndaugs =3;
756 daugs[0]=EvtPDL::getId(std::string("K_S0"));
757 daugs[1]=EvtPDL::getId(std::string("K-"));
758 daugs[2]=EvtPDL::getId(std::string("pi+"));
759 }else if(mode ==11 ){
760 _ndaugs =3;
761 daugs[0]=EvtPDL::getId(std::string("K+"));
762 daugs[1]=EvtPDL::getId(std::string("K+"));
763 daugs[2]=EvtPDL::getId(std::string("eta"));
764 }else if(mode ==12 ){
765 _ndaugs =4;
766 daugs[0]=EvtPDL::getId(std::string("pi+"));
767 daugs[1]=EvtPDL::getId(std::string("pi-"));
768 daugs[2]=EvtPDL::getId(std::string("pi+"));
769 daugs[3]=EvtPDL::getId(std::string("pi-"));
770 }else if(mode ==13 ){
771 _ndaugs =4;
772 daugs[0]=EvtPDL::getId(std::string("pi+"));
773 daugs[1]=EvtPDL::getId(std::string("pi-"));
774 daugs[2]=EvtPDL::getId(std::string("pi0"));
775 daugs[3]=EvtPDL::getId(std::string("pi0"));
776 }else if(mode ==14 ){
777 _ndaugs =4;
778 daugs[0]=EvtPDL::getId(std::string("K+"));
779 daugs[1]=EvtPDL::getId(std::string("K-"));
780 daugs[2]=EvtPDL::getId(std::string("pi+"));
781 daugs[3]=EvtPDL::getId(std::string("pi-"));
782 }else if(mode ==15 ){
783 _ndaugs =4;
784 daugs[0]=EvtPDL::getId(std::string("K+"));
785 daugs[1]=EvtPDL::getId(std::string("K-"));
786 daugs[2]=EvtPDL::getId(std::string("pi0"));
787 daugs[3]=EvtPDL::getId(std::string("pi0"));
788 }else if(mode ==16 ){
789 _ndaugs =4;
790 daugs[0]=EvtPDL::getId(std::string("K+"));
791 daugs[1]=EvtPDL::getId(std::string("K-"));
792 daugs[2]=EvtPDL::getId(std::string("K+"));
793 daugs[3]=EvtPDL::getId(std::string("K-"));
794 }else if(mode ==17 ){
795 _ndaugs =5;
796 daugs[0]=EvtPDL::getId(std::string("pi+"));
797 daugs[1]=EvtPDL::getId(std::string("pi-"));
798 daugs[2]=EvtPDL::getId(std::string("pi+"));
799 daugs[3]=EvtPDL::getId(std::string("pi-"));
800 daugs[4]=EvtPDL::getId(std::string("pi0"));
801 }else if(mode ==18 ){
802 _ndaugs =5;
803 daugs[0]=EvtPDL::getId(std::string("pi+"));
804 daugs[1]=EvtPDL::getId(std::string("pi-"));
805 daugs[2]=EvtPDL::getId(std::string("pi+"));
806 daugs[3]=EvtPDL::getId(std::string("pi-"));
807 daugs[4]=EvtPDL::getId(std::string("eta"));
808 }else if(mode ==19 ){
809 _ndaugs =5;
810 daugs[0]=EvtPDL::getId(std::string("K+"));
811 daugs[1]=EvtPDL::getId(std::string("K-"));
812 daugs[2]=EvtPDL::getId(std::string("pi+"));
813 daugs[3]=EvtPDL::getId(std::string("pi-"));
814 daugs[4]=EvtPDL::getId(std::string("pi0"));
815 }else if(mode ==20 ){
816 _ndaugs =5;
817 daugs[0]=EvtPDL::getId(std::string("K+"));
818 daugs[1]=EvtPDL::getId(std::string("K-"));
819 daugs[2]=EvtPDL::getId(std::string("pi+"));
820 daugs[3]=EvtPDL::getId(std::string("pi-"));
821 daugs[4]=EvtPDL::getId(std::string("eta"));
822 }else if(mode ==21 ){
823 _ndaugs =6;
824 daugs[0]=EvtPDL::getId(std::string("pi+"));
825 daugs[1]=EvtPDL::getId(std::string("pi-"));
826 daugs[2]=EvtPDL::getId(std::string("pi+"));
827 daugs[3]=EvtPDL::getId(std::string("pi-"));
828 daugs[4]=EvtPDL::getId(std::string("pi+"));
829 daugs[5]=EvtPDL::getId(std::string("pi-"));
830 }else if(mode ==22 ){
831 _ndaugs =6;
832 daugs[0]=EvtPDL::getId(std::string("pi+"));
833 daugs[1]=EvtPDL::getId(std::string("pi-"));
834 daugs[2]=EvtPDL::getId(std::string("pi+"));
835 daugs[3]=EvtPDL::getId(std::string("pi-"));
836 daugs[4]=EvtPDL::getId(std::string("pi0"));
837 daugs[5]=EvtPDL::getId(std::string("pi0"));
838 }else if(mode == 23){
839 _ndaugs =2;
840 daugs[0]=EvtPDL::getId(std::string("eta"));
841 daugs[1]=EvtPDL::getId(std::string("phi"));
842 }else if(mode == 24){
843 _ndaugs =2;
844 daugs[0]=EvtPDL::getId(std::string("phi"));
845 daugs[1]=EvtPDL::getId(std::string("pi0"));
846 }else if(mode == 25){
847 _ndaugs =2;
848 daugs[0]=EvtPDL::getId(std::string("K+"));
849 daugs[1]=EvtPDL::getId(std::string("K*-"));
850 }else if(mode == 26){
851 _ndaugs =2;
852 daugs[0]=EvtPDL::getId(std::string("K-"));
853 daugs[1]=EvtPDL::getId(std::string("K*+"));
854 }else if(mode == 27){
855 _ndaugs =2;
856 daugs[0]=EvtPDL::getId(std::string("K_S0"));
857 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
858 }else if(mode == 28){
859 _ndaugs =3;
860 daugs[0]=EvtPDL::getId(std::string("K*0"));
861 daugs[1]=EvtPDL::getId(std::string("K+"));
862 daugs[2]=EvtPDL::getId(std::string("pi-"));
863 }else if(mode == 29){
864 _ndaugs =3;
865 daugs[0]=EvtPDL::getId(std::string("K*0"));
866 daugs[1]=EvtPDL::getId(std::string("K-"));
867 daugs[2]=EvtPDL::getId(std::string("pi+"));
868 }else if(mode == 30){
869 _ndaugs =3;
870 daugs[0]=EvtPDL::getId(std::string("K*+"));
871 daugs[1]=EvtPDL::getId(std::string("K-"));
872 daugs[2]=EvtPDL::getId(std::string("pi0"));
873 }else if(mode == 31){
874 _ndaugs =3;
875 daugs[0]=EvtPDL::getId(std::string("K*-"));
876 daugs[1]=EvtPDL::getId(std::string("K+"));
877 daugs[2]=EvtPDL::getId(std::string("pi0"));
878 }else if(mode == 32){
879 _ndaugs =3;
880 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
881 daugs[1]=EvtPDL::getId(std::string("K+"));
882 daugs[2]=EvtPDL::getId(std::string("pi-"));
883 }else if(mode == 33){
884 _ndaugs =3;
885 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
886 daugs[1]=EvtPDL::getId(std::string("K-"));
887 daugs[2]=EvtPDL::getId(std::string("pi+"));
888 }else if(mode == 34){
889 _ndaugs =3;
890 daugs[0]=EvtPDL::getId(std::string("K+"));
891 daugs[1]=EvtPDL::getId(std::string("K-"));
892 daugs[2]=EvtPDL::getId(std::string("rho0"));
893 }else if(mode == 35){
894 _ndaugs =3;
895 daugs[0]=EvtPDL::getId(std::string("phi"));
896 daugs[1]=EvtPDL::getId(std::string("pi-"));
897 daugs[2]=EvtPDL::getId(std::string("pi+"));
898 }else if(mode == 36){
899 _ndaugs =2;
900 daugs[0]=EvtPDL::getId(std::string("phi"));
901 daugs[1]=EvtPDL::getId(std::string("f_0"));
902 }else if(mode == 37){
903 _ndaugs =3;
904 daugs[0]=EvtPDL::getId(std::string("eta"));
905 daugs[1]=EvtPDL::getId(std::string("pi+"));
906 daugs[2]=EvtPDL::getId(std::string("pi-"));
907 }else if(mode == 38){
908 _ndaugs =3;
909 daugs[0]=EvtPDL::getId(std::string("omega"));
910 daugs[1]=EvtPDL::getId(std::string("pi+"));
911 daugs[2]=EvtPDL::getId(std::string("pi-"));
912 }else if(mode == 39){
913 _ndaugs =2;
914 daugs[0]=EvtPDL::getId(std::string("omega"));
915 daugs[1]=EvtPDL::getId(std::string("f_0"));
916 }else if(mode == 40){
917 _ndaugs =3;
918 daugs[0]=EvtPDL::getId(std::string("eta'"));
919 daugs[1]=EvtPDL::getId(std::string("pi+"));
920 daugs[2]=EvtPDL::getId(std::string("pi-"));
921 }else if(mode == 41){
922 _ndaugs =3;
923 daugs[0]=EvtPDL::getId(std::string("f_1"));
924 daugs[1]=EvtPDL::getId(std::string("pi+"));
925 daugs[2]=EvtPDL::getId(std::string("pi-"));
926 }else if(mode == 42){
927 _ndaugs =3;
928 daugs[0]=EvtPDL::getId(std::string("omega"));
929 daugs[1]=EvtPDL::getId(std::string("K+"));
930 daugs[2]=EvtPDL::getId(std::string("K-"));
931 }else if(mode == 43){
932 _ndaugs =4;
933 daugs[0]=EvtPDL::getId(std::string("omega"));
934 daugs[1]=EvtPDL::getId(std::string("pi+"));
935 daugs[2]=EvtPDL::getId(std::string("pi-"));
936 daugs[3]=EvtPDL::getId(std::string("pi0"));
937 }else if(mode == 44){
938 _ndaugs =2;
939 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
940 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
941 }else if(mode == 45){
942 _ndaugs =2;
943 daugs[0]=EvtPDL::getId(std::string("K+"));
944 daugs[1]=EvtPDL::getId(std::string("K-"));
945 }else if(mode == 46){
946 _ndaugs =2;
947 daugs[0]=EvtPDL::getId(std::string("K_S0"));
948 daugs[1]=EvtPDL::getId(std::string("K_L0"));
949 }else if(mode == 47){
950 _ndaugs =2;
951 daugs[0]=EvtPDL::getId(std::string("omega"));
952 daugs[1]=EvtPDL::getId(std::string("eta"));
953 }else if(mode == 48){
954 _ndaugs =2;
955 daugs[0]=EvtPDL::getId(std::string("phi"));
956 daugs[1]=EvtPDL::getId(std::string("eta'"));
957 }else if(mode == 49){
958 _ndaugs =3;
959 daugs[0]=EvtPDL::getId(std::string("phi"));
960 daugs[1]=EvtPDL::getId(std::string("K+"));
961 daugs[2]=EvtPDL::getId(std::string("K-"));
962 }else if(mode == 50){
963 _ndaugs =3;
964 daugs[0]=EvtPDL::getId(std::string("p+"));
965 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
966 daugs[2]=EvtPDL::getId(std::string("pi0"));
967 }else if(mode == 51){
968 _ndaugs =3;
969 daugs[0]=EvtPDL::getId(std::string("p+"));
970 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
971 daugs[2]=EvtPDL::getId(std::string("eta"));
972 }else if(mode == 52){
973 _ndaugs =2;
974 daugs[0]=EvtPDL::getId(std::string("omega"));
975 daugs[1]=EvtPDL::getId(std::string("pi0"));
976
977 }else if(mode == 59){
978 _ndaugs =3;
979 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
980 daugs[1]=EvtPDL::getId(std::string("D0"));
981 daugs[2]=EvtPDL::getId(std::string("pi0"));
982 }else if(mode == 60){
983 _ndaugs =3;
984 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
985 daugs[1]=EvtPDL::getId(std::string("D*0"));
986 daugs[2]=EvtPDL::getId(std::string("pi0"));
987 }else if(mode == 61){
988 _ndaugs =3;
989 daugs[0]=EvtPDL::getId(std::string("D0"));
990 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
991 daugs[2]=EvtPDL::getId(std::string("pi0"));
992 }else if(mode == 62){
993 _ndaugs =3;
994 daugs[0]=EvtPDL::getId(std::string("D+"));
995 daugs[1]=EvtPDL::getId(std::string("D*-"));
996 daugs[2]=EvtPDL::getId(std::string("pi0"));
997 }else if(mode == 63){
998 _ndaugs =3;
999 daugs[0]=EvtPDL::getId(std::string("D-"));
1000 daugs[1]=EvtPDL::getId(std::string("D+"));
1001 daugs[2]=EvtPDL::getId(std::string("pi0"));
1002 }else if(mode == 64){
1003 _ndaugs =3;
1004 daugs[0]=EvtPDL::getId(std::string("D-"));
1005 daugs[1]=EvtPDL::getId(std::string("D*+"));
1006 daugs[2]=EvtPDL::getId(std::string("pi0"));
1007 }else if(mode == 65){
1008 _ndaugs =3;
1009 daugs[0]=EvtPDL::getId(std::string("D-"));
1010 daugs[1]=EvtPDL::getId(std::string("D*0"));
1011 daugs[2]=EvtPDL::getId(std::string("pi+"));
1012 }else if(mode == 66){
1013 _ndaugs =3;
1014 daugs[0]=EvtPDL::getId(std::string("D+"));
1015 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1016 daugs[2]=EvtPDL::getId(std::string("pi-"));
1017 }else if(mode == 67){
1018 _ndaugs =2;
1019 daugs[0]=EvtPDL::getId(std::string("D*0"));
1020 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1021 }else if(mode == 68){
1022 _ndaugs =2;
1023 daugs[0]=EvtPDL::getId(std::string("D0"));
1024 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1025
1026 }else if(mode == 69){
1027 _ndaugs =2;
1028 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1029 daugs[1]=EvtPDL::getId(std::string("D*0"));
1030
1031 }else if(mode == 70){
1032 _ndaugs =2;
1033 daugs[0]=EvtPDL::getId(std::string("D0"));
1034 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
1035
1036}else if(mode == 71){
1037 _ndaugs =2;
1038 daugs[0]=EvtPDL::getId(std::string("D+"));
1039 daugs[1]=EvtPDL::getId(std::string("D-"));
1040 }else if(mode == 72){
1041 _ndaugs =2;
1042 daugs[0]=EvtPDL::getId(std::string("D+"));
1043 daugs[1]=EvtPDL::getId(std::string("D*-"));
1044
1045 }else if(mode == 73){
1046 _ndaugs =2;
1047 daugs[0]=EvtPDL::getId(std::string("D-"));
1048 daugs[1]=EvtPDL::getId(std::string("D*+"));
1049
1050 }else if(mode == 74){
1051 _ndaugs =2;
1052 daugs[0]=EvtPDL::getId(std::string("D*+"));
1053 daugs[1]=EvtPDL::getId(std::string("D*-"));
1054
1055 }else if(mode == 75){
1056 _ndaugs =3;
1057 daugs[0]=EvtPDL::getId(std::string("D0"));
1058 daugs[1]=EvtPDL::getId(std::string("D-"));
1059 daugs[2]=EvtPDL::getId(std::string("pi+"));
1060 }else if(mode == 76){
1061 _ndaugs =3;
1062 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1063 daugs[1]=EvtPDL::getId(std::string("D+"));
1064 daugs[2]=EvtPDL::getId(std::string("pi-"));
1065 }else if(mode == 77){
1066 _ndaugs =3;
1067 daugs[0]=EvtPDL::getId(std::string("D0"));
1068 daugs[1]=EvtPDL::getId(std::string("D*-"));
1069 daugs[2]=EvtPDL::getId(std::string("pi+"));
1070 }else if(mode == 78){
1071 _ndaugs =3;
1072 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1073 daugs[1]=EvtPDL::getId(std::string("D*+"));
1074 daugs[2]=EvtPDL::getId(std::string("pi-"));
1075 }else if(mode == 79){
1076 _ndaugs =3;
1077 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1078 daugs[1]=EvtPDL::getId(std::string("pi0"));
1079 daugs[2]=EvtPDL::getId(std::string("pi0"));
1080 }else if(mode == 80){
1081 _ndaugs =2;
1082 daugs[0]=EvtPDL::getId(std::string("eta"));
1083 daugs[1]=EvtPDL::getId(std::string("J/psi"));
1084 }else if(mode == 81){
1085 _ndaugs =3;
1086 daugs[0]=EvtPDL::getId(std::string("h_c"));
1087 daugs[1]=EvtPDL::getId(std::string("pi+"));
1088 daugs[2]=EvtPDL::getId(std::string("pi-"));
1089 }else if(mode == 82){
1090 _ndaugs =3;
1091 daugs[0]=EvtPDL::getId(std::string("h_c"));
1092 daugs[1]=EvtPDL::getId(std::string("pi0"));
1093 daugs[2]=EvtPDL::getId(std::string("pi0"));
1094 }else if(mode == 83){
1095 _ndaugs =3;
1096 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1097 daugs[1]=EvtPDL::getId(std::string("K+"));
1098 daugs[2]=EvtPDL::getId(std::string("K-"));
1099 }else if(mode == 84){
1100 _ndaugs =3;
1101 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1102 daugs[1]=EvtPDL::getId(std::string("K_S0"));
1103 daugs[2]=EvtPDL::getId(std::string("K_S0"));
1104 }else if(mode == 90){
1105 _ndaugs =3;
1106 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1107 daugs[1]=EvtPDL::getId(std::string("pi+"));
1108 daugs[2]=EvtPDL::getId(std::string("pi-"));
1109 }else if(mode == 91){
1110 _ndaugs =3;
1111 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1112 daugs[1]=EvtPDL::getId(std::string("pi+"));
1113 daugs[2]=EvtPDL::getId(std::string("pi-"));
1114 }else if(mode == 92){
1115 _ndaugs =3;
1116 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1117 daugs[1]=EvtPDL::getId(std::string("K+"));
1118 daugs[2]=EvtPDL::getId(std::string("K-"));
1119 }else if(mode == 93){
1120 _ndaugs =2;
1121 daugs[0]=EvtPDL::getId(std::string("D_s+"));
1122 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1123 }else if(mode == 94){
1124 _ndaugs =2;
1125 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
1126 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1127 }else if(mode == 95){
1128 _ndaugs =2;
1129 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
1130 daugs[1]=EvtPDL::getId(std::string("D_s+"));
1131 }else if(mode == 96){
1132 _ndaugs =2;
1133 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
1134 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
1135 }else if(mode == 74100){
1136 _ndaugs =1;
1137 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1138 }else if(mode == 74101){
1139 _ndaugs =1;
1140 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1141 }else if(mode == 74102){
1142 _ndaugs =1;
1143 daugs[0]=EvtPDL::getId(std::string("psi(3770)"));
1144 }else if(mode == 74103){
1145 _ndaugs =1;
1146 daugs[0]=EvtPDL::getId(std::string("phi"));
1147 }else if(mode == 74104){
1148 _ndaugs =1;
1149 daugs[0]=EvtPDL::getId(std::string("omega"));
1150 }else if(mode == 74105){
1151 _ndaugs =1;
1152 daugs[0]=EvtPDL::getId(std::string("rho0"));
1153 }else if(mode == 74106){
1154 _ndaugs =1;
1155 daugs[0]=EvtPDL::getId(std::string("rho(3S)0"));
1156 }else if(mode == 74107){
1157 _ndaugs =1;
1158 daugs[0]=EvtPDL::getId(std::string("omega(2S)"));
1159 }else if(mode == 74110){
1160 _ndaugs =1;
1161 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1162 EvtId myvpho = EvtPDL::getId(std::string("gamma*"));
1163 EvtPDL::reSetMass(myvpho,mhds*0.9); //for calculating maxXS, it will be resent when this mode is selected
1164 _modeFlag.clear();
1165 _modeFlag.push_back(74110); //R-value generator tag
1166 _modeFlag.push_back(0); //to contain the mode selected
1167 }else if(mode == -1){
1168 _ndaugs = nson;
1169 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1170 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
1171 }else if(mode == -2){
1172 _ndaugs = nson;
1173 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1174 }else if(mode == -100){
1175 _ndaugs =1;
1176 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1177 _modeFlag.clear();
1178 _modeFlag.push_back(-100); //R-value generator tag
1179 _modeFlag.push_back(0); //to contain the mode selected
1180 return;
1181 }else if(mode == 10000){//use for check ISR
1182 _ndaugs =2;
1183 daugs[0]=EvtPDL::getId(std::string("pi+"));
1184 daugs[1]=EvtPDL::getId(std::string("pi-"));
1185 }else {
1186 std::cout<<"Bad mode_index number " <<mode<<", please refer to the manual."<<std::endl;
1187 ::abort();
1188 }
1189 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1190
1191 if(VISR){
1192 _ndaugs=2;
1193 daugs[0]=getDaugs()[0];//EvtPDL::getId(std::string("gamma"));
1194 daugs[1]=getDaugs()[1];//EvtPDL::getId(std::string("gamma*"));
1195 }
1196
1197 double fmth=0;
1198 for(int i=0;i<_ndaugs;i++){
1199 double xmi=EvtPDL::getMinMass(daugs[i]);
1200 fmth +=xmi;
1201 }
1202 myxsection = new EvtXsection (mode);
1203 Mthr=myxsection->getXlw();
1204 if(!(_mode==74110 || _mode==-100) ){
1205 if(Mthr<fmth) {std::cout<<"Low end of cross section: ("<<Mthr<<" < (mass of final state)"<<fmth<<") GeV."<< std::endl; abort();}
1206 }
1207}
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition: EvtPDL.hh:51

Referenced by decay(), EvtRexc::decay(), init(), and sumExc().

◆ InitPars()

void EvtConExc::InitPars ( )

Definition at line 3611 of file EvtConExc.cc.

3611 {
3612 threshold=0;
3613 beamEnergySpread=0;
3614 std::string pattern="=";
3615 for(int i=0;i<ncommand;i++){
3616 std::vector<std::string> result=split(commands[i],pattern);
3617 if(result[0]=="threshold") { threshold=atof(result[1].c_str());}else
3618 if(result[0]=="beamEnergySpread"){ beamEnergySpread=atof(result[1].c_str());}
3619 else{
3620 std::cout<<"Your parameter in decay card \""<<result[0]<<"\" incorect"<<std::endl;
3621 std::cout<<"Possible words: threshold, beamEnergySpread"<<std::endl;
3622 abort();
3623 }
3624 }
3625
3626}
EvtConExc::split
std::vector< std::string > split(std::string str, std::string pattern)
Definition: EvtConExc.cc:3590
EvtCyclic3::c_str
char * c_str(Index i)
Definition: EvtCyclic3.cc:252

Referenced by init().

◆ initProbMax()

void EvtConExc::initProbMax ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 1716 of file EvtConExc.cc.

1716 {
1717
1718 noProbMax();
1719
1720}
EvtDecayBase::noProbMax
void noProbMax()
Definition: EvtDecayBase.cc:304

◆ islgr()

bool EvtConExc::islgr ( double *  x,
double *  y,
int  n,
double  t 
)

Definition at line 2792 of file EvtConExc.cc.

2793{ int n0=-1;
2794 double z;
2795 for(int i=0;i<n-1;i++){
2796 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
2797 }
2798 double xstotal=y[599];
2799 if(n0==-1) {return 0;}else{
2800 double p1=y[n0]/xstotal;
2801 double p2=y[n0+1]/xstotal;
2802 double pm= EvtRandom::Flat(0.,1);
2803 if(p1<pm && pm<=p2) {return 1;}else{return 0;}
2804 }
2805}
t
int t()
Definition: t.c:1

◆ ISR_ang_integrate()

double EvtConExc::ISR_ang_integrate ( double  x,
double  theta 
)

Definition at line 2851 of file EvtConExc.cc.

2851 {
2852 //see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
2853 double costheta=cos(theta);
2854 double eb=_cms/2;
2855 double cos2=costheta*costheta;
2856 double sin2=1-cos2;
2857 double me=0.51*0.001;
2858 double L=2*log(_cms/me);
2859 double meE2= me*me/eb/eb;
2860 double hpi=L-1;
2861 double hq1= meE2/2*costheta/(sin2+meE2*cos2);
2862 double hq2= -0.5*log((1+sqrt(1-meE2)*costheta)/(1-sqrt(1-meE2)*costheta));
2863 double hq3=x*x*costheta/2/(x*x-2*x+2)*(1-meE2/(sin2+meE2*cos2));
2864 double hq=(L-1)/2.+hq1+hq2+hq3;
2865 hq /= hpi; //Eq (11) in arXif:hep-ph/9910523
2866 return hq;
2867}
me
const double me
Definition: PipiJpsi.cxx:48

Referenced by ISR_ang_sampling().

◆ ISR_ang_sampling()

double EvtConExc::ISR_ang_sampling ( double  x)

Definition at line 2869 of file EvtConExc.cc.

2869 {
2870 double xx[180],yy[180];
2871 double dgr2rad=1./180.*3.1415926;
2872 xx[0]=0;
2873 xx[1]=5*dgr2rad; //first bin is 5 degree
2874 int nc=2;
2875 for(int i=6;i<=175;i++){ //last bin is 5 degree
2876 xx[nc]=i*dgr2rad;
2877 nc++;
2878 }
2879 xx[nc]=180*dgr2rad;
2880 for(int j=0;j<=nc;j++){
2881 yy[j]=ISR_ang_integrate(x,xx[j]);
2882 }
2883 double pm= EvtRandom::Flat(0.,1);
2884 int mypt=1;
2885 for(int j=1;j<=nc;j++){
2886 if(yy[j-1]<pm && pm<=yy[j]){mypt=j;break;}
2887 }
2888 pm= EvtRandom::Flat(0.,1);
2889 double mytheta= xx[mypt-1]+ (xx[mypt]-xx[mypt-1])*pm;
2890 return mytheta; //theta in rad unit
2891}
EvtConExc::ISR_ang_integrate
double ISR_ang_integrate(double x, double theta)
Definition: EvtConExc.cc:2851

Referenced by decay().

◆ lgr()

double EvtConExc::lgr ( double *  x,
double *  y,
int  n,
double  t 
)

Definition at line 2777 of file EvtConExc.cc.

2778{ int n0=-1;
2779 double z;
2780 for(int i=0;i<n-1;i++){
2781 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
2782 }
2783 if(n0==-1) {return 0.0;}else{
2784 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
2785 z= y[n0+1]+k*(t-x[n0+1]);
2786 }
2787 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
2788 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
2789 return z;
2790}

◆ Li2()

double EvtConExc::Li2 ( double  x)

Definition at line 2771 of file EvtConExc.cc.

2771 {
2772 double li2= -x +x*x/4. - x*x*x/9;
2773 return li2;
2774}

Referenced by SoftPhoton_xs().

◆ LLr()

double EvtConExc::LLr ( double *  x,
double *  y,
int  n,
double  t 
)

Definition at line 2807 of file EvtConExc.cc.

2808{ int n0=-1;
2809 double z;
2810 if( t==x[n-1] ) return y[n-1];
2811 for(int i=0;i<n-1;i++){
2812 if(x[i]<=t && t< x[i+1]) {n0=i;break;}
2813 }
2814 if(n0==-1) {return 0.0;}else{
2815 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
2816 z= y[n0+1]+k*(t-x[n0+1]);
2817 }
2818 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
2819 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
2820 return z;
2821}

◆ meson_sampling()

bool EvtConExc::meson_sampling ( EvtVector4R  pcm,
EvtVector4R  pi 
)

Definition at line 2505 of file EvtConExc.cc.

2505 {
2506 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2507 double theta=angles.getHelAng(1);
2508 double phi =angles.getHelAng(2);
2509 double costheta=cos(theta); //using helicity angles in parent system
2510
2511 double pm= EvtRandom::Flat(0.,1);
2512 double ang = 1 - costheta*costheta;
2513 if(pm< ang/1.) {return true;}else{return false;}
2514}

Referenced by hadron_angle_sampling().

◆ Mhad_sampling()

double EvtConExc::Mhad_sampling ( double *  x,
double *  y 
)

Definition at line 2823 of file EvtConExc.cc.

2823 {//sample ISR hadrons, return Mhadron
2824 //x=MH: array contains the Mhad
2825 //y=AF: array contains the Mhad*Xsection
2826 //n: specify how many bins in the hadron sampling
2827 int n=598; //AF[599] is the total xs including Ecms point
2828 double pm= EvtRandom::Flat(0.,1);
2829 double xrt=1;
2830 int mybin=1;
2831 double xst=y[n];
2832 for(int i=0;i<n;i++){
2833 if((y[i]/xst)<pm && pm<=(y[i+1]/xst)){
2834 mybin=i+1;
2835 break;
2836 }
2837 }
2838 pm= EvtRandom::Flat(0.,1);
2839 if(pm>1){std::cout<<"random larger than 1: "<<pm<<std::endl;}
2840 double mhd=x[mybin-1]+(x[mybin]-x[mybin-1])*pm;
2841
2842 if((mhd - _cms)>0 ){std::cout<<"selected mhd larger than Ecms: "<<mhd<<" > "<<x[mybin] <<std::endl;abort();}
2843 if(mhd<_mhdL){
2844 std::cout<<"the sample mhassrons is less than the low end of XS"<<mhd<<" <"<<_mhdL<<std::endl;
2845 for(int i=0;i<598;i++){std::cout<<i<<" "<<x[i]<<" "<<y[i]<<std::endl;}
2846 abort();
2847 }
2848 return mhd;
2849}

Referenced by decay().

◆ mk_VXS()

void EvtConExc::mk_VXS ( double  Esig,
double  Egamcut,
double  EgamH,
int  midx 
)

Definition at line 3362 of file EvtConExc.cc.

3362 {
3363 //the mode index is put into vmode
3364 //initial
3365 //midx: mode index in vmode[midx] and VXS[midx][bin]
3366 int mode=vmode[midx];
3367 double uscale;
3368
3369 EvtXsection* myxsection = new EvtXsection (mode);
3370 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3371 double s = _cms*_cms;
3372 double mlow=myxsection->getXlw();
3373 if(_cms <= mlow){
3374 for(int i=0;i<600;i++){VXS[midx][i]=0;}
3375 return;
3376 }
3377
3378 int nint=50;
3379 int nmc=800;
3380 //double myxs0 =uscale*gamHXSection(s,Esig,Egamcut,nmc);
3381 double myxs0 =uscale*trapezoid(s,Esig,Egamcut,nint,myxsection);
3382 double xb= 2*Esig/_cms;
3383 double sig_xs = uscale*SoftPhoton_xs(s,xb)*(myxsection->getXsection(_cms));
3384 myxs0 += sig_xs;
3385
3386 int Nstp=598;
3387 double stp=(EgamH - Egamcut)/Nstp;
3388 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
3389 double Eg0=EgamH - i*stp;
3390 double Eg1=EgamH - (i+1)*stp;
3391 //double xsi= uscale*gamHXSection(s,Eg1,Eg0,nmc);
3392 double xsi= uscale*trapezoid(s,Eg1,Eg0,nint,myxsection);
3393 if(i==0) VXS[midx][0]=xsi;
3394 if(i>=1) VXS[midx][i]=VXS[midx][i-1]+xsi;
3395 }
3396 VXS[midx][598]=VXS[midx][597];
3397 VXS[midx][599]=VXS[midx][598]+ myxs0;
3398 //std::cout<<"mode "<<mode<<" "<<uscale<<" "<<VXS[midx][599]<<std::endl;
3399}

Referenced by init().

◆ narrowRXS()

double EvtConExc::narrowRXS ( double  mxL,
double  mxH 
)

Definition at line 3250 of file EvtConExc.cc.

3250 {
3251 //br for ee
3252 double psippee,psipee,jsiee,phiee,omegaee,rhoee;
3253 double kev2Gev=0.000001;
3254 psippee=0.262*kev2Gev;
3255 psipee =2.36*kev2Gev;
3256 jsiee =5.55*kev2Gev;
3257 phiee=4.266*0.001*2.954*0.0001;
3258 omegaee=0.6*kev2Gev;
3259 rhoee=7.04*kev2Gev;
3260 double s=_cms*_cms;
3261 double sigv=0;
3262 double mv=0;
3263 double widee=0;
3264 double xpi=12*3.1415926*3.1415926;
3265 if(mxL<=3.0957 && 3.0957<=mxH || mxL<=3.0983 && 3.0983<=mxH) {
3266 widee = jsiee;
3267 mv = 3.096916;
3268 }else if(mxL<=3.6847 && 3.6847<=mxH || mxL<=3.6873 && 3.6873<=mxH){
3269 widee = psipee;
3270 mv = 3.686109;
3271 }else if(mxL<=1.01906 && 1.01906<=mxH || mxL<= 1.01986 && 1.01986<=mxH){
3272 widee = phiee;
3273 mv = 1.01946;
3274 }else{return -1.0;}
3275
3276 if(_cms<=mv) return -1.;
3277 double xv = 1-mv*mv/s;
3278 sigv += xpi*widee/mv/s*Rad2(s,xv);
3279 double unic = 3.89379304*100000; //in unit nb
3280 return sigv*unic; //vaccum factor has included in the Rad2
3281}
phiee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rh2ee DOUBLE PRECISION omepi DOUBLE PRECISION om3ee DOUBLE PRECISION phiee
Definition: RRes.h:37
rhoee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rhoee
Definition: RRes.h:34
EvtConExc::Rad2
double Rad2(double s, double x)
Definition: EvtConExc.cc:2539

◆ OutStaISR()

void EvtConExc::OutStaISR ( )

Definition at line 3698 of file EvtConExc.cc.

3698 {
3699 int ntot=myFisr.size();
3700 double mymu=0;
3701 for(int i=0;i<ntot;i++){mymu += myFisr[i];}
3702 mymu /= ntot;
3703 double mysig=0;
3704 for(int i=0;i<ntot;i++){ mysig += (myFisr[i]-mymu)*(myFisr[i]-mymu);}
3705 mysig /=ntot;
3706 mysig = sqrt(mysig);
3707 std::cout<<"========= Iteration times over ISR factor: "<<ntot<<std::endl;
3708 std::cout<<" ISR factor * VP factor= observedXS/BornXS: "<<mymu<<" + "<<mysig<<std::endl;
3709}

Referenced by ~EvtConExc().

◆ photonSampling()

bool EvtConExc::photonSampling ( EvtParticle *  part)

Definition at line 3219 of file EvtConExc.cc.

3219 {
3220 bool tagp,tagk;
3221 tagk=0;
3222 tagp=0;
3223 int nds = par->getNDaug();
3224 for(int i=0;i<par->getNDaug();i++){
3225 EvtId di=par->getDaug(i)->getId();
3226 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
3227 int pdgcode =EvtPDL::getStdHep( di );
3228 double alpha=0;
3229
3230 if(pdgcode==111 ||pdgcode==221 || pdgcode==331 ){//for photon
3231 alpha = 1;
3232 }else continue;
3233
3234 double angmax = 1+alpha;
3235 double costheta = cos(p4i.theta());
3236 double ang=1+alpha*costheta*costheta;
3237 double xratio = ang/angmax;
3238 double xi=EvtRandom::Flat(0.,1);
3239 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
3240 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
3241 if(xi>xratio) return false;
3242 }//loop over duaghters
3243 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
3244 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
3245
3246 return true;
3247}

Referenced by decay().

◆ Rad1()

double EvtConExc::Rad1 ( double  s,
double  x 
)

Definition at line 2527 of file EvtConExc.cc.

2527 { //first order correction
2528 //radiator correction upto second order, see Int. J. of Moder.Phys. A
2529 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
2530 double alpha = 1./137.;
2531 double pi=3.1415926;
2532 double me = 0.5*0.001; //e mass
2533 double sme = sqrt(s)/me;
2534 double L = 2*log(sme);
2535 double wsx = 2*alpha/pi/x*(L-1)*(1 - x + x*x/2 );
2536 return wsx;
2537}

Referenced by Rad1difXs().

◆ Rad1difXs()

double EvtConExc::Rad1difXs ( EvtParticle *  p)

Definition at line 2617 of file EvtConExc.cc.

2617 {// // first order xsection
2618 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2619 double summass = p->getDaug(0)->getP4().mass();
2620 for(int i=1;i<_ndaugs;i++){
2621 summass += p->getDaug(i)->getP4().mass();
2622 }
2623
2624 double cms = p->getP4().mass();
2625 double mrg = cms - summass;
2626 double pm= EvtRandom::Flat(0.,1);
2627 double mhs = pm*mrg + summass;
2628
2629 double s = cms * cms;
2630 double x = 1 - mhs*mhs/s;
2631 double wsx = Rad1(s,x);
2632
2633 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2634
2635 double xs = myxsection->getXsection(mhs);
2636 if(xs>XS_max){XS_max = xs;}
2637 double difxs = 2*mhs/cms/cms * wsx *xs;
2638
2639 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2640 differ *= mrg; //Jacobi factor for variable m_hds
2641 difxs *= mrg;
2642 return difxs;
2643}
EvtConExc::Rad1
double Rad1(double s, double x)
Definition: EvtConExc.cc:2527

Referenced by gamHXSection().

◆ Rad2()

double EvtConExc::Rad2 ( double  s,
double  x 
)

Definition at line 2539 of file EvtConExc.cc.

2539 {
2540 //radiator correction upto second order, see Int. J. of Moder.Phys. A, hep-ph/9910523, Eq(8)
2541 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
2542 double alpha = 1./137.;
2543 double pi=3.1415926;
2544 double me = 0.5*0.001; //e mass
2545 double xi2 = 1.64493407;
2546 double xi3=1.2020569;
2547 double sme = sqrt(s)/me;
2548 double L = 2*log(sme);
2549 double beta = 2*alpha/pi*(L-1);
2550 double delta2 = (9./8.-2*xi2)*L*L - (45./16.-5.5*xi2-3*xi3)*L-6./5.*xi2*xi2-4.5*xi3-6*xi2*log(2.)+3./8.*xi2+57./12.;
2551 double ap = alpha/pi;
2552 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
2553 double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
2554 double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
2555 wsx = wsx + beta*beta/8. *wsx2;
2556 double mymx = sqrt(s*(1-x));
2557 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
2558 return wsx*getVP(mymx);//vph is vaccum polarization factor
2559}

Referenced by narrowRXS(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), and Ros_xs().

◆ Rad2difXs() [1/3]

double EvtConExc::Rad2difXs ( double  s,
double  x 
)

Definition at line 2594 of file EvtConExc.cc.

2594 {// leading second order xsection
2595 double wsx = Rad2(s,x);
2596 double mhs = sqrt(s*(1-x));
2597 double xs = myxsection->getXsection(mhs);
2598 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
2599 if(xs>XS_max){XS_max = xs;}
2600 double difxs = wsx *xs;
2601 differ2 = wsx *(myxsection->getErr(mhs));
2602 return difxs;
2603}

◆ Rad2difXs() [2/3]

double EvtConExc::Rad2difXs ( double  s,
double  x,
EvtXsection *  myxsection 
)

Definition at line 2605 of file EvtConExc.cc.

2605 {// leading second order xsection
2606 double wsx = Rad2(s,x);
2607 double mhs = sqrt(s*(1-x));
2608 double xs = myxsection->getXsection(mhs);
2609 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
2610 if(xs>XS_max){XS_max = xs;}
2611 double difxs = wsx *xs;
2612 differ2 = wsx *(myxsection->getErr(mhs));
2613 return difxs;
2614}

◆ Rad2difXs() [3/3]

double EvtConExc::Rad2difXs ( EvtParticle *  p)

Definition at line 2563 of file EvtConExc.cc.

2563 {// leading second order xsection
2564 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2565 double summass = p->getDaug(0)->getP4().mass();
2566 EvtVector4R v4p=p->getDaug(0)->getP4();
2567 for(int i=1;i<_ndaugs;i++){
2568 summass += p->getDaug(i)->getP4().mass();
2569 v4p += p->getDaug(i)->getP4();
2570 }
2571
2572 double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
2573 double cms = p->getP4().mass();
2574 double mrg = cms - summass;
2575 double pm= EvtRandom::Flat(0.,1);
2576 double mhs = pm*mrg + summass;
2577
2578 double s = cms * cms;
2579 double x = 2*Egam/cms;
2580 //double mhs = sqrt( s*(1-x) );
2581 double wsx = Rad2(s,x);
2582
2583 // std::cout<<"random : "<<pm<<std::endl;
2584
2585 double xs = myxsection->getXsection(mhs);
2586 if(xs>XS_max){XS_max = xs;}
2587 double difxs = 2*mhs/cms/cms * wsx *xs;
2588 differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2589 differ *= mrg; //Jacobi factor for variable m_hds
2590 difxs *= mrg;
2591 return difxs;
2592}

Referenced by gamHXSection(), and trapezoid().

◆ ReadVP()

void EvtConExc::ReadVP ( )

Definition at line 3339 of file EvtConExc.cc.

3339 {
3340 vpx.clear();vpr.clear();vpi.clear();
3341 int vpsize=4001;
3342 vpx.resize(vpsize);
3343 vpr.resize(vpsize);
3344 vpi.resize(vpsize);
3345 std::string locvp=getenv("BESEVTGENROOT");
3346 locvp +="/share/vp.dat";
3347 ifstream m_inputFile;
3348 m_inputFile.open(locvp.c_str());
3349 double xx,r1,i1;
3350 double x1,y1,y2;
3351 xx=0;
3352 int mytag=1;
3353 for(int i=0;i<4001;i++){
3354 m_inputFile>>vpx[i]>>vpr[i]>>vpi[i];
3355 }
3356
3357}

Referenced by init().

◆ resetResMass()

void EvtConExc::resetResMass ( )

Definition at line 3041 of file EvtConExc.cc.

3041 {
3042 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
3043 EvtId myres = EvtPDL::getId(std::string("J/psi"));
3044 EvtPDL::reSetMass(myres,mjsi);
3045 //EvtPDL::reSetWidth(myres,wjsi);
3046
3047 myres = EvtPDL::getId(std::string("psi(2S)"));
3048 EvtPDL::reSetMass(myres,mpsip);
3049 //EvtPDL::reSetWidth(myres,wpsip);
3050
3051 myres = EvtPDL::getId(std::string("psi(3770)"));
3052 EvtPDL::reSetMass(myres,mpsipp);
3053 //EvtPDL::reSetWidth(myres,wpsipp);
3054
3055 myres = EvtPDL::getId(std::string("phi"));
3056 EvtPDL::reSetMass(myres,mphi);
3057 //EvtPDL::reSetWidth(myres,wphi);
3058
3059 myres = EvtPDL::getId(std::string("omega"));
3060 EvtPDL::reSetMass(myres,momega);
3061 //EvtPDL::reSetWidth(myres,womega);
3062
3063 myres = EvtPDL::getId(std::string("rho0"));
3064 EvtPDL::reSetMass(myres,mrho0);
3065 //EvtPDL::reSetWidth(myres,wrho0);
3066
3067 myres = EvtPDL::getId(std::string("rho(3S)0"));
3068 EvtPDL::reSetMass(myres,mrho3s);
3069 //EvtPDL::reSetWidth(myres,wrho3s);
3070
3071 myres = EvtPDL::getId(std::string("omega(2S)"));
3072 EvtPDL::reSetMass(myres,momega2s);
3073 //EvtPDL::reSetWidth(myres,womega2s);
3074}

Referenced by decay().

◆ Ros_xs()

double EvtConExc::Ros_xs ( double  mx,
double  bree,
EvtId  pid 
)

Definition at line 2677 of file EvtConExc.cc.

2677 {// in unit nb
2678 double pi=3.1415926;
2679 double s= mx*mx;
2680 double width = EvtPDL::getWidth(pid);
2681 double mass = EvtPDL::getMeanMass(pid);
2682 double xv = 1-mass*mass/s;
2683 double sigma = 12*pi*pi*bree*width/mass/s;
2684 sigma *= Rad2(s, xv);
2685 double nbar = 3.89379304*100000; // ! GeV-2 = 3.89379304*10^5 nbar
2686 return sigma*nbar;
2687}
mass
double mass
Definition: CosmicGenerator.cxx:138

Referenced by init().

◆ selectMass()

double EvtConExc::selectMass ( )

Definition at line 3293 of file EvtConExc.cc.

3293 {
3294 double pm,mhdL,mhds;
3295 pm = EvtRandom::Flat(0.,1);
3296 mhdL =_mhdL;
3297 mhds = pm*(_cms - mhdL)+mhdL; //non narrow resonance
3298 std::vector<double> sxs;
3299 for(int i=0;i<ISRID.size();i++){
3300 double ri=ISRXS[i]/AF[598];
3301 sxs.push_back(ri);
3302 }
3303 for(int j=0;j<sxs.size();j++){
3304 if(j>0) sxs[j] += sxs[j-1];
3305 }
3306 pm = EvtRandom::Flat(0.,1);
3307 if(pm<sxs[0]) {
3308 mhds = EvtPDL::getMass(ISRID[0]); //narrow resonance
3309 }
3310 for(int j=1;j<sxs.size();j++){
3311 double x0 = sxs[j-1];
3312 double x1 = sxs[j];
3313 if(x0<pm && pm<=x1) mhds = EvtPDL::getMass(ISRID[j]); //narrow resonance
3314 }
3315 return mhds;
3316}

◆ selectMode()

int EvtConExc::selectMode ( std::vector< int >  vmod,
double  mhds 
)

Definition at line 2961 of file EvtConExc.cc.

2961 {
2962 //first get xsection for each mode in vmod array
2963 //std::cout<<"vmod.size, mhds "<<vmod.size()<<" "<<mhds<<std::endl;
2964 std::vector<int>excmod;
2965 excmod.push_back(0);
2966 excmod.push_back(1);
2967 excmod.push_back(2);
2968 excmod.push_back(6);
2969 excmod.push_back(7);
2970 excmod.push_back(12);
2971 excmod.push_back(13);
2972 excmod.push_back(45);
2973 excmod.push_back(46);
2974 double uscale = 1;
2975 std::vector<double> vxs;vxs.clear();
2976 for (int i=0;i<vmod.size();i++){
2977 int imod = vmod[i];
2978 delete myxsection;
2979 myxsection =new EvtXsection (imod);
2980 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
2981 //if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
2982
2983 double myxs=uscale*myxsection->getXsection(mhds);
2984
2985 if(i==0) {vxs.push_back(myxs);}
2986 else if(imod==74110){//for continuum process
2987 double xcon = myxs - vxs[i-1];
2988 if(xcon<0) {xcon=vxs[i-1];}else{xcon=myxs;}
2989 if(mhds<2.0) xcon=vxs[i-1]; //continuum cut: _cms energy less than 1.1, sampling phi,rho0 and omega resonance, see parp(2)=1.08 at Pythia.F
2990 vxs.push_back(xcon);
2991 }else{
2992 vxs.push_back(myxs+vxs[i-1]);
2993 }
2994 }//for_i_block
2995 //degugging
2996 // for(int i=0;i<vxs.size();i++){std::cout<<"Mhds="<<mhds<<" Mode="<<vmod[i]<<" xs_i= "<<vxs[i]<<" totalxs "<< vxs[vxs.size()-1]<<std::endl;}
2997
2998
2999 double totxs = vxs[vxs.size()-1];
3000 //if(totxs==0){std::cout<<"totalXs=0, vxs.size()= "<<vxs.size()<<std::endl;}
3001 int icount=0;
3002 int themode;
3003 mode_iter:
3004 double pm= EvtRandom::Flat(0.,1); //std::cout<<"pm= "<<pm<<" mhds= "<<mhds<<std::endl;
3005 if(pm <= vxs[0]/totxs) {
3006 themode= vmod[0];
3007 std::vector<EvtId> theDaug=get_mode(themode);
3008 EvtId daugs[100];
3009 for(int i=0;i<theDaug.size();i++){daugs[i]=theDaug[i];}
3010 int exmode=EvtDecayTable::inChannelList(theparent->getId(),theDaug.size(),daugs);
3011 if(_mode!=74110){return themode;}
3012 else if(exmode==-1){ return themode;}else{goto mycount;}
3013 }
3014
3015 for(int j=1;j<vxs.size();j++){
3016 double x0 = vxs[j-1]/totxs;
3017 double x1 = vxs[j]/totxs; //std::cout<<"j,x0,x1 "<<j<<" "<<x0<<" "<<x1<<std::endl;
3018 if(x0<pm && pm<=x1){
3019 themode=vmod[j];
3020 std::vector<EvtId> theDaug=get_mode(themode);
3021 EvtId daugs[100];
3022 for(int i=0;i<theDaug.size();i++){daugs[i]=theDaug[i];}
3023 int exmode=EvtDecayTable::inChannelList(theparent->getId(),theDaug.size(),daugs);
3024 if(_mode!=74110){return themode;}
3025 else if(exmode==-1){ return themode;} else{ break;}
3026 }
3027 }
3028
3029 mycount:
3030 icount++;
3031 if(icount<20000 ) goto mode_iter;
3032 //debugging
3033 //for(int i=0;i<vxs.size();i++){std::cout<<"Random="<<pm<<" Mode,Mhad "<<vmod[i]<<", "<<mhds<<" xs_i "<<vxs[i]<<" xsi/totalxs "<<vxs[i]/totxs<<std::endl;}
3034 std::cout<<"EvtConExc::selectMode can not find a mode with 20000 iteration for Mhadrons="<<mhds<<std::endl;
3035 return -10;
3036 //abort();
3037
3038}
EvtConExc::get_mode
std::vector< EvtId > get_mode(int mode)
Definition: EvtConExc.cc:1210
EvtDecayTable::inChannelList
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition: EvtDecayTable.cc:795

Referenced by decay().

◆ SetP4()

void EvtConExc::SetP4 ( EvtParticle *  part,
double  mhdr,
double  xeng,
double  theta 
)

Definition at line 2893 of file EvtConExc.cc.

2893 { //set the gamma and gamma* momentum according sampled results
2894 double pm= EvtRandom::Flat(0.,1);
2895 double phi=2*3.1415926*pm;
2896 double gamE = _cms/2*xeng;
2897 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2898 double px = gamE*sin(theta)*cos(phi);
2899 double py = gamE*sin(theta)*sin(phi);
2900 double pz = gamE*cos(theta);
2901 EvtVector4R p4[2];
2902 p4[0].set(gamE,px,py,pz);//ISR photon
2903 p4[1].set(hdrE,-px,-py,-pz); //mhdraons
2904
2905 for(int i=0;i<2;i++){
2906 EvtId myid = p->getDaug(i)->getId();
2907 p->getDaug(i)->init(myid,p4[i]);
2908 if(EvtPDL::name(myid)=="gamma*"){
2909 if( (_cms-mhdr)<0.002 ){
2910 EvtVector4R PG(_cms,0,0,0);
2911 p->getDaug(i)->init(myid,PG);
2912 }
2913 }
2914 }
2915}
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210

Referenced by decay().

◆ SetP4Rvalue()

void EvtConExc::SetP4Rvalue ( EvtParticle *  part,
double  mhdr,
double  xeng,
double  theta 
)

Definition at line 2917 of file EvtConExc.cc.

2917 { //set the gamma and gamma* momentum according sampled results
2918 double pm= EvtRandom::Flat(0.,1);
2919 double phi=2*3.1415926*pm;
2920 double gamE = _cms/2*xeng;
2921 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2922 double px = gamE*sin(theta)*cos(phi);
2923 double py = gamE*sin(theta)*sin(phi);
2924 double pz = gamE*cos(theta);
2925 EvtVector4R p4[2];
2926 p4[0].set(hdrE,px,py,pz);//mhdraons
2927 p4[1].set(gamE,-px,-py,-pz); //ISR photon
2928 for(int i=0;i<2;i++){
2929 EvtId myid = p->getDaug(i)->getId();
2930 p->getDaug(i)->init(myid,p4[i]);
2931 }
2932}

Referenced by decay().

◆ showResMass()

void EvtConExc::showResMass ( )

Definition at line 3112 of file EvtConExc.cc.

3112 {
3113 std::cout<<"J/psi: "<<mjsi<<" "<<wjsi<<std::endl;
3114 std::cout<<"psipr: "<<mpsip<<" "<<wpsip<<std::endl;
3115 std::cout<<"psipp: "<<mpsipp<<" "<<wpsipp<<std::endl;
3116 std::cout<<"phi : "<<mphi<<" "<<wphi<<std::endl;
3117 std::cout<<"omega: "<<momega<<" "<<womega<<std::endl;
3118 std::cout<<"rho0 : "<<mrho0<<" "<<wrho0<<std::endl;
3119 std::cout<<"rho(3S)0: "<<mrho3s<<" "<<wrho3s<<std::endl;
3120 std::cout<<"omega(2S): "<<momega2s<<" "<<womega2s<<std::endl;
3121}

◆ SoftPhoton_xs()

double EvtConExc::SoftPhoton_xs ( double  s,
double  b 
)

Definition at line 2746 of file EvtConExc.cc.

2746 {
2747 double alpha = 1./137.;
2748 double pi=3.1415926;
2749 double me = 0.5*0.001; //e mass
2750 double xi2 = 1.64493407;
2751 double xi3=1.2020569;
2752 double sme = sqrt(s)/me;
2753 double L = 2*log(sme);
2754 double beta = 2*alpha/pi*(L-1);
2755 double delta2 = (9./8.-2*xi2)*L*L - (45./16.-5.5*xi2-3*xi3)*L-6./5.*xi2*xi2-4.5*xi3-6*xi2*log(2.)+3./8.*xi2+57./12.;
2756 double ap = alpha/pi;
2757 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
2758
2759 double beta2 = beta*beta;
2760 double b2 = b*b;
2761
2762 double xs=(-32*b*beta + 8*pow(b,2)*beta - 10*b*pow(beta,2) + pow(b,2)*pow(beta,2) + 32*pow(b,beta)*Delta -
2763 6*(3 - 4*b + pow(b,2))*pow(beta,2)*log(1 - b) - 32*b*pow(beta,2)*log(b) + 8*pow(b,2)*pow(beta,2)*log(b) +
2764 16*pow(beta,2)*Li2(b))/32. ;
2765 double mymx = sqrt(s*(1-b));
2766 //return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
2767 return xs*getVP(mymx); //vph :the vaccum polarzation factor
2768
2769}
EvtConExc::Li2
double Li2(double x)
Definition: EvtConExc.cc:2771

Referenced by calAF(), init(), and mk_VXS().

◆ split()

std::vector< std::string > EvtConExc::split ( std::string  str,
std::string  pattern 
)

Definition at line 3590 of file EvtConExc.cc.

3591{
3592 std::string::size_type pos;
3593 std::vector<std::string> result;
3594 str+=pattern;
3595 int size=str.size();
3596
3597 for(int i=0; i<size; i++)
3598 {
3599 pos=str.find(pattern,i);
3600 if(pos<size)
3601 {
3602 std::string s=str.substr(i,pos-i);
3603 result.push_back(s);
3604 i=pos+pattern.size()-1;
3605 }
3606 }
3607 return result;
3608}

Referenced by InitPars().

◆ sumExc()

double EvtConExc::sumExc ( double  mx)

Definition at line 3137 of file EvtConExc.cc.

3137 {//summation all cross section of exclusive decays
3138 std::vector<int> vmod; vmod.clear();
3139 int mn[82]={0,1,2,3,4,5,6,7,8,9,10,11,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46, // 12,14, 30,31,39,42 and 43 is removed
3140 50,51,
3141 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
3142 90,91,92,93,94,95,96,
3143 74100,74101,74102,74103,74104,74105,74110};
3144 int mn2[83]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46,//mode 14,30,31,42, 43 are removed
3145 50,51,
3146 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
3147 90,91,92,93,94,95,96,
3148 74100,74101,74102,74103,74104,74105,74110};
3149
3150 if(_cms > 2.5){
3151 for(int i=0;i<82;i++){vmod.push_back(mn[i]);}
3152 }else{
3153 for(int i=0;i<83;i++){vmod.push_back(mn2[i]);}
3154 }
3155
3156 // for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
3157
3158 double xsum=0;
3159 double uscale = 1;
3160 for(int i=0;i<vmod.size();i++){
3161 int mymode = vmod[i];
3162 if(mymode ==74110) continue;
3163 delete myxsection;
3164 myxsection =new EvtXsection (mymode);
3165 init_mode(mymode);
3166 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3167 xsum += uscale*myxsection->getXsection(mx);
3168 //if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
3169 }
3170 return xsum;
3171}

Referenced by init().

◆ trapezoid() [1/2]

double EvtConExc::trapezoid ( double  s,
double  a,
double  b,
int  n 
)

Definition at line 3712 of file EvtConExc.cc.

3713{
3714 double xL=2*El/sqrt(s);
3715 double xH=2*Eh/sqrt(s);
3716 double sum = 0.0;
3717 double gaps = (xH-xL)/double(n); //interval
3718 for (int i = 0; i < n; i++)
3719 {
3720 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, staxsection) + Rad2difXs(s,xL + (i+1)*gaps,staxsection) );
3721 }
3722 return sum;
3723}

Referenced by calAF(), init(), and mk_VXS().

◆ trapezoid() [2/2]

double EvtConExc::trapezoid ( double  s,
double  El,
double  Eh,
int  n,
EvtXsection *  myxs 
)

Definition at line 3725 of file EvtConExc.cc.

3726{
3727 double xL=2*El/sqrt(s);
3728 double xH=2*Eh/sqrt(s);
3729 double sum = 0.0;
3730 double gaps = (xH-xL)/double(n); //interval
3731 for (int i = 0; i < n; i++)
3732 {
3733 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, myxs) + Rad2difXs(s,xL + (i+1)*gaps,myxs) );
3734 }
3735 return sum;
3736}

◆ VP_sampling()

bool EvtConExc::VP_sampling ( EvtVector4R  pcm,
EvtVector4R  pi 
)

Definition at line 2516 of file EvtConExc.cc.

2516 {
2517 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2518 double theta=angles.getHelAng(1);
2519 double phi =angles.getHelAng(2);
2520 double costheta=cos(theta); //using helicity angles in parent system
2521
2522 double pm= EvtRandom::Flat(0.,1);
2523 double ang = 1 + costheta*costheta;
2524 if(pm< ang/2.) {return true;}else{return false;}
2525}

Referenced by hadron_angle_sampling().

◆ writeDecayTabel()

void EvtConExc::writeDecayTabel ( )

Definition at line 3429 of file EvtConExc.cc.

3429 {
3430 ofstream oa,ob;
3431 oa.open("_pkhr.dec");
3432 ob.open("obsxs.dat");
3433 //
3434 int im=getModeIndex(74110);
3435
3436 double xscon=VXS[im][599];
3437 //std::cout<<"tsize "<<tsize<<" "<<xscon<<" "<<VXS[70][599]<<std::endl;
3438 std::vector<int> Vmode;
3439 Vmode.push_back(6);//1:pi+pi-
3440 Vmode.push_back(13);//2:pi+pi-2pi0
3441 Vmode.push_back(12);//3:2(pi+pi-)
3442 Vmode.push_back(0);//4:ppbar
3443 Vmode.push_back(1);//5:nnbar
3444 Vmode.push_back(45);//6:K+K-
3445 Vmode.push_back(46);//7:K0K0bar
3446 Vmode.push_back(7);//8:pi+pi-pi0
3447 Vmode.push_back(2);//9:Lambda anti-Lambda
3448 Vmode.push_back(37);//10: pi+pi- eta
3449 std::vector<std::string> vmdl;
3450 vmdl.push_back("PHOKHARA_pipi");
3451 vmdl.push_back("PHOKHARA_pi0pi0pipi");
3452 vmdl.push_back("PHOKHARA_4pi");
3453 vmdl.push_back("PHOKHARA_ppbar");
3454 vmdl.push_back("PHOKHARA_nnbar");
3455 vmdl.push_back("PHOKHARA_KK");
3456 vmdl.push_back("PHOKHARA_K0K0");
3457 vmdl.push_back("PHOKHARA_pipipi0");
3458 vmdl.push_back("PHOKHARA_LLB");
3459 vmdl.push_back("PHOKHARA_pipieta");
3460
3461 ob<<"mode_index observed cross /nb"<<std::endl;
3462 for(int i=0;i<vmode.size();i++){
3463 ob<<vmode[i]<<" "<<VXS[getModeIndex(vmode[i])][599]<<std::endl;
3464 }
3465
3466 for(int i=0;i<Vmode.size();i++){
3467 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
3468 xscon -= VXS[ getModeIndex(Vmode[i]) ][599];
3469 }
3470 //debugging
3471 for(int i=0;i<Vmode.size();i++){
3472 // std::cout<<Vmode[i]<<"-th mode: "<<VXS[getModeIndex(Vmode[i])][599]<<" "<<std::endl;
3473 }
3474
3475 oa<<"LundaPar PARJ(11)=0.611798"<<std::endl;
3476 oa<<"LundaPar PARJ(12)=7.92527E-12"<<std::endl;
3477 oa<<"LundaPar PARJ(14)=0.244495"<<std::endl;
3478 oa<<"LundaPar PARJ(15)=1.16573E-15"<<std::endl;
3479 oa<<"LundaPar PARJ(16)=0.436516"<<std::endl;
3480 oa<<"LundaPar PARJ(17)=0.530517"<<std::endl;
3481 oa<<"LundaPar PARJ(1)=0.0651577"<<std::endl;
3482 oa<<"LundaPar PARJ(2)=0.260378"<<std::endl;
3483 oa<<"LundaPar PARJ(21)=0.0664835"<<std::endl;
3484 oa<<"LundaPar RALPA(15)=0.576687"<<std::endl;
3485 oa<<"LundaPar RALPA(16)=0.364796"<<std::endl;
3486 oa<<"LundaPar RALPA(17)=3.19486E-06"<<std::endl;
3487 oa<<"noPhotos"<<std::endl;
3488 oa<<"Particle vpho "<<_cms<<" 0"<<std::endl;
3489 oa<<"Decay vpho"<<std::endl;
3490 oa<<"0 pi+ pi- PHSP ;"<<std::endl;
3491 oa<<"0 pi0 pi0 pi- pi+ PHSP ;"<<std::endl;
3492 oa<<"0 pi+ pi- pi- pi+ PHSP ;"<<std::endl;
3493 oa<<"0 anti-p- p+ PHSP ;"<<std::endl;
3494 oa<<"0 anti-n0 n0 PHSP ;"<<std::endl;
3495 oa<<"0 K+ K- PHSP ;"<<std::endl;
3496 oa<<"0 K_S0 K_L0 PHSP ;"<<std::endl;
3497 oa<<"0 pi+ pi- pi0 PHSP ;"<<std::endl;
3498 oa<<"0 Lambda0 anti-Lambda0 PHSP ;"<<std::endl;
3499 oa<<"0 pi+ pi- eta PHSP ;"<<std::endl;
3500 oa<<"0 gamma phi PHSP;"<<std::endl;
3501 oa<<"0 gamma rho0 PHSP;"<<std::endl;
3502 oa<<"0 gamma omega PHSP;"<<std::endl;
3503 oa<<xscon<<" ConExc 74110;"<<std::endl;
3504 for(int i=0;i<Vmode.size();i++){
3505 oa<<VXS[ getModeIndex(Vmode[i]) ][599]<<" "<<vmdl[i]<<" ;"<<std::endl;
3506 }
3507 oa<<"Enddecay"<<std::endl;
3508 oa<<"End"<<std::endl;
3509}

Referenced by init().

◆ xs_sampling() [1/2]

bool EvtConExc::xs_sampling ( double  xs)

Definition at line 2476 of file EvtConExc.cc.

2476 {
2477 double pm= EvtRandom::Flat(0.,1);
2478 //std::cout<<"Rdm= "<<pm<<std::endl;
2479 if(pm <xs/XS_max) {return true;} else {return false;}
2480}

◆ xs_sampling() [2/2]

bool EvtConExc::xs_sampling ( double  xs,
double  xs1 
)

Definition at line 2482 of file EvtConExc.cc.

2482 {
2483 double pm= EvtRandom::Flat(0.,1);
2484 // std::cout<<"Rdm= "<<pm<<std::endl;
2485 if(pm <xs/(xs0*1.1)) {return true;} else {return false;}
2486}

Member Data Documentation

◆ _cms

double EvtConExc::_cms
static

Definition at line 140 of file EvtConExc.hh.

Referenced by baryon_sampling(), calAF(), checkEvtRatio(), decay(), difgamXs(), Egam2Mhds(), findMaxXS(), getObsXsection(), init(), ISR_ang_integrate(), Mhad_sampling(), mk_VXS(), narrowRXS(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMass(), SetP4(), SetP4Rvalue(), sumExc(), and writeDecayTabel().

◆ conexcmode

int EvtConExc::conexcmode =-1
static

Definition at line 161 of file EvtConExc.hh.

Referenced by decay(), and EvtRexc::decay().

◆ getMode

int EvtConExc::getMode
static

Definition at line 144 of file EvtConExc.hh.

Referenced by init().

◆ mlow

double EvtConExc::mlow =0
static

Definition at line 184 of file EvtConExc.hh.

Referenced by EvtDecay::energySpread(), init(), and mk_VXS().

◆ multimode

int EvtConExc::multimode =0
static

Definition at line 161 of file EvtConExc.hh.

Referenced by EvtRexc::decay(), and init().

◆ mup

double EvtConExc::mup =0
static

Definition at line 184 of file EvtConExc.hh.

Referenced by EvtDecay::energySpread(), and init().

◆ myxsection

EvtXsection * EvtConExc::myxsection
static

Definition at line 139 of file EvtConExc.hh.

Referenced by decay(), difgamXs(), findMaxXS(), init(), init_mode(), mk_VXS(), Rad1difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMode(), sumExc(), and ~EvtConExc().

◆ SetMthr

double EvtConExc::SetMthr =0
static

Definition at line 142 of file EvtConExc.hh.

Referenced by decay(), init(), and EvtDecay::initialize().

◆ staxsection

EvtXsection * EvtConExc::staxsection
static

Definition at line 139 of file EvtConExc.hh.

Referenced by calAF(), energySpread(), init(), trapezoid(), and ~EvtConExc().

◆ XS_max

double EvtConExc::XS_max
static

Definition at line 141 of file EvtConExc.hh.

Referenced by init(), Rad1difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), and xs_sampling().

The documentation for this class was generated from the following files: