BOSS 7.0.1
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Static Public Attributes | List of all members
EvtConExc Class Reference

#include <EvtConExc.hh>

+ Inheritance diagram for EvtConExc:

Public Member Functions

 EvtConExc ()
 
virtual ~EvtConExc ()
 
void getName (std::string &name)
 
EvtDecayBase * clone ()
 
void initProbMax ()
 
void init ()
 
void init_Br_ee ()
 
void decay (EvtParticle *p)
 
void init_mode (int mode)
 
double gamHXSection (EvtParticle *p, double El, double Eh, int nmc=100000)
 
double gamHXSection (double s, double El, double Eh, int nmc=100000)
 
double gamHXSection (double El, double Eh)
 
double gamHXSection_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 ()
 
EvtId * getDaugId ()
 
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 EvtDecayBase * clone ()=0
 
virtual void init ()
 
virtual void initProbMax ()
 
virtual std::string commandName ()
 
virtual void command (std::string cmd)
 
double getProbMax (double prob)
 
double resetProbMax (double prob)
 
 EvtDecayBase ()
 
virtual ~EvtDecayBase ()
 
virtual bool matchingDecay (const EvtDecayBase &other) const
 
EvtId getParentId ()
 
double getBranchingFraction ()
 
void disableCheckQ ()
 
void checkQ ()
 
int getNDaug ()
 
EvtId * getDaugs ()
 
EvtId getDaug (int i)
 
int getNArg ()
 
int getPHOTOS ()
 
void setPHOTOS ()
 
void setVerbose ()
 
void setSummary ()
 
double * getArgs ()
 
std::string * getArgsStr ()
 
double getArg (int j)
 
std::string getArgStr (int j)
 
std::string getModelName ()
 
int getDSum ()
 
int summary ()
 
int verbose ()
 
void saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
 
void printSummary ()
 
void setProbMax (double prbmx)
 
void noProbMax ()
 
void checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
 
void checkNDaug (int d1, int d2=-1)
 
void checkSpinParent (EvtSpinType::spintype sp)
 
void checkSpinDaughter (int d1, EvtSpinType::spintype sp)
 
virtual int nRealDaughters ()
 

Static Public Attributes

static EvtXsection * myxsection
 
static 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}
EvtConExc::myxsection
static EvtXsection * myxsection
Definition: EvtConExc.hh:136
EvtParticle::deleteTree
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}
alpha
const double alpha
Definition: FastVertexFit.cxx:4
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
EvtId
Definition: EvtId.hh:27
EvtId::getId
int getId() const
Definition: EvtId.hh:41
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:73
EvtVector4R::theta
double theta()
Definition: EvtVector4R.cc:249

Referenced by decay().

◆ baryon_sampling()

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

Definition at line 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}
EvtConExc::baryonAng
double baryonAng(double mx)
Definition: EvtConExc.cc:2225
EvtConExc::_cms
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}
mp
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}
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle::getNDaug
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:127

Referenced by decay().

◆ 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:
1127 EvtGlobalSet::ConExcPythia=1;
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}
EvtConExc::checkdecay
bool checkdecay(EvtParticle *p)
Definition: EvtConExc.cc:2143
EvtConExc::findMaxXS
void findMaxXS(EvtParticle *p)
Definition: EvtConExc.cc:1466
EvtConExc::conexcmode
static int conexcmode
Definition: EvtConExc.hh:157
EvtConExc::SetP4Rvalue
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:1970
EvtConExc::photonSampling
bool photonSampling(EvtParticle *part)
Definition: EvtConExc.cc:2235
EvtConExc::selectMode
int selectMode(std::vector< int > vmod, double mhds)
Definition: EvtConExc.cc:2014
EvtConExc::init_mode
void init_mode(int mode)
Definition: EvtConExc.cc:576
EvtConExc::angularSampling
bool angularSampling(EvtParticle *part)
Definition: EvtConExc.cc:2189
EvtConExc::ISR_ang_sampling
double ISR_ang_sampling(double x)
Definition: EvtConExc.cc:1929
EvtConExc::Mhad_sampling
double Mhad_sampling(double *x, double *y)
Definition: EvtConExc.cc:1888
EvtConExc::resetResMass
void resetResMass()
Definition: EvtConExc.cc:2061
EvtConExc::hadron_angle_sampling
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition: EvtConExc.cc:1322
EvtConExc::gam_sampling
bool gam_sampling(EvtParticle *p)
Definition: EvtConExc.cc:1538
EvtConExc::SetP4
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:1953
EvtDecayBase::getArg
double getArg(int j)
Definition: EvtDecayBase.cc:564
EvtGlobalSet::ConExcPythia
static bool ConExcPythia
Definition: EvtGlobalSet.hh:20
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45
EvtPDL::reSetWidth
static void reSetWidth(EvtId i, double width)
Definition: EvtPDL.hh:69
EvtPDL::reSetMass
static void reSetMass(EvtId i, double mass)
Definition: EvtPDL.hh:68
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
EvtPDL::getMass
static double getMass(EvtId i)
Definition: EvtPDL.hh:46
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
EvtParticle::setMass
void setMass(double m)
Definition: EvtParticle.hh:372
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1177
EvtParticle::setIntFlag
void setIntFlag(std::vector< int > vi)
Definition: EvtParticle.hh:151
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
EvtParticle::printTree
void printTree() const
Definition: EvtParticle.cc:897
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1071
EvtVector4R::mass
double mass() const
Definition: EvtVector4R.cc:39
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:179
EvtVector4R::d3mag
double d3mag() const
Definition: EvtVector4R.cc:186

◆ difgamXs() [1/2]

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

Definition at line 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}
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
EvtXsection::getXsection
double getXsection(double mx)
Definition: EvtXsection.cc:11959

◆ 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}
EvtXsection::getErr
double getErr(double mx)
Definition: EvtXsection.cc:11965

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}
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition: KarFin.h:34
EvtConExc::difgamXs
double difgamXs(EvtParticle *p)
Definition: EvtConExc.cc:1511
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition: EvtParticleFactory.cc:46

Referenced by decay(), and init().

◆ gam_sampling() [1/2]

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

Definition at line 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}
Rad2difXs2
double Rad2difXs2(double *mhs)
Definition: EvtConExc.cc:1782
dgauss_
double dgauss_(double(*fun)(double *), double *, double *, double *)
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
eps
EvtTensor3C eps(const EvtVector3R &v)
Definition: EvtTensor3C.cc:307
EvtConExc::Rad2difXs
double Rad2difXs(EvtParticle *p)
Definition: EvtConExc.cc: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}
s
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}
EvtConExc::Rad1difXs
double Rad1difXs(EvtParticle *p)
Definition: EvtConExc.cc:1682
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
EvtParticle::deleteDaughters
void deleteDaughters(bool keepChannel=false)
Definition: EvtParticle.cc:540

Referenced by init().

◆ gamHXSection_er()

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

Definition at line 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}
Rad2difXs_er2
double Rad2difXs_er2(double *mhs)
Definition: EvtConExc.cc:1797
Rad2difXs_er
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}
EvtPDL::getWidth
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}
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition: EvtDiracSpinor.cc:75
EvtConExc::VP_sampling
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1592
EvtConExc::meson_sampling
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1581
EvtConExc::baryon_sampling
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1564
EvtVector4R::set
void set(int i, double d)
Definition: EvtVector4R.hh:183

Referenced by decay().

◆ init()

void EvtConExc::init ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 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);
252 EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
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}
mstp
int mstp[200]
Definition: EvtPycont.cc:54
EvtConExc::getMode
static int getMode
Definition: EvtConExc.hh:140
EvtConExc::addNarrowRXS
double addNarrowRXS(double mhi, double binwidth)
Definition: EvtConExc.cc:2300
EvtConExc::init_Br_ee
void init_Br_ee()
Definition: EvtConExc.cc:1710
EvtConExc::ReadVP
void ReadVP()
Definition: EvtConExc.cc:2355
EvtConExc::Ros_xs
double Ros_xs(double mx, double bree, EvtId pid)
Definition: EvtConExc.cc:1742
EvtConExc::gamHXSection
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition: EvtConExc.cc:1340
EvtConExc::sumExc
double sumExc(double mx)
Definition: EvtConExc.cc:2157
EvtConExc::SoftPhoton_xs
double SoftPhoton_xs(double s, double b)
Definition: EvtConExc.cc:1811
EvtConExc::getResMass
void getResMass()
Definition: EvtConExc.cc:2096
EvtConExc::getVP
double getVP(double cms)
Definition: EvtConExc.cc:2334
EvtConExc::XS_max
static double XS_max
Definition: EvtConExc.hh:138
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtPDL::evtIdFromStdHep
static EvtId evtIdFromStdHep(int stdhep)
Definition: EvtPDL.cc:244
EvtXsection::getXlw
double getXlw()
Definition: EvtXsection.hh:153
EvtXsection::getMsg
std::string getMsg()
Definition: EvtXsection.hh:154
EvtXsection::getEr
std::vector< double > getEr()
Definition: EvtXsection.hh:151
EvtXsection::getXup
double getXup()
Definition: EvtXsection.hh:152
EvtXsection::getUnit
std::string getUnit()
Definition: EvtXsection.hh:147
EvtXsection::getYY
std::vector< double > getYY()
Definition: EvtXsection.hh:150
EvtXsection::setBW
void setBW(int pdg)
Definition: EvtXsection.cc:12067
EvtXsection::getXX
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}
EvtDecayBase::noProbMax
void noProbMax()
Definition: EvtDecayBase.cc:304

◆ 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}
n
const Int_t n
Definition: DataBase/tau_mode.c:65
t
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}
me
const double me
Definition: PipiJpsi.cxx:48

Referenced by ISR_ang_sampling().

◆ ISR_ang_sampling()

double EvtConExc::ISR_ang_sampling ( double  x)

Definition at line 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}
EvtConExc::ISR_ang_integrate
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}
phiee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rh2ee DOUBLE PRECISION omepi DOUBLE PRECISION om3ee DOUBLE PRECISION phiee
Definition: RRes.h:37
rhoee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rhoee
Definition: RRes.h:34
EvtConExc::Rad2
double Rad2(double s, double x)
Definition: EvtConExc.cc: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}
EvtConExc::Rad1
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}
mass
double mass
Definition: CosmicGenerator.cxx:138

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}
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210

Referenced by decay().

◆ SetP4Rvalue()

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

Definition at line 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}
EvtConExc::Li2
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

double EvtConExc::_cms
static

Definition at line 137 of file EvtConExc.hh.

Referenced by baryon_sampling(), decay(), difgamXs(), findMaxXS(), init(), ISR_ang_integrate(), Mhad_sampling(), narrowRXS(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), Rad2difXs_er2(), selectMass(), SetP4(), SetP4Rvalue(), and sumExc().

◆ 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

EvtXsection * EvtConExc::myxsection
static

Definition at line 136 of file EvtConExc.hh.

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

◆ 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: