BOSS 7.0.5
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtConExc Class Reference

#include <EvtConExc.hh>

+ Inheritance diagram for EvtConExc:

Public Member Functions

 EvtConExc ()
 
virtual ~EvtConExc ()
 
void getName (std::string &name)
 
EvtDecayBaseclone ()
 
void initProbMax ()
 
void init ()
 
void init_Br_ee ()
 
void decay (EvtParticle *p)
 
void init_mode (int mode)
 
double gamHXSection (EvtParticle *p, double El, double Eh, int nmc=100000)
 
double gamHXSection (double s, double El, double Eh, int nmc=100000)
 
double gamHXSection (double El, double Eh)
 
double gamHXSection (double El, double Eh, int mode)
 
double gamHXSection_er (double El, double Eh)
 
void findMaxXS (EvtParticle *p)
 
double difgamXs (EvtParticle *p)
 
double difgamXs (double mhds, double sintheta)
 
double Mhad_sampling (double *x, double *y)
 
double ISR_ang_integrate (double x, double theta)
 
double ISR_ang_sampling (double x)
 
bool hadron_angle_sampling (EvtVector4R ppi, EvtVector4R pcm)
 
void SetP4 (EvtParticle *part, double mhdr, double xeng, double theta)
 
void SetP4Rvalue (EvtParticle *part, double mhdr, double xeng, double theta)
 
bool gam_sampling (EvtParticle *p)
 
bool xs_sampling (double xs)
 
bool xs_sampling (double xs, double xs1)
 
bool baryon_sampling (EvtVector4R pcm, EvtVector4R pi)
 
bool meson_sampling (EvtVector4R pcm, EvtVector4R pi)
 
bool VP_sampling (EvtVector4R pcm, EvtVector4R pi)
 
bool angularSampling (EvtParticle *part)
 
bool photonSampling (EvtParticle *part)
 
double baryonAng (double mx)
 
double Rad1 (double s, double x)
 
double Rad2 (double s, double x)
 
double Rad2difXs (EvtParticle *p)
 
double Rad2difXs (double s, double x)
 
double Rad2difXs (double s, double x, EvtXsection *myxsection)
 
double Rad1difXs (EvtParticle *p)
 
double Ros_xs (double mx, double bree, EvtId pid)
 
double Li2 (double x)
 
double SoftPhoton_xs (double s, double b)
 
double lgr (double *x, double *y, int n, double t)
 
bool islgr (double *x, double *y, int n, double t)
 
double LLr (double *x, double *y, int n, double t)
 
int selectMode (std::vector< int > vmod, double mhds)
 
void findMaxXS (double mhds)
 
bool gam_sampling (double mhds, double sintheta)
 
void resetResMass ()
 
void getResMass ()
 
bool checkdecay (EvtParticle *p)
 
double sumExc (double mx)
 
void PrintXS (double mx)
 
void showResMass ()
 
int getNdaugs ()
 
EvtIdgetDaugId ()
 
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< EvtIdget_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 EvtDecayBaseclone ()=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 ()
 
EvtIdgetDaugs ()
 
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 EvtXsectionmyxsection
 
static EvtXsectionstaxsection
 
static double _cms
 
static double XS_max
 
static double SetMthr =0
 
static int getMode
 
static int conexcmode =-1
 
static int multimode =0
 
static double mlow =0
 
static double mup =0
 

Additional Inherited Members

- Static Public Member Functions inherited from EvtDecayBase
static void findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
 
static void findMass (EvtParticle *p)
 
static double findMaxMass (EvtParticle *p)
 
- Protected Member Functions inherited from EvtDecayBase
bool daugsDecayedByParentModel ()
 
- Protected Attributes inherited from EvtDecayBase
bool _daugsDecayedByParentModel
 

Detailed Description

Definition at line 138 of file EvtConExc.hh.

Constructor & Destructor Documentation

◆ EvtConExc()

EvtConExc::EvtConExc ( )
inline

Definition at line 142 of file EvtConExc.hh.

142 {
143 } //constructor

Referenced by clone().

◆ ~EvtConExc()

EvtConExc::~EvtConExc ( )
virtual

Definition at line 207 of file EvtConExc.cc.

207 {
208// if(_mode==74110)checkEvtRatio();
209 if(myFisr.size())OutStaISR();
210 if(mydbg){
211 // xs->Write();
212 myfile->Write();
213 }
214
215 delete myxsection;
216 delete staxsection;
217 gamH->deleteTree();
218
219 //the deletion of commands is really uggly!
220
221 int i;
222 if (nconexcdecays==0) {
223 delete [] commands;
224 commands=0;
225 return;
226 }
227
228 for(i=0;i<nconexcdecays;i++){
229 if (conexcdecays[i]==this){
230 conexcdecays[i]=conexcdecays[nconexcdecays-1];
231 nconexcdecays--;
232 if (nconexcdecays==0) {
233 delete [] commands;
234 commands=0;
235 }
236 return;
237 }
238 }
239
240
241}
void OutStaISR()
Definition: EvtConExc.cc:3803
static EvtXsection * staxsection
Definition: EvtConExc.hh:189
static EvtXsection * myxsection
Definition: EvtConExc.hh:189
void deleteTree()
Definition: EvtParticle.cc:557

Member Function Documentation

◆ addNarrowRXS()

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

Definition at line 3391 of file EvtConExc.cc.

3391 {
3392 double myxs = 0;
3393 for(int i=0;i<ISRXS.size();i++){ // std::cout<<"ISRM "<<ISRM[i]<<std::endl;
3394 if( (ISRM[i]-mhi)<binwidth && ISRFLAG[i]){ISRFLAG[i]=0; myxs = ISRXS[i];}
3395 }
3396 //std::cout<<mhi<<" "<<binwidth<<" "<<myxs<<std::endl;
3397 return myxs;
3398}

◆ angularSampling()

bool EvtConExc::angularSampling ( EvtParticle part)

Definition at line 3280 of file EvtConExc.cc.

3280 {
3281 bool tagp,tagk;
3282 tagk=0;
3283 tagp=0;
3284 int nds = par->getNDaug();
3285 for(int i=0;i<par->getNDaug();i++){
3286 EvtId di=par->getDaug(i)->getId();
3287 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
3288 int pdgcode =EvtPDL::getStdHep( di );
3289 double alpha=1;
3290 /*
3291 if(fabs(pdgcode)==2212){//proton and anti-proton
3292 alpha = baryonAng(p4i.mass());
3293 cosp = cos(p4i.theta());
3294 tagp=1;
3295 }else if(fabs(pdgcode)==321||fabs(pdgcode)==211){
3296 alpha=1;
3297 cosk = cos(p4i.theta());
3298 tagk=1;
3299 }else continue;
3300 */
3301 double angmax = 1+alpha;
3302 double costheta = cos(p4i.theta());
3303 double ang=1+alpha*costheta*costheta;
3304 double xratio = ang/angmax;
3305 double xi=EvtRandom::Flat(0.,1);
3306 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
3307 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
3308 if(xi>xratio) return false;
3309 }//loop over duaghters
3310 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
3311 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
3312
3313 return true;
3314}
const double alpha
double cos(const BesAngle a)
Definition: EvtId.hh:27
int getId() const
Definition: EvtId.hh:41
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
static double Flat()
Definition: EvtRandom.cc:74
double theta()
Definition: EvtVector4R.cc:249
float costheta

Referenced by decay().

◆ baryon_sampling()

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

Definition at line 2549 of file EvtConExc.cc.

2549 {
2550 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2551 double theta=angles.getHelAng(1);
2552 double phi =angles.getHelAng(2);
2553 double costheta=cos(theta); //using helicity angles in parent system
2554
2555 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
2556 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
2557 double alpha = baryonAng(_cms);
2558 if(_mode ==96){alpha=-0.34;}
2559 double pm= EvtRandom::Flat(0.,1);
2560 double ang = 1 + alpha*costheta*costheta;
2561 double myang;
2562 if(alpha>=0){myang=1+alpha;}else{myang=1;}
2563 if(pm< ang/myang) {return true;}else{return false;}
2564}
double baryonAng(double mx)
Definition: EvtConExc.cc:3316
static double _cms
Definition: EvtConExc.hh:190

Referenced by hadron_angle_sampling().

◆ baryonAng()

double EvtConExc::baryonAng ( double  mx)

Definition at line 3316 of file EvtConExc.cc.

3316 {
3317 double mp=0.938;
3318 double u = 0.938/mx;
3319 u = u*u;
3320 double u2 = (1+u)*(1+u);
3321 double uu = u*(1+6*u);
3322 double alpha = (u2-uu)/(u2+uu);
3323 return alpha;
3324}
const double mp
Definition: incllambda.cxx:45

Referenced by baryon_sampling().

◆ calAF()

void EvtConExc::calAF ( double  myecms)

Definition at line 3752 of file EvtConExc.cc.

3752 {
3753
3754
3755 double _cms=myecms;
3756 double Esig = 0.0001; //sigularity engy
3757 double x = 2*Egamcut/_cms;
3758 double s = _cms*_cms;
3759
3760 double mhdL=staxsection->getXlw();
3761 double EgamH = (s-mhdL*mhdL)/2/sqrt(s);
3762
3763 int nint=50;
3764 int nmc= 1000;
3765 _xs0 = trapezoid(s,Esig,Egamcut,nint); // std::cout<<"Int: _xs0= "<<_xs0<<std::endl;
3766 //--- sigularity point x from 0 to 0.0001GeV
3767 double xb= 2*Esig/_cms;
3768 double sig_xs = SoftPhoton_xs(s,xb)*(staxsection->getXsection(_cms));
3769 _xs0 += sig_xs;
3770
3771 int Nstp=598;
3772 double stp=(EgamH - Egamcut)/Nstp;
3773 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
3774 double Eg0=EgamH - i*stp;
3775 double Eg1=EgamH - (i+1)*stp;
3776 double xsi= trapezoid(s,Eg1,Eg0,nint); // std::cout<<"Int xsi= " <<xsi<<std::endl;
3777
3778 AA[i]=(Eg1+Eg0)/2;
3779 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
3780 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
3781 double binwidth = mh2-mhi;
3782 if(i==0) AF[0]=xsi;
3783 if(i>=1) AF[i]=AF[i-1]+xsi;
3784 RadXS[i]=xsi/stp;
3785 }
3786 //add the interval 0~Esig
3787 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
3788 AF[598]=AF[597];
3789 AF[599]=AF[598]+ _xs0;
3790 //--
3791 for(int i=0;i<600;i++){
3792 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
3793 }
3794 //--debugging
3795 double bornXS = staxsection->getXsection(_cms);
3796 if(bornXS==0){std::cout<<"EvtConExc::calAF: bad BornXS at "<<_cms<<" GeV"<<std::endl;abort();}
3797 double fisr=AF[599]/bornXS;
3798 myFisr.push_back(fisr);
3799 //std::cout<<"fisr= "<<fisr<<std::endl;
3800}
Double_t x[10]
XmlRpcServer s
Definition: HelloServer.cpp:11
double SoftPhoton_xs(double s, double b)
Definition: EvtConExc.cc:2807
double trapezoid(double s, double a, double b, int n)
Definition: EvtConExc.cc:3817
double getXlw()
Definition: EvtXsection.hh:181
double getXsection(double mx)
Definition: EvtXsection.cc:254

Referenced by decay().

◆ checkdecay()

bool EvtConExc::checkdecay ( EvtParticle p)

Definition at line 3185 of file EvtConExc.cc.

3185 {
3186 int nson=p->getNDaug();
3187 double massflag=1;
3188 for(int i=0;i<nson;i++){
3189 if( EvtPDL::getStdHep(p->getDaug(i)->getId())==22 ) continue;
3190 massflag *= p->getDaug(i)->mass();
3191 }
3192 if(massflag==0){
3193 std::cout<<"Zero mass!"<<std::endl;
3194 return 0;
3195 }else{return 1;}
3196}
EvtId getId() const
Definition: EvtParticle.cc:113
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
double mass() const
Definition: EvtParticle.cc:127

Referenced by decay().

◆ checkEvtRatio()

void EvtConExc::checkEvtRatio ( )

Definition at line 3618 of file EvtConExc.cc.

3618 {
3619 ofstream oa;
3620 oa.open("_evtR.dat");
3621 //
3622 int im=getModeIndex(74110);
3623 double xscon=VXS[im][599];
3624 double xscon0=xscon;
3625 oa<<"Ecms= "<<_cms<<" GeV"<<std::endl;
3626 oa<<"The total observed xs: "<<xscon<<" nb"<<std::endl;
3627 oa<<"=== mode =========== ratio ==========="<<std::endl;
3628 for(int i=0;i<vmode.size();i++){
3629 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
3630 int themode=getModeIndex(vmode[i]);
3631 if(vmode[i]==74110) continue;
3632 xscon -= VXS[themode ][599];
3633 }
3634 if(xscon<0) xscon=0;
3635 //debugging
3636 for(int i=0;i<vmode.size();i++){
3637 int themode=getModeIndex(vmode[i]);
3638 if(vmode[i]==74110) continue;
3639 oa<<vmode[i]<<"-th mode: "<<100*VXS[themode][599]/xscon0<<" % "<<std::endl;
3640 }
3641 oa<<"74110-th mode: "<<100*xscon/xscon0<<" % "<<std::endl;
3642
3643
3644}
std::ofstream oa
Definition: EvtConExc.cc:201
int getModeIndex(int m)
Definition: EvtConExc.cc:3646

◆ clone()

EvtDecayBase * EvtConExc::clone ( )
virtual

Implements EvtDecayBase.

Definition at line 249 of file EvtConExc.cc.

249 {
250
251 return new EvtConExc;
252
253}

◆ command()

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

Reimplemented from EvtDecayBase.

Definition at line 3660 of file EvtConExc.cc.

3660 {
3661 if (ncommand==lcommand){
3662 lcommand=10+2*lcommand;
3663 std::string* newcommands=new std::string[lcommand];
3664 int i;
3665 for(i=0;i<ncommand;i++){
3666 newcommands[i]=commands[i];
3667 }
3668 delete [] commands;
3669 commands=newcommands;
3670 }
3671 commands[ncommand]=cmd;
3672 ncommand++;
3673}

◆ commandName()

std::string EvtConExc::commandName ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 3654 of file EvtConExc.cc.

3654 {
3655
3656 return std::string("ConExcPar");
3657
3658}

◆ decay()

void EvtConExc::decay ( EvtParticle p)
virtual

Implements EvtDecayBase.

Definition at line 1783 of file EvtConExc.cc.

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

◆ difgamXs() [1/2]

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

Definition at line 3012 of file EvtConExc.cc.

3012 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
3013 double x = 1-(mhds/_cms)*(mhds/_cms);
3014 double sin2theta = sintheta*sintheta;
3015 double alpha = 1./137.;
3016 double pi = 3.1415926;
3017 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
3018 double xs = myxsection->getXsection(mhds);
3019 double difxs = 2*mhds/_cms/_cms * wsx *xs;
3020 return difxs;
3021}

◆ difgamXs() [2/2]

double EvtConExc::difgamXs ( EvtParticle p)

Definition at line 2496 of file EvtConExc.cc.

2496 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2497 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2498 EvtVector4R p0 = p->getDaug(0)->getP4();
2499 for(int i=1;i<_ndaugs;i++){
2500 p0 += p->getDaug(i)->getP4();
2501 }
2502 EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
2503 double mhs = p0.mass();
2504 double egam = pgam.get(0);
2505 double sin2theta = pgam.get(3)/pgam.d3mag();
2506 sin2theta = 1-sin2theta*sin2theta;
2507
2508 double cms = p->getP4().mass();
2509 double alpha = 1./137.;
2510 double pi = 3.1415926;
2511 double x = 2*egam/cms;
2512 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2513 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2514
2515 double xs = myxsection->getXsection(mhs);
2516 double difxs = 2*mhs/cms/cms * wsx *xs;
2517 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2518 return difxs;
2519
2520}
double getErr(double mx)
Definition: EvtXsection.cc:260

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

◆ Egam2Mhds()

double EvtConExc::Egam2Mhds ( double  Egam)

Definition at line 3527 of file EvtConExc.cc.

3527 {
3528 double s=_cms*_cms;
3529 double mhds = sqrt( s - 2*Egam*_cms);
3530 return mhds;
3531}

◆ energySpread()

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

Definition at line 3734 of file EvtConExc.cc.

3734 {
3735 //mu: mean value in Gaussian
3736 //sigma: variance in Gaussian
3737 rloop:
3738 int n=12;
3739 double ri=0;
3740 for(int i=0;i<n;i++){
3741 double pm= EvtRandom::Flat(0.,1);
3742 ri += pm;
3743 }
3744 double eta=sqrt(n*12.0)*(ri/12-0.5);
3745 double xsig= sigma*eta+mu;//smapling Gaussion value
3746 double mx0=staxsection->getXlw();
3747 double mx1=staxsection->getXup();
3748 if(xsig<mx0 || xsig>mx1) goto rloop;
3749 return xsig;
3750}
TTree * sigma
const Int_t n
double getXup()
Definition: EvtXsection.hh:180

◆ findMaxXS() [1/2]

void EvtConExc::findMaxXS ( double  mhds)

Definition at line 2997 of file EvtConExc.cc.

2997 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2998 maxXS=-1;
2999 for(int i=0;i<90000;i++){
3000 double x = 1-(mhds/_cms)*(mhds/_cms);
3001 double cos = EvtRandom::Flat(0.0006,0.9994); //set cut on ISR gamma 2^o
3002 double sin2theta =sqrt(1-cos*cos);
3003 double alpha = 1./137.;
3004 double pi = 3.1415926;
3005 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
3006 double xs = myxsection->getXsection(mhds);
3007 double difxs = 2*mhds/_cms/_cms * wsx *xs;
3008 if(difxs>maxXS)maxXS=difxs;
3009 }
3010}

◆ findMaxXS() [2/2]

void EvtConExc::findMaxXS ( EvtParticle p)

Definition at line 2451 of file EvtConExc.cc.

2451 {
2452 //std::cout<<"nmc= "<<nmc<<std::endl;
2453 gamId = EvtPDL::getId(std::string("gamma"));
2454 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
2455 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
2456 double totxs = 0;
2457 maxXS=-1;
2458 _er1=0;
2459 Rad2Xs =0;
2460 int nmc = 50000;
2461 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
2462 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
2463 int gamdaugs = _ndaugs+1;
2464
2465 EvtId theDaugs[10];
2466 for(int i=0;i<_ndaugs;i++){
2467 theDaugs[i] = daugs[i];
2468 }
2469 theDaugs[_ndaugs]=gamId;
2470
2471 gamH->makeDaughters(gamdaugs,theDaugs);
2472 loop:
2473 double totm=0;
2474 for(int i=0;i<gamdaugs;i++){
2475 EvtParticle* di=gamH->getDaug(i);
2476 double xmi=EvtPDL::getMass(di->getId());
2477 di->setMass(xmi);
2478 totm += xmi;
2479 }
2480
2481 //std::cout<<totm<<" "<<p_init.get(0)<<std::endl;
2482 if(totm >= p_init.get(0) ) goto loop;
2483
2484 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
2485 double thexs = difgamXs(gamH); //prepare the photon angular distribution
2486 EvtVector4R p4gam = gamH->getDaug(_ndaugs)->getP4();
2487 double costheta = p4gam.get(3)/p4gam.d3mag();
2488 double sintheta = sqrt(1-costheta*costheta);
2489 bool acut=(sintheta>0.04) && (sintheta<0.96); //set photon auglar cut 2^o
2490 if(thexs>maxXS && acut ) {maxXS=thexs;}
2491 //gamH->deleteTree();
2492 }
2493
2494}
*********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
double difgamXs(EvtParticle *p)
Definition: EvtConExc.cc:2496
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)

Referenced by decay(), and init().

◆ gam_sampling() [1/2]

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

Definition at line 2530 of file EvtConExc.cc.

2530 {
2531 double pm= EvtRandom::Flat(0.,1);
2532 double xs = difgamXs( mhds,sintheta );
2533 double xsratio = xs/maxXS;
2534 if(pm<xsratio){return true;}else{return false;}
2535}

◆ gam_sampling() [2/2]

bool EvtConExc::gam_sampling ( EvtParticle p)

Definition at line 2523 of file EvtConExc.cc.

2523 {//photon angular distribution sampling
2524 double pm= EvtRandom::Flat(0.,1);
2525 double xs = difgamXs( p );
2526 double xsratio = xs/maxXS;
2527 if(pm<xsratio){return true;}else{return false;}
2528}

Referenced by decay().

◆ gamHXSection() [1/4]

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

Definition at line 2397 of file EvtConExc.cc.

2397 {// using Gaussian integration subroutine from Cern lib
2398 std::cout<<" "<<std::endl;
2399 extern double Rad2difXs(double *x);
2400 extern double Rad2difXs2(double *x);
2401 double eps = 0.01;
2402 double Rad2Xs;
2403 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2404 if(Rad2Xs==0){
2405 for(int iter=0;iter<10;iter++){
2406 eps += 0.01;
2407 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2408 if(!Rad2Xs) return Rad2Xs;
2409 }
2410 }
2411 return Rad2Xs; // the leading second order correction
2412}
double Rad2difXs2(double *mhs)
Definition: EvtConExc.cc:2778
double dgauss_(double(*fun)(double *), double *, double *, double *)
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
EvtTensor3C eps(const EvtVector3R &v)
Definition: EvtTensor3C.cc:307
double Rad2difXs(EvtParticle *p)
Definition: EvtConExc.cc:2624

◆ gamHXSection() [2/4]

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

Definition at line 2415 of file EvtConExc.cc.

2415 {// using Gaussian integration subroutine from Cern lib
2416 std::cout<<" "<<std::endl;
2417 extern double Rad2difXs(double *x);
2418 extern double Rad2difXs2(double *x);
2419 double eps = 0.01;
2420 double Rad2Xs;
2421 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2422 if(Rad2Xs==0){
2423 for(int iter=0;iter<10;iter++){
2424 eps += 0.01;
2425 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2426 if(!Rad2Xs) return Rad2Xs;
2427 }
2428 }
2429 return Rad2Xs; // the leading second order correction
2430}

◆ gamHXSection() [3/4]

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

Definition at line 2365 of file EvtConExc.cc.

2365 {//El, Eh : the lower and higher energy of radiative photons
2366 //--for narrow resonance psi(2S),J/psi, phi, which can not calculated with the MC integration
2367 //double mxL = sqrt( s-2*Eh*sqrt(s) ); //low mass
2368 //double mxH = sqrt( s-2*El*sqrt(s) ); //high mass
2369 //double sigv = narrowRXS(mxL,mxH);
2370 //--
2371
2372 double totxs = 0;
2373 maxXS=-1;
2374 _er1=0;
2375 Rad2Xs =0;
2376 double xL=2*El/sqrt(s);
2377 double xH=2*Eh/sqrt(s);
2378 for(int i=0;i<nmc;i++){//It is found that the narrow resonance can not included in this MC integartion
2379 double rdm = EvtRandom::Flat(0.,1.);// set angular cut 1^o
2380 double x= xL+ (xH-xL)*rdm;
2381 Rad2Xs += Rad2difXs(s,x);
2382 _er1 += differ2; //stored when calculate Rad2Xs
2383 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
2384 }
2385 _er1/=nmc;
2386 _er1*=(xH-xL);
2387 //std::cout<<"_er1= "<<_er1<<std::endl;
2388 Rad2Xs/=nmc; // the leading second order correction
2389 double thexs= Rad2Xs*(xH-xL); //xH-xL is the Jacobi factor
2390 //std::cout<<"thexs= "<<thexs<<std::endl;
2391 //if(sigv != -1) thexs += sigv; //add narrow resonance ISR cross section
2392 return thexs;
2393}

◆ gamHXSection() [4/4]

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

Definition at line 2306 of file EvtConExc.cc.

2306 {//El, Eh : the lower and higher energy of radiative photons
2307 //std::cout<<"nmc= "<<nmc<<std::endl;
2308 gamId = EvtPDL::getId(std::string("gamma"));
2309 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
2310 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
2311 double totxs = 0;
2312 maxXS=-1;
2313 _er1=0;
2314 Rad2Xs =0;
2315 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
2316 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
2317 int gamdaugs = _ndaugs+1;
2318
2319 EvtId theDaugs[10];
2320 for(int i=0;i<_ndaugs;i++){
2321 theDaugs[i] = daugs[i];
2322 }
2323 theDaugs[_ndaugs]=gamId;
2324
2325 gamH->makeDaughters(gamdaugs,theDaugs);
2326
2327 for(int i=0;i<gamdaugs;i++){
2328 EvtParticle* di=gamH->getDaug(i);
2329 double xmi=EvtPDL::getMass(di->getId());
2330 di->setMass(xmi);
2331 }
2332 loop:
2333 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
2334 //-- slection:theta_gamma > 1 degree
2335 EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
2336 double pmag = pgam.d3mag();
2337 double pz = pgam.get(3);
2338 double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
2339 double onedegree = 1./180.*3.1415926;
2340 //if(sintheta < onedegree) {goto loop;}
2341 if(pmag <El || pmag >Eh) {goto loop;}
2342 //--------
2343 // std::cout<<"pmag= "<<pmag<<std::endl;
2344
2345 double thexs = difgamXs(gamH); //prepare the photon angular distribution
2346 Rad2Xs += Rad2difXs( gamH );
2347 if(thexs>maxXS) {maxXS=thexs;}
2348 double s = p_init.mass2();
2349 double x = 2*pgam.get(0)/sqrt(s);
2350 double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
2351 totxs += rad1xs;
2352 _er1 += differ;
2353 gamH->deleteDaughters();
2354 } //for int_i block
2355 _er1/=nmc;
2356
2357 Rad2Xs/=nmc; // the leading second order correction
2358 totxs/=nmc; // first order correction XS
2359
2360 // return totxs; // first order correction XS
2361 return Rad2Xs; // the leading second order correction
2362}
double Rad1difXs(EvtParticle *p)
Definition: EvtConExc.cc:2678
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
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 2433 of file EvtConExc.cc.

2433 {// using Gaussian integration subroutine from Cern lib
2434 std::cout<<" "<<std::endl;
2435 extern double Rad2difXs_er(double *x);
2436 extern double Rad2difXs_er2(double *x);
2437 double eps = 0.01;
2438 double Rad2Xs;
2439 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
2440 if(Rad2Xs==0){
2441 for(int iter=0;iter<10;iter++){
2442 eps += 0.01;
2443 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
2444 if(!Rad2Xs) return Rad2Xs;;
2445 }
2446 }
2447 return Rad2Xs; // the leading second order correction
2448}
double Rad2difXs_er2(double *mhs)
Definition: EvtConExc.cc:2793
double Rad2difXs_er(double *m)
Definition: EvtConExc.cc:2765

◆ get_mode()

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

Definition at line 1261 of file EvtConExc.cc.

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

Referenced by selectMode().

◆ get_mode_index()

int EvtConExc::get_mode_index ( int  mode)

Definition at line 3505 of file EvtConExc.cc.

3505 {
3506 int i=0;
3507 for(int j=0;i<vmode.size();j++){
3508 if(mode==vmode[j]) return j;
3509 }
3510 std::cout<<" EvtConExc::get_mode_index: no index is found in vmode"<<std::endl;
3511 abort();
3512}

Referenced by getObsXsection().

◆ getDaugId()

EvtId * EvtConExc::getDaugId ( )
inline

Definition at line 210 of file EvtConExc.hh.

210{return daugs;}

Referenced by EvtRexc::decay().

◆ getModeIndex()

int EvtConExc::getModeIndex ( int  m)

Definition at line 3646 of file EvtConExc.cc.

3646 {
3647 for (int i=0;i<vmode.size();i++){
3648 if(m==vmode[i]) return i;
3649 }
3650 std::cout<<"EvtConExc::getModeIndex: no index in vmode found "<<std::endl;
3651 abort();
3652}

Referenced by checkEvtRatio(), and writeDecayTabel().

◆ getName()

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

Implements EvtDecayBase.

Definition at line 243 of file EvtConExc.cc.

243 {
244
245 model_name="ConExc";
246
247}

◆ getNdaugs()

int EvtConExc::getNdaugs ( )
inline

Definition at line 209 of file EvtConExc.hh.

209{return _ndaugs;}

Referenced by EvtRexc::decay().

◆ getObsXsection()

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

Definition at line 3514 of file EvtConExc.cc.

3514 {
3515 double hwid=(AA[0]-AA[1])/2.;
3516 double s=_cms*_cms;
3517 int idx=get_mode_index(mode);
3518 double Egam=(s-mhds*mhds)/2/_cms;
3519 for(int i=0;i<600;i++){
3520 if( (AA[i]-hwid)<=Egam && Egam<(AA[i]+hwid) ) return VXS[idx][i];
3521 }
3522
3523 std:cout<<"EvtConExc::getObsXsection : no observed xs is found "<<endl;
3524 abort();
3525}
int get_mode_index(int mode)
Definition: EvtConExc.cc:3505

◆ getResMass()

void EvtConExc::getResMass ( )

Definition at line 3138 of file EvtConExc.cc.

3138 {
3139 EvtId myres = EvtPDL::getId(std::string("J/psi"));
3140 mjsi = EvtPDL::getMeanMass(myres);
3141 wjsi = EvtPDL::getWidth(myres);
3142
3143 myres = EvtPDL::getId(std::string("psi(2S)"));
3144 mpsip= EvtPDL::getMeanMass(myres);
3145 wpsip= EvtPDL::getWidth(myres);
3146
3147 myres = EvtPDL::getId(std::string("psi(3770)"));
3148 mpsipp= EvtPDL::getMeanMass(myres);
3149 wpsipp= EvtPDL::getWidth(myres);
3150
3151 myres = EvtPDL::getId(std::string("phi"));
3152 mphi = EvtPDL::getMeanMass(myres);
3153 wphi = EvtPDL::getWidth(myres);
3154
3155 myres = EvtPDL::getId(std::string("omega"));
3156 momega= EvtPDL::getMeanMass(myres);
3157 womega= EvtPDL::getWidth(myres);
3158
3159 myres = EvtPDL::getId(std::string("rho0"));
3160 mrho0 = EvtPDL::getMeanMass(myres);
3161 wrho0 = EvtPDL::getWidth(myres);
3162
3163 myres = EvtPDL::getId(std::string("rho(3S)0"));
3164 mrho3s= EvtPDL::getMeanMass(myres);
3165 wrho3s= EvtPDL::getWidth(myres);
3166
3167
3168 myres = EvtPDL::getId(std::string("omega(2S)"));
3169 momega2s= EvtPDL::getMeanMass(myres);
3170 womega2s= EvtPDL::getWidth(myres);
3171
3172}
static double getWidth(EvtId i)
Definition: EvtPDL.hh:54

Referenced by init().

◆ getSelectedMode()

int EvtConExc::getSelectedMode ( )
inline

Definition at line 211 of file EvtConExc.hh.

211{return _selectedMode;}

◆ getVP()

double EvtConExc::getVP ( double  cms)

Definition at line 3425 of file EvtConExc.cc.

3425 {
3426 double xx,r1,i1;
3427 double x1,y1,y2;
3428 xx=0;
3429 int mytag=1;
3430 for(int i=0;i<4001;i++){
3431 x1=vpx[i];
3432 y1=vpr[i];
3433 y2=vpi[i];
3434 if(xx<=mx && mx <x1){xx=x1; r1=y1; i1=y2; mytag=0; break;}
3435 xx=x1; r1=y1; i1=y2;
3436 }
3437 if(mytag==1){std::cout<<"No vacuum polarization value found"<<std::endl;abort();}
3438 EvtComplex cvp(r1,i1);
3439 double thevp=abs(1./(1-cvp));
3440 if(3.0933<mx && mx<3.1035) return 1.; //J/psi peak excluded
3441 if(3.6810<mx && mx<3.6913) return 1.; //psi(2S) peak excluded
3442 return thevp*thevp;
3443}

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

◆ hadron_angle_sampling()

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

Definition at line 2278 of file EvtConExc.cc.

2278 {
2279 EvtVector4R pbst=-1*pcm;
2280 pbst.set(0,pcm.get(0));
2281 EvtVector4R p4i = boostTo(ppi,pbst);
2282 if((_mode>=0 && _mode<=5) || _mode==44 || _mode==96){//ee->two baryon mode, VP, SP, mode=-2 is excluded
2283 bool byn_ang = baryon_sampling(pcm, p4i);
2284 return byn_ang;
2285 }else if(_mode==6 || _mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
2286 bool msn_ang = meson_sampling(pcm,p4i);
2287 return msn_ang;
2288 }else if(_mode==23 || _mode==24 || _mode==25 || _mode==26 || _mode==27 || _mode==36 || _mode==39 || _mode==47 || _mode==48 || _mode==52 || _mode==68 || _mode==69 || _mode==72 || _mode==73 || _mode==80 || _mode==94 || _mode==95){
2289 bool msn_ang = VP_sampling(pcm,p4i);
2290 return msn_ang;
2291 }else if(_mode==-2){
2294 if(type0==EvtSpinType::SCALAR &&type1==EvtSpinType::SCALAR){
2295 bool msn_ang = meson_sampling(pcm,p4i);
2296 return msn_ang;
2297 }else if(type0==EvtSpinType::VECTOR &&type1==EvtSpinType::SCALAR || type1==EvtSpinType::VECTOR &&type0==EvtSpinType::SCALAR){
2298 bool msn_ang = VP_sampling(pcm,p4i);
2299 return msn_ang;
2300 }
2301 }
2302 return true;
2303}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2577
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2566
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2549
static EvtSpinType::spintype getSpinType(EvtId i)
Definition: EvtPDL.hh:61
void set(int i, double d)
Definition: EvtVector4R.hh:183

Referenced by decay().

◆ init()

void EvtConExc::init ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 256 of file EvtConExc.cc.

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

◆ init_Br_ee()

void EvtConExc::init_Br_ee ( )

Definition at line 2706 of file EvtConExc.cc.

2706 {
2707 // double psipp_ee=9.6E-06;
2708 double psip_ee =7.73E-03;
2709 double jsi_ee =5.94*0.01;
2710 double phi_ee =2.954E-04;
2711 // double omega_ee=7.28E-05;
2712 // double rho_ee = 4.72E-05;
2713 EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
2714 EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
2715 EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
2716 EvtId phiId=EvtPDL::getId(std::string("phi"));
2717 EvtId omegaId=EvtPDL::getId(std::string("omega"));
2718 EvtId rhoId=EvtPDL::getId(std::string("rho0"));
2719 BR_ee.clear();
2720 ResId.clear();
2721
2722 //BR_ee.push_back(rho_ee);
2723 //BR_ee.push_back(omega_ee);
2724 BR_ee.push_back(phi_ee);
2725 BR_ee.push_back(jsi_ee);
2726 BR_ee.push_back(psip_ee);
2727 //BR_ee.push_back(psipp_ee);
2728
2729 //ResId.push_back(rhoId);
2730 //ResId.push_back(omegaId);
2731 ResId.push_back(phiId);
2732 ResId.push_back(jsiId);
2733 ResId.push_back(pspId);
2734 //ResId.push_back(psppId);
2735
2736}

Referenced by init().

◆ init_mode()

void EvtConExc::init_mode ( int  mode)

Definition at line 751 of file EvtConExc.cc.

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

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

◆ InitPars()

void EvtConExc::InitPars ( )

Definition at line 3716 of file EvtConExc.cc.

3716 {
3717 threshold=0;
3718 beamEnergySpread=0;
3719 std::string pattern="=";
3720 for(int i=0;i<ncommand;i++){
3721 std::vector<std::string> result=split(commands[i],pattern);
3722 if(result[0]=="threshold") { threshold=atof(result[1].c_str());}else
3723 if(result[0]=="beamEnergySpread"){ beamEnergySpread=atof(result[1].c_str());}
3724 else{
3725 std::cout<<"Your parameter in decay card \""<<result[0]<<"\" incorect"<<std::endl;
3726 std::cout<<"Possible words: threshold, beamEnergySpread"<<std::endl;
3727 abort();
3728 }
3729 }
3730
3731}
std::vector< std::string > split(std::string str, std::string pattern)
Definition: EvtConExc.cc:3695
char * c_str(Index i)
Definition: EvtCyclic3.cc:252

Referenced by init().

◆ initProbMax()

void EvtConExc::initProbMax ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 1777 of file EvtConExc.cc.

1777 {
1778
1779 noProbMax();
1780
1781}
void noProbMax()

◆ islgr()

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

Definition at line 2854 of file EvtConExc.cc.

2855{ int n0=-1;
2856 double z;
2857 for(int i=0;i<n-1;i++){
2858 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
2859 }
2860 double xstotal=y[599];
2861 if(n0==-1) {return 0;}else{
2862 double p1=y[n0]/xstotal;
2863 double p2=y[n0+1]/xstotal;
2864 double pm= EvtRandom::Flat(0.,1);
2865 if(p1<pm && pm<=p2) {return 1;}else{return 0;}
2866 }
2867}
TTree * t
Definition: binning.cxx:23
double y[1000]

◆ ISR_ang_integrate()

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

Definition at line 2913 of file EvtConExc.cc.

2913 {
2914 //see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
2915 double costheta=cos(theta);
2916 double eb=_cms/2;
2917 double cos2=costheta*costheta;
2918 double sin2=1-cos2;
2919 double me=0.51*0.001;
2920 double L=2*log(_cms/me);
2921 double meE2= me*me/eb/eb;
2922 double hpi=L-1;
2923 double hq1= meE2/2*costheta/(sin2+meE2*cos2);
2924 double hq2= -0.5*log((1+sqrt(1-meE2)*costheta)/(1-sqrt(1-meE2)*costheta));
2925 double hq3=x*x*costheta/2/(x*x-2*x+2)*(1-meE2/(sin2+meE2*cos2));
2926 double hq=(L-1)/2.+hq1+hq2+hq3;
2927 hq /= hpi; //Eq (11) in arXif:hep-ph/9910523
2928 return hq;
2929}
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 2931 of file EvtConExc.cc.

2931 {
2932 double xx[180],yy[180];
2933 double dgr2rad=1./180.*3.1415926;
2934 xx[0]=0;
2935 xx[1]=5*dgr2rad; //first bin is 5 degree
2936 int nc=2;
2937 for(int i=6;i<=175;i++){ //last bin is 5 degree
2938 xx[nc]=i*dgr2rad;
2939 nc++;
2940 }
2941 xx[nc]=180*dgr2rad;
2942 for(int j=0;j<=nc;j++){
2943 yy[j]=ISR_ang_integrate(x,xx[j]);
2944 }
2945 double pm= EvtRandom::Flat(0.,1);
2946 int mypt=1;
2947 for(int j=1;j<=nc;j++){
2948 if(yy[j-1]<pm && pm<=yy[j]){mypt=j;break;}
2949 }
2950 pm= EvtRandom::Flat(0.,1);
2951 double mytheta= xx[mypt-1]+ (xx[mypt]-xx[mypt-1])*pm;
2952 return mytheta; //theta in rad unit
2953}
double ISR_ang_integrate(double x, double theta)
Definition: EvtConExc.cc:2913
const int nc
Definition: histgen.cxx:26

Referenced by decay().

◆ lgr()

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

Definition at line 2839 of file EvtConExc.cc.

2840{ int n0=-1;
2841 double z;
2842 for(int i=0;i<n-1;i++){
2843 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
2844 }
2845 if(n0==-1) {return 0.0;}else{
2846 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
2847 z= y[n0+1]+k*(t-x[n0+1]);
2848 }
2849 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
2850 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
2851 return z;
2852}

◆ Li2()

double EvtConExc::Li2 ( double  x)

Definition at line 2832 of file EvtConExc.cc.

2832 {
2833// double li2= -x +x*x/4. - x*x*x/9; // wangwp: may be not correct!
2834 double li2= +x +x*x/4. + x*x*x/9; // corrected by wangwp
2835 return li2;
2836}

Referenced by SoftPhoton_xs().

◆ LLr()

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

Definition at line 2869 of file EvtConExc.cc.

2870{ int n0=-1;
2871 double z;
2872 if( t==x[n-1] ) return y[n-1];
2873 for(int i=0;i<n-1;i++){
2874 if(x[i]<=t && t< x[i+1]) {n0=i;break;}
2875 }
2876 if(n0==-1) {return 0.0;}else{
2877 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
2878 z= y[n0+1]+k*(t-x[n0+1]);
2879 }
2880 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
2881 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
2882 return z;
2883}

◆ meson_sampling()

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

Definition at line 2566 of file EvtConExc.cc.

2566 {
2567 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2568 double theta=angles.getHelAng(1);
2569 double phi =angles.getHelAng(2);
2570 double costheta=cos(theta); //using helicity angles in parent system
2571
2572 double pm= EvtRandom::Flat(0.,1);
2573 double ang = 1 - costheta*costheta;
2574 if(pm< ang/1.) {return true;}else{return false;}
2575}

Referenced by hadron_angle_sampling().

◆ Mhad_sampling()

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

Definition at line 2885 of file EvtConExc.cc.

2885 {//sample ISR hadrons, return Mhadron
2886 //x=MH: array contains the Mhad
2887 //y=AF: array contains the Mhad*Xsection
2888 //n: specify how many bins in the hadron sampling
2889 int n=598; //AF[599] is the total xs including Ecms point
2890 double pm= EvtRandom::Flat(0.,1);
2891 double xrt=1;
2892 int mybin=1;
2893 double xst=y[n];
2894 for(int i=0;i<n;i++){
2895 if((y[i]/xst)<pm && pm<=(y[i+1]/xst)){
2896 mybin=i+1;
2897 break;
2898 }
2899 }
2900 pm= EvtRandom::Flat(0.,1);
2901 if(pm>1){std::cout<<"random larger than 1: "<<pm<<std::endl;}
2902 double mhd=x[mybin-1]+(x[mybin]-x[mybin-1])*pm;
2903
2904 if((mhd - _cms)>0){std::cout<<"selected mhd larger than Ecms: "<<mhd<<" > "<<x[mybin] <<std::endl;abort();}
2905 if(mhd<_mhdL){
2906 std::cout<<"the sample mhassrons is less than the low end of XS"<<mhd<<" <"<<_mhdL<<std::endl;
2907 for(int i=0;i<598;i++){std::cout<<i<<" "<<x[i]<<" "<<y[i]<<std::endl;}
2908 abort();
2909 }
2910 return mhd;
2911}

Referenced by decay().

◆ mk_VXS()

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

Definition at line 3466 of file EvtConExc.cc.

3466 {
3467 //the mode index is put into vmode
3468 //initial
3469 //midx: mode index in vmode[midx] and VXS[midx][bin]
3470 int mode=vmode[midx];
3471 double uscale;
3472
3473 EvtXsection* myxsection = new EvtXsection (mode);
3474 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3475 double s = _cms*_cms;
3476 double mlow=myxsection->getXlw();
3477 if(_cms <= mlow){
3478 for(int i=0;i<600;i++){VXS[midx][i]=0;}
3479 return;
3480 }
3481
3482 int nint=50;
3483 int nmc=800;
3484 //double myxs0 =uscale*gamHXSection(s,Esig,Egamcut,nmc);
3485 double myxs0 = uscale*trapezoid(s,Esig,Egamcut,nint,myxsection);
3486 double xb= 2*Esig/_cms;
3487 double sig_xs = uscale*SoftPhoton_xs(s,xb)*(myxsection->getXsection(_cms));
3488 myxs0 += sig_xs;
3489
3490 int Nstp=598;
3491 double stp=(EgamH - Egamcut)/Nstp;
3492 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
3493 double Eg0=EgamH - i*stp;
3494 double Eg1=EgamH - (i+1)*stp;
3495 //double xsi= uscale*gamHXSection(s,Eg1,Eg0,nmc);
3496 double xsi= uscale*trapezoid(s,Eg1,Eg0,nint,myxsection);
3497 if(i==0) VXS[midx][0]=xsi;
3498 if(i>=1) VXS[midx][i]=VXS[midx][i-1]+xsi;
3499 }
3500 VXS[midx][598]=VXS[midx][597];
3501 VXS[midx][599]=VXS[midx][598]+ myxs0;
3502 //std::cout<<"mode "<<mode<<" "<<uscale<<" "<<VXS[midx][599]<<std::endl;
3503}

Referenced by init().

◆ narrowRXS()

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

Definition at line 3357 of file EvtConExc.cc.

3357 {
3358 //br for ee
3359 double psippee,psipee,jsiee,phiee,omegaee,rhoee;
3360 double kev2Gev=0.000001;
3361 psippee=0.262*kev2Gev;
3362 psipee =2.36*kev2Gev;
3363 jsiee =5.55*kev2Gev;
3364 phiee=4.266*0.001*2.954*0.0001;
3365 omegaee=0.6*kev2Gev;
3366 rhoee=7.04*kev2Gev;
3367 double s=_cms*_cms;
3368 double sigv=0;
3369 double mv=0;
3370 double widee=0;
3371 double xpi=12*3.1415926*3.1415926;
3372 if(mxL<=3.0957 && 3.0957<=mxH || mxL<=3.0983 && 3.0983<=mxH) {
3373 widee = jsiee;
3374 mv = 3.096916;
3375 }else if(mxL<=3.6847 && 3.6847<=mxH || mxL<=3.6873 && 3.6873<=mxH){
3376 widee = psipee;
3377 mv = 3.686109;
3378 }else if(mxL<=1.01906 && 1.01906<=mxH || mxL<= 1.01986 && 1.01986<=mxH){
3379 widee = phiee;
3380 mv = 1.01946;
3381 }else{return -1.0;}
3382
3383 if(_cms<=mv) return -1.;
3384 double xv = 1-mv*mv/s;
3385 sigv += xpi*widee/mv/s*Rad2(s,xv);
3386 double unic = 3.89379304*100000; //in unit nb
3387 return sigv*unic; //vaccum factor has included in the Rad2
3388}
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rh2ee DOUBLE PRECISION omepi DOUBLE PRECISION om3ee DOUBLE PRECISION phiee
Definition: RRes.h:37
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rhoee
Definition: RRes.h:34
double Rad2(double s, double x)
Definition: EvtConExc.cc:2600

◆ OutStaISR()

void EvtConExc::OutStaISR ( )

Definition at line 3803 of file EvtConExc.cc.

3803 {
3804 int ntot=myFisr.size();
3805 double mymu=0;
3806 for(int i=0;i<ntot;i++){mymu += myFisr[i];}
3807 mymu /= ntot;
3808 double mysig=0;
3809 for(int i=0;i<ntot;i++){ mysig += (myFisr[i]-mymu)*(myFisr[i]-mymu);}
3810 mysig /=ntot;
3811 mysig = sqrt(mysig);
3812 std::cout<<"========= Iteration times over ISR factor: "<<ntot<<std::endl;
3813 std::cout<<" ISR factor * VP factor= observedXS/BornXS: "<<mymu<<" + "<<mysig<<std::endl;
3814}

Referenced by ~EvtConExc().

◆ photonSampling()

bool EvtConExc::photonSampling ( EvtParticle part)

Definition at line 3326 of file EvtConExc.cc.

3326 {
3327 bool tagp,tagk;
3328 tagk=0;
3329 tagp=0;
3330 int nds = par->getNDaug();
3331 for(int i=0;i<par->getNDaug();i++){
3332 EvtId di=par->getDaug(i)->getId();
3333 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
3334 int pdgcode =EvtPDL::getStdHep( di );
3335 double alpha=0;
3336
3337 if(pdgcode==111 ||pdgcode==221 || pdgcode==331){//for photon
3338 alpha = 1;
3339 }else continue;
3340
3341 double angmax = 1+alpha;
3342 double costheta = cos(p4i.theta());
3343 double ang=1+alpha*costheta*costheta;
3344 double xratio = ang/angmax;
3345 double xi=EvtRandom::Flat(0.,1);
3346 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
3347 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
3348 if(xi>xratio) return false;
3349 }//loop over duaghters
3350 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
3351 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
3352
3353 return true;
3354}

Referenced by decay().

◆ PrintXS()

void EvtConExc::PrintXS ( double  mx)

Definition at line 3236 of file EvtConExc.cc.

3236 {//print xsection mode by mode at mx energy point
3237 std::vector<int> vmod; vmod.clear();
3238 int mn[83]={0,1,2,3,4,5,6,7,8,9,10,11,13,15,16,17,18,19,20,21,22,23,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
3239 50,51,
3240 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,85,
3241 90,91,92,93,94,95,96,
3242 74100,74101,74102,74103,74104,74105,74110};
3243 int mn2[84]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20,21,22,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
3244 50,51,
3245 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,85,
3246 90,91,92,93,94,95,96,
3247 74100,74101,74102,74103,74104,74105,74110};
3248
3249 if(_cms > 2.5){
3250 for(int i=0;i<83;i++){vmod.push_back(mn[i]);}
3251 }else{
3252 for(int i=0;i<84;i++){vmod.push_back(mn2[i]);}
3253 }
3254
3255 // for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
3256
3257 double xsum=0;
3258 double uscale = 1;
3259 double hxs=0;
3260 for(int i=0;i<vmod.size();i++){
3261 int mymode = vmod[i];
3262 delete myxsection;
3263 myxsection =new EvtXsection (mymode);
3264 init_mode(mymode);
3265 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3266 double mixs = uscale*myxsection->getXsection(mx);
3267 if(mymode ==74110){ hxs=mixs; continue;}
3268 xsum += mixs;
3269 //if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
3270 std::cout<<setw(5)<<mymode<<" : ";
3271 std::cout<<setw(10)<<mixs<<" nb ";
3272 for(int ii=0;ii<_ndaugs;ii++){std::cout<<setw(10)<<EvtPDL::name(daugs[ii])<<" ";}
3273 std::cout<<std::endl;
3274 }
3275 std::cout<<setw(6)<<"LUARLW "<<setw(10)<<hxs-xsum<<" nb "<<setw(10)<<" anything "<<std::endl;
3276 std::cout<<"Total hadron cross section here is "<<hxs<<" nb"<<std::endl;
3277}

Referenced by init().

◆ Rad1()

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

Definition at line 2588 of file EvtConExc.cc.

2588 { //first order correction
2589 //radiator correction upto second order, see Int. J. of Moder.Phys. A
2590 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
2591 double alpha = 1./137.;
2592 double pi=3.1415926;
2593 double me = 0.5*0.001; //e mass
2594 double sme = sqrt(s)/me;
2595 double L = 2*log(sme);
2596 double wsx = 2*alpha/pi/x*(L-1)*(1 - x + x*x/2 );
2597 return wsx;
2598}

Referenced by Rad1difXs().

◆ Rad1difXs()

double EvtConExc::Rad1difXs ( EvtParticle p)

Definition at line 2678 of file EvtConExc.cc.

2678 {// // first order xsection
2679 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2680 double summass = p->getDaug(0)->getP4().mass();
2681 for(int i=1;i<_ndaugs;i++){
2682 summass += p->getDaug(i)->getP4().mass();
2683 }
2684
2685 double cms = p->getP4().mass();
2686 double mrg = cms - summass;
2687 double pm= EvtRandom::Flat(0.,1);
2688 double mhs = pm*mrg + summass;
2689
2690 double s = cms * cms;
2691 double x = 1 - mhs*mhs/s;
2692 double wsx = Rad1(s,x);
2693
2694 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2695
2696 double xs = myxsection->getXsection(mhs);
2697 if(xs>XS_max){XS_max = xs;}
2698 double difxs = 2*mhs/cms/cms * wsx *xs;
2699
2700 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2701 differ *= mrg; //Jacobi factor for variable m_hds
2702 difxs *= mrg;
2703 return difxs;
2704}
double Rad1(double s, double x)
Definition: EvtConExc.cc:2588

Referenced by gamHXSection().

◆ Rad2()

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

Definition at line 2600 of file EvtConExc.cc.

2600 {
2601 //radiator correction upto second order, see Int. J. of Moder.Phys. A, hep-ph/9910523, Eq(8)
2602 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
2603 double alpha = 1./137.;
2604 double pi=3.1415926;
2605 double me = 0.5*0.001; //e mass
2606 double xi2 = 1.64493407;
2607 double xi3=1.2020569;
2608 double sme = sqrt(s)/me;
2609 double L = 2*log(sme);
2610 double beta = 2*alpha/pi*(L-1);
2611 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.;
2612 double ap = alpha/pi;
2613 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
2614 double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
2615 double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
2616 wsx = wsx + beta*beta/8. *wsx2;
2617 double mymx = sqrt(s*(1-x));
2618 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
2619 return wsx*getVP(mymx);//vph is vaccum polarization factor
2620}

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 2655 of file EvtConExc.cc.

2655 {// leading second order xsection
2656 double wsx = Rad2(s,x);
2657 double mhs = sqrt(s*(1-x));
2658 double xs = myxsection->getXsection(mhs);
2659 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
2660 if(xs>XS_max){XS_max = xs;}
2661 double difxs = wsx *xs;
2662 differ2 = wsx *(myxsection->getErr(mhs));
2663 return difxs;
2664}

◆ Rad2difXs() [2/3]

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

Definition at line 2666 of file EvtConExc.cc.

2666 {// leading second order xsection
2667 double wsx = Rad2(s,x);
2668 double mhs = sqrt(s*(1-x));
2669 double xs = myxsection->getXsection(mhs);
2670 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
2671 if(xs>XS_max){XS_max = xs;}
2672 double difxs = wsx *xs;
2673 differ2 = wsx *(myxsection->getErr(mhs));
2674 return difxs;
2675}

◆ Rad2difXs() [3/3]

double EvtConExc::Rad2difXs ( EvtParticle p)

Definition at line 2624 of file EvtConExc.cc.

2624 {// leading second order xsection
2625 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2626 double summass = p->getDaug(0)->getP4().mass();
2627 EvtVector4R v4p=p->getDaug(0)->getP4();
2628 for(int i=1;i<_ndaugs;i++){
2629 summass += p->getDaug(i)->getP4().mass();
2630 v4p += p->getDaug(i)->getP4();
2631 }
2632
2633 double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
2634 double cms = p->getP4().mass();
2635 double mrg = cms - summass;
2636 double pm= EvtRandom::Flat(0.,1);
2637 double mhs = pm*mrg + summass;
2638
2639 double s = cms * cms;
2640 double x = 2*Egam/cms;
2641 //double mhs = sqrt( s*(1-x) );
2642 double wsx = Rad2(s,x);
2643
2644 // std::cout<<"random : "<<pm<<std::endl;
2645
2646 double xs = myxsection->getXsection(mhs);
2647 if(xs>XS_max){XS_max = xs;}
2648 double difxs = 2*mhs/cms/cms * wsx *xs;
2649 differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2650 differ *= mrg; //Jacobi factor for variable m_hds
2651 difxs *= mrg;
2652 return difxs;
2653}

Referenced by gamHXSection(), and trapezoid().

◆ ReadVP()

void EvtConExc::ReadVP ( )

Definition at line 3446 of file EvtConExc.cc.

3446 {
3447 vpx.clear();vpr.clear();vpi.clear();
3448 int vpsize=4001;
3449 vpx.resize(vpsize);
3450 vpr.resize(vpsize);
3451 vpi.resize(vpsize);
3452 std::string locvp=getenv("BESEVTGENROOT");
3453 locvp +="/share/vp.dat";
3454 ifstream m_inputFile;
3455 m_inputFile.open(locvp.c_str());
3456 double xx,r1,i1;
3457 double x1,y1,y2;
3458 xx=0;
3459 int mytag=1;
3460 for(int i=0;i<4001;i++){
3461 m_inputFile>>vpx[i]>>vpr[i]>>vpi[i];
3462 }
3463
3464}

Referenced by init().

◆ resetResMass()

void EvtConExc::resetResMass ( )

Definition at line 3103 of file EvtConExc.cc.

3103 {
3104 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
3105 EvtId myres = EvtPDL::getId(std::string("J/psi"));
3106 EvtPDL::reSetMass(myres,mjsi);
3107 //EvtPDL::reSetWidth(myres,wjsi);
3108
3109 myres = EvtPDL::getId(std::string("psi(2S)"));
3110 EvtPDL::reSetMass(myres,mpsip);
3111 //EvtPDL::reSetWidth(myres,wpsip);
3112
3113 myres = EvtPDL::getId(std::string("psi(3770)"));
3114 EvtPDL::reSetMass(myres,mpsipp);
3115 //EvtPDL::reSetWidth(myres,wpsipp);
3116
3117 myres = EvtPDL::getId(std::string("phi"));
3118 EvtPDL::reSetMass(myres,mphi);
3119 //EvtPDL::reSetWidth(myres,wphi);
3120
3121 myres = EvtPDL::getId(std::string("omega"));
3122 EvtPDL::reSetMass(myres,momega);
3123 //EvtPDL::reSetWidth(myres,womega);
3124
3125 myres = EvtPDL::getId(std::string("rho0"));
3126 EvtPDL::reSetMass(myres,mrho0);
3127 //EvtPDL::reSetWidth(myres,wrho0);
3128
3129 myres = EvtPDL::getId(std::string("rho(3S)0"));
3130 EvtPDL::reSetMass(myres,mrho3s);
3131 //EvtPDL::reSetWidth(myres,wrho3s);
3132
3133 myres = EvtPDL::getId(std::string("omega(2S)"));
3134 EvtPDL::reSetMass(myres,momega2s);
3135 //EvtPDL::reSetWidth(myres,womega2s);
3136}

Referenced by decay().

◆ Ros_xs()

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

Definition at line 2738 of file EvtConExc.cc.

2738 {// in unit nb
2739 double pi=3.1415926;
2740 double s= mx*mx;
2741 double width = EvtPDL::getWidth(pid);
2742 double mass = EvtPDL::getMeanMass(pid);
2743 double xv = 1-mass*mass/s;
2744 double sigma = 12*pi*pi*bree*width/mass/s;
2745 sigma *= Rad2(s, xv);
2746 double nbar = 3.89379304*100000; // ! GeV-2 = 3.89379304*10^5 nbar
2747 return sigma*nbar;
2748}
double mass

Referenced by init().

◆ selectMass()

double EvtConExc::selectMass ( )

Definition at line 3400 of file EvtConExc.cc.

3400 {
3401 double pm,mhdL,mhds;
3402 pm = EvtRandom::Flat(0.,1);
3403 mhdL =_mhdL;
3404 mhds = pm*(_cms - mhdL)+mhdL; //non narrow resonance
3405 std::vector<double> sxs;
3406 for(int i=0;i<ISRID.size();i++){
3407 double ri=ISRXS[i]/AF[598];
3408 sxs.push_back(ri);
3409 }
3410 for(int j=0;j<sxs.size();j++){
3411 if(j>0) sxs[j] += sxs[j-1];
3412 }
3413 pm = EvtRandom::Flat(0.,1);
3414 if(pm<sxs[0]) {
3415 mhds = EvtPDL::getMass(ISRID[0]); //narrow resonance
3416 }
3417 for(int j=1;j<sxs.size();j++){
3418 double x0 = sxs[j-1];
3419 double x1 = sxs[j];
3420 if(x0<pm && pm<=x1) mhds = EvtPDL::getMass(ISRID[j]); //narrow resonance
3421 }
3422 return mhds;
3423}

◆ selectMode()

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

Definition at line 3023 of file EvtConExc.cc.

3023 {
3024 //first get xsection for each mode in vmod array
3025 //std::cout<<"vmod.size, mhds "<<vmod.size()<<" "<<mhds<<std::endl;
3026 std::vector<int>excmod;
3027 excmod.push_back(0);
3028 excmod.push_back(1);
3029 excmod.push_back(2);
3030 excmod.push_back(6);
3031 excmod.push_back(7);
3032 excmod.push_back(12);
3033 excmod.push_back(13);
3034 excmod.push_back(45);
3035 excmod.push_back(46);
3036 double uscale = 1;
3037 std::vector<double> vxs;vxs.clear();
3038 for (int i=0;i<vmod.size();i++){
3039 int imod = vmod[i];
3040 delete myxsection;
3041 myxsection =new EvtXsection (imod);
3042 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3043 //if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
3044
3045 double myxs=uscale*myxsection->getXsection(mhds);
3046
3047 if(i==0) {vxs.push_back(myxs);}
3048 else if(imod==74110){//for continuum process
3049 double xcon = myxs - vxs[i-1];
3050 if(xcon<0) {xcon=vxs[i-1];}else{xcon=myxs;}
3051 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
3052 vxs.push_back(xcon);
3053 }else{
3054 vxs.push_back(myxs+vxs[i-1]);
3055 }
3056 }//for_i_block
3057 //degugging
3058 // 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;}
3059
3060
3061 double totxs = vxs[vxs.size()-1];
3062 //if(totxs==0){std::cout<<"totalXs=0, vxs.size()= "<<vxs.size()<<std::endl;}
3063 int icount=0;
3064 int themode;
3065 mode_iter:
3066 double pm= EvtRandom::Flat(0.,1); //std::cout<<"pm= "<<pm<<" mhds= "<<mhds<<std::endl;
3067 if(pm <= vxs[0]/totxs) {
3068 themode= vmod[0];
3069 std::vector<EvtId> theDaug=get_mode(themode);
3070 EvtId daugs[100];
3071 for(int i=0;i<theDaug.size();i++){daugs[i]=theDaug[i];}
3072 int exmode=EvtDecayTable::inChannelList(theparent->getId(),theDaug.size(),daugs);
3073 if(_mode!=74110){return themode;}
3074 else if(exmode==-1){ return themode;}else{goto mycount;}
3075 }
3076
3077 for(int j=1;j<vxs.size();j++){
3078 double x0 = vxs[j-1]/totxs;
3079 double x1 = vxs[j]/totxs; //std::cout<<"j,x0,x1 "<<j<<" "<<x0<<" "<<x1<<std::endl;
3080 if(x0<pm && pm<=x1){
3081 themode=vmod[j];
3082 std::vector<EvtId> theDaug=get_mode(themode);
3083 EvtId daugs[100];
3084 for(int i=0;i<theDaug.size();i++){daugs[i]=theDaug[i];}
3085 int exmode=EvtDecayTable::inChannelList(theparent->getId(),theDaug.size(),daugs);
3086 if(_mode!=74110){return themode;}
3087 else if(exmode==-1){ return themode;} else{ break;}
3088 }
3089 }
3090
3091 mycount:
3092 icount++;
3093 if(icount<20000 ) goto mode_iter;
3094 //debugging
3095 //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;}
3096 std::cout<<"EvtConExc::selectMode can not find a mode with 20000 iteration for Mhadrons="<<mhds<<std::endl;
3097 return -10;
3098 //abort();
3099
3100}
std::vector< EvtId > get_mode(int mode)
Definition: EvtConExc.cc:1261
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)

Referenced by decay().

◆ SetP4()

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

Definition at line 2955 of file EvtConExc.cc.

2955 { //set the gamma and gamma* momentum according sampled results
2956 double pm= EvtRandom::Flat(0.,1);
2957 double phi=2*3.1415926*pm;
2958 double gamE = _cms/2*xeng;
2959 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2960 double px = gamE*sin(theta)*cos(phi);
2961 double py = gamE*sin(theta)*sin(phi);
2962 double pz = gamE*cos(theta);
2963 EvtVector4R p4[2];
2964 p4[0].set(gamE,px,py,pz);//ISR photon
2965 p4[1].set(hdrE,-px,-py,-pz); //mhdraons
2966
2967 for(int i=0;i<2;i++){
2968 EvtId myid = p->getDaug(i)->getId();
2969 p->getDaug(i)->init(myid,p4[i]);
2970 if(EvtPDL::name(myid)=="gamma*"){
2971 if( (_cms-mhdr)<0.002){
2972 EvtVector4R PG(_cms,0,0,0);
2973 p->getDaug(i)->init(myid,PG);
2974 }
2975 }
2976 }
2977}
double sin(const BesAngle a)

Referenced by decay().

◆ SetP4Rvalue()

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

Definition at line 2979 of file EvtConExc.cc.

2979 { //set the gamma and gamma* momentum according sampled results
2980 double pm= EvtRandom::Flat(0.,1);
2981 double phi=2*3.1415926*pm;
2982 double gamE = _cms/2*xeng;
2983 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2984 double px = gamE*sin(theta)*cos(phi);
2985 double py = gamE*sin(theta)*sin(phi);
2986 double pz = gamE*cos(theta);
2987 EvtVector4R p4[2];
2988 p4[0].set(hdrE,px,py,pz);//mhdraons
2989 p4[1].set(gamE,-px,-py,-pz); //ISR photon
2990 for(int i=0;i<2;i++){
2991 EvtId myid = p->getDaug(i)->getId();
2992 p->getDaug(i)->init(myid,p4[i]);
2993 }
2994}

Referenced by decay().

◆ showResMass()

void EvtConExc::showResMass ( )

Definition at line 3174 of file EvtConExc.cc.

3174 {
3175 std::cout<<"J/psi: "<<mjsi<<" "<<wjsi<<std::endl;
3176 std::cout<<"psipr: "<<mpsip<<" "<<wpsip<<std::endl;
3177 std::cout<<"psipp: "<<mpsipp<<" "<<wpsipp<<std::endl;
3178 std::cout<<"phi : "<<mphi<<" "<<wphi<<std::endl;
3179 std::cout<<"omega: "<<momega<<" "<<womega<<std::endl;
3180 std::cout<<"rho0 : "<<mrho0<<" "<<wrho0<<std::endl;
3181 std::cout<<"rho(3S)0: "<<mrho3s<<" "<<wrho3s<<std::endl;
3182 std::cout<<"omega(2S): "<<momega2s<<" "<<womega2s<<std::endl;
3183}

◆ SoftPhoton_xs()

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

Definition at line 2807 of file EvtConExc.cc.

2807 {
2808 double alpha = 1./137.;
2809 double pi=3.1415926;
2810 double me = 0.5*0.001; //e mass
2811 double xi2 = 1.64493407;
2812 double xi3=1.2020569;
2813 double sme = sqrt(s)/me;
2814 double L = 2*log(sme);
2815 double beta = 2*alpha/pi*(L-1);
2816 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.;
2817 double ap = alpha/pi;
2818 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
2819
2820 double beta2 = beta*beta;
2821 double b2 = b*b;
2822
2823 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 -
2824 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) +
2825 16*pow(beta,2)*Li2(b))/32. ;
2826 double mymx = sqrt(s*(1-b));
2827 //return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
2828 return xs*getVP(mymx); //vph :the vaccum polarzation factor
2829
2830}
double Li2(double x)
Definition: EvtConExc.cc:2832
const double b
Definition: slope.cxx:9

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

◆ split()

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

Definition at line 3695 of file EvtConExc.cc.

3696{
3697 std::string::size_type pos;
3698 std::vector<std::string> result;
3699 str+=pattern;
3700 int size=str.size();
3701
3702 for(int i=0; i<size; i++)
3703 {
3704 pos=str.find(pattern,i);
3705 if(pos<size)
3706 {
3707 std::string s=str.substr(i,pos-i);
3708 result.push_back(s);
3709 i=pos+pattern.size()-1;
3710 }
3711 }
3712 return result;
3713}

Referenced by InitPars().

◆ sumExc()

double EvtConExc::sumExc ( double  mx)

Definition at line 3199 of file EvtConExc.cc.

3199 {//summation all cross section of exclusive decays
3200 std::vector<int> vmod; vmod.clear();
3201 int mn[83]={0,1,2,3,4,5,6,7,8,9,10,11,13,15,16,17,18,19,20,21,22,23,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
3202 50,51,
3203 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,85,
3204 90,91,92,93,94,95,96,
3205 74100,74101,74102,74103,74104,74105,74110};
3206 int mn2[84]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20,21,22,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
3207 50,51,
3208 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,85,
3209 90,91,92,93,94,95,96,
3210 74100,74101,74102,74103,74104,74105,74110};
3211
3212 if(_cms > 2.5){
3213 for(int i=0;i<83;i++){vmod.push_back(mn[i]);}
3214 }else{
3215 for(int i=0;i<84;i++){vmod.push_back(mn2[i]);}
3216 }
3217
3218 // for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
3219
3220 double xsum=0;
3221 double uscale = 1;
3222 for(int i=0;i<vmod.size();i++){
3223 int mymode = vmod[i];
3224 if(mymode ==74110) continue;
3225 delete myxsection;
3226 myxsection =new EvtXsection (mymode);
3227 init_mode(mymode);
3228 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3229 xsum += uscale*myxsection->getXsection(mx);
3230 //if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
3231 }
3232 return xsum;
3233}

Referenced by init().

◆ trapezoid() [1/2]

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

Definition at line 3817 of file EvtConExc.cc.

3818{
3819 double xL=2*El/sqrt(s);
3820 double xH=2*Eh/sqrt(s);
3821 double sum = 0.0;
3822 double gaps = (xH-xL)/double(n); //interval
3823 for (int i = 0; i < n; i++)
3824 {
3825 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, staxsection) + Rad2difXs(s,xL + (i+1)*gaps,staxsection) );
3826 }
3827 return sum;
3828}

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 3830 of file EvtConExc.cc.

3831{
3832 double xL=2*El/sqrt(s);
3833 double xH=2*Eh/sqrt(s);
3834 double sum = 0.0;
3835 double gaps = (xH-xL)/double(n); //interval
3836 for (int i = 0; i < n; i++)
3837 {
3838 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, myxs) + Rad2difXs(s,xL + (i+1)*gaps,myxs) );
3839 }
3840 return sum;
3841}

◆ VP_sampling()

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

Definition at line 2577 of file EvtConExc.cc.

2577 {
2578 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2579 double theta=angles.getHelAng(1);
2580 double phi =angles.getHelAng(2);
2581 double costheta=cos(theta); //using helicity angles in parent system
2582
2583 double pm= EvtRandom::Flat(0.,1);
2584 double ang = 1 + costheta*costheta;
2585 if(pm< ang/2.) {return true;}else{return false;}
2586}

Referenced by hadron_angle_sampling().

◆ writeDecayTabel()

void EvtConExc::writeDecayTabel ( )

Definition at line 3533 of file EvtConExc.cc.

3533 {
3534 ofstream oa,ob;
3535 oa.open("_pkhr.dec");
3536 ob.open("obsxs.dat");
3537 //
3538 int im=getModeIndex(74110);
3539
3540 double xscon=VXS[im][599];
3541 //std::cout<<"tsize "<<tsize<<" "<<xscon<<" "<<VXS[70][599]<<std::endl;
3542 std::vector<int> Vmode;
3543 Vmode.push_back(6);//1:pi+pi-
3544 Vmode.push_back(13);//2:pi+pi-2pi0
3545 Vmode.push_back(12);//3:2(pi+pi-)
3546 Vmode.push_back(0);//4:ppbar
3547 Vmode.push_back(1);//5:nnbar
3548 Vmode.push_back(45);//6:K+K-
3549 Vmode.push_back(46);//7:K0K0bar
3550 Vmode.push_back(7);//8:pi+pi-pi0
3551 Vmode.push_back(2);//9:Lambda anti-Lambda
3552 Vmode.push_back(37);//10: pi+pi- eta
3553 std::vector<std::string> vmdl;
3554 vmdl.push_back("PHOKHARA_pipi");
3555 vmdl.push_back("PHOKHARA_pi0pi0pipi");
3556 vmdl.push_back("PHOKHARA_4pi");
3557 vmdl.push_back("PHOKHARA_ppbar");
3558 vmdl.push_back("PHOKHARA_nnbar");
3559 vmdl.push_back("PHOKHARA_KK");
3560 vmdl.push_back("PHOKHARA_K0K0");
3561 vmdl.push_back("PHOKHARA_pipipi0");
3562 vmdl.push_back("PHOKHARA_LLB");
3563 vmdl.push_back("PHOKHARA_pipieta");
3564
3565 ob<<"mode_index observed cross /nb"<<std::endl;
3566 for(int i=0;i<vmode.size();i++){
3567 ob<<vmode[i]<<" "<<VXS[getModeIndex(vmode[i])][599]<<std::endl;
3568//cout<<"vmode["<<i<<"] = "<<vmode[i]<<", VXS["<<getModeIndex(vmode[i])<<"][599] = "<<XS[getModeIndex(vmode[i])][599]<<std::endl; // wangwp
3569 }
3570
3571 for(int i=0;i<Vmode.size();i++){
3572 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
3573 xscon -= VXS[ getModeIndex(Vmode[i]) ][599];
3574 }
3575 //debugging
3576 for(int i=0;i<Vmode.size();i++){
3577 // std::cout<<Vmode[i]<<"-th mode: "<<VXS[getModeIndex(Vmode[i])][599]<<" "<<std::endl;
3578 }
3579
3580 oa<<"LundaPar PARJ(11)=0.611798"<<std::endl;
3581 oa<<"LundaPar PARJ(12)=7.92527E-12"<<std::endl;
3582 oa<<"LundaPar PARJ(14)=0.244495"<<std::endl;
3583 oa<<"LundaPar PARJ(15)=1.16573E-15"<<std::endl;
3584 oa<<"LundaPar PARJ(16)=0.436516"<<std::endl;
3585 oa<<"LundaPar PARJ(17)=0.530517"<<std::endl;
3586 oa<<"LundaPar PARJ(1)=0.0651577"<<std::endl;
3587 oa<<"LundaPar PARJ(2)=0.260378"<<std::endl;
3588 oa<<"LundaPar PARJ(21)=0.0664835"<<std::endl;
3589 oa<<"LundaPar RALPA(15)=0.576687"<<std::endl;
3590 oa<<"LundaPar RALPA(16)=0.364796"<<std::endl;
3591 oa<<"LundaPar RALPA(17)=3.19486E-06"<<std::endl;
3592 oa<<"noPhotos"<<std::endl;
3593 oa<<"Particle vpho "<<_cms<<" 0"<<std::endl;
3594 oa<<"Decay vpho"<<std::endl;
3595 oa<<"0 pi+ pi- PHSP ;"<<std::endl;
3596 oa<<"0 pi0 pi0 pi- pi+ PHSP ;"<<std::endl;
3597 oa<<"0 pi+ pi- pi- pi+ PHSP ;"<<std::endl;
3598 oa<<"0 anti-p- p+ PHSP ;"<<std::endl;
3599 oa<<"0 anti-n0 n0 PHSP ;"<<std::endl;
3600 oa<<"0 K+ K- PHSP ;"<<std::endl;
3601 oa<<"0 K_S0 K_L0 PHSP ;"<<std::endl;
3602 oa<<"0 pi+ pi- pi0 PHSP ;"<<std::endl;
3603 oa<<"0 Lambda0 anti-Lambda0 PHSP ;"<<std::endl;
3604 oa<<"0 pi+ pi- eta PHSP ;"<<std::endl;
3605 oa<<"0 gamma phi PHSP;"<<std::endl;
3606 oa<<"0 gamma rho0 PHSP;"<<std::endl;
3607 oa<<"0 gamma omega PHSP;"<<std::endl;
3608 oa<<xscon<<" ConExc 74110;"<<std::endl;
3609 for(int i=0;i<Vmode.size();i++){
3610 oa<<VXS[ getModeIndex(Vmode[i]) ][599]<<" "<<vmdl[i]<<" ;"<<std::endl;
3611 }
3612 oa<<"Enddecay"<<std::endl;
3613 oa<<"End"<<std::endl;
3614}

Referenced by init().

◆ xs_sampling() [1/2]

bool EvtConExc::xs_sampling ( double  xs)

Definition at line 2537 of file EvtConExc.cc.

2537 {
2538 double pm= EvtRandom::Flat(0.,1);
2539 //std::cout<<"Rdm= "<<pm<<std::endl;
2540 if(pm <xs/XS_max) {return true;} else {return false;}
2541}

◆ xs_sampling() [2/2]

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

Definition at line 2543 of file EvtConExc.cc.

2543 {
2544 double pm= EvtRandom::Flat(0.,1);
2545 // std::cout<<"Rdm= "<<pm<<std::endl;
2546 if(pm <xs/(xs0*1.1)) {return true;} else {return false;}
2547}

Member Data Documentation

◆ _cms

◆ conexcmode

int EvtConExc::conexcmode =-1
static

Definition at line 212 of file EvtConExc.hh.

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

◆ getMode

int EvtConExc::getMode
static

Definition at line 194 of file EvtConExc.hh.

Referenced by init().

◆ mlow

double EvtConExc::mlow =0
static

Definition at line 235 of file EvtConExc.hh.

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

◆ multimode

int EvtConExc::multimode =0
static

Definition at line 212 of file EvtConExc.hh.

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

◆ mup

double EvtConExc::mup =0
static

Definition at line 235 of file EvtConExc.hh.

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

◆ myxsection

◆ SetMthr

double EvtConExc::SetMthr =0
static

Definition at line 192 of file EvtConExc.hh.

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

◆ staxsection

EvtXsection * EvtConExc::staxsection
static

Definition at line 189 of file EvtConExc.hh.

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

◆ XS_max

double EvtConExc::XS_max
static

Definition at line 191 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: