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

#include <EvtConExc.hh>

+ Inheritance diagram for EvtConExc:

Public Member Functions

 EvtConExc ()
 
virtual ~EvtConExc ()
 
void getName (std::string &name)
 
EvtDecayBaseclone ()
 
void initProbMax ()
 
void init ()
 
void init_Br_ee ()
 
void decay (EvtParticle *p)
 
void init_mode (int mode)
 
double gamHXSection (EvtParticle *p, double El, double Eh, int nmc=100000)
 
double gamHXSection (double s, double El, double Eh, int nmc=100000)
 
double gamHXSection (double El, double Eh)
 
double gamHXSection_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 Rad1difXs (EvtParticle *p)
 
double Ros_xs (double mx, double bree, EvtId pid)
 
double Li2 (double x)
 
double SoftPhoton_xs (double s, double b)
 
double lgr (double *x, double *y, int n, double t)
 
bool islgr (double *x, double *y, int n, double t)
 
double LLr (double *x, double *y, int n, double t)
 
int selectMode (std::vector< int > vmod, double mhds)
 
void findMaxXS (double mhds)
 
bool gam_sampling (double mhds, double sintheta)
 
void resetResMass ()
 
void getResMass ()
 
bool checkdecay (EvtParticle *p)
 
double sumExc (double mx)
 
void showResMass ()
 
int getNdaugs ()
 
EvtIdgetDaugId ()
 
int getSelectedMode ()
 
double narrowRXS (double mxL, double mxH)
 
double selectMass ()
 
double addNarrowRXS (double mhi, double binwidth)
 
void ReadVP ()
 
double getVP (double cms)
 
- Public Member Functions inherited from EvtDecayIncoherent
void makeDecay (EvtParticle *p)
 
virtual ~EvtDecayIncoherent ()
 
void setDaughterSpinDensity (int daughter)
 
int isDaughterSpinDensitySet (int daughter)
 
- Public Member Functions inherited from EvtDecayBase
virtual void getName (std::string &name)=0
 
virtual void decay (EvtParticle *p)=0
 
virtual void makeDecay (EvtParticle *p)=0
 
virtual EvtDecayBaseclone ()=0
 
virtual void init ()
 
virtual void initProbMax ()
 
virtual std::string commandName ()
 
virtual void command (std::string cmd)
 
double getProbMax (double prob)
 
double resetProbMax (double prob)
 
 EvtDecayBase ()
 
virtual ~EvtDecayBase ()
 
virtual bool matchingDecay (const EvtDecayBase &other) const
 
EvtId getParentId ()
 
double getBranchingFraction ()
 
void disableCheckQ ()
 
void checkQ ()
 
int getNDaug ()
 
EvtIdgetDaugs ()
 
EvtId getDaug (int i)
 
int getNArg ()
 
int getPHOTOS ()
 
void setPHOTOS ()
 
void setVerbose ()
 
void setSummary ()
 
double * getArgs ()
 
std::string * getArgsStr ()
 
double getArg (int j)
 
std::string getArgStr (int j)
 
std::string getModelName ()
 
int getDSum ()
 
int summary ()
 
int verbose ()
 
void saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
 
void printSummary ()
 
void setProbMax (double prbmx)
 
void noProbMax ()
 
void checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
 
void checkNDaug (int d1, int d2=-1)
 
void checkSpinParent (EvtSpinType::spintype sp)
 
void checkSpinDaughter (int d1, EvtSpinType::spintype sp)
 
virtual int nRealDaughters ()
 

Static Public Attributes

static EvtXsectionmyxsection
 
static double _cms
 
static double XS_max
 
static int getMode
 
static int conexcmode =-1
 

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 87 of file EvtConExc.hh.

Constructor & Destructor Documentation

◆ EvtConExc()

EvtConExc::EvtConExc ( )
inline

Definition at line 91 of file EvtConExc.hh.

91 {
92 } //constructor

Referenced by clone().

◆ ~EvtConExc()

EvtConExc::~EvtConExc ( )
virtual

Definition at line 165 of file EvtConExc.cc.

165 {
166 if(mydbg){
167 // xs->Write();
168 myfile->Write();
169 }
170 delete myxsection;
171 gamH->deleteTree();
172}
static EvtXsection * myxsection
Definition: EvtConExc.hh:136
void deleteTree()
Definition: EvtParticle.cc:557

Member Function Documentation

◆ addNarrowRXS()

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

Definition at line 2300 of file EvtConExc.cc.

2300 {
2301 double myxs = 0;
2302 for(int i=0;i<ISRXS.size();i++){ // std::cout<<"ISRM "<<ISRM[i]<<std::endl;
2303 if( (ISRM[i]-mhi)<binwidth && ISRFLAG[i]){ISRFLAG[i]=0; myxs = ISRXS[i];}
2304 }
2305 //std::cout<<mhi<<" "<<binwidth<<" "<<myxs<<std::endl;
2306 return myxs;
2307}

Referenced by init().

◆ angularSampling()

bool EvtConExc::angularSampling ( EvtParticle part)

Definition at line 2189 of file EvtConExc.cc.

2189 {
2190 bool tagp,tagk;
2191 tagk=0;
2192 tagp=0;
2193 int nds = par->getNDaug();
2194 for(int i=0;i<par->getNDaug();i++){
2195 EvtId di=par->getDaug(i)->getId();
2196 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
2197 int pdgcode =EvtPDL::getStdHep( di );
2198 double alpha=1;
2199 /*
2200 if(fabs(pdgcode)==2212){//proton and anti-proton
2201 alpha = baryonAng(p4i.mass());
2202 cosp = cos(p4i.theta());
2203 tagp=1;
2204 }else if(fabs(pdgcode)==321||fabs(pdgcode)==211 ){
2205 alpha=1;
2206 cosk = cos(p4i.theta());
2207 tagk=1;
2208 }else continue;
2209 */
2210 double angmax = 1+alpha;
2211 double costheta = cos(p4i.theta());
2212 double ang=1+alpha*costheta*costheta;
2213 double xratio = ang/angmax;
2214 double xi=EvtRandom::Flat(0.,1);
2215 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
2216 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
2217 if(xi>xratio) return false;
2218 }//loop over duaghters
2219 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
2220 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
2221
2222 return true;
2223}
const double alpha
double cos(const BesAngle a)
Definition: EvtId.hh:27
int getId() const
Definition: EvtId.hh:41
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
static double Flat()
Definition: EvtRandom.cc:73
double theta()
Definition: EvtVector4R.cc:249

Referenced by decay().

◆ baryon_sampling()

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

Definition at line 1564 of file EvtConExc.cc.

1564 {
1565 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1566 double theta=angles.getHelAng(1);
1567 double phi =angles.getHelAng(2);
1568 double costheta=cos(theta); //using helicity angles in parent system
1569
1570 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
1571 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
1572 double alpha = baryonAng(_cms);
1573 if(_mode ==96){alpha=-0.34;}
1574 double pm= EvtRandom::Flat(0.,1);
1575 double ang = 1 + alpha*costheta*costheta;
1576 double myang;
1577 if(alpha>=0){myang=1+alpha;}else{myang=1;}
1578 if(pm< ang/myang) {return true;}else{return false;}
1579}
double baryonAng(double mx)
Definition: EvtConExc.cc:2225
static double _cms
Definition: EvtConExc.hh:137

Referenced by hadron_angle_sampling().

◆ baryonAng()

double EvtConExc::baryonAng ( double  mx)

Definition at line 2225 of file EvtConExc.cc.

2225 {
2226 double mp=0.938;
2227 double u = 0.938/mx;
2228 u = u*u;
2229 double u2 = (1+u)*(1+u);
2230 double uu = u*(1+6*u);
2231 double alpha = (u2-uu)/(u2+uu);
2232 return alpha;
2233}
const double mp
Definition: incllambda.cxx:45

Referenced by baryon_sampling().

◆ checkdecay()

bool EvtConExc::checkdecay ( EvtParticle p)

Definition at line 2143 of file EvtConExc.cc.

2143 {
2144 int nson=p->getNDaug();
2145 double massflag=1;
2146 for(int i=0;i<nson;i++){
2147 if( EvtPDL::getStdHep(p->getDaug(i)->getId())==22 ) continue;
2148 massflag *= p->getDaug(i)->mass();
2149 }
2150 if(massflag==0){
2151 std::cout<<"Zero mass!"<<std::endl;
2152 return 0;
2153 }else{return 1;}
2154}
EvtId getId() const
Definition: EvtParticle.cc:113
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
double mass() const
Definition: EvtParticle.cc:127

Referenced by decay().

◆ clone()

EvtDecayBase * EvtConExc::clone ( )
virtual

Implements EvtDecayBase.

Definition at line 180 of file EvtConExc.cc.

180 {
181
182 return new EvtConExc;
183
184}

◆ decay()

void EvtConExc::decay ( EvtParticle p)
virtual

Implements EvtDecayBase.

Definition at line 1030 of file EvtConExc.cc.

1030 {
1031 //std::cout<<"_cms= "<<_cms<<std::endl;
1032 EvtId myvpho = EvtPDL::getId(std::string("vpho"));
1033 if(myvpho != p->getId()){
1034 std::cout<<"Parent needs to be vpho, but found "<<EvtPDL::name(p->getId())<<std::endl;
1035 abort();
1036 }
1037 //for fill root tree
1038 EvtVector4R vgam,hd1,hd2,vhds;
1039
1040 //first for Rscan generator with _mode=74110
1041 if(_mode==74110){
1042 //preparation of mode sampling
1043 std::vector<int> vmod; vmod.clear();
1044 int mn[78]={0,1,2,3,4,5,6,7,8,9,10,11,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, // 12 and 43 is removed
1045 50,51,
1046 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1047 90,91,92,93,94,95,96,
1048 74100,74101,74102,74103,74104,74105,74110};
1049 int mn2[79]={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,//mode 43 are removed
1050 50,51,
1051 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1052 90,91,92,93,94,95,96,
1053 74100,74101,74102,74103,74104,74105,74110};
1054 double mx = p->getP4().mass();
1055 if(_cms>2.5){
1056 for(int i=0;i<78;i++){vmod.push_back(mn[i]);}
1057 }else{
1058 for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
1059 }
1060 int mymode;
1061 double pm= EvtRandom::Flat(0.,1);
1062
1063 //std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
1064 //--
1065 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross section within Egam=(0.025-Egamma^max GeV)
1066 if(pm <_xs0/(_xs0 + _xs1) ){//without ISR photon for mode=74110
1067
1068 mhds=_cms;
1069 mymode=selectMode(vmod,mhds);
1070 _selectedMode = mymode;
1071 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1072 delete myxsection;
1073 myxsection =new EvtXsection (mymode);
1074 init_mode(mymode);
1075 resetResMass();
1076
1077
1078 if(_ndaugs>1){//for e+e- -> exclusive decays
1079 checkA:
1080 p->makeDaughters(_ndaugs,daugs);
1081 double totMass=0;
1082 for(int i=0;i< _ndaugs;i++){
1083 EvtParticle* di=p->getDaug(i);
1084 double xmi=EvtPDL::getMass(di->getId());
1085 di->setMass(xmi);
1086 totMass+=xmi;
1087 }
1088 if(totMass>p->mass()) goto checkA;
1089 p->initializePhaseSpace( _ndaugs , daugs);
1090 if(!checkdecay(p)) goto checkA;
1091 vhds = p->getP4();
1092 if(_cms>2.5 && !angularSampling(p)) goto checkA;
1093 if(_cms<2.5 && !photonSampling(p)) goto checkA;
1094 }else{// for e+e- -> vector resonace, add a very soft photon
1095 EvtId mydaugs[2];
1096 mydaugs[0]=daugs[0];
1097 EvtPDL::reSetMass(mydaugs[0],mhds*0.9999);
1098 //EvtPDL::reSetWidth(mydaugs[0],0);
1099 mydaugs[1]=gamId; //ISR photon
1100 p->makeDaughters(2,mydaugs);
1101 checkB:
1102 double totMass=0;
1103 for(int i=0;i< 2;i++){
1104 EvtParticle* di=p->getDaug(i);
1105 double xmi=EvtPDL::getMass(di->getId());
1106 di->setMass(xmi);
1107 totMass+=xmi;
1108 }
1109 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
1110 if(totMass>p->mass()) goto checkB;
1111 p->initializePhaseSpace(2,mydaugs);
1112
1113 if(!checkdecay(p)) goto checkB;
1114 vhds = p->getDaug(0)->getP4();;
1115 vgam = p->getDaug(1)->getP4();
1116 }
1117 //-----
1118 }else{// with ISR photon for mode=74110
1119 mhds=Mhad_sampling(MH,AF);
1120 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
1121 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1122
1123 if(mydbg) mass2=mhds;
1124
1125 //generating events
1126 ModeSelection:
1128 mymode=selectMode(vmod,mhds);
1129 conexcmode = mymode; //save for model REXC ( see EvtRexc.cc) decay
1130 if(mhds<2.3 && mymode==74110) {goto ModeSelection;} // E<2.3 GeV, fully exclusive production
1131 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1132 _selectedMode = mymode;
1133 delete myxsection;
1134 myxsection =new EvtXsection (mymode);
1135 init_mode(mymode);
1136 EvtId myvpho = EvtPDL::getId(std::string("vgam"));
1137 EvtPDL::reSetMass(myvpho,mhds); //set to continum cms energy
1138 EvtPDL::reSetWidth(myvpho,0);
1139
1140 //debugging
1141 //for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<" "<<EvtPDL::name(daugs[i])<<std::endl;}
1142
1143 //decay e+e- ->gamma_ISR gamma*
1144 EvtId mydaugs[2];
1145 if(_ndaugs>1) { //gamma* -> exclusive decays
1146 resetResMass();
1147 mydaugs[0]=myvpho;
1148 }else{// vector resonance decays
1149 resetResMass();
1150 mydaugs[0]=daugs[0];
1151 EvtPDL::reSetMass(mydaugs[0],mhds);
1152 //EvtPDL::reSetWidth(mydaugs[0],0);
1153 }
1154 mydaugs[1]=gamId; //ISR photon
1155 int maxflag=0;
1156 int trycount=0;
1157 ISRphotonAng_sampling:
1158 double totMass=0;
1159 p->makeDaughters(2,mydaugs);
1160 for(int i=0;i< 2;i++){
1161 EvtParticle* di=p->getDaug(i);
1162 double xmi=EvtPDL::getMass(di->getId());
1163 di->setMass(xmi);
1164 totMass+=xmi;
1165 }
1166 if(totMass>p->mass()) goto ISRphotonAng_sampling;
1167 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
1168 double weight1 = p->initializePhaseSpace(2,mydaugs);
1169 if(!checkdecay(p)) goto ISRphotonAng_sampling;
1170 //set the ISR photon angular distribution
1171 SetP4Rvalue(p,mhds,xeng,theta); //end of decay e+e- -> vpho gamma_ISR
1172
1173 if(maxflag==0) {findMaxXS(mhds); maxflag=1;}
1174 vhds = p->getDaug(0)->getP4();
1175 vgam=p->getDaug(1)->getP4();
1176 double gx=vgam.get(1);
1177 double gy=vgam.get(2);
1178 double sintheta= sqrt(gx*gx+gy*gy)/vgam.d3mag();
1179 bool gam_ang=gam_sampling(mhds,sintheta);
1180 trycount++;
1181
1182 } // with and without ISR sampling block
1183 // finish event generation
1184 // for debugging
1185 if(mydbg){
1186 pgam[0]=vgam.get(0);
1187 pgam[1]=vgam.get(1);
1188 pgam[2]=vgam.get(2);
1189 pgam[3]=vgam.get(3);
1190
1191 phds[0]=vhds.get(0);
1192 phds[1]=vhds.get(1);
1193 phds[2]=vhds.get(2);
1194 phds[3]=vhds.get(3);
1195 costheta = vgam.get(3)/vgam.d3mag();
1196 selectmode = mymode;
1197 xs->Fill();
1198 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
1199 }
1200 _modeFlag[1]= _selectedMode;
1201 p->setIntFlag(_modeFlag);
1202 return;
1203 }//end block if(_mode==74110)
1204
1205 //=======================================================
1206 // e+e- -> gamma_ISR gamma*
1207 //=======================================================
1208 if(VISR){
1209 EvtId gid=EvtPDL::getId("gamma*");
1210 double pm= EvtRandom::Flat(0.,1);
1211
1212 if(pm <_xs0/(_xs0 + _xs1) ){//with a very soft ISR photon
1213 EvtPDL::reSetMass(gid,_cms*0.995);
1214 mhds = _cms*0.995;
1215 }else{
1216 mhds=Mhad_sampling(MH,AF);
1217 EvtPDL::reSetMass(gid,mhds);
1218 }
1219 //debugging
1220 std::cout<<"generatedMass "<<EvtPDL::getMeanMass(gid)<<std::endl;
1221 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
1222 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1223 p->makeDaughters(2,daugs);
1224 for(int i=0;i< 2;i++){
1225 EvtParticle* di=p->getDaug(i);
1226 double xmi=EvtPDL::getMeanMass(di->getId());
1227 di->setMass(xmi);
1228 }
1229 p->initializePhaseSpace(2,daugs);
1230 SetP4(p,mhds,xeng,theta);
1231 return;
1232 }
1233
1234
1235 //========================================================
1236 //=== for other mode generation with _mode != 74110
1237 //========================================================
1238
1239 if((_xs0 + _xs1)==0) {std::cout<<"EvtConExc::zero total xsection"<<std::endl;::abort();}
1240 double pm= EvtRandom::Flat(0.,1);
1241 double xsratio = _xs0/(_xs0 + _xs1);
1242 int iflag=2; //flag = 0 for ee->hadrons, 1 for ee->gamma_ISR hadrons, 2: mix 0 and 1 case
1243 if(getArg(0)!= -2){// exclude DIY case
1244 if(getNArg()==3 && radflag==1){iflag=1;xsratio=0;} // only generating gamma hadrons mode
1245 else if(getNArg()==3 && radflag==0) {iflag=0;xsratio=1;}
1246 }
1247
1248 // std::cout<<"xsratio= "<<xsratio<<std::endl;
1249
1250 if(pm<xsratio ){// no ISR photon
1251
1252 p->makeDaughters(_ndaugs,daugs);
1253 for(int i=0;i< _ndaugs;i++){
1254 EvtParticle* di=p->getDaug(i);
1255 double xmi=EvtPDL::getMass(di->getId());
1256 di->setMass(xmi);
1257 }
1258 angular_hadron:
1259 p->initializePhaseSpace(_ndaugs,daugs);
1260 bool byn_ang;
1261 EvtVector4R pd1 = p->getDaug(0)->getP4();
1262 EvtVector4R pcm(_cms,0,0,0);
1263 if(_ndaugs==2){//angular distribution for the hadrons only for two-body decays
1264 byn_ang=hadron_angle_sampling(pd1, pcm);
1265 if(!byn_ang) goto angular_hadron;
1266 }
1267
1268 //for histogram
1269 cosp = pd1.get(3)/pd1.d3mag();
1270 mhds = _cms;
1271 //std::cout<<"cosp, mhds "<<cosp<<" "<<mhds<<std::endl;
1272 //p->printTree();
1273
1274 }else{// sampling m_hadrons and decay e+e- ->gamma gamma*
1275 double mhdr=Mhad_sampling(MH,AF);
1276 double xeng=1-mhdr*mhdr/(_cms*_cms); //photon energy ratio
1277 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1278 EvtId gid =EvtPDL::getId("vhdr");
1279 EvtPDL::reSetMass(gid,mhdr);
1280 int ndaugs =2;
1281 EvtId mydaugs[2];
1282 mydaugs[0] =EvtPDL::getId(std::string("gamma"));
1283 mydaugs[1]=EvtPDL::getId(std::string("vhdr"));
1284
1285 //for histogram
1286 mhds = mhdr;
1287 costheta = cos(theta);
1288 //--
1289
1290 p->makeDaughters(2,mydaugs);
1291 for(int i=0;i< 2;i++){
1292 EvtParticle* di=p->getDaug(i);
1293 double xmi=EvtPDL::getMass(di->getId());
1294 di->setMass(xmi);
1295 }
1296 p->initializePhaseSpace(2,mydaugs);
1297 SetP4(p,mhdr,xeng,theta); //end of decay e+e- -> gamma_ISR gamma*
1298 //p->printTree();
1299
1300 //----
1301 }//end of gamma gamma*, gamma*->hadrons generation
1302
1303
1304 //-----------------------------------------
1305 // End of decays for all topology
1306 //----------------------------------------
1307 //================================= event analysis
1308 if(_nevt ==0 ){
1309 std::cout<<"The decay chain: "<<std::endl;
1310 p->printTree();
1311 }
1312 _nevt++;
1313 //--- for debuggint to make root file
1314 if(mydbg){
1315 xs->Fill();
1316 }
1317
1318 //----
1319 return ;
1320}
bool checkdecay(EvtParticle *p)
Definition: EvtConExc.cc:2143
void findMaxXS(EvtParticle *p)
Definition: EvtConExc.cc:1466
static int conexcmode
Definition: EvtConExc.hh:157
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:1970
bool photonSampling(EvtParticle *part)
Definition: EvtConExc.cc:2235
int selectMode(std::vector< int > vmod, double mhds)
Definition: EvtConExc.cc:2014
void init_mode(int mode)
Definition: EvtConExc.cc:576
bool angularSampling(EvtParticle *part)
Definition: EvtConExc.cc:2189
double ISR_ang_sampling(double x)
Definition: EvtConExc.cc:1929
double Mhad_sampling(double *x, double *y)
Definition: EvtConExc.cc:1888
void resetResMass()
Definition: EvtConExc.cc:2061
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition: EvtConExc.cc:1322
bool gam_sampling(EvtParticle *p)
Definition: EvtConExc.cc:1538
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:1953
double getArg(int j)
static bool ConExcPythia
Definition: EvtGlobalSet.hh:20
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45
static void reSetWidth(EvtId i, double width)
Definition: EvtPDL.hh:69
static void reSetMass(EvtId i, double mass)
Definition: EvtPDL.hh:68
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
static double getMass(EvtId i)
Definition: EvtPDL.hh:46
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
void setMass(double m)
Definition: EvtParticle.hh:372
void makeDaughters(int ndaug, EvtId *id)
void setIntFlag(std::vector< int > vi)
Definition: EvtParticle.hh:151
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
void printTree() const
Definition: EvtParticle.cc:897
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double mass() const
Definition: EvtVector4R.cc:39
double get(int i) const
Definition: EvtVector4R.hh:179
double d3mag() const
Definition: EvtVector4R.cc:186

◆ difgamXs() [1/2]

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

Definition at line 2003 of file EvtConExc.cc.

2003 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2004 double x = 1-(mhds/_cms)*(mhds/_cms);
2005 double sin2theta = sintheta*sintheta;
2006 double alpha = 1./137.;
2007 double pi = 3.1415926;
2008 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2009 double xs = myxsection->getXsection(mhds);
2010 double difxs = 2*mhds/_cms/_cms * wsx *xs;
2011 return difxs;
2012}
Double_t x[10]
double getXsection(double mx)

◆ difgamXs() [2/2]

double EvtConExc::difgamXs ( EvtParticle p)

Definition at line 1511 of file EvtConExc.cc.

1511 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
1512 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1513 EvtVector4R p0 = p->getDaug(0)->getP4();
1514 for(int i=1;i<_ndaugs;i++){
1515 p0 += p->getDaug(i)->getP4();
1516 }
1517 EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
1518 double mhs = p0.mass();
1519 double egam = pgam.get(0);
1520 double sin2theta = pgam.get(3)/pgam.d3mag();
1521 sin2theta = 1-sin2theta*sin2theta;
1522
1523 double cms = p->getP4().mass();
1524 double alpha = 1./137.;
1525 double pi = 3.1415926;
1526 double x = 2*egam/cms;
1527 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
1528 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
1529
1530 double xs = myxsection->getXsection(mhs);
1531 double difxs = 2*mhs/cms/cms * wsx *xs;
1532 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1533 return difxs;
1534
1535}
double getErr(double mx)

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

◆ findMaxXS() [1/2]

void EvtConExc::findMaxXS ( double  mhds)

Definition at line 1988 of file EvtConExc.cc.

1988 {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
1989 maxXS=-1;
1990 for(int i=0;i<90000;i++){
1991 double x = 1-(mhds/_cms)*(mhds/_cms);
1992 double cos = EvtRandom::Flat(0.0006,0.9994); //set cut on ISR gamma 2^o
1993 double sin2theta =sqrt(1-cos*cos);
1994 double alpha = 1./137.;
1995 double pi = 3.1415926;
1996 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
1997 double xs = myxsection->getXsection(mhds);
1998 double difxs = 2*mhds/_cms/_cms * wsx *xs;
1999 if(difxs>maxXS)maxXS=difxs;
2000 }
2001}

◆ findMaxXS() [2/2]

void EvtConExc::findMaxXS ( EvtParticle p)

Definition at line 1466 of file EvtConExc.cc.

1466 {
1467 //std::cout<<"nmc= "<<nmc<<std::endl;
1468 gamId = EvtPDL::getId(std::string("gamma"));
1469 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
1470 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
1471 double totxs = 0;
1472 maxXS=-1;
1473 _er1=0;
1474 Rad2Xs =0;
1475 int nmc = 50000;
1476 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
1477 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
1478 int gamdaugs = _ndaugs+1;
1479
1480 EvtId theDaugs[10];
1481 for(int i=0;i<_ndaugs;i++){
1482 theDaugs[i] = daugs[i];
1483 }
1484 theDaugs[_ndaugs]=gamId;
1485
1486 gamH->makeDaughters(gamdaugs,theDaugs);
1487 loop:
1488 double totm=0;
1489 for(int i=0;i<gamdaugs;i++){
1490 EvtParticle* di=gamH->getDaug(i);
1491 double xmi=EvtPDL::getMass(di->getId());
1492 di->setMass(xmi);
1493 totm += xmi;
1494 }
1495
1496 //std::cout<<totm<<" "<<p_init.get(0)<<std::endl;
1497 if(totm >= p_init.get(0) ) goto loop;
1498
1499 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
1500 double thexs = difgamXs(gamH); //prepare the photon angular distribution
1501 EvtVector4R p4gam = gamH->getDaug(_ndaugs)->getP4();
1502 double costheta = p4gam.get(3)/p4gam.d3mag();
1503 double sintheta = sqrt(1-costheta*costheta);
1504 bool acut=(sintheta>0.04) && (sintheta<0.96); //set photon auglar cut 2^o
1505 if(thexs>maxXS && acut ) {maxXS=thexs;}
1506 //gamH->deleteTree();
1507 }
1508
1509}
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition: KarFin.h:34
double difgamXs(EvtParticle *p)
Definition: EvtConExc.cc:1511
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)

Referenced by decay(), and init().

◆ gam_sampling() [1/2]

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

Definition at line 1545 of file EvtConExc.cc.

1545 {
1546 double pm= EvtRandom::Flat(0.,1);
1547 double xs = difgamXs( mhds,sintheta );
1548 double xsratio = xs/maxXS;
1549 if(pm<xsratio){return true;}else{return false;}
1550}

◆ gam_sampling() [2/2]

bool EvtConExc::gam_sampling ( EvtParticle p)

Definition at line 1538 of file EvtConExc.cc.

1538 {//photon angular distribution sampling
1539 double pm= EvtRandom::Flat(0.,1);
1540 double xs = difgamXs( p );
1541 double xsratio = xs/maxXS;
1542 if(pm<xsratio){return true;}else{return false;}
1543}

Referenced by decay().

◆ gamHXSection() [1/3]

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

Definition at line 1431 of file EvtConExc.cc.

1431 {// using Gaussian integration subroutine from Cern lib
1432 std::cout<<" "<<std::endl;
1433 extern double Rad2difXs(double *x);
1434 extern double Rad2difXs2(double *x);
1435 double eps = 0.01;
1436 double Rad2Xs;
1437 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
1438 if(Rad2Xs==0){
1439 for(int iter=0;iter<10;iter++){
1440 eps += 0.01;
1441 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
1442 if(!Rad2Xs) return Rad2Xs;
1443 }
1444 }
1445 return Rad2Xs; // the leading second order correction
1446}
double Rad2difXs2(double *mhs)
Definition: EvtConExc.cc:1782
double dgauss_(double(*fun)(double *), double *, double *, double *)
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
EvtTensor3C eps(const EvtVector3R &v)
Definition: EvtTensor3C.cc:307
double Rad2difXs(EvtParticle *p)
Definition: EvtConExc.cc:1639

◆ gamHXSection() [2/3]

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

Definition at line 1399 of file EvtConExc.cc.

1399 {//El, Eh : the lower and higher energy of radiative photons
1400 //--for narrow resonance psi(2S),J/psi, phi, which can not calculated with the MC integration
1401 //double mxL = sqrt( s-2*Eh*sqrt(s) ); //low mass
1402 //double mxH = sqrt( s-2*El*sqrt(s) ); //high mass
1403 //double sigv = narrowRXS(mxL,mxH);
1404 //--
1405
1406 double totxs = 0;
1407 maxXS=-1;
1408 _er1=0;
1409 Rad2Xs =0;
1410 double xL=2*El/sqrt(s);
1411 double xH=2*Eh/sqrt(s);
1412 for(int i=0;i<nmc;i++){//It is found that the narrow resonance can not included in this MC integartion
1413 double rdm = EvtRandom::Flat(0.,1.);// set angular cut 1^o
1414 double x= xL+ (xH-xL)*rdm;
1415 Rad2Xs += Rad2difXs(s,x);
1416 _er1 += differ2; //stored when calculate Rad2Xs
1417 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
1418 }
1419 _er1/=nmc;
1420 _er1*=(xH-xL);
1421 //std::cout<<"_er1= "<<_er1<<std::endl;
1422 Rad2Xs/=nmc; // the leading second order correction
1423 double thexs= Rad2Xs*(xH-xL); //xH-xL is the Jacobi factor
1424 //std::cout<<"thexs= "<<thexs<<std::endl;
1425 //if(sigv != -1) thexs += sigv; //add narrow resonance ISR cross section
1426 return thexs;
1427}
XmlRpcServer s
Definition: HelloServer.cpp:11

◆ gamHXSection() [3/3]

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

Definition at line 1340 of file EvtConExc.cc.

1340 {//El, Eh : the lower and higher energy of radiative photons
1341 //std::cout<<"nmc= "<<nmc<<std::endl;
1342 gamId = EvtPDL::getId(std::string("gamma"));
1343 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
1344 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
1345 double totxs = 0;
1346 maxXS=-1;
1347 _er1=0;
1348 Rad2Xs =0;
1349 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
1350 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
1351 int gamdaugs = _ndaugs+1;
1352
1353 EvtId theDaugs[10];
1354 for(int i=0;i<_ndaugs;i++){
1355 theDaugs[i] = daugs[i];
1356 }
1357 theDaugs[_ndaugs]=gamId;
1358
1359 gamH->makeDaughters(gamdaugs,theDaugs);
1360
1361 for(int i=0;i<gamdaugs;i++){
1362 EvtParticle* di=gamH->getDaug(i);
1363 double xmi=EvtPDL::getMass(di->getId());
1364 di->setMass(xmi);
1365 }
1366 loop:
1367 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
1368 //-- slection:theta_gamma > 1 degree
1369 EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
1370 double pmag = pgam.d3mag();
1371 double pz = pgam.get(3);
1372 double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
1373 double onedegree = 1./180.*3.1415926;
1374 //if(sintheta < onedegree) {goto loop;}
1375 if(pmag <El || pmag >Eh) {goto loop;}
1376 //--------
1377 // std::cout<<"pmag= "<<pmag<<std::endl;
1378
1379 double thexs = difgamXs(gamH); //prepare the photon angular distribution
1380 Rad2Xs += Rad2difXs( gamH );
1381 if(thexs>maxXS) {maxXS=thexs;}
1382 double s = p_init.mass2();
1383 double x = 2*pgam.get(0)/sqrt(s);
1384 double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
1385 totxs += rad1xs;
1386 _er1 += differ;
1387 gamH->deleteDaughters();
1388 } //for int_i block
1389 _er1/=nmc;
1390
1391 Rad2Xs/=nmc; // the leading second order correction
1392 totxs/=nmc; // first order correction XS
1393
1394 // return totxs; // first order correction XS
1395 return Rad2Xs; // the leading second order correction
1396}
double Rad1difXs(EvtParticle *p)
Definition: EvtConExc.cc:1682
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
void deleteDaughters(bool keepChannel=false)
Definition: EvtParticle.cc:540

Referenced by init().

◆ gamHXSection_er()

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

Definition at line 1448 of file EvtConExc.cc.

1448 {// using Gaussian integration subroutine from Cern lib
1449 std::cout<<" "<<std::endl;
1450 extern double Rad2difXs_er(double *x);
1451 extern double Rad2difXs_er2(double *x);
1452 double eps = 0.01;
1453 double Rad2Xs;
1454 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
1455 if(Rad2Xs==0){
1456 for(int iter=0;iter<10;iter++){
1457 eps += 0.01;
1458 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
1459 if(!Rad2Xs) return Rad2Xs;;
1460 }
1461 }
1462 return Rad2Xs; // the leading second order correction
1463}
double Rad2difXs_er2(double *mhs)
Definition: EvtConExc.cc:1797
double Rad2difXs_er(double *m)
Definition: EvtConExc.cc:1769

◆ getDaugId()

EvtId * EvtConExc::getDaugId ( )
inline

Definition at line 155 of file EvtConExc.hh.

155{return daugs;}

Referenced by EvtRexc::decay().

◆ getName()

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

Implements EvtDecayBase.

Definition at line 174 of file EvtConExc.cc.

174 {
175
176 model_name="ConExc";
177
178}

◆ getNdaugs()

int EvtConExc::getNdaugs ( )
inline

Definition at line 154 of file EvtConExc.hh.

154{return _ndaugs;}

Referenced by EvtRexc::decay().

◆ getResMass()

void EvtConExc::getResMass ( )

Definition at line 2096 of file EvtConExc.cc.

2096 {
2097 EvtId myres = EvtPDL::getId(std::string("J/psi"));
2098 mjsi = EvtPDL::getMeanMass(myres);
2099 wjsi = EvtPDL::getWidth(myres);
2100
2101 myres = EvtPDL::getId(std::string("psi(2S)"));
2102 mpsip= EvtPDL::getMeanMass(myres);
2103 wpsip= EvtPDL::getWidth(myres);
2104
2105 myres = EvtPDL::getId(std::string("psi(3770)"));
2106 mpsipp= EvtPDL::getMeanMass(myres);
2107 wpsipp= EvtPDL::getWidth(myres);
2108
2109 myres = EvtPDL::getId(std::string("phi"));
2110 mphi = EvtPDL::getMeanMass(myres);
2111 wphi = EvtPDL::getWidth(myres);
2112
2113 myres = EvtPDL::getId(std::string("omega"));
2114 momega= EvtPDL::getMeanMass(myres);
2115 womega= EvtPDL::getWidth(myres);
2116
2117 myres = EvtPDL::getId(std::string("rho0"));
2118 mrho0 = EvtPDL::getMeanMass(myres);
2119 wrho0 = EvtPDL::getWidth(myres);
2120
2121 myres = EvtPDL::getId(std::string("rho(3S)0"));
2122 mrho3s= EvtPDL::getMeanMass(myres);
2123 wrho3s= EvtPDL::getWidth(myres);
2124
2125
2126 myres = EvtPDL::getId(std::string("omega(2S)"));
2127 momega2s= EvtPDL::getMeanMass(myres);
2128 womega2s= EvtPDL::getWidth(myres);
2129
2130}
static double getWidth(EvtId i)
Definition: EvtPDL.hh:54

Referenced by init().

◆ getSelectedMode()

int EvtConExc::getSelectedMode ( )
inline

Definition at line 156 of file EvtConExc.hh.

156{return _selectedMode;}

◆ getVP()

double EvtConExc::getVP ( double  cms)

Definition at line 2334 of file EvtConExc.cc.

2334 {
2335 double xx,r1,i1;
2336 double x1,y1,y2;
2337 xx=0;
2338 int mytag=1;
2339 for(int i=0;i<4001;i++){
2340 x1=vpx[i];
2341 y1=vpr[i];
2342 y2=vpi[i];
2343 if(xx<=mx && mx <x1){xx=x1; r1=y1; i1=y2; mytag=0; break;}
2344 xx=x1; r1=y1; i1=y2;
2345 }
2346 if(mytag==1){std::cout<<"No vacuum polarization value found"<<std::endl;abort();}
2347 EvtComplex cvp(r1,i1);
2348 double thevp=abs(1./(1-cvp));
2349 if(3.0933<mx && mx<3.1035) return 1.; //J/psi peak excluded
2350 if(3.6810<mx && mx<3.6913) return 1.; //psi(2S) peak excluded
2351 return thevp*thevp;
2352}

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

◆ hadron_angle_sampling()

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

Definition at line 1322 of file EvtConExc.cc.

1322 {
1323 EvtVector4R pbst=-1*pcm;
1324 pbst.set(0,pcm.get(0));
1325 EvtVector4R p4i = boostTo(ppi,pbst);
1326 if(_mode<=5 || _mode ==44|| _mode ==47|| _mode ==48 ||_mode ==72 ||_mode == 73||_mode==94||_mode==95 ||_mode ==96 || _mode>=23 &&_mode<=27 || _mode==36|| _mode ==39 || _mode == 80 ){//ee->two baryon mode, VP, SP
1327 bool byn_ang = baryon_sampling(pcm, p4i);
1328 return byn_ang;
1329 }else if(_mode ==6||_mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
1330 bool msn_ang = meson_sampling(pcm,p4i);
1331 return msn_ang;
1332 }else if(_mode==23 || _mode==24 || _mode==25 || _mode==26 || _mode==27 || _mode==36||_mode==47||_mode==48||_mode==68||_mode==69||_mode==72||_mode==73||_mode==80||_mode==94||_mode==95){
1333 bool msn_ang = VP_sampling(pcm,p4i);
1334 return msn_ang;
1335 }
1336 return true;
1337}
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1592
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1581
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1564
void set(int i, double d)
Definition: EvtVector4R.hh:183

Referenced by decay().

◆ init()

void EvtConExc::init ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 187 of file EvtConExc.cc.

187 {
188 ReadVP();
189 //----------------
190 // check that there are 0 arguments
191 // checkNArg(1);
192 //Vector ISR process
193 VISR=0;
194 if(getNDaug()==2){
195 if(getDaugs()[0]==EvtPDL::getId("gamma") &&getDaugs()[1]==EvtPDL::getId("gamma*") || getDaugs()[0]==EvtPDL::getId("gamma*") &&getDaugs()[1]==EvtPDL::getId("gamma") ) VISR=1;
196 }else if(getNDaug()>2){std::cout<<"Bad daughter specified"<<std::endl;abort();}
197
198 cmsspread=0.0015; //CMC energy spread
199 testflag=0;
200 getResMass();
201 if(getArg(0) == -1){
202 radflag=0;mydbg=false;
203 _mode = getArg(0);
204 pdgcode = getArg(1);
205 pid = EvtPDL::evtIdFromStdHep(pdgcode );
206 nson = getNArg()-2;
207 std::cout<<"The decay daughters are "<<std::endl;
208 for(int i=2;i<getNArg();i++ ){son[i-2]=EvtPDL::evtIdFromStdHep(getArg(i));std::cout<<EvtPDL::name(son[i-2])<<" ";}
209 std::cout<<std::endl;
210 }else if(getArg(0)==-2){
211 radflag=0;mydbg=false;
212 _mode = getArg(0);
213 nson = getNArg()-1;
214 for(int i=1;i<getNArg();i++ ){son[i-1]=EvtPDL::evtIdFromStdHep(getArg(i));}
215 }
216 else if(getNArg()==1){ _mode = getArg(0);radflag=0;mydbg=false;}
217 else if(getNArg()==2){ _mode = getArg(0); radflag=getArg(1);mydbg=false;}
218 else if(getNArg()==3){ _mode = getArg(0); radflag=getArg(1);mydbg=getArg(2);}
219 else{std::cout<<"ConExc:umber of parameter should be 1,2 or 3, but you set "<<getNArg()<<std::endl;::abort(); }
220 //--- debugging
221 //std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<" "<<EvtPDL::name(son[1])<<std::endl;
222
223 gamId = EvtPDL::getId(std::string("gamma"));
224 init_mode(_mode);
225 XS_max=-1;
226 //-----debugging to make root file
227 if(mydbg){
228 myfile = new TFile("mydbg.root","RECREATE");
229 xs=new TTree ("xs","momentum of rad. gamma and hadrons");
230 xs->Branch("imode",&imode,"imode/I");
231 xs->Branch("pgam",pgam,"pgam[4]/D");
232 xs->Branch("phds",phds,"phds[4]/D");
233 xs->Branch("ph1",ph1,"ph1[4]/D");
234 xs->Branch("ph2",ph2,"ph2[4]/D");
235 xs->Branch("mhds",&mhds,"mhds/D");
236 xs->Branch("mass1",&mass1,"mass1/D");
237 xs->Branch("mass2",&mass2,"mass2/D");
238 xs->Branch("costheta",&costheta,"costheta/D");
239 xs->Branch("cosp",&cosp,"cosp/D");
240 xs->Branch("cosk",&cosk,"cosk/D");
241 xs->Branch("sumxs",&sumxs,"sumxs/D");
242 xs->Branch("selectmode",&selectmode,"selectmode/D");
243 }
244 //--- prepare the output information
245 EvtId parentId =EvtPDL::getId(std::string("vpho"));
246 EvtPDL::reSetWidth(parentId,0.0);
247 double parentMass = EvtPDL::getMass(parentId);
248 std::cout<<"parent mass = "<<parentMass<<std::endl;
249
250
251 EvtVector4R p_init(parentMass,0.0,0.0,0.0);
253 myxsection =new EvtXsection (_mode);
254 double mth=0;
255 double mx = EvtPDL::getMeanMass(parentId); //p->getP4().mass();
256 if(_mode ==-1){
257 myxsection->setBW(pdgcode);
258 for(int i=0;i<nson;i++){mth +=EvtPDL::getMeanMass(son[i]);}
259 if(mx<mth){std::cout<<"The CMS energy is lower than the threshold of the final states"<<std::endl;abort();}
260 }else if(_mode == -2){
261 mth=myxsection->getXlw();
262 }else{ mth=myxsection->getXlw();}
263 _cms = mx;
264 _unit = myxsection->getUnit();
265
266 std::cout<<"The specified mode is "<<_mode<<std::endl;
267 EvtConExc::getMode = _mode;
268 //std::cout<<"getXlw= "<<mth<<std::endl;
269
270 Egamcut = 0.001; //set photon energy cut according to the BESIII detector
271 double Esig = 0.0001; //sigularity engy
272 double x = 2*Egamcut/_cms;
273 double s = _cms*_cms;
274 double mhscut = sqrt(s*(1-x));
275
276 //for vaccum polarization
277 float fe,fst2,fder,ferrder,fdeg,ferrdeg;
278 fe=_cms;
279 //using the vacuum pol. value as given by http://www-com.physik.hu-berlin.de/~fjeger/software.html
280 vph=getVP(_cms);
281 if(3.0943<_cms && _cms<3.102) vph=1;// For J/psi, the vacuum factor is included in the resonance
282 std::cout<<"sum of leptonic, hadronic contributions(u,d,s,c,b) to vacuum polarization: "<<vph<<" @"<<fe<<"GeV"<<std::endl;
283
284 //caculate the xs for ISR to psiprime, J/psi and phi narrow resonance.
285 double totxsIRS=0;
286 init_Br_ee();
287 double mthrld = EvtPDL::getMeanMass(daugs[0]); //threshold mass of hadrons
288 for(int jj=1;jj<_ndaugs;jj++){
289 mthrld += EvtPDL::getMeanMass(daugs[jj]);
290 }
291
292 ISRXS.clear();ISRM.clear();ISRFLAG.clear();
293 for(int ii=0;ii<3;ii++){
294 double mR = EvtPDL::getMeanMass(ResId[ii]);
295 if(mx<mR || _mode !=74110) continue;
296 double narRxs=Ros_xs(mx,BR_ee[ii],ResId[ii]);
297 std::cout<<"The corss section for gamma "<<EvtPDL::name(ResId[ii]).c_str()<<" is: "<< narRxs<<"*Br (nb)."<<std::endl;
298 ISRXS.push_back(narRxs);
299 ISRM.push_back(mR);
300 ISRFLAG.push_back(1.);
301 ISRID.push_back(ResId[ii]);
302 totxsIRS += narRxs;
303 }
304 std::cout<<"==========================================================================="<<std::endl;
305
306 //-- calculated with MC integration method
307 double mhdL=myxsection->getXlw();
308 EgamH = (s-mhdL*mhdL)/2/sqrt(s);
309 int nmc=1000000;
310 _xs0 = gamHXSection(s,Esig,Egamcut,nmc);
311 _er0 = _er1; // stored when calculate _xs0
312 std::cout<<"_er0= "<<_er0<<std::endl;
313 _xs1 = gamHXSection(s,Egamcut,EgamH,nmc);
314 double xs1_err = _er1;
315
316
317 //--- sigularity point x from 0 to 0.0001GeV
318 double xb= 2*Esig/_cms;
319 double sig_xs = SoftPhoton_xs(s,xb)*(myxsection->getXsection(mx));
320 _xs0 += sig_xs;
321
322 //prepare the observed cross section with interpotation function divdif in CERNLIB
323 testflag=1;
324 int Nstp=598;
325 double stp=(EgamH - Egamcut)/Nstp;
326 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
327 double Eg0=EgamH - i*stp;
328 double Eg1=EgamH - (i+1)*stp;
329 int nmc=20000;
330 double xsi= gamHXSection(s,Eg1,Eg0,nmc);
331 AA[i]=(Eg1+Eg0)/2;
332 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
333 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
334 double binwidth = mh2-mhi;
335 if(_mode==74110) xsi += addNarrowRXS(mhi,binwidth);
336 if(i==0) AF[0]=xsi;
337 if(i>=1) AF[i]=AF[i-1]+xsi;
338 RadXS[i]=xsi/stp;
339 }
340 //add the interval 0~Esig
341 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
342 AF[598]=AF[597];
343 AF[599]=AF[598]+ _xs0;
344 RadXS[599]=_xs0;
345 std::cout<<"mhadL= "<<mhdL<<" ecms "<<_cms<<std::endl;
346 for(int i=0;i<600;i++){
347 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
348 //std::cout<<i<<" Egamma "<<AA[i]<<" Mhadons "<<MH[i]<<std::endl;
349 }
350 //for debugging
351 //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;}
352 std::cout<<"EgamH="<<EgamH<<" "<<_xs0+_xs1<<" AF[599]="<<AF[599]<<" AF[598] "<<AF[598]<<std::endl;
353 //double Etst=0.0241838;
354 //std::cout<<"Etst="<<Etst<<" "<<gamHXSection(s,Etst,EgamH,nmc)<<" "<< lgr(AA,AF,600,Etst)<<std::endl; abort();
355
356 //for information output
357 double bornXS = myxsection->getXsection(mx); // std::cout<<"BornXS= "<<bornXS<<std::endl;
358 double bornXS_er=myxsection->getErr(mx);
359 double obsvXS = _xs0 + _xs1; //gamHXSection(mth ,mx);
360 double obsvXS_er= _er0 + xs1_err;
361 double corr = obsvXS/bornXS;
362 double corr_er =corr*sqrt(bornXS_er*bornXS_er/bornXS/bornXS + obsvXS_er*obsvXS_er/obsvXS/obsvXS);
363
364
365 if(bornXS==0){std::cout<<"EvtConExc::init : The Born cross section at this point is zero!"<<std::endl;abort();}
366 if(obsvXS==0){std::cout<<"EvtConExc::init : The observed cross section at this point is zero!"<<std::endl;abort();}
367
368 //_xs0 += bornXS;// events with very soft photon (Egam<0.025) are taken as the born process
369 //_er0 = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
370
371 std::cout<<""<<std::endl;
372 std::cout<<"========= Generation using cross section (Egamcut="<<Egamcut<<" GeV) =============="<<std::endl;
373 std::cout<<" sqrt(s)= "<<mx<< " GeV "<<std::endl;
374 if(_mode>=0 || _mode ==-2 ){
375 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<"+/-"<<_er0<<" in Unit "<<_unit.c_str()<<std::endl;
376 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"+/-"<<xs1_err<<" in Unit "<<_unit.c_str()<<std::endl;
377 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
378 std::cout<<"The radiative correction factor f_ISR*f_vacuum= sigma_obs/sigma_born(s) = "<<corr<<"+/-"<< fabs(corr_er)<<","<<std::endl;
379 std::cout<<"which is calcualted with these quantities:"<<std::endl;
380 std::cout<<"the observed cross section is "<<obsvXS<<"+/-"<< obsvXS_er<<_unit.c_str()<<std::endl;
381 std::cout<<"and the Born cross section (s) is "<<bornXS<<"+/-"<<bornXS_er<<_unit.c_str()<<std::endl;
382 std::cout<<"and the vaccum polarziation factor (lep. + hadr.) 1/|1-Pi|^2= "<<vph<<std::endl;
383 std::cout<<"Within the range "<<myxsection->getXlw()<<" ~"<<myxsection->getXup()<<" GeV, "<<myxsection->getMsg().c_str()<<std::endl;
384 std::cout<<"==========================================================================="<<std::endl;
385 }else if(_mode==-1){
386 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<" *Br_ee"<<" in Unit "<<_unit.c_str()<<std::endl;
387 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"*Br_ee in Unit "<<_unit.c_str()<<std::endl;
388 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
389 std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<std::endl;
390 std::cout<<"==========================================================================="<<std::endl;
391 }
392 std::cout<<std::endl;
393 std::cout<<std::endl;
394
395 findMaxXS(p);
396 _mhdL=myxsection->getXlw();
397 //--debugging
398 //std::cout<<"maxXS= "<<maxXS<<std::endl;
399
400 if(_xs0 == 0 && _xs1==0){std::cout<<"EvtConExc::ini() has zero cross section"<<std::endl;abort();}
401
402 std::cout<<std::endl;
403 std::cout<<std::endl;
404
405 //--- for debugging
406 if(mydbg && _mode!=74110){
407 int nd = myxsection->getYY().size();
408 double xx[10000],yy[10000],xer[10000],yer[10000];
409 for(int i=0;i<nd;i++){
410 xx[i]=myxsection->getXX()[i];
411 yy[i]=myxsection->getYY()[i];
412 yer[i]=myxsection->getEr()[i];
413 xer[i]=0.;
414 //std::cout<<"yy "<<yy[i]<<std::endl;
415 }
416 myth=new TH1F("myth","Exp.data",200,xx[0],xx[nd]);
417 for(int i=0;i<nd;i++){
418 myth->Fill(xx[i],yy[i]);
419 }
420 myth->SetError(yer);
421 myth->Write();
422
423 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
424 }else
425
426 if(mydbg && _mode==74110 ){
427 int nd = myxsection->getYY().size();
428 double xx[10000],yy[10000],xer[10000],yer[10000],ysum[10000],yysum[10000];
429 for(int i=0;i<nd;i++){
430 xx[i]=myxsection->getXX()[i];
431 yy[i]=myxsection->getYY()[i];
432 yer[i]=myxsection->getEr()[i];
433 xer[i]=0.;
434 //std::cout<<"yy "<<yy[i]<<std::endl;
435 }
436#include "sty.C"
437 double xhigh=5.0;
438 double xlow = myxsection->getXlw();
439 TGraphErrors *Gdt = new TGraphErrors(nd,xx,yy,xer,yer);
440
441 myth=new TH1F("myth","Exp.data",600,xlow,xhigh);
442 Xsum=new TH1F("sumxs","sum of exclusive xs",600,xlow,xhigh);
443 double mstp=(xhigh-xlow)/600;
444 for(int i=0;i<600;i++){
445 double mx=i*mstp+xlow;
446 double xsi = myxsection->getXsection(mx);
447 if(xsi==0) continue;
448 myth->Fill(mx,xsi);
449 //std::cout<<mx<<" "<<xsi<<std::endl;
450 }
451
452 for(int i=0;i<600;i++){
453 double mx=i*mstp+xlow;
454 ysum[i]=sumExc(mx);
455 if(ysum[i]==0) continue;
456 Xsum->Fill(mx,ysum[i]);
457 //std::cout<<mx<<" "<<ysum[i]<<std::endl;
458 }
459
460 for(int i=0;i<nd;i++){
461 yysum[i]=sumExc(xx[i]);
462 }
463
464 myth->SetError(yer);
465 myth->Write();
466 Xsum->Write();
467
468 TGraphErrors *Gsum = new TGraphErrors(nd,xx,yysum,xer,yer);
469 //draw graph
470 double ex[600]={600*0};
471 double ey[600],ta[600];
472 double exz[600]={600*0};
473 double eyz[600]={600*0};
474 for(int i=0;i<600;i++){
475 ey[i]=AF[i]*0.001;
476 }
477 TGraphErrors *gr0 = new TGraphErrors(600,MH,AF,ex,ey);
478 TGraphErrors *gr1 = new TGraphErrors(600,MH,RadXS,exz,eyz);
479 TPostScript *ps = new TPostScript("xsobs.ps",111);
480 TCanvas *c1 = new TCanvas("c1","TGraphs for data",200,10,600,400);
481 gPad-> SetLogy(1);
482 // c1->SetLogy(1);
483 gStyle->SetCanvasColor(0);
484 gStyle->SetStatBorderSize(1);
485 c1->Divide(1,1);
486
487 c1->Update();
488 ps->NewPage();
489 c1->Draw();
490 c1->cd(1);
491 c1->SetLogy(1);
492 gr0->SetMarkerStyle(10);
493 gr0->Draw("AP");
494 gr0->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
495 gr0->GetYaxis()->SetTitle("observed cross section (nb)");
496 gr0->Write();
497
498 c1->Update();
499 ps->NewPage();
500 c1->Draw();
501 c1->cd(1);
502 c1->SetLogy(1);
503 gr1->SetMarkerStyle(10);
504 gr1->Draw("AP");
505 gr1->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
506 gr1->GetYaxis()->SetTitle("Rad*BornXS");
507 gr1->Write();
508
509 //--------- imposing graphs to a pad
510 TMultiGraph *mg = new TMultiGraph();
511 mg->SetTitle("Exclusion graphs");
512 Gdt->SetMarkerStyle(2);
513 Gdt->SetMarkerSize(1);
514 Gsum->SetLineColor(2);
515 Gsum->SetLineWidth(1);
516 mg->Add(Gdt);
517 mg->Add(Gsum);
518
519 c1->Update();
520 ps->NewPage();
521 c1->Draw();
522 c1->cd(1);
523 gPad-> SetLogy(1);
524 mg->Draw("APL");
525 mg->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
526 mg->GetYaxis()->SetTitle("observed cross section (nb)");
527 mg->Write();
528 //-------
529
530 c1->Update();
531 ps->NewPage();
532 c1->Draw();
533 c1->cd(1);
534 Gdt->SetMarkerStyle(2);
535 Gdt->Draw("AP");
536 Gdt->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
537 Gdt->GetYaxis()->SetTitle("observed cross section (nb)");
538 Gdt->Write();
539
540 c1->Update();
541 ps->NewPage();
542 c1->Draw();
543 c1->cd(1);
544 Gsum->SetMarkerStyle(2);
545 Gsum->SetMarkerSize(1);
546 Gsum->Draw("AP");
547 Gsum->SetLineWidth(0);
548 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
549 Gsum->GetYaxis()->SetTitle("observed cross section (nb)");
550 Gsum->Write();
551
552 c1->Update();
553 ps->NewPage();
554 c1->Draw();
555 c1->cd(1);
556 // gPad-> SetLogx(1);
557 Gdt->SetMarkerStyle(2);
558 Gdt->SetMarkerSize(1);
559 Gdt->SetMaximum(8000.0);
560 Gsum->SetLineColor(2);
561 Gsum->SetLineWidth(1.5);
562 Gsum->Draw("AL"); //A: draw axis
563 Gdt->Draw("P"); // superpose to the Gsum, without A-option
564 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
565 Gsum->GetYaxis()->SetTitle("cross section (nb)");
566 Gsum->Write();
567
568 ps->Close();
569 } //if(mydbg)_block
570 //-------------------------
571
572}
int mstp[200]
Definition: EvtPycont.cc:54
static int getMode
Definition: EvtConExc.hh:140
double addNarrowRXS(double mhi, double binwidth)
Definition: EvtConExc.cc:2300
void init_Br_ee()
Definition: EvtConExc.cc:1710
void ReadVP()
Definition: EvtConExc.cc:2355
double Ros_xs(double mx, double bree, EvtId pid)
Definition: EvtConExc.cc:1742
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition: EvtConExc.cc:1340
double sumExc(double mx)
Definition: EvtConExc.cc:2157
double SoftPhoton_xs(double s, double b)
Definition: EvtConExc.cc:1811
void getResMass()
Definition: EvtConExc.cc:2096
double getVP(double cms)
Definition: EvtConExc.cc:2334
static double XS_max
Definition: EvtConExc.hh:138
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
static EvtId evtIdFromStdHep(int stdhep)
Definition: EvtPDL.cc:244
double getXlw()
Definition: EvtXsection.hh:153
std::string getMsg()
Definition: EvtXsection.hh:154
std::vector< double > getEr()
Definition: EvtXsection.hh:151
double getXup()
Definition: EvtXsection.hh:152
std::string getUnit()
Definition: EvtXsection.hh:147
std::vector< double > getYY()
Definition: EvtXsection.hh:150
void setBW(int pdg)
std::vector< double > getXX()
Definition: EvtXsection.hh:149

◆ init_Br_ee()

void EvtConExc::init_Br_ee ( )

Definition at line 1710 of file EvtConExc.cc.

1710 {
1711 // double psipp_ee=9.6E-06;
1712 double psip_ee =7.73E-03;
1713 double jsi_ee =5.94*0.01;
1714 double phi_ee =2.954E-04;
1715 // double omega_ee=7.28E-05;
1716 // double rho_ee = 4.72E-05;
1717 EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
1718 EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
1719 EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
1720 EvtId phiId=EvtPDL::getId(std::string("phi"));
1721 EvtId omegaId=EvtPDL::getId(std::string("omega"));
1722 EvtId rhoId=EvtPDL::getId(std::string("rho0"));
1723 BR_ee.clear();
1724 ResId.clear();
1725
1726 //BR_ee.push_back(rho_ee);
1727 //BR_ee.push_back(omega_ee);
1728 BR_ee.push_back(phi_ee);
1729 BR_ee.push_back(jsi_ee);
1730 BR_ee.push_back(psip_ee);
1731 //BR_ee.push_back(psipp_ee);
1732
1733 //ResId.push_back(rhoId);
1734 //ResId.push_back(omegaId);
1735 ResId.push_back(phiId);
1736 ResId.push_back(jsiId);
1737 ResId.push_back(pspId);
1738 //ResId.push_back(psppId);
1739
1740}

Referenced by init().

◆ init_mode()

void EvtConExc::init_mode ( int  mode)

Definition at line 576 of file EvtConExc.cc.

576 {
577 if(mode ==0 ){
578 _ndaugs =2;
579 daugs[0]=EvtPDL::getId(std::string("p+"));
580 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
581 }else if(mode ==1 ){
582 _ndaugs =2;
583 daugs[0]=EvtPDL::getId(std::string("n0"));
584 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
585 }else if(mode ==2 ){
586 _ndaugs =2;
587 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
588 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
589 }else if(mode ==3 ){
590 _ndaugs =2;
591 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
592 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
593 }else if(mode ==4 ){
594 _ndaugs =2;
595 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
596 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
597 }else if(mode ==5 ){
598 _ndaugs =2;
599 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
600 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
601 } else if(mode ==6 ){
602 _ndaugs =2;
603 daugs[0]=EvtPDL::getId(std::string("pi+"));
604 daugs[1]=EvtPDL::getId(std::string("pi-"));
605 }else if(mode ==7 ){
606 _ndaugs =3;
607 daugs[0]=EvtPDL::getId(std::string("pi+"));
608 daugs[1]=EvtPDL::getId(std::string("pi-"));
609 daugs[2]=EvtPDL::getId(std::string("pi0"));
610 }else if(mode ==8 ){
611 _ndaugs =3;
612 daugs[0]=EvtPDL::getId(std::string("K+"));
613 daugs[1]=EvtPDL::getId(std::string("K-"));
614 daugs[2]=EvtPDL::getId(std::string("pi0"));
615 }else if(mode ==9 ){
616 _ndaugs =3;
617 daugs[0]=EvtPDL::getId(std::string("K_S0"));
618 daugs[1]=EvtPDL::getId(std::string("K+"));
619 daugs[2]=EvtPDL::getId(std::string("pi-"));
620 }else if(mode ==10 ){
621 _ndaugs =3;
622 daugs[0]=EvtPDL::getId(std::string("K_S0"));
623 daugs[1]=EvtPDL::getId(std::string("K-"));
624 daugs[2]=EvtPDL::getId(std::string("pi+"));
625 }else if(mode ==11 ){
626 _ndaugs =3;
627 daugs[0]=EvtPDL::getId(std::string("K+"));
628 daugs[1]=EvtPDL::getId(std::string("K+"));
629 daugs[2]=EvtPDL::getId(std::string("eta"));
630 }else if(mode ==12 ){
631 _ndaugs =4;
632 daugs[0]=EvtPDL::getId(std::string("pi+"));
633 daugs[1]=EvtPDL::getId(std::string("pi-"));
634 daugs[2]=EvtPDL::getId(std::string("pi+"));
635 daugs[3]=EvtPDL::getId(std::string("pi-"));
636 }else if(mode ==13 ){
637 _ndaugs =4;
638 daugs[0]=EvtPDL::getId(std::string("pi+"));
639 daugs[1]=EvtPDL::getId(std::string("pi-"));
640 daugs[2]=EvtPDL::getId(std::string("pi0"));
641 daugs[3]=EvtPDL::getId(std::string("pi0"));
642 }else if(mode ==14 ){
643 _ndaugs =4;
644 daugs[0]=EvtPDL::getId(std::string("K+"));
645 daugs[1]=EvtPDL::getId(std::string("K-"));
646 daugs[2]=EvtPDL::getId(std::string("pi+"));
647 daugs[3]=EvtPDL::getId(std::string("pi-"));
648 }else if(mode ==15 ){
649 _ndaugs =4;
650 daugs[0]=EvtPDL::getId(std::string("K+"));
651 daugs[1]=EvtPDL::getId(std::string("K-"));
652 daugs[2]=EvtPDL::getId(std::string("pi0"));
653 daugs[3]=EvtPDL::getId(std::string("pi0"));
654 }else if(mode ==16 ){
655 _ndaugs =4;
656 daugs[0]=EvtPDL::getId(std::string("K+"));
657 daugs[1]=EvtPDL::getId(std::string("K-"));
658 daugs[2]=EvtPDL::getId(std::string("K+"));
659 daugs[3]=EvtPDL::getId(std::string("K-"));
660 }else if(mode ==17 ){
661 _ndaugs =5;
662 daugs[0]=EvtPDL::getId(std::string("pi+"));
663 daugs[1]=EvtPDL::getId(std::string("pi-"));
664 daugs[2]=EvtPDL::getId(std::string("pi+"));
665 daugs[3]=EvtPDL::getId(std::string("pi-"));
666 daugs[4]=EvtPDL::getId(std::string("pi0"));
667 }else if(mode ==18 ){
668 _ndaugs =5;
669 daugs[0]=EvtPDL::getId(std::string("pi+"));
670 daugs[1]=EvtPDL::getId(std::string("pi-"));
671 daugs[2]=EvtPDL::getId(std::string("pi+"));
672 daugs[3]=EvtPDL::getId(std::string("pi-"));
673 daugs[4]=EvtPDL::getId(std::string("eta"));
674 }else if(mode ==19 ){
675 _ndaugs =5;
676 daugs[0]=EvtPDL::getId(std::string("K+"));
677 daugs[1]=EvtPDL::getId(std::string("K-"));
678 daugs[2]=EvtPDL::getId(std::string("pi+"));
679 daugs[3]=EvtPDL::getId(std::string("pi-"));
680 daugs[4]=EvtPDL::getId(std::string("pi0"));
681 }else if(mode ==20 ){
682 _ndaugs =5;
683 daugs[0]=EvtPDL::getId(std::string("K+"));
684 daugs[1]=EvtPDL::getId(std::string("K-"));
685 daugs[2]=EvtPDL::getId(std::string("pi+"));
686 daugs[3]=EvtPDL::getId(std::string("pi-"));
687 daugs[4]=EvtPDL::getId(std::string("eta"));
688 }else if(mode ==21 ){
689 _ndaugs =6;
690 daugs[0]=EvtPDL::getId(std::string("pi+"));
691 daugs[1]=EvtPDL::getId(std::string("pi-"));
692 daugs[2]=EvtPDL::getId(std::string("pi+"));
693 daugs[3]=EvtPDL::getId(std::string("pi-"));
694 daugs[4]=EvtPDL::getId(std::string("pi+"));
695 daugs[5]=EvtPDL::getId(std::string("pi-"));
696 }else if(mode ==22 ){
697 _ndaugs =6;
698 daugs[0]=EvtPDL::getId(std::string("pi+"));
699 daugs[1]=EvtPDL::getId(std::string("pi-"));
700 daugs[2]=EvtPDL::getId(std::string("pi+"));
701 daugs[3]=EvtPDL::getId(std::string("pi-"));
702 daugs[4]=EvtPDL::getId(std::string("pi0"));
703 daugs[5]=EvtPDL::getId(std::string("pi0"));
704 }else if(mode == 23){
705 _ndaugs =2;
706 daugs[0]=EvtPDL::getId(std::string("eta"));
707 daugs[1]=EvtPDL::getId(std::string("phi"));
708 }else if(mode == 24){
709 _ndaugs =2;
710 daugs[0]=EvtPDL::getId(std::string("phi"));
711 daugs[1]=EvtPDL::getId(std::string("pi0"));
712 }else if(mode == 25){
713 _ndaugs =2;
714 daugs[0]=EvtPDL::getId(std::string("K+"));
715 daugs[1]=EvtPDL::getId(std::string("K*-"));
716 }else if(mode == 26){
717 _ndaugs =2;
718 daugs[0]=EvtPDL::getId(std::string("K-"));
719 daugs[1]=EvtPDL::getId(std::string("K*+"));
720 }else if(mode == 27){
721 _ndaugs =2;
722 daugs[0]=EvtPDL::getId(std::string("K_S0"));
723 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
724 }else if(mode == 28){
725 _ndaugs =3;
726 daugs[0]=EvtPDL::getId(std::string("K*0"));
727 daugs[1]=EvtPDL::getId(std::string("K+"));
728 daugs[2]=EvtPDL::getId(std::string("pi-"));
729 }else if(mode == 29){
730 _ndaugs =3;
731 daugs[0]=EvtPDL::getId(std::string("K*0"));
732 daugs[1]=EvtPDL::getId(std::string("K-"));
733 daugs[2]=EvtPDL::getId(std::string("pi+"));
734 }else if(mode == 30){
735 _ndaugs =3;
736 daugs[0]=EvtPDL::getId(std::string("K*+"));
737 daugs[1]=EvtPDL::getId(std::string("K-"));
738 daugs[2]=EvtPDL::getId(std::string("pi0"));
739 }else if(mode == 31){
740 _ndaugs =3;
741 daugs[0]=EvtPDL::getId(std::string("K*-"));
742 daugs[1]=EvtPDL::getId(std::string("K+"));
743 daugs[2]=EvtPDL::getId(std::string("pi0"));
744 }else if(mode == 32){
745 _ndaugs =3;
746 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
747 daugs[1]=EvtPDL::getId(std::string("K+"));
748 daugs[2]=EvtPDL::getId(std::string("pi-"));
749 }else if(mode == 33){
750 _ndaugs =3;
751 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
752 daugs[1]=EvtPDL::getId(std::string("K-"));
753 daugs[2]=EvtPDL::getId(std::string("pi+"));
754 }else if(mode == 34){
755 _ndaugs =3;
756 daugs[0]=EvtPDL::getId(std::string("K+"));
757 daugs[1]=EvtPDL::getId(std::string("K-"));
758 daugs[2]=EvtPDL::getId(std::string("rho0"));
759 }else if(mode == 35){
760 _ndaugs =3;
761 daugs[0]=EvtPDL::getId(std::string("phi"));
762 daugs[1]=EvtPDL::getId(std::string("pi-"));
763 daugs[2]=EvtPDL::getId(std::string("pi+"));
764 }else if(mode == 36){
765 _ndaugs =2;
766 daugs[0]=EvtPDL::getId(std::string("phi"));
767 daugs[1]=EvtPDL::getId(std::string("f_0"));
768 }else if(mode == 37){
769 _ndaugs =3;
770 daugs[0]=EvtPDL::getId(std::string("eta"));
771 daugs[1]=EvtPDL::getId(std::string("pi+"));
772 daugs[2]=EvtPDL::getId(std::string("pi-"));
773 }else if(mode == 38){
774 _ndaugs =3;
775 daugs[0]=EvtPDL::getId(std::string("omega"));
776 daugs[1]=EvtPDL::getId(std::string("pi+"));
777 daugs[2]=EvtPDL::getId(std::string("pi-"));
778 }else if(mode == 39){
779 _ndaugs =2;
780 daugs[0]=EvtPDL::getId(std::string("omega"));
781 daugs[1]=EvtPDL::getId(std::string("f_0"));
782 }else if(mode == 40){
783 _ndaugs =3;
784 daugs[0]=EvtPDL::getId(std::string("eta'"));
785 daugs[1]=EvtPDL::getId(std::string("pi+"));
786 daugs[2]=EvtPDL::getId(std::string("pi-"));
787 }else if(mode == 41){
788 _ndaugs =3;
789 daugs[0]=EvtPDL::getId(std::string("f_1"));
790 daugs[1]=EvtPDL::getId(std::string("pi+"));
791 daugs[2]=EvtPDL::getId(std::string("pi-"));
792 }else if(mode == 42){
793 _ndaugs =3;
794 daugs[0]=EvtPDL::getId(std::string("omega"));
795 daugs[1]=EvtPDL::getId(std::string("K+"));
796 daugs[2]=EvtPDL::getId(std::string("K-"));
797 }else if(mode == 43){
798 _ndaugs =4;
799 daugs[0]=EvtPDL::getId(std::string("omega"));
800 daugs[1]=EvtPDL::getId(std::string("pi+"));
801 daugs[2]=EvtPDL::getId(std::string("pi-"));
802 daugs[3]=EvtPDL::getId(std::string("pi0"));
803 }else if(mode == 44){
804 _ndaugs =2;
805 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
806 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
807 }else if(mode == 45){
808 _ndaugs =2;
809 daugs[0]=EvtPDL::getId(std::string("K+"));
810 daugs[1]=EvtPDL::getId(std::string("K-"));
811 }else if(mode == 46){
812 _ndaugs =2;
813 daugs[0]=EvtPDL::getId(std::string("K_S0"));
814 daugs[1]=EvtPDL::getId(std::string("K_L0"));
815 }else if(mode == 47){
816 _ndaugs =2;
817 daugs[0]=EvtPDL::getId(std::string("omega"));
818 daugs[1]=EvtPDL::getId(std::string("eta"));
819 }else if(mode == 48){
820 _ndaugs =2;
821 daugs[0]=EvtPDL::getId(std::string("phi"));
822 daugs[1]=EvtPDL::getId(std::string("eta'"));
823 }else if(mode == 49){
824 _ndaugs =3;
825 daugs[0]=EvtPDL::getId(std::string("phi"));
826 daugs[1]=EvtPDL::getId(std::string("K+"));
827 daugs[2]=EvtPDL::getId(std::string("K-"));
828 }else if(mode == 50){
829 _ndaugs =3;
830 daugs[0]=EvtPDL::getId(std::string("p+"));
831 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
832 daugs[2]=EvtPDL::getId(std::string("pi0"));
833 }else if(mode == 51){
834 _ndaugs =3;
835 daugs[0]=EvtPDL::getId(std::string("p+"));
836 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
837 daugs[2]=EvtPDL::getId(std::string("eta"));
838 }else if(mode == 65){
839 _ndaugs =3;
840 daugs[0]=EvtPDL::getId(std::string("D-"));
841 daugs[1]=EvtPDL::getId(std::string("D*0"));
842 daugs[2]=EvtPDL::getId(std::string("pi+"));
843 }else if(mode == 66){
844 _ndaugs =3;
845 daugs[0]=EvtPDL::getId(std::string("D+"));
846 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
847 daugs[2]=EvtPDL::getId(std::string("pi-"));
848 }else if(mode == 67){
849 _ndaugs =2;
850 daugs[0]=EvtPDL::getId(std::string("D*0"));
851 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
852 }else if(mode == 68){
853 _ndaugs =2;
854 daugs[0]=EvtPDL::getId(std::string("D0"));
855 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
856
857 }else if(mode == 69){
858 _ndaugs =2;
859 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
860 daugs[1]=EvtPDL::getId(std::string("D*0"));
861
862 }else if(mode == 70){
863 _ndaugs =2;
864 daugs[0]=EvtPDL::getId(std::string("D0"));
865 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
866
867}else if(mode == 71){
868 _ndaugs =2;
869 daugs[0]=EvtPDL::getId(std::string("D+"));
870 daugs[1]=EvtPDL::getId(std::string("D-"));
871 }else if(mode == 72){
872 _ndaugs =2;
873 daugs[0]=EvtPDL::getId(std::string("D+"));
874 daugs[1]=EvtPDL::getId(std::string("D*-"));
875
876 }else if(mode == 73){
877 _ndaugs =2;
878 daugs[0]=EvtPDL::getId(std::string("D-"));
879 daugs[1]=EvtPDL::getId(std::string("D*+"));
880
881 }else if(mode == 74){
882 _ndaugs =2;
883 daugs[0]=EvtPDL::getId(std::string("D*+"));
884 daugs[1]=EvtPDL::getId(std::string("D*-"));
885
886 }else if(mode == 75){
887 _ndaugs =3;
888 daugs[0]=EvtPDL::getId(std::string("D0"));
889 daugs[1]=EvtPDL::getId(std::string("D-"));
890 daugs[2]=EvtPDL::getId(std::string("pi+"));
891 }else if(mode == 76){
892 _ndaugs =3;
893 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
894 daugs[1]=EvtPDL::getId(std::string("D+"));
895 daugs[2]=EvtPDL::getId(std::string("pi-"));
896 }else if(mode == 77){
897 _ndaugs =3;
898 daugs[0]=EvtPDL::getId(std::string("D0"));
899 daugs[1]=EvtPDL::getId(std::string("D*-"));
900 daugs[2]=EvtPDL::getId(std::string("pi+"));
901 }else if(mode == 78){
902 _ndaugs =3;
903 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
904 daugs[1]=EvtPDL::getId(std::string("D*+"));
905 daugs[2]=EvtPDL::getId(std::string("pi-"));
906 }else if(mode == 79){
907 _ndaugs =3;
908 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
909 daugs[1]=EvtPDL::getId(std::string("pi0"));
910 daugs[2]=EvtPDL::getId(std::string("pi0"));
911 }else if(mode == 80){
912 _ndaugs =2;
913 daugs[0]=EvtPDL::getId(std::string("eta"));
914 daugs[1]=EvtPDL::getId(std::string("J/psi"));
915 }else if(mode == 81){
916 _ndaugs =3;
917 daugs[0]=EvtPDL::getId(std::string("h_c"));
918 daugs[1]=EvtPDL::getId(std::string("pi+"));
919 daugs[2]=EvtPDL::getId(std::string("pi-"));
920 }else if(mode == 82){
921 _ndaugs =3;
922 daugs[0]=EvtPDL::getId(std::string("h_c"));
923 daugs[1]=EvtPDL::getId(std::string("pi0"));
924 daugs[2]=EvtPDL::getId(std::string("pi0"));
925 }else if(mode == 83){
926 _ndaugs =3;
927 daugs[0]=EvtPDL::getId(std::string("J/psi"));
928 daugs[1]=EvtPDL::getId(std::string("K+"));
929 daugs[2]=EvtPDL::getId(std::string("K-"));
930 }else if(mode == 84){
931 _ndaugs =3;
932 daugs[0]=EvtPDL::getId(std::string("J/psi"));
933 daugs[1]=EvtPDL::getId(std::string("K_S0"));
934 daugs[2]=EvtPDL::getId(std::string("K_S0"));
935 }else if(mode == 90){
936 _ndaugs =3;
937 daugs[0]=EvtPDL::getId(std::string("J/psi"));
938 daugs[1]=EvtPDL::getId(std::string("pi+"));
939 daugs[2]=EvtPDL::getId(std::string("pi-"));
940 }else if(mode == 91){
941 _ndaugs =3;
942 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
943 daugs[1]=EvtPDL::getId(std::string("pi+"));
944 daugs[2]=EvtPDL::getId(std::string("pi-"));
945 }else if(mode == 92){
946 _ndaugs =3;
947 daugs[0]=EvtPDL::getId(std::string("J/psi"));
948 daugs[1]=EvtPDL::getId(std::string("K+"));
949 daugs[2]=EvtPDL::getId(std::string("K-"));
950 }else if(mode == 93){
951 _ndaugs =2;
952 daugs[0]=EvtPDL::getId(std::string("D_s+"));
953 daugs[1]=EvtPDL::getId(std::string("D_s-"));
954 }else if(mode == 94){
955 _ndaugs =2;
956 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
957 daugs[1]=EvtPDL::getId(std::string("D_s-"));
958 }else if(mode == 95){
959 _ndaugs =2;
960 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
961 daugs[1]=EvtPDL::getId(std::string("D_s+"));
962 }else if(mode == 96){
963 _ndaugs =2;
964 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
965 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
966 }else if(mode == 74100){
967 _ndaugs =1;
968 daugs[0]=EvtPDL::getId(std::string("J/psi"));
969 }else if(mode == 74101){
970 _ndaugs =1;
971 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
972 }else if(mode == 74102){
973 _ndaugs =1;
974 daugs[0]=EvtPDL::getId(std::string("psi(3770)"));
975 }else if(mode == 74103){
976 _ndaugs =1;
977 daugs[0]=EvtPDL::getId(std::string("phi"));
978 }else if(mode == 74104){
979 _ndaugs =1;
980 daugs[0]=EvtPDL::getId(std::string("omega"));
981 }else if(mode == 74105){
982 _ndaugs =1;
983 daugs[0]=EvtPDL::getId(std::string("rho0"));
984 }else if(mode == 74106){
985 _ndaugs =1;
986 daugs[0]=EvtPDL::getId(std::string("rho(3S)0"));
987 }else if(mode == 74107){
988 _ndaugs =1;
989 daugs[0]=EvtPDL::getId(std::string("omega(2S)"));
990 }else if(mode == 74110){
991 _ndaugs =1;
992 daugs[0]=EvtPDL::getId(std::string("gamma*"));
993 EvtId myvpho = EvtPDL::getId(std::string("gamma*"));
994 EvtPDL::reSetMass(myvpho,mhds*0.9); //for calculating maxXS, it will be resent when this mode is selected
995 _modeFlag.clear();
996 _modeFlag.push_back(74110); //R-value generator tag
997 _modeFlag.push_back(0); //to contain the mode selected
998 }else if(mode == -1){
999 _ndaugs = nson;
1000 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1001 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
1002 }else if(mode == -2){
1003 _ndaugs = nson;
1004 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1005 }else if(mode == 10000){//use for check ISR
1006 _ndaugs =2;
1007 daugs[0]=EvtPDL::getId(std::string("pi+"));
1008 daugs[1]=EvtPDL::getId(std::string("pi-"));
1009 }else {
1010 std::cout<<"Bad mode_index number " <<mode<<", please refer to the manual."<<std::endl;
1011 ::abort();
1012 }
1013 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1014
1015 if(VISR){
1016 _ndaugs=2;
1017 daugs[0]=getDaugs()[0];//EvtPDL::getId(std::string("gamma"));
1018 daugs[1]=getDaugs()[1];//EvtPDL::getId(std::string("gamma*"));
1019 }
1020
1021}

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

◆ initProbMax()

void EvtConExc::initProbMax ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 1024 of file EvtConExc.cc.

1024 {
1025
1026 noProbMax();
1027
1028}
void noProbMax()

◆ islgr()

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

Definition at line 1857 of file EvtConExc.cc.

1858{ int n0=-1;
1859 double z;
1860 for(int i=0;i<n-1;i++){
1861 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
1862 }
1863 double xstotal=y[599];
1864 if(n0==-1) {return 0;}else{
1865 double p1=y[n0]/xstotal;
1866 double p2=y[n0+1]/xstotal;
1867 double pm= EvtRandom::Flat(0.,1);
1868 if(p1<pm && pm<=p2) {return 1;}else{return 0;}
1869 }
1870}
const Int_t n
int t()
Definition: t.c:1

◆ ISR_ang_integrate()

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

Definition at line 1911 of file EvtConExc.cc.

1911 {
1912 //see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
1913 double costheta=cos(theta);
1914 double eb=_cms/2;
1915 double cos2=costheta*costheta;
1916 double sin2=1-cos2;
1917 double me=0.51*0.001;
1918 double L=2*log(_cms/me);
1919 double meE2= me*me/eb/eb;
1920 double hpi=L-1;
1921 double hq1= meE2/2*costheta/(sin2+meE2*cos2);
1922 double hq2= -0.5*log((1+sqrt(1-meE2)*costheta)/(1-sqrt(1-meE2)*costheta));
1923 double hq3=x*x*costheta/2/(x*x-2*x+2)*(1-meE2/(sin2+meE2*cos2));
1924 double hq=(L-1)/2.+hq1+hq2+hq3;
1925 hq /= hpi; //Eq (11) in arXif:hep-ph/9910523
1926 return hq;
1927}
const double me
Definition: PipiJpsi.cxx:48

Referenced by ISR_ang_sampling().

◆ ISR_ang_sampling()

double EvtConExc::ISR_ang_sampling ( double  x)

Definition at line 1929 of file EvtConExc.cc.

1929 {
1930 double xx[180],yy[180];
1931 double dgr2rad=1./180.*3.1415926;
1932 xx[0]=0;
1933 xx[1]=5*dgr2rad; //first bin is 5 degree
1934 int nc=2;
1935 for(int i=6;i<=175;i++){ //last bin is 5 degree
1936 xx[nc]=i*dgr2rad;
1937 nc++;
1938 }
1939 xx[nc]=180*dgr2rad;
1940 for(int j=0;j<=nc;j++){
1941 yy[j]=ISR_ang_integrate(x,xx[j]);
1942 }
1943 double pm= EvtRandom::Flat(0.,1);
1944 int mypt=1;
1945 for(int j=1;j<=nc;j++){
1946 if(yy[j-1]<pm && pm<=yy[j]){mypt=j;break;}
1947 }
1948 pm= EvtRandom::Flat(0.,1);
1949 double mytheta= xx[mypt-1]+ (xx[mypt]-xx[mypt-1])*pm;
1950 return mytheta; //theta in rad unit
1951}
double ISR_ang_integrate(double x, double theta)
Definition: EvtConExc.cc:1911

Referenced by decay().

◆ lgr()

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

Definition at line 1842 of file EvtConExc.cc.

1843{ int n0=-1;
1844 double z;
1845 for(int i=0;i<n-1;i++){
1846 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
1847 }
1848 if(n0==-1) {return 0.0;}else{
1849 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
1850 z= y[n0+1]+k*(t-x[n0+1]);
1851 }
1852 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
1853 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
1854 return z;
1855}

◆ Li2()

double EvtConExc::Li2 ( double  x)

Definition at line 1836 of file EvtConExc.cc.

1836 {
1837 double li2= -x +x*x/4. - x*x*x/9;
1838 return li2;
1839}

Referenced by SoftPhoton_xs().

◆ LLr()

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

Definition at line 1872 of file EvtConExc.cc.

1873{ int n0=-1;
1874 double z;
1875 if( t==x[n-1] ) return y[n-1];
1876 for(int i=0;i<n-1;i++){
1877 if(x[i]<=t && t< x[i+1]) {n0=i;break;}
1878 }
1879 if(n0==-1) {return 0.0;}else{
1880 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
1881 z= y[n0+1]+k*(t-x[n0+1]);
1882 }
1883 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
1884 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
1885 return z;
1886}

◆ meson_sampling()

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

Definition at line 1581 of file EvtConExc.cc.

1581 {
1582 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1583 double theta=angles.getHelAng(1);
1584 double phi =angles.getHelAng(2);
1585 double costheta=cos(theta); //using helicity angles in parent system
1586
1587 double pm= EvtRandom::Flat(0.,1);
1588 double ang = 1 - costheta*costheta;
1589 if(pm< ang/1.) {return true;}else{return false;}
1590}

Referenced by hadron_angle_sampling().

◆ Mhad_sampling()

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

Definition at line 1888 of file EvtConExc.cc.

1888 {//sample ISR hadrons, return Mhadron
1889 //x=MH: array contains the Mhad
1890 //y=AF: array contains the Mhad*Xsection
1891 //n: specify how many bins in the hadron sampling
1892 int n=598; //AF[599] is the total xs including Ecms point
1893 double pm= EvtRandom::Flat(0.,1);
1894 int mybin=1;
1895 double xst=y[n];
1896 for(int i=0;i<n;i++){
1897 if((y[i]/xst)<pm && pm<=(y[i+1]/xst)){mybin=i+1;break;}
1898 }
1899 pm= EvtRandom::Flat(0.,1);
1900 double mhd=x[mybin-1]+(x[mybin]-x[mybin-1])*pm;
1901
1902 if(mhd>_cms) {std::cout<<"selected mhd larger than Ecms "<<mhd<<" > "<<x[mybin]<<std::endl;abort();}
1903 if(mhd<_mhdL){
1904 std::cout<<"the sample mhassrons is less than the low end of XS"<<mhd<<" <"<<_mhdL<<std::endl;
1905 for(int i=0;i<598;i++){std::cout<<i<<" "<<x[i]<<" "<<y[i]<<std::endl;}
1906 abort();
1907 }
1908 return mhd;
1909}

Referenced by decay().

◆ narrowRXS()

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

Definition at line 2266 of file EvtConExc.cc.

2266 {
2267 //br for ee
2268 double psippee,psipee,jsiee,phiee,omegaee,rhoee;
2269 double kev2Gev=0.000001;
2270 psippee=0.262*kev2Gev;
2271 psipee =2.36*kev2Gev;
2272 jsiee =5.55*kev2Gev;
2273 phiee=4.266*0.001*2.954*0.0001;
2274 omegaee=0.6*kev2Gev;
2275 rhoee=7.04*kev2Gev;
2276 double s=_cms*_cms;
2277 double sigv=0;
2278 double mv=0;
2279 double widee=0;
2280 double xpi=12*3.1415926*3.1415926;
2281 if(mxL<=3.0957 && 3.0957<=mxH || mxL<=3.0983 && 3.0983<=mxH) {
2282 widee = jsiee;
2283 mv = 3.096916;
2284 }else if(mxL<=3.6847 && 3.6847<=mxH || mxL<=3.6873 && 3.6873<=mxH){
2285 widee = psipee;
2286 mv = 3.686109;
2287 }else if(mxL<=1.01906 && 1.01906<=mxH || mxL<= 1.01986 && 1.01986<=mxH){
2288 widee = phiee;
2289 mv = 1.01946;
2290 }else{return -1.0;}
2291
2292 if(_cms<=mv) return -1.;
2293 double xv = 1-mv*mv/s;
2294 sigv += xpi*widee/mv/s*Rad2(s,xv);
2295 double unic = 3.89379304*100000; //in unit nb
2296 return sigv*unic; //vaccum factor has included in the Rad2
2297}
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rh2ee DOUBLE PRECISION omepi DOUBLE PRECISION om3ee DOUBLE PRECISION phiee
Definition: RRes.h:37
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rhoee
Definition: RRes.h:34
double Rad2(double s, double x)
Definition: EvtConExc.cc:1615

◆ photonSampling()

bool EvtConExc::photonSampling ( EvtParticle part)

Definition at line 2235 of file EvtConExc.cc.

2235 {
2236 bool tagp,tagk;
2237 tagk=0;
2238 tagp=0;
2239 int nds = par->getNDaug();
2240 for(int i=0;i<par->getNDaug();i++){
2241 EvtId di=par->getDaug(i)->getId();
2242 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
2243 int pdgcode =EvtPDL::getStdHep( di );
2244 double alpha=0;
2245
2246 if(pdgcode==111 ||pdgcode==221 || pdgcode==331 ){//for photon
2247 alpha = 1;
2248 }else continue;
2249
2250 double angmax = 1+alpha;
2251 double costheta = cos(p4i.theta());
2252 double ang=1+alpha*costheta*costheta;
2253 double xratio = ang/angmax;
2254 double xi=EvtRandom::Flat(0.,1);
2255 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
2256 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
2257 if(xi>xratio) return false;
2258 }//loop over duaghters
2259 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
2260 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
2261
2262 return true;
2263}

Referenced by decay().

◆ Rad1()

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

Definition at line 1603 of file EvtConExc.cc.

1603 { //first order correction
1604 //radiator correction upto second order, see Int. J. of Moder.Phys. A
1605 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
1606 double alpha = 1./137.;
1607 double pi=3.1415926;
1608 double me = 0.5*0.001; //e mass
1609 double sme = sqrt(s)/me;
1610 double L = 2*log(sme);
1611 double wsx = 2*alpha/pi/x*(L-1)*(1 - x + x*x/2 );
1612 return wsx;
1613}

Referenced by Rad1difXs().

◆ Rad1difXs()

double EvtConExc::Rad1difXs ( EvtParticle p)

Definition at line 1682 of file EvtConExc.cc.

1682 {// // first order xsection
1683 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1684 double summass = p->getDaug(0)->getP4().mass();
1685 for(int i=1;i<_ndaugs;i++){
1686 summass += p->getDaug(i)->getP4().mass();
1687 }
1688
1689 double cms = p->getP4().mass();
1690 double mrg = cms - summass;
1691 double pm= EvtRandom::Flat(0.,1);
1692 double mhs = pm*mrg + summass;
1693
1694 double s = cms * cms;
1695 double x = 1 - mhs*mhs/s;
1696 double wsx = Rad1(s,x);
1697
1698 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
1699
1700 double xs = myxsection->getXsection(mhs);
1701 if(xs>XS_max){XS_max = xs;}
1702 double difxs = 2*mhs/cms/cms * wsx *xs;
1703
1704 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1705 differ *= mrg; //Jacobi factor for variable m_hds
1706 difxs *= mrg;
1707 return difxs;
1708}
double Rad1(double s, double x)
Definition: EvtConExc.cc:1603

Referenced by gamHXSection().

◆ Rad2()

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

Definition at line 1615 of file EvtConExc.cc.

1615 {
1616 //radiator correction upto second order, see Int. J. of Moder.Phys. A, hep-ph/9910523, Eq(8)
1617 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
1618 double alpha = 1./137.;
1619 double pi=3.1415926;
1620 double me = 0.5*0.001; //e mass
1621 double xi2 = 1.64493407;
1622 double xi3=1.2020569;
1623 double sme = sqrt(s)/me;
1624 double L = 2*log(sme);
1625 double beta = 2*alpha/pi*(L-1);
1626 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.;
1627 double ap = alpha/pi;
1628 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
1629 double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
1630 double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
1631 wsx = wsx + beta*beta/8. *wsx2;
1632 double mymx = sqrt(s*(1-x));
1633 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
1634 return wsx*getVP(mymx);//vph is vaccum polarization factor
1635}

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

◆ Rad2difXs() [1/2]

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

Definition at line 1670 of file EvtConExc.cc.

1670 {// leading second order xsection
1671 double wsx = Rad2(s,x);
1672 double mhs = sqrt(s*(1-x));
1673 double xs = myxsection->getXsection(mhs);
1674 //if(testflag==1)std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
1675 if(xs>XS_max){XS_max = xs;}
1676 double difxs = wsx *xs;
1677 differ2 = wsx *(myxsection->getErr(mhs));
1678 return difxs;
1679}

◆ Rad2difXs() [2/2]

double EvtConExc::Rad2difXs ( EvtParticle p)

Definition at line 1639 of file EvtConExc.cc.

1639 {// leading second order xsection
1640 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1641 double summass = p->getDaug(0)->getP4().mass();
1642 EvtVector4R v4p=p->getDaug(0)->getP4();
1643 for(int i=1;i<_ndaugs;i++){
1644 summass += p->getDaug(i)->getP4().mass();
1645 v4p += p->getDaug(i)->getP4();
1646 }
1647
1648 double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
1649 double cms = p->getP4().mass();
1650 double mrg = cms - summass;
1651 double pm= EvtRandom::Flat(0.,1);
1652 double mhs = pm*mrg + summass;
1653
1654 double s = cms * cms;
1655 double x = 2*Egam/cms;
1656 //double mhs = sqrt( s*(1-x) );
1657 double wsx = Rad2(s,x);
1658
1659 // std::cout<<"random : "<<pm<<std::endl;
1660
1661 double xs = myxsection->getXsection(mhs);
1662 if(xs>XS_max){XS_max = xs;}
1663 double difxs = 2*mhs/cms/cms * wsx *xs;
1664 differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1665 differ *= mrg; //Jacobi factor for variable m_hds
1666 difxs *= mrg;
1667 return difxs;
1668}

Referenced by gamHXSection().

◆ ReadVP()

void EvtConExc::ReadVP ( )

Definition at line 2355 of file EvtConExc.cc.

2355 {
2356 vpx.clear();vpr.clear();vpi.clear();
2357 int vpsize=4001;
2358 vpx.resize(vpsize);
2359 vpr.resize(vpsize);
2360 vpi.resize(vpsize);
2361 std::string locvp=getenv("BESEVTGENROOT");
2362 locvp +="/share/vp.dat";
2363 ifstream m_inputFile;
2364 m_inputFile.open(locvp.c_str());
2365 double xx,r1,i1;
2366 double x1,y1,y2;
2367 xx=0;
2368 int mytag=1;
2369 for(int i=0;i<4001;i++){
2370 m_inputFile>>vpx[i]>>vpr[i]>>vpi[i];
2371 }
2372
2373}

Referenced by init().

◆ resetResMass()

void EvtConExc::resetResMass ( )

Definition at line 2061 of file EvtConExc.cc.

2061 {
2062 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
2063 EvtId myres = EvtPDL::getId(std::string("J/psi"));
2064 EvtPDL::reSetMass(myres,mjsi);
2065 //EvtPDL::reSetWidth(myres,wjsi);
2066
2067 myres = EvtPDL::getId(std::string("psi(2S)"));
2068 EvtPDL::reSetMass(myres,mpsip);
2069 //EvtPDL::reSetWidth(myres,wpsip);
2070
2071 myres = EvtPDL::getId(std::string("psi(3770)"));
2072 EvtPDL::reSetMass(myres,mpsipp);
2073 //EvtPDL::reSetWidth(myres,wpsipp);
2074
2075 myres = EvtPDL::getId(std::string("phi"));
2076 EvtPDL::reSetMass(myres,mphi);
2077 //EvtPDL::reSetWidth(myres,wphi);
2078
2079 myres = EvtPDL::getId(std::string("omega"));
2080 EvtPDL::reSetMass(myres,momega);
2081 //EvtPDL::reSetWidth(myres,womega);
2082
2083 myres = EvtPDL::getId(std::string("rho0"));
2084 EvtPDL::reSetMass(myres,mrho0);
2085 //EvtPDL::reSetWidth(myres,wrho0);
2086
2087 myres = EvtPDL::getId(std::string("rho(3S)0"));
2088 EvtPDL::reSetMass(myres,mrho3s);
2089 //EvtPDL::reSetWidth(myres,wrho3s);
2090
2091 myres = EvtPDL::getId(std::string("omega(2S)"));
2092 EvtPDL::reSetMass(myres,momega2s);
2093 //EvtPDL::reSetWidth(myres,womega2s);
2094}

Referenced by decay().

◆ Ros_xs()

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

Definition at line 1742 of file EvtConExc.cc.

1742 {// in unit nb
1743 double pi=3.1415926;
1744 double s= mx*mx;
1745 double width = EvtPDL::getWidth(pid);
1746 double mass = EvtPDL::getMeanMass(pid);
1747 double xv = 1-mass*mass/s;
1748 double sigma = 12*pi*pi*bree*width/mass/s;
1749 sigma *= Rad2(s, xv);
1750 double nbar = 3.89379304*100000; // ! GeV-2 = 3.89379304*10^5 nbar
1751 return sigma*nbar;
1752}
double mass

Referenced by init().

◆ selectMass()

double EvtConExc::selectMass ( )

Definition at line 2309 of file EvtConExc.cc.

2309 {
2310 double pm,mhdL,mhds;
2311 pm = EvtRandom::Flat(0.,1);
2312 mhdL =_mhdL;
2313 mhds = pm*(_cms - mhdL)+mhdL; //non narrow resonance
2314 std::vector<double> sxs;
2315 for(int i=0;i<ISRID.size();i++){
2316 double ri=ISRXS[i]/AF[598];
2317 sxs.push_back(ri);
2318 }
2319 for(int j=0;j<sxs.size();j++){
2320 if(j>0) sxs[j] += sxs[j-1];
2321 }
2322 pm = EvtRandom::Flat(0.,1);
2323 if(pm<sxs[0]) {
2324 mhds = EvtPDL::getMass(ISRID[0]); //narrow resonance
2325 }
2326 for(int j=1;j<sxs.size();j++){
2327 double x0 = sxs[j-1];
2328 double x1 = sxs[j];
2329 if(x0<pm && pm<=x1) mhds = EvtPDL::getMass(ISRID[j]); //narrow resonance
2330 }
2331 return mhds;
2332}

◆ selectMode()

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

Definition at line 2014 of file EvtConExc.cc.

2014 {
2015 //first get xsection for each mode in vmod array
2016 double uscale = 1;
2017 std::vector<double> vxs;vxs.clear();
2018 for (int i=0;i<vmod.size();i++){
2019 int imod = vmod[i];
2020 delete myxsection;
2021 myxsection =new EvtXsection (imod);
2022 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
2023 //if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
2024
2025 double myxs=uscale*myxsection->getXsection(mhds);
2026 if(imod==0) {vxs.push_back(myxs);}
2027 else if(imod==74110){//for continuum process
2028 double xcon = myxs - vxs[i-1];
2029 if(xcon<0) {xcon=vxs[i-1];}else{xcon=myxs;}
2030 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
2031 vxs.push_back(xcon);
2032 }else{
2033 vxs.push_back(myxs+vxs[i-1]);
2034 }
2035 }//for_i_block
2036 //degugging
2037 // 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;}
2038
2039 double totxs = vxs[vxs.size()-1];
2040 int icount=0;
2041 mode_iter:
2042 double pm= EvtRandom::Flat(0.,1);
2043 if(pm <= vxs[0]/totxs) return vmod[0];
2044 for(int j=1;j<vxs.size();j++){
2045 double x0 = vxs[j-1]/totxs;
2046 double x1 = vxs[j]/totxs;
2047 if(x0<pm && pm<=x1) return vmod[j];
2048 }
2049
2050 icount++;
2051 if(icount<10000) goto mode_iter;
2052
2053 //debugging
2054 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;}
2055
2056 std::cout<<"EvtConExc::selectMode can not find a mode with 100 iteration"<<std::endl;
2057 abort();
2058}

Referenced by decay().

◆ SetP4()

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

Definition at line 1953 of file EvtConExc.cc.

1953 { //set the gamma and gamma* momentum according sampled results
1954 double pm= EvtRandom::Flat(0.,1);
1955 double phi=2*3.1415926*pm;
1956 double gamE = _cms/2*xeng;
1957 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
1958 double px = gamE*sin(theta)*cos(phi);
1959 double py = gamE*sin(theta)*sin(phi);
1960 double pz = gamE*cos(theta);
1961 EvtVector4R p4[2];
1962 p4[0].set(gamE,px,py,pz);//ISR photon
1963 p4[1].set(hdrE,-px,-py,-pz); //mhdraons
1964 for(int i=0;i<2;i++){
1965 EvtId myid = p->getDaug(i)->getId();
1966 p->getDaug(i)->init(myid,p4[i]);
1967 }
1968}
double sin(const BesAngle a)

Referenced by decay().

◆ SetP4Rvalue()

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

Definition at line 1970 of file EvtConExc.cc.

1970 { //set the gamma and gamma* momentum according sampled results
1971 double pm= EvtRandom::Flat(0.,1);
1972 double phi=2*3.1415926*pm;
1973 double gamE = _cms/2*xeng;
1974 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
1975 double px = gamE*sin(theta)*cos(phi);
1976 double py = gamE*sin(theta)*sin(phi);
1977 double pz = gamE*cos(theta);
1978 EvtVector4R p4[2];
1979 p4[0].set(hdrE,px,py,pz);//mhdraons
1980 p4[1].set(gamE,-px,-py,-pz); //ISR photon
1981 for(int i=0;i<2;i++){
1982 EvtId myid = p->getDaug(i)->getId();
1983 p->getDaug(i)->init(myid,p4[i]);
1984 }
1985}

Referenced by decay().

◆ showResMass()

void EvtConExc::showResMass ( )

Definition at line 2132 of file EvtConExc.cc.

2132 {
2133 std::cout<<"J/psi: "<<mjsi<<" "<<wjsi<<std::endl;
2134 std::cout<<"psipr: "<<mpsip<<" "<<wpsip<<std::endl;
2135 std::cout<<"psipp: "<<mpsipp<<" "<<wpsipp<<std::endl;
2136 std::cout<<"phi : "<<mphi<<" "<<wphi<<std::endl;
2137 std::cout<<"omega: "<<momega<<" "<<womega<<std::endl;
2138 std::cout<<"rho0 : "<<mrho0<<" "<<wrho0<<std::endl;
2139 std::cout<<"rho(3S)0: "<<mrho3s<<" "<<wrho3s<<std::endl;
2140 std::cout<<"omega(2S): "<<momega2s<<" "<<womega2s<<std::endl;
2141}

◆ SoftPhoton_xs()

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

Definition at line 1811 of file EvtConExc.cc.

1811 {
1812 double alpha = 1./137.;
1813 double pi=3.1415926;
1814 double me = 0.5*0.001; //e mass
1815 double xi2 = 1.64493407;
1816 double xi3=1.2020569;
1817 double sme = sqrt(s)/me;
1818 double L = 2*log(sme);
1819 double beta = 2*alpha/pi*(L-1);
1820 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.;
1821 double ap = alpha/pi;
1822 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
1823
1824 double beta2 = beta*beta;
1825 double b2 = b*b;
1826
1827 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 -
1828 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) +
1829 16*pow(beta,2)*Li2(b))/32. ;
1830 double mymx = sqrt(s*(1-b));
1831 //return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
1832 return xs*getVP(mymx); //vph :the vaccum polarzation factor
1833
1834}
double Li2(double x)
Definition: EvtConExc.cc:1836

Referenced by init().

◆ sumExc()

double EvtConExc::sumExc ( double  mx)

Definition at line 2157 of file EvtConExc.cc.

2157 {//summation all cross section of exclusive decays
2158 std::vector<int> vmod; vmod.clear();
2159 int mn[78]={0,1,2,3,4,5,6,7,8,9,10,11,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,//mode 12 and 43 are removed
2160 50,51,
2161 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
2162 90,91,92,93,94,95,96,
2163 74100,74101,74102,74103,74104,74105,74110};
2164 int mn2[79]={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,//mode 12 and 43 are removed
2165 50,51,
2166 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
2167 90,91,92,93,94,95,96,
2168 74100,74101,74102,74103,74104,74105,74110};
2169 if(_cms > 2.5){
2170 for(int i=0;i<78;i++){vmod.push_back(mn[i]);}
2171 }else{
2172 for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
2173 }
2174 double xsum=0;
2175 double uscale = 1;
2176 for(int i=0;i<vmod.size();i++){
2177 int mymode = vmod[i];
2178 if(mymode ==74110) continue;
2179 delete myxsection;
2180 myxsection =new EvtXsection (mymode);
2181 init_mode(mymode);
2182 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
2183 xsum += uscale*myxsection->getXsection(mx);
2184 //if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
2185 }
2186 return xsum;
2187}

Referenced by init().

◆ VP_sampling()

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

Definition at line 1592 of file EvtConExc.cc.

1592 {
1593 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1594 double theta=angles.getHelAng(1);
1595 double phi =angles.getHelAng(2);
1596 double costheta=cos(theta); //using helicity angles in parent system
1597
1598 double pm= EvtRandom::Flat(0.,1);
1599 double ang = 1 + costheta*costheta;
1600 if(pm< ang/2.) {return true;}else{return false;}
1601}

Referenced by hadron_angle_sampling().

◆ xs_sampling() [1/2]

bool EvtConExc::xs_sampling ( double  xs)

Definition at line 1552 of file EvtConExc.cc.

1552 {
1553 double pm= EvtRandom::Flat(0.,1);
1554 //std::cout<<"Rdm= "<<pm<<std::endl;
1555 if(pm <xs/XS_max) {return true;} else {return false;}
1556}

◆ xs_sampling() [2/2]

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

Definition at line 1558 of file EvtConExc.cc.

1558 {
1559 double pm= EvtRandom::Flat(0.,1);
1560 // std::cout<<"Rdm= "<<pm<<std::endl;
1561 if(pm <xs/(xs0*1.1)) {return true;} else {return false;}
1562}

Member Data Documentation

◆ _cms

◆ conexcmode

int EvtConExc::conexcmode =-1
static

Definition at line 157 of file EvtConExc.hh.

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

◆ getMode

int EvtConExc::getMode
static

Definition at line 140 of file EvtConExc.hh.

Referenced by init().

◆ myxsection

◆ XS_max

double EvtConExc::XS_max
static

Definition at line 138 of file EvtConExc.hh.

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


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