Loading...
Searching...
No Matches
EvtConExc Class Reference
#include <EvtConExc.hh>
Inheritance diagram for EvtConExc:
Public Member Functions | |
| EvtConExc () | |
| virtual | ~EvtConExc () |
| void | getName (std::string &name) |
| EvtDecayBase * | clone () |
| void | initProbMax () |
| void | init () |
| void | init_Br_ee () |
| void | decay (EvtParticle *p) |
| void | init_mode (int mode) |
| double | gamHXSection (EvtParticle *p, double El, double Eh, int nmc=100000) |
| double | gamHXSection (double s, double El, double Eh, int nmc=100000) |
| double | gamHXSection (double El, double Eh) |
| double | gamHXSection (double El, double Eh, int mode) |
| double | gamHXSection_er (double El, double Eh) |
| void | findMaxXS (EvtParticle *p) |
| double | difgamXs (EvtParticle *p) |
| double | difgamXs (double mhds, double sintheta) |
| double | Mhad_sampling (double *x, double *y) |
| double | ISR_ang_integrate (double x, double theta) |
| double | ISR_ang_sampling (double x) |
| bool | hadron_angle_sampling (EvtVector4R ppi, EvtVector4R pcm) |
| void | SetP4 (EvtParticle *part, double mhdr, double xeng, double theta) |
| void | SetP4Rvalue (EvtParticle *part, double mhdr, double xeng, double theta) |
| bool | gam_sampling (EvtParticle *p) |
| bool | xs_sampling (double xs) |
| bool | xs_sampling (double xs, double xs1) |
| bool | baryon_sampling (EvtVector4R pcm, EvtVector4R pi) |
| bool | meson_sampling (EvtVector4R pcm, EvtVector4R pi) |
| bool | VP_sampling (EvtVector4R pcm, EvtVector4R pi) |
| bool | angularSampling (EvtParticle *part) |
| bool | photonSampling (EvtParticle *part) |
| double | baryonAng (double mx) |
| double | Rad1 (double s, double x) |
| double | Rad2 (double s, double x) |
| double | Rad2difXs (EvtParticle *p) |
| double | Rad2difXs (double s, double x) |
| double | Rad2difXs (double s, double x, EvtXsection *myxsection) |
| double | Rad1difXs (EvtParticle *p) |
| double | Ros_xs (double mx, double bree, EvtId pid) |
| double | Li2 (double x) |
| double | SoftPhoton_xs (double s, double b) |
| double | lgr (double *x, double *y, int n, double t) |
| bool | islgr (double *x, double *y, int n, double t) |
| double | LLr (double *x, double *y, int n, double t) |
| int | selectMode (std::vector< int > vmod, double mhds) |
| void | findMaxXS (double mhds) |
| bool | gam_sampling (double mhds, double sintheta) |
| void | resetResMass () |
| void | getResMass () |
| bool | checkdecay (EvtParticle *p) |
| double | sumExc (double mx) |
| void | PrintXS (double mx) |
| void | showResMass () |
| int | getNdaugs () |
| EvtId * | getDaugId () |
| int | getSelectedMode () |
| double | narrowRXS (double mxL, double mxH) |
| double | selectMass () |
| double | addNarrowRXS (double mhi, double binwidth) |
| void | ReadVP () |
| double | getVP (double cms) |
| void | mk_VXS (double Esig, double Egamcut, double EgamH, int midx) |
| int | get_mode_index (int mode) |
| double | getObsXsection (double mhds, int mode) |
| double | Egam2Mhds (double Egam) |
| std::vector< EvtId > | get_mode (int mode) |
| void | writeDecayTabel () |
| void | checkEvtRatio () |
| int | getModeIndex (int m) |
| std::string | commandName () |
| void | command (std::string cmd) |
| std::vector< std::string > | split (std::string str, std::string pattern) |
| void | InitPars () |
| double | energySpread (double mu, double sigma) |
| void | calAF (double myecms) |
| void | OutStaISR () |
| double | trapezoid (double s, double a, double b, int n) |
| double | trapezoid (double s, double El, double Eh, int n, EvtXsection *myxs) |
| Public Member Functions inherited from EvtDecayIncoherent | |
| void | makeDecay (EvtParticle *p) |
| virtual | ~EvtDecayIncoherent () |
| void | setDaughterSpinDensity (int daughter) |
| int | isDaughterSpinDensitySet (int daughter) |
| Public Member Functions inherited from EvtDecayBase | |
| double | getProbMax (double prob) |
| double | resetProbMax (double prob) |
| EvtDecayBase () | |
| virtual | ~EvtDecayBase () |
| virtual bool | matchingDecay (const EvtDecayBase &other) const |
| EvtId | getParentId () |
| double | getBranchingFraction () |
| void | disableCheckQ () |
| void | checkQ () |
| int | getNDaug () |
| EvtId * | getDaugs () |
| EvtId | getDaug (int i) |
| int | getNArg () |
| int | getPHOTOS () |
| void | setPHOTOS () |
| void | setVerbose () |
| void | setSummary () |
| double * | getArgs () |
| std::string * | getArgsStr () |
| double | getArg (int j) |
| std::string | getArgStr (int j) |
| std::string | getModelName () |
| int | getDSum () |
| int | summary () |
| int | verbose () |
| void | saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr) |
| void | printSummary () |
| void | setProbMax (double prbmx) |
| void | noProbMax () |
| void | checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1) |
| void | checkNDaug (int d1, int d2=-1) |
| void | checkSpinParent (EvtSpinType::spintype sp) |
| void | checkSpinDaughter (int d1, EvtSpinType::spintype sp) |
| virtual int | nRealDaughters () |
Static Public Attributes | |
| static EvtXsection * | myxsection |
| static EvtXsection * | staxsection |
| static double | _cms |
| static double | XS_max |
| static double | SetMthr =0 |
| static int | getMode |
| static int | conexcmode =-1 |
| static int | multimode =0 |
| static double | mlow =0 |
| static double | mup =0 |
Additional Inherited Members | |
| Static Public Member Functions inherited from EvtDecayBase | |
| static void | findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10]) |
| static void | findMass (EvtParticle *p) |
| static double | findMaxMass (EvtParticle *p) |
| Protected Member Functions inherited from EvtDecayBase | |
| bool | daugsDecayedByParentModel () |
| Protected Attributes inherited from EvtDecayBase | |
| bool | _daugsDecayedByParentModel |
Detailed Description
Definition at line 138 of file EvtConExc.hh.
Constructor & Destructor Documentation
◆ EvtConExc()
|
inline |
◆ ~EvtConExc()
|
virtual |
Definition at line 207 of file EvtConExc.cc.
207 {
208// if(_mode==74110)checkEvtRatio();
210 if(mydbg){
211 // xs->Write();
212 myfile->Write();
213 }
214
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}
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();
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 */
3304 double xratio = ang/angmax;
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}
Referenced by decay().
◆ baryon_sampling()
| bool EvtConExc::baryon_sampling | ( | EvtVector4R | pcm, |
| EvtVector4R | pi ) |
Definition at line 2549 of file EvtConExc.cc.
2549 {
2551 double theta=angles.getHelAng(1);
2552 double phi =angles.getHelAng(2);
2554
2555 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
2556 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
2561 double myang;
2563 if(pm< ang/myang) {return true;}else{return false;}
2564}
Referenced by hadron_angle_sampling().
◆ baryonAng()
| double EvtConExc::baryonAng | ( | double | mx | ) |
Definition at line 3316 of file EvtConExc.cc.
3316 {
3318 double u = 0.938/mx;
3319 u = u*u;
3320 double u2 = (1+u)*(1+u);
3321 double uu = u*(1+6*u);
3324}
Referenced by baryon_sampling().
◆ calAF()
| void EvtConExc::calAF | ( | double | myecms | ) |
Definition at line 3752 of file EvtConExc.cc.
3752 {
3753
3754
3756 double Esig = 0.0001; //sigularity engy
3759
3762
3763 int nint=50;
3764 int nmc= 1000;
3766 //--- sigularity point x from 0 to 0.0001GeV
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;
3777
3778 AA[i]=(Eg1+Eg0)/2;
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++){
3793 }
3794 //--debugging
3797 double fisr=AF[599]/bornXS;
3798 myFisr.push_back(fisr);
3799 //std::cout<<"fisr= "<<fisr<<std::endl;
3800}
double trapezoid(double s, double a, double b, int n)
Definition EvtConExc.cc:3817
Referenced by decay().
◆ checkdecay()
| bool EvtConExc::checkdecay | ( | EvtParticle * | p | ) |
Definition at line 3185 of file EvtConExc.cc.
3185 {
3187 double massflag=1;
3188 for(int i=0;i<nson;i++){
3191 }
3192 if(massflag==0){
3193 std::cout<<"Zero mass!"<<std::endl;
3194 return 0;
3195 }else{return 1;}
3196}
Referenced by decay().
◆ checkEvtRatio()
| void EvtConExc::checkEvtRatio | ( | ) |
Definition at line 3618 of file EvtConExc.cc.
3618 {
3621 //
3623 double xscon=VXS[im][599];
3624 double xscon0=xscon;
3628 for(int i=0;i<vmode.size();i++){
3629 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
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++){
3638 if(vmode[i]==74110) continue;
3640 }
3642
3643
3644}
◆ clone()
|
virtual |
◆ command()
|
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()
|
virtual |
Reimplemented from EvtDecayBase.
Definition at line 3654 of file EvtConExc.cc.
3654 {
3655
3656 return std::string("ConExcPar");
3657
3658}
◆ decay()
|
virtual |
Implements EvtDecayBase.
Definition at line 1783 of file EvtConExc.cc.
1783 {
1784 //std::cout<<"_cms= "<<_cms<<" mode= "<<_mode<<std::endl;
1788 abort();
1789 }
1790 //-- to read _cms from database or to consider beam energy spread
1792 bool beamflag=0;
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
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
1837
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;
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
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);
1870 di->setMass(xmi);
1871 totMass+=xmi;
1872 }
1874 p->initializePhaseSpace( _ndaugs , daugs);
1876 vhds = p->getP4();
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);
1891 di->setMass(xmi);
1892 totMass+=xmi;
1893 }
1894 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
1896 p->initializePhaseSpace(2,mydaugs);
1897
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;
1910
1911 if(mydbg) mass2=mhds;
1912
1913 //generating events
1914 ModeSelection:
1915 EvtGlobalSet::ConExcPythia=1;
1916 mymode=selectMode(vmod,mhds);
1917 if(mymode==-10) goto Sampling_mhds;
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;
1924 init_mode(mymode);
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);
1952 di->setMass(xmi);
1953 totMass+=xmi;
1954 }
1956 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
1959 //set the ISR photon angular distribution
1961
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);
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;
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
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);
2022 di->setMass(xmi);
2023 totMass+=xmi;
2024 }
2026 p->initializePhaseSpace( _ndaugs , daugs);
2028 vhds = p->getP4();
2029 bool byn_ang;
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);
2050 di->setMass(xmi);
2051 totMass+=xmi;
2052 }
2053 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
2055 p->initializePhaseSpace(2,mydaugs);
2056
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;
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;
2077 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
2078 _selectedMode = mymode;
2081 init_mode(mymode);
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);
2109 di->setMass(xmi);
2110 totMass+=xmi;
2111 }
2113 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2116 //set the ISR photon angular distribution
2118
2119 //--debugging
2120 //p->printTree();
2121
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);
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){
2161
2162 if(pm <_xs0/(_xs0 + _xs1)){//with a very soft ISR photon
2164 mhds = _cms-0.0001;
2165 }else{
2166 mhds=Mhad_sampling(MH,AF);
2167 EvtPDL::reSetMass(gid,mhds);
2168 }
2169 //debugging
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();
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();}
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
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);
2208 di->setMass(xmi);
2209 summassx += xmi;
2210 //std::cout<<"PID and mass "<<di->getId()<<" "<<xmi<<std::endl;
2211 }
2213 angular_hadron:
2214 p->initializePhaseSpace(_ndaugs,daugs);
2215 bool byn_ang;
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
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*
2234 EvtPDL::reSetMass(gid,mhdr);
2235 int ndaugs =2;
2236 EvtId mydaugs[2];
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);
2249 //if(EvtPDL::name(di->getId())=="vhdr") di->setVectorSpinDensity();
2250 di->setMass(xmi);
2251 }
2252 p->initializePhaseSpace(2,mydaugs);
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}
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition EvtConExc.cc:2979
int selectMode(std::vector< int > vmod, double mhds)
Definition EvtConExc.cc:3023
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition EvtConExc.cc:2278
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition EvtConExc.cc:2955
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1070
◆ 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)
3014 double sin2theta = sintheta*sintheta;
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
2499 for(int i=1;i<_ndaugs;i++){
2501 }
2504 double egam = pgam.get(0);
2505 double sin2theta = pgam.get(3)/pgam.d3mag();
2506 sin2theta = 1-sin2theta*sin2theta;
2507
2513 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2514
2516 double difxs = 2*mhs/cms/cms * wsx *xs;
2517 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2518 return difxs;
2519
2520}
Referenced by findMaxXS(), gam_sampling(), gam_sampling(), and gamHXSection().
◆ Egam2Mhds()
| double EvtConExc::Egam2Mhds | ( | double | Egam | ) |
◆ energySpread()
| double EvtConExc::energySpread | ( | double | mu, |
| double | sigma ) |
Definition at line 3734 of file EvtConExc.cc.
3734 {
3735 //mu: mean value in Gaussian
3736 //sigma: variance in Gaussian
3737 rloop:
3739 double ri=0;
3742 ri += pm;
3743 }
3748 if(xsig<mx0 || xsig>mx1) goto rloop;
3749 return xsig;
3750}
◆ 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++){
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;
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);
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
2486 EvtVector4R p4gam = gamH->getDaug(_ndaugs)->getP4();
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
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition EvtParticleFactory.cc:46
◆ gam_sampling() [1/2]
| bool EvtConExc::gam_sampling | ( | double | mhds, |
| double | sintheta ) |
Definition at line 2530 of file EvtConExc.cc.
2530 {
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
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;
2402 double Rad2Xs;
2403 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2404 if(Rad2Xs==0){
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 dgauss_(double(*fun)(double *), double *, double *, double *)
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:100
◆ 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;
2420 double Rad2Xs;
2421 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2422 if(Rad2Xs==0){
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;
2378 for(int i=0;i<nmc;i++){//It is found that the narrow resonance can not included in this MC integartion
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;
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);
2330 di->setMass(xmi);
2331 }
2332 loop:
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
2346 Rad2Xs += Rad2difXs( gamH );
2347 if(thexs>maxXS) {maxXS=thexs;}
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}
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;
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){
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}
◆ 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;
1268 }else if(mode==1){
1269 _ndaugs =2;
1272 }else if(mode==2){
1273 _ndaugs =2;
1276 }else if(mode==3){
1277 _ndaugs =2;
1280 }else if(mode==4){
1281 _ndaugs =2;
1284 }else if(mode==5){
1285 _ndaugs =2;
1288 }else if(mode ==6){
1289 _ndaugs =2;
1292 }else if(mode==7){
1293 _ndaugs =3;
1297 }else if(mode==8){
1298 _ndaugs =3;
1302 }else if(mode==9){
1303 _ndaugs =3;
1307 }else if(mode==10){
1308 _ndaugs =3;
1312 }else if(mode==11){
1313 _ndaugs =3;
1317 }else if(mode==12){
1318 _ndaugs =4;
1323 }else if(mode==13){
1324 _ndaugs =4;
1329 }else if(mode==14){
1330 _ndaugs =4;
1335 }else if(mode==15){
1336 _ndaugs =4;
1341 }else if(mode==16){
1342 _ndaugs =4;
1347 }else if(mode==17){
1348 _ndaugs =5;
1354 }else if(mode==18){
1355 _ndaugs =5;
1361 }else if(mode==19){
1362 _ndaugs =5;
1368 }else if(mode==20){
1369 _ndaugs =5;
1375 }else if(mode==21){
1376 _ndaugs =6;
1383 }else if(mode==22){
1384 _ndaugs =6;
1391 }else if(mode==23){
1392 _ndaugs =2;
1395 }else if(mode==24){
1396 _ndaugs =2;
1399 }else if(mode==25){
1400 _ndaugs =2;
1403 }else if(mode==26){
1404 _ndaugs =2;
1407 }else if(mode==27){
1408 _ndaugs =2;
1411 }else if(mode==28){
1412 _ndaugs =3;
1416 }else if(mode==29){
1417 _ndaugs =3;
1421 }else if(mode==30){
1422 _ndaugs =3;
1426 }else if(mode==31){
1427 _ndaugs =3;
1431 }else if(mode==32){
1432 _ndaugs =3;
1436 }else if(mode==33){
1437 _ndaugs =3;
1441 }else if(mode==34){
1442 _ndaugs =3;
1446 }else if(mode==35){
1447 _ndaugs =3;
1451 }else if(mode==36){
1452 _ndaugs =2;
1455 }else if(mode==37){
1456 _ndaugs =3;
1460 }else if(mode==38){
1461 _ndaugs =3;
1465 }else if(mode==39){
1466 _ndaugs =2;
1469 }else if(mode==40){
1470 _ndaugs =3;
1474 }else if(mode==41){
1475 _ndaugs =3;
1479 }else if(mode==42){
1480 _ndaugs =3;
1484 }else if(mode==43){
1485 _ndaugs =4;
1490 }else if(mode==44){
1491 _ndaugs =2;
1494 }else if(mode==45){
1495 _ndaugs =2;
1498 }else if(mode==46){
1499 _ndaugs =2;
1502 }else if(mode==47){
1503 _ndaugs =2;
1506 }else if(mode==48){
1507 _ndaugs =2;
1510 }else if(mode==49){
1511 _ndaugs =3;
1515 }else if(mode==50){
1516 _ndaugs =3;
1520 }else if(mode==51){
1521 _ndaugs =3;
1525 }else if(mode==52){
1526 _ndaugs =2;
1529 }else if(mode==59){
1530 _ndaugs =3;
1534 }else if(mode==60){
1535 _ndaugs =3;
1539 }else if(mode==61){
1540 _ndaugs =3;
1544 }else if(mode==62){
1545 _ndaugs =3;
1549 }else if(mode==63){
1550 _ndaugs =3;
1554 }else if(mode==64){
1555 _ndaugs =3;
1559 }else if(mode==65){
1560 _ndaugs =3;
1564 }else if(mode==66){
1565 _ndaugs =3;
1569 }else if(mode==67){
1570 _ndaugs =2;
1573 }else if(mode==68){
1574 _ndaugs =2;
1577 }else if(mode==69){
1578 _ndaugs =2;
1581 }else if(mode==70){
1582 _ndaugs =2;
1585 }else if(mode==71){
1586 _ndaugs =2;
1589 }else if(mode==72){
1590 _ndaugs =2;
1593 }else if(mode==73){
1594 _ndaugs =2;
1597 }else if(mode==74){
1598 _ndaugs =2;
1601 }else if(mode==75){
1602 _ndaugs =3;
1606 }else if(mode==76){
1607 _ndaugs =3;
1611 }else if(mode==77){
1612 _ndaugs =3;
1616 }else if(mode==78){
1617 _ndaugs =3;
1621 }else if(mode==79){
1622 _ndaugs =3;
1626 }else if(mode==80){
1627 _ndaugs =2;
1630 }else if(mode==81){
1631 _ndaugs =3;
1635 }else if(mode==82){
1636 _ndaugs =3;
1640 }else if(mode==83){
1641 _ndaugs =3;
1645 }else if(mode==84){
1646 _ndaugs =3;
1650 }else if(mode==85){
1651 _ndaugs =2;
1654 }else if(mode==90){
1655 _ndaugs =3;
1659 }else if(mode==91){
1660 _ndaugs =3;
1664 }else if(mode==92){
1665 _ndaugs =3;
1669 }else if(mode==93){
1670 _ndaugs =2;
1673 }else if(mode==94){
1674 _ndaugs =2;
1677 }else if(mode==95){
1678 _ndaugs =2;
1681 }else if(mode==96){
1682 _ndaugs =2;
1685 }else if(mode==97){
1686 _ndaugs =3;
1690 }else if(mode==111){
1691 _ndaugs =2;
1694 }else if(mode==112){
1695 _ndaugs =2;
1698 }else if(mode==74100){
1699 _ndaugs =2;
1702 }else if(mode==74101){
1703 _ndaugs =2;
1706 }else if(mode==74102){
1707 _ndaugs =2;
1710 }else if(mode==74103){
1711 _ndaugs =2;
1714 }else if(mode==74104){
1715 _ndaugs =2;
1718 }else if(mode==74105){
1719 _ndaugs =2;
1722 }else if(mode==74106){
1723 _ndaugs =2;
1726 }else if(mode==74107){
1727 _ndaugs =2;
1730 }else if(mode==74110){
1731 _ndaugs =2;
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]; }
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;
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;
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;
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}
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()
|
inline |
Definition at line 210 of file EvtConExc.hh.
210{return daugs;}
◆ 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()
|
virtual |
Implements EvtDecayBase.
Definition at line 243 of file EvtConExc.cc.
243 {
244
245 model_name="ConExc";
246
247}
◆ getNdaugs()
|
inline |
Definition at line 209 of file EvtConExc.hh.
209{return _ndaugs;}
◆ 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.;
3519 for(int i=0;i<600;i++){
3520 if( (AA[i]-hwid)<=Egam && Egam<(AA[i]+hwid) ) return VXS[idx][i];
3521 }
3522
3524 abort();
3525}
◆ getResMass()
| void EvtConExc::getResMass | ( | ) |
Definition at line 3138 of file EvtConExc.cc.
3138 {
3140 mjsi = EvtPDL::getMeanMass(myres);
3141 wjsi = EvtPDL::getWidth(myres);
3142
3144 mpsip= EvtPDL::getMeanMass(myres);
3145 wpsip= EvtPDL::getWidth(myres);
3146
3148 mpsipp= EvtPDL::getMeanMass(myres);
3149 wpsipp= EvtPDL::getWidth(myres);
3150
3152 mphi = EvtPDL::getMeanMass(myres);
3153 wphi = EvtPDL::getWidth(myres);
3154
3156 momega= EvtPDL::getMeanMass(myres);
3157 womega= EvtPDL::getWidth(myres);
3158
3160 mrho0 = EvtPDL::getMeanMass(myres);
3161 wrho0 = EvtPDL::getWidth(myres);
3162
3164 mrho3s= EvtPDL::getMeanMass(myres);
3165 wrho3s= EvtPDL::getWidth(myres);
3166
3167
3169 momega2s= EvtPDL::getMeanMass(myres);
3170 womega2s= EvtPDL::getWidth(myres);
3171
3172}
Referenced by init().
◆ 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);
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;
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
2284 return byn_ang;
2285 }else if(_mode==6 || _mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
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){
2290 return msn_ang;
2291 }else if(_mode==-2){
2296 return msn_ang;
2297 }else if(type0==EvtSpinType::VECTOR &&type1==EvtSpinType::SCALAR || type1==EvtSpinType::VECTOR &&type0==EvtSpinType::SCALAR){
2299 return msn_ang;
2300 }
2301 }
2302 return true;
2303}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:75
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:2566
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:2549
Referenced by decay().
◆ 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;
293 if(getDaugs()[0]==EvtPDL::getId("gamma")&&getDaugs()[1]==EvtPDL::getId("gamma*") || getDaugs()[0]==EvtPDL::getId("gamma*")&&getDaugs()[1]==EvtPDL::getId("gamma")) VISR=1;
295 }
296
297 // cmsspread=0.0015; //CMC energy spread
298 testflag=0;
299 getResMass();
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;
311 radflag=0;mydbg=false;
312 _mode = getArg(0);
313 nson = getNArg()-1;
316 _mode = getArg(0);
317 multimode=1;
318 vmd.clear();
320 std::cout<<" multi-exclusive mode "<<std::endl;
322 std::cout<<std::endl;
323 if(vmd[vmd.size()-1]==74110){std::cout<<"74110 is not allowd for multimode simulation"<<std::endl;abort();}
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
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
355 EvtPDL::reSetWidth(parentId,0.0);
357 std::cout<<"parent mass = "<<parentMass<<std::endl;
358
359
360 EvtVector4R p_init(parentMass,0.0,0.0,0.0);
361 EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
362 theparent = p;
365 if(_mode==-100) {
366 myxsection->setModes(vmd);
367 myxsection->ini_data_multimode();
368 staxsection->setModes(vmd);
369 staxsection->ini_data_multimode();
370 }
371 double mth=0;
373 if(_mode ==-1){
374 myxsection->setBW(pdgcode);
375 staxsection->setBW(pdgcode);
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();
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
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
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();
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++){
413 if(mx<mR || _mode !=74110) continue;
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
429 int nmc=1000000;
431 _er0 = _er1; // stored when calculate _xs0
432 std::cout<<"_er0= "<<_er0<<std::endl;
434 std::cout<<"_xs1= "<<_xs0<<std::endl;
435 double xs1_err = _er1;
436
437
438 //--- sigularity point x from 0 to 0.0001GeV
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);
455 AA[i]=(Eg1+Eg0)/2;
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
473 for(int i=0;i<vmode.size();i++){
475 }
477 //after mk_VXS, restore the myxsection to _mode
478 if(_mode !=-100){
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++){
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
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){
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){
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;
614 TGraphErrors *Gdt = new TGraphErrors(nd,xx,yy,xer,yer);
615
619 for(int i=0;i<600;i++){
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++){
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 gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition EvtConExc.cc:2306
void mk_VXS(double Esig, double Egamcut, double EgamH, int midx)
Definition EvtConExc.cc:3466
◆ 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;
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;
756 }else if(mode==1){
757 _ndaugs =2;
760 }else if(mode==2){
761 _ndaugs =2;
764 }else if(mode==3){
765 _ndaugs =2;
768 }else if(mode==4){
769 _ndaugs =2;
772 }else if(mode==5){
773 _ndaugs =2;
776 }else if(mode==6){
777 _ndaugs =2;
780 }else if(mode==7){
781 _ndaugs =3;
785 }else if(mode==8){
786 _ndaugs =3;
790 }else if(mode==9){
791 _ndaugs =3;
795 }else if(mode==10){
796 _ndaugs =3;
800 }else if(mode==11){
801 _ndaugs =3;
805 }else if(mode==12){
806 _ndaugs =4;
811 }else if(mode==13){
812 _ndaugs =4;
817 }else if(mode==14){
818 _ndaugs =4;
823 }else if(mode==15){
824 _ndaugs =4;
829 }else if(mode==16){
830 _ndaugs =4;
835 }else if(mode==17){
836 _ndaugs =5;
842 }else if(mode==18){
843 _ndaugs =5;
849 }else if(mode==19){
850 _ndaugs =5;
856 }else if(mode==20){
857 _ndaugs =5;
863 }else if(mode==21){
864 _ndaugs =6;
871 }else if(mode==22){
872 _ndaugs =6;
879 }else if(mode==23){
880 _ndaugs =2;
883 }else if(mode==24){
884 _ndaugs =2;
887 }else if(mode==25){
888 _ndaugs =2;
891 }else if(mode==26){
892 _ndaugs =2;
895 }else if(mode==27){
896 _ndaugs =2;
899 }else if(mode==28){
900 _ndaugs =3;
904 }else if(mode==29){
905 _ndaugs =3;
909 }else if(mode==30){
910 _ndaugs =3;
914 }else if(mode==31){
915 _ndaugs =3;
919 }else if(mode==32){
920 _ndaugs =3;
924 }else if(mode==33){
925 _ndaugs =3;
929 }else if(mode==34){
930 _ndaugs =3;
934 }else if(mode==35){
935 _ndaugs =3;
939 }else if(mode==36){
940 _ndaugs =2;
943 }else if(mode==37){
944 _ndaugs =3;
948 }else if(mode==38){
949 _ndaugs =3;
953 }else if(mode==39){
954 _ndaugs =2;
957 }else if(mode==40){
958 _ndaugs =3;
962 }else if(mode==41){
963 _ndaugs =3;
967 }else if(mode==42){
968 _ndaugs =3;
972 }else if(mode==43){
973 _ndaugs =4;
978 }else if(mode==44){
979 _ndaugs =2;
982 }else if(mode==45){
983 _ndaugs =2;
986 }else if(mode==46){
987 _ndaugs =2;
990 }else if(mode==47){
991 _ndaugs =2;
994 }else if(mode==48){
995 _ndaugs =2;
998 }else if(mode==49){
999 _ndaugs =3;
1003 }else if(mode==50){
1004 _ndaugs =3;
1008 }else if(mode==51){
1009 _ndaugs =3;
1013 }else if(mode==52){
1014 _ndaugs =2;
1017 }else if(mode==59){
1018 _ndaugs =3;
1022 }else if(mode==60){
1023 _ndaugs =3;
1027 }else if(mode==61){
1028 _ndaugs =3;
1032 }else if(mode==62){
1033 _ndaugs =3;
1037 }else if(mode==63){
1038 _ndaugs =3;
1042 }else if(mode==64){
1043 _ndaugs =3;
1047 }else if(mode==65){
1048 _ndaugs =3;
1052 }else if(mode==66){
1053 _ndaugs =3;
1057 }else if(mode==67){
1058 _ndaugs =2;
1061 }else if(mode==68){
1062 _ndaugs =2;
1065 }else if(mode==69){
1066 _ndaugs =2;
1069 }else if(mode==70){
1070 _ndaugs =2;
1073 }else if(mode==71){
1074 _ndaugs =2;
1077 }else if(mode==72){
1078 _ndaugs =2;
1081 }else if(mode==73){
1082 _ndaugs =2;
1085 }else if(mode==74){
1086 _ndaugs =2;
1089 }else if(mode==75){
1090 _ndaugs =3;
1094 }else if(mode==76){
1095 _ndaugs =3;
1099 }else if(mode==77){
1100 _ndaugs =3;
1104 }else if(mode==78){
1105 _ndaugs =3;
1109 }else if(mode==79){
1110 _ndaugs =3;
1114 }else if(mode==80){
1115 _ndaugs =2;
1118 }else if(mode==81){
1119 _ndaugs =3;
1123 }else if(mode==82){
1124 _ndaugs =3;
1128 }else if(mode==83){
1129 _ndaugs =3;
1133 }else if(mode==84){
1134 _ndaugs =3;
1138 }else if(mode==85){
1139 _ndaugs =2;
1142 }else if(mode==90){
1143 _ndaugs =3;
1147 }else if(mode==91){
1148 _ndaugs =3;
1152 }else if(mode==92){
1153 _ndaugs =3;
1157 }else if(mode==93){
1158 _ndaugs =2;
1161 }else if(mode==94){
1162 _ndaugs =2;
1165 }else if(mode==95){
1166 _ndaugs =2;
1169 }else if(mode==96){
1170 _ndaugs =2;
1173 }else if(mode==97){
1174 _ndaugs =3;
1178 }else if(mode==111){
1179 _ndaugs =2;
1182 }else if(mode==112){
1183 _ndaugs =2;
1186 }else if(mode==74100){
1187 _ndaugs =1;
1189 }else if(mode==74101){
1190 _ndaugs =1;
1192 }else if(mode==74102){
1193 _ndaugs =1;
1195 }else if(mode==74103){
1196 _ndaugs =1;
1198 }else if(mode==74104){
1199 _ndaugs =1;
1201 }else if(mode==74105){
1202 _ndaugs =1;
1204 }else if(mode==74106){
1205 _ndaugs =1;
1207 }else if(mode==74107){
1208 _ndaugs =1;
1210 }else if(mode==74110){
1211 _ndaugs =1;
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]; }
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;
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;
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;
1246 }
1247
1248 double fmth=0;
1249 for(int i=0;i<_ndaugs;i++){
1251 fmth +=xmi;
1252 }
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}
◆ 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);
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
Referenced by init().
◆ initProbMax()
|
virtual |
Reimplemented from EvtDecayBase.
Definition at line 1777 of file EvtConExc.cc.
◆ 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;
2859 }
2861 if(n0==-1) {return 0;}else{
2866 }
2867}
◆ 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)
2917 double cos2=costheta*costheta;
2918 double sin2=1-cos2;
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));
2926 double hq=(L-1)/2.+hq1+hq2+hq3;
2927 hq /= hpi; //Eq (11) in arXif:hep-ph/9910523
2928 return hq;
2929}
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
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;
2943 yy[j]=ISR_ang_integrate(x,xx[j]);
2944 }
2946 int mypt=1;
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
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;
2844 }
2845 if(n0==-1) {return 0.0;}else{
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!
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;
2875 }
2876 if(n0==-1) {return 0.0;}else{
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 {
2568 double theta=angles.getHelAng(1);
2569 double phi =angles.getHelAng(2);
2571
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
2891 double xrt=1;
2892 int mybin=1;
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;}
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;
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
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);
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);
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
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;
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
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
◆ 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();
3336
3337 if(pdgcode==111 ||pdgcode==221 || pdgcode==331){//for photon
3338 alpha = 1;
3339 }else continue;
3340
3344 double xratio = ang/angmax;
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
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];
3264 init_mode(mymode);
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 ";
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
2595 double L = 2*log(sme);
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
2681 for(int i=1;i<_ndaugs;i++){
2683 }
2684
2686 double mrg = cms - summass;
2688 double mhs = pm*mrg + summass;
2689
2693
2694 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2695
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}
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
2606 double xi2 = 1.64493407;
2607 double xi3=1.2020569;
2609 double L = 2*log(sme);
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.;
2616 wsx = wsx + beta*beta/8. *wsx2;
2618 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
2620}
Referenced by narrowRXS(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), and Ros_xs().
◆ Rad2difXs() [1/3]
| double EvtConExc::Rad2difXs | ( | double | s, |
| double | x ) |
Definition at line 2655 of file EvtConExc.cc.
2655 {// leading second order xsection
2659 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
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
2670 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
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
2628 for(int i=1;i<_ndaugs;i++){
2631 }
2632
2635 double mrg = cms - summass;
2637 double mhs = pm*mrg + summass;
2638
2641 //double mhs = sqrt( s*(1-x) );
2643
2644 // std::cout<<"random : "<<pm<<std::endl;
2645
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(), gamHXSection(), gamHXSection(), gamHXSection(), trapezoid(), 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.
3106 EvtPDL::reSetMass(myres,mjsi);
3107 //EvtPDL::reSetWidth(myres,wjsi);
3108
3110 EvtPDL::reSetMass(myres,mpsip);
3111 //EvtPDL::reSetWidth(myres,wpsip);
3112
3114 EvtPDL::reSetMass(myres,mpsipp);
3115 //EvtPDL::reSetWidth(myres,wpsipp);
3116
3118 EvtPDL::reSetMass(myres,mphi);
3119 //EvtPDL::reSetWidth(myres,wphi);
3120
3122 EvtPDL::reSetMass(myres,momega);
3123 //EvtPDL::reSetWidth(myres,womega);
3124
3126 EvtPDL::reSetMass(myres,mrho0);
3127 //EvtPDL::reSetWidth(myres,wrho0);
3128
3130 EvtPDL::reSetMass(myres,mrho3s);
3131 //EvtPDL::reSetWidth(myres,wrho3s);
3132
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
2746 double nbar = 3.89379304*100000; // ! GeV-2 = 3.89379304*10^5 nbar
2748}
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;
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]) {
3416 }
3417 for(int j=1;j<sxs.size();j++){
3418 double x0 = sxs[j-1];
3419 double x1 = sxs[j];
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];
3043 //if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
3044
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:
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];}
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];}
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}
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition EvtDecayTable.cc:795
Referenced by decay().
◆ SetP4()
| void EvtConExc::SetP4 | ( | EvtParticle * | part, |
| double | mhdr, | ||
| double | xeng, | ||
| double | theta ) |
Definition at line 2955 of file EvtConExc.cc.
2955 { //set the gamma and gamma* momentum according sampled results
2957 double phi=2*3.1415926*pm;
2959 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2963 EvtVector4R p4[2];
2966
2967 for(int i=0;i<2;i++){
2968 EvtId myid = p->getDaug(i)->getId();
2969 p->getDaug(i)->init(myid,p4[i]);
2972 EvtVector4R PG(_cms,0,0,0);
2973 p->getDaug(i)->init(myid,PG);
2974 }
2975 }
2976 }
2977}
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
2981 double phi=2*3.1415926*pm;
2983 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2987 EvtVector4R p4[2];
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 {
2811 double xi2 = 1.64493407;
2812 double xi3=1.2020569;
2814 double L = 2*log(sme);
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.;
2819
2820 double beta2 = beta*beta;
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) +
2827 //return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
2829
2830}
◆ 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
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;
3227 init_mode(mymode);
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{
3821 double sum = 0.0;
3824 {
3825 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, staxsection) + Rad2difXs(s,xL + (i+1)*gaps,staxsection) );
3826 }
3827 return sum;
3828}
◆ 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{
3834 double sum = 0.0;
3837 {
3839 }
3840 return sum;
3841}
◆ VP_sampling()
| bool EvtConExc::VP_sampling | ( | EvtVector4R | pcm, |
| EvtVector4R | pi ) |
Definition at line 2577 of file EvtConExc.cc.
2577 {
2579 double theta=angles.getHelAng(1);
2580 double phi =angles.getHelAng(2);
2582
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 {
3536 ob.open("obsxs.dat");
3537 //
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++){
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
3609 for(int i=0;i<Vmode.size();i++){
3611 }
3614}
Referenced by init().
◆ xs_sampling() [1/2]
| bool EvtConExc::xs_sampling | ( | double | xs | ) |
Definition at line 2537 of file EvtConExc.cc.
2537 {
2539 //std::cout<<"Rdm= "<<pm<<std::endl;
2541}
◆ xs_sampling() [2/2]
| bool EvtConExc::xs_sampling | ( | double | xs, |
| double | xs1 ) |
Definition at line 2543 of file EvtConExc.cc.
2543 {
2545 // std::cout<<"Rdm= "<<pm<<std::endl;
2546 if(pm <xs/(xs0*1.1)) {return true;} else {return false;}
2547}
Member Data Documentation
◆ _cms
|
static |
Definition at line 190 of file EvtConExc.hh.
Referenced by baryon_sampling(), calAF(), checkEvtRatio(), decay(), difgamXs(), Egam2Mhds(), findMaxXS(), getObsXsection(), init(), ISR_ang_integrate(), Mhad_sampling(), mk_VXS(), narrowRXS(), PrintXS(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMass(), SetP4(), SetP4Rvalue(), sumExc(), and writeDecayTabel().
◆ conexcmode
|
static |
Definition at line 212 of file EvtConExc.hh.
Referenced by decay(), and EvtRexc::decay().
◆ getMode
|
static |
Definition at line 194 of file EvtConExc.hh.
Referenced by init().
◆ mlow
|
static |
Definition at line 235 of file EvtConExc.hh.
Referenced by EvtDecay::energySpread(), init(), and mk_VXS().
◆ multimode
|
static |
Definition at line 212 of file EvtConExc.hh.
Referenced by EvtRexc::decay(), and init().
◆ mup
|
static |
Definition at line 235 of file EvtConExc.hh.
Referenced by EvtDecay::energySpread(), and init().
◆ myxsection
|
static |
Definition at line 189 of file EvtConExc.hh.
Referenced by decay(), difgamXs(), difgamXs(), findMaxXS(), init(), init_mode(), mk_VXS(), PrintXS(), Rad1difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMode(), sumExc(), and ~EvtConExc().
◆ SetMthr
|
static |
Definition at line 192 of file EvtConExc.hh.
Referenced by decay(), init(), and EvtDecay::initialize().
◆ staxsection
|
static |
Definition at line 189 of file EvtConExc.hh.
Referenced by calAF(), energySpread(), init(), trapezoid(), and ~EvtConExc().
◆ XS_max
|
static |
Definition at line 191 of file EvtConExc.hh.
Referenced by init(), Rad1difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs(), Rad2difXs2(), and xs_sampling().
The documentation for this class was generated from the following files:
