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)
bool	gam_sampling (EvtParticle *p)
bool	xs_sampling (double xs)
bool	baryon_sampling (EvtVector4R pcm, EvtVector4R pi)
bool	meson_sampling (EvtVector4R pcm, EvtVector4R pi)
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)
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

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.

              {
  } //constructor

Referenced by clone().

~EvtConExc()

EvtConExc::~EvtConExc ( )

virtual

Definition at line 136 of file EvtConExc.cc.

                      {
  if(mydbg){
    // xs->Write();
    myfile->Write();
  }
  delete myxsection;
  gamH->deleteTree();
}

Member Function Documentation

baryon_sampling()

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

Definition at line 1011 of file EvtConExc.cc.

                                                              {
  EvtHelSys angles(pcm,pi);   //using helicity sys.angles
  double theta=angles.getHelAng(1);
  double phi  =angles.getHelAng(2);
  double costheta=cos(theta);  //using helicity angles in parent system
 
  // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
  //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
 
  double pm= EvtRandom::Flat(0.,1);
  double ang = 1 + costheta*costheta;
  if(pm< ang/2.) {return true;}else{return false;}
}

Referenced by decay().

clone()

EvtDecayBase * EvtConExc::clone ( )

virtual

Implements EvtDecayBase.

Definition at line 151 of file EvtConExc.cc.

                              {
 
  return new EvtConExc;
 
}

decay()

void EvtConExc::decay ( EvtParticle *  p )

virtual

Implements EvtDecayBase.

Definition at line 717 of file EvtConExc.cc.

                                     {
  //--- debugging
  //  std::cout<<"_mode= "<<_mode<<", XS_max= "<<XS_max<<std::endl;
  //------
 
  if((_xs0 + _xs1)==0) {std::cout<<"EvtConExc::zero total xsection"<<std::endl;::abort();}
  double pm= EvtRandom::Flat(0.,1);
  double xsratio = _xs0/(_xs0 + _xs1);
  int iflag; //flag = 0 for ee->hadrons, 1 for ee->gamma hadrons
  if(pm<xsratio){iflag = 0;}else{iflag=1;}
 
  if(radflag==1){xsratio=1;iflag=1;} // only generating gamma hadrons mode
 
   daugs[_ndaugs]=gamId;
   p->makeDaughters(_ndaugs+iflag,daugs); 
 
  for(int i=0;i< _ndaugs+iflag;i++){
    EvtParticle* di=p->getDaug(i);
    double xmi=EvtPDL::getMass(di->getId());
    di->setMass(xmi);
   }
 angular_sampling:
  double weight1 = p->initializePhaseSpace(  _ndaugs+iflag , daugs);
  bool gam_ang,byn_ang,msn_ang,xs_ang;
  EvtVector4R vhds=p->getDaug(0)->getP4();
  EvtVector4R hd1=vhds;
  EvtVector4R hd2 =p->getDaug(1)->getP4();;
  mass1 = hd1.mass();
  mass2 = hd2.mass();
  for(int i=1;i<_ndaugs;i++){
    vhds += p->getDaug(i)->getP4();
  }
  double xm = vhds.mass();
  double mxs=myxsection->getXsection(xm);
  xs_ang = xs_sampling(mxs);
  //---debugging
  //std::cout<<"mhds= "<<xm<<", xs= "<<mxs<<" XS_max="<<XS_max<<std::endl;
 
  // if(!xs_ang && iflag==1){ goto angular_sampling;}
 
  if(iflag ==1 ){
    gam_ang = gam_sampling(p);
    if(!gam_ang || !xs_ang ){ goto angular_sampling;}
  }
 
  EvtVector4R ppi = p->getDaug(0)->getP4();
  EvtVector4R pcm = p->getDaug(1)->getP4() + ppi;
  EvtVector4R pbst=-1*pcm;
  pbst.set(0,pcm.get(0));
 
  EvtVector4R pi = boostTo(ppi,pbst);
  //--debugging
  //EvtVector4R d2 = boostTo(p->getDaug(1)->getP4(),pbst);
  //std::cout<<"pi,d2"<<pi<<d2<<std::endl;
 
  if(_mode<=5 || _mode ==44 ||_mode ==72 ||_mode == 73||_mode==94||_mode==95 ||_mode ==79 || _mode>=23 &&_mode<=27 || _mode==36|| _mode ==39 || _mode == 80){//ee->two baryon mode, VP, SP
    byn_ang = baryon_sampling(pcm, pi);
    if(!byn_ang) { goto angular_sampling;}
  }else if(_mode ==6||_mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
    msn_ang = meson_sampling(pcm,pi);
    if(!msn_ang){goto angular_sampling;}
  }
 
  if(_nevt ==0 ){
    std::cout<<"The decay chain: "<<std::endl;
    p->printTree();
  }
  _nevt++;
  //--- for debuggint to make root file
  if(mydbg){
  imode = iflag;
  if(iflag==1){
    EvtVector4R vgam = p->getDaug(_ndaugs)->getP4();
    pgam[0]=vgam.get(0);
    pgam[1]=vgam.get(1);
    pgam[2]=vgam.get(2);
    pgam[3]=vgam.get(3);   
    }
 
    ph1[0]=hd1.get(0);
    ph1[1]=hd1.get(1);
    ph1[2]=hd1.get(2);
    ph1[3]=hd1.get(3);
 
    ph2[0]=hd2.get(0);
    ph2[1]=hd2.get(1);
    ph2[2]=hd2.get(2);
    ph2[3]=hd2.get(3);
 
    phds[0]=vhds.get(0);
    phds[1]=vhds.get(1);    
    phds[2]=vhds.get(2);
    phds[3]=vhds.get(3);
    mhds = vhds.mass();
    xs->Fill();
    //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
  }
  //----
  return ;
}

difgamXs()

double EvtConExc::difgamXs

(

EvtParticle *

p

)

Definition at line 973 of file EvtConExc.cc.

                                        {// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
  // where p -> hdarons(0 ~ _ndaugs-1) + gamma
  EvtVector4R p0 = p->getDaug(0)->getP4();
  for(int i=1;i<_ndaugs;i++){
  p0 += p->getDaug(i)->getP4();
  }
  EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
  double mhs = p0.mass();
  double egam = pgam.get(0);
  double sin2theta =  pgam.get(3)/pgam.d3mag();
  sin2theta = 1-sin2theta*sin2theta;
 
  double cms = p->getP4().mass();
  double alpha = 1./137.;
  double pi = 3.1415926;
  double x = 2*egam/cms;
  double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
  //  std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<<  mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
 
  double xs = myxsection->getXsection(mhs);
  double difxs = 2*mhs/cms/cms * wsx *xs;
  differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
  return difxs;
}

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

findMaxXS()

void EvtConExc::findMaxXS

(

EvtParticle *

p

)

Definition at line 938 of file EvtConExc.cc.

                                        {
  //std::cout<<"nmc= "<<nmc<<std::endl;
  gamId = EvtPDL::getId(std::string("gamma"));
  EvtId xvp = EvtPDL::getId(std::string("xvpho"));
  EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
  double totxs = 0;
  maxXS=-1;
  _er1=0;
  Rad2Xs =0;
  int nmc =  50000;
  for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
  gamH=EvtParticleFactory::particleFactory(xvp, p_init);
  int gamdaugs = _ndaugs+1;
 
  EvtId theDaugs[10];
  for(int i=0;i<_ndaugs;i++){
    theDaugs[i] = daugs[i];
  }
  theDaugs[_ndaugs]=gamId;
 
  gamH->makeDaughters(gamdaugs,theDaugs); 
 
  for(int i=0;i<gamdaugs;i++){
    EvtParticle* di=gamH->getDaug(i);
    double xmi=EvtPDL::getMass(di->getId());
    di->setMass(xmi);
   }
 
  double weight = gamH->initializePhaseSpace(  gamdaugs , theDaugs);
  double thexs = difgamXs(gamH); //prepare the photon angular distribution
  if(thexs>maxXS) {maxXS=thexs;}
  }
 
}

Referenced by init().

gam_sampling()

bool EvtConExc::gam_sampling

(

EvtParticle *

p

)

Definition at line 999 of file EvtConExc.cc.

                                          {
  double pm= EvtRandom::Flat(0.,1);
  double xs = difgamXs( p );
  double xsratio = xs/maxXS;
  if(pm<xsratio){return true;}else{return false;}
}

Referenced by decay().

gamHXSection() [1/3]

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

Definition at line 903 of file EvtConExc.cc.

                                                 {// using Gaussian integration subroutine from Cern lib
  std::cout<<" "<<std::endl;
  extern double Rad2difXs(double *x);
  extern double Rad2difXs2(double *x);
  double eps = 0.01;
  double Rad2Xs;
  if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
  if(Rad2Xs==0){
    for(int iter=0;iter<10;iter++){
      eps += 0.01;
      if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
      if(!Rad2Xs)  return Rad2Xs;
    }
  }
  return Rad2Xs;  // the leading second order correction 
}

gamHXSection() [2/3]

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

Definition at line 878 of file EvtConExc.cc.

                                                                   {//El, Eh : the lower and higher energy of radiative photons
  double totxs = 0;
  maxXS=-1;
  _er1=0;
  Rad2Xs =0;
  double xL=2*El/sqrt(s);
  double xH=2*Eh/sqrt(s);
  for(int i=0;i<nmc;i++){
    double rdm = EvtRandom::Flat(0.,1);// std::cout<<"rdm= "<<rdm<<std::endl;
    double x= xL+ (xH-xL)*rdm;
    Rad2Xs += Rad2difXs(s,x);
    _er1 += differ2; //stored when calculate Rad2Xs
    //  std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
  }
  _er1/=nmc;
  _er1*=(xH-xL);
  //std::cout<<"_er1= "<<_er1<<std::endl;
  Rad2Xs/=nmc; // the leading second order correction 
  double thexs= Rad2Xs*(xH-xL);  //xH-xL is the Jacobi factor 
  //std::cout<<"thexs= "<<thexs<<std::endl;
  return thexs;
}

gamHXSection() [3/3]

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

Definition at line 819 of file EvtConExc.cc.

                                                                        {//El, Eh : the lower and higher energy of radiative photons
  //std::cout<<"nmc= "<<nmc<<std::endl;
  gamId = EvtPDL::getId(std::string("gamma"));
  EvtId xvp = EvtPDL::getId(std::string("xvpho"));
  EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
  double totxs = 0;
  maxXS=-1;
  _er1=0;
  Rad2Xs =0;
  for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
  gamH=EvtParticleFactory::particleFactory(xvp, p_init);
  int gamdaugs = _ndaugs+1;
 
  EvtId theDaugs[10];
  for(int i=0;i<_ndaugs;i++){
    theDaugs[i] = daugs[i];
  }
  theDaugs[_ndaugs]=gamId;
 
  gamH->makeDaughters(gamdaugs,theDaugs); 
 
  for(int i=0;i<gamdaugs;i++){
    EvtParticle* di=gamH->getDaug(i);
    double xmi=EvtPDL::getMass(di->getId());
    di->setMass(xmi);
   }
  loop:
  double weight = gamH->initializePhaseSpace(  gamdaugs , theDaugs);
  //-- slection:theta_gamma > 1 degree
  EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
  double pmag = pgam.d3mag();
  double pz = pgam.get(3);
  double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
  double onedegree = 1./180.*3.1415926;
  //if(sintheta < onedegree) {goto loop;}
  if(pmag <El || pmag >Eh) {goto loop;}
  //--------
  // std::cout<<"pmag= "<<pmag<<std::endl;
  
  double thexs = difgamXs(gamH); //prepare the photon angular distribution
  Rad2Xs += Rad2difXs( gamH );
  if(thexs>maxXS) {maxXS=thexs;}
  double s = p_init.mass2();
  double x = 2*pgam.get(0)/sqrt(s);
  double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
  totxs += rad1xs;
  _er1 += differ;
  gamH->deleteDaughters();
  }  //for int_i block
  _er1/=nmc;
 
  Rad2Xs/=nmc; // the leading second order correction 
  totxs/=nmc;  // first order correction XS
 
  // return totxs;  // first order correction XS
   return Rad2Xs;  // the leading second order correction 
}

Referenced by init().

gamHXSection_er()

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

Definition at line 920 of file EvtConExc.cc.

                                                    {// using Gaussian integration subroutine from Cern lib
  std::cout<<"  "<<std::endl;
  extern double Rad2difXs_er(double *x);
  extern double Rad2difXs_er2(double *x);
  double eps = 0.01;
  double Rad2Xs;
  if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
  if(Rad2Xs==0){
    for(int iter=0;iter<10;iter++){
      eps += 0.01;
      if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
      if(!Rad2Xs) return Rad2Xs;;
    }
  }
  return Rad2Xs;  // the leading second order correction 
}

getName()

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

virtual

Implements EvtDecayBase.

Definition at line 145 of file EvtConExc.cc.

                                            {
 
  model_name="ConExc";     
 
}

init()

void EvtConExc::init ( )

virtual

Reimplemented from EvtDecayBase.

Definition at line 158 of file EvtConExc.cc.

                    {
 
  // check that there are 0 arguments
  // checkNArg(1);
  if(getArg(0) == -1){
    radflag=0;mydbg=false;
    _mode   = getArg(0);
    pdgcode = getArg(1);
    pid = EvtPDL::evtIdFromStdHep(pdgcode );
    nson = getNArg()-2;
    std::cout<<"The decay daughters are "<<std::endl;
    for(int i=2;i<getNArg();i++ ){son[i-2]=EvtPDL::evtIdFromStdHep(getArg(i));std::cout<<EvtPDL::name(son[i-2])<<" ";}
    std::cout<<std::endl;
  }else if(getArg(0)==-2){
    radflag=0;mydbg=false;
    _mode   = getArg(0);
    nson = getNArg()-1;
    for(int i=1;i<getNArg();i++ ){son[i-1]=EvtPDL::evtIdFromStdHep(getArg(i));}
  }
  else if(getNArg()==1){ _mode = getArg(0);radflag=0;mydbg=false;} 
  else if(getNArg()==2){_mode = getArg(0); radflag=getArg(1);mydbg=false;}
  else if(getNArg()==3){ _mode = getArg(0); radflag=getArg(1);mydbg=getArg(2);}
  else{std::cout<<"ConExc:umber of parameter should be 1,2 or 3, but you set "<<getNArg()<<std::endl;::abort(); }
  //--- debugging
  //std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<" "<<EvtPDL::name(son[1])<<std::endl;
  
  gamId = EvtPDL::getId(std::string("gamma"));
  init_mode(_mode);
  XS_max=-1;
  //-----debugging to make root file
    if(mydbg){
      myfile = new TFile("mydbg.root","RECREATE");
      xs=new TTree ("xs","momentum of rad. gamma and hadrons");
      xs->Branch("imode",&imode,"imode/I");
      xs->Branch("pgam",pgam,"pgam[4]/D");
      xs->Branch("phds",phds,"phds[4]/D");
      xs->Branch("ph1",ph1,"ph1[4]/D");
      xs->Branch("ph2",ph2,"ph2[4]/D");
      xs->Branch("mhds",&mhds,"mhds/D");
      xs->Branch("mass1",&mass1,"mass1/D");
      xs->Branch("mass2",&mass2,"mass2/D");
    }
    //--- prepare the output information
    EvtId parentId =EvtPDL::getId(std::string("vpho"));
    double parentMass = EvtPDL::getMass(parentId);
    // std::cout<<"parent mass = "<<parentMass<<std::endl;
    EvtVector4R p_init(parentMass,0.0,0.0,0.0);
    EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
    myxsection =new EvtXsection (_mode);
    double mth=0;
    double mx = p->getP4().mass();
    if(_mode ==-1){
      myxsection->setBW(pdgcode);
      for(int i=0;i<nson;i++){mth +=EvtPDL::getMeanMass(son[i]);}
      if(mx<mth){std::cout<<"The CMS energy is lower than the threshold of the final states"<<std::endl;abort();}
    }else if(_mode == -2){
      mth=myxsection->getXlw();
    }else{ mth=myxsection->getXlw();}
    _cms = mx;
    _unit = myxsection->getUnit();
 
    std::cout<<"The specified mode is "<<_mode<<std::endl;
    EvtConExc::getMode = _mode;
    //std::cout<<"getXlw= "<<mth<<std::endl;
 
    double Egamcut = 0.025; //set photon energy cut according to the BESIII detector
    double Esig = 0.0001;  //sigularity engy 
    double x = 2*Egamcut/_cms;
    double s = _cms*_cms;
    double mhscut = sqrt(s*(1-x));
 
    // std::cout<<"mhscut= "<<mhscut<<" _cms= "<<_cms<<" mth= "<<mth<<std::endl;
    //--- calculated with Gaussian integration
    /*
    _xs0 = gamHXSection(mhscut,_cms); 
    _er0 = gamHXSection_er(mhscut,_cms);
    _xs1 = gamHXSection(mth,mhscut);
    double xs1_err = gamHXSection_er(mth,mhscut);
    */
    //-- calculated with MC integration method
    
    double mhdL=myxsection->getXlw();
    double EgamH = (s-mhdL*mhdL)/2/sqrt(s);
    int nmc=100000; 
    _xs0 = gamHXSection(s,Esig,Egamcut,nmc); 
    _er0 = _er1; // stored when calculate _xs0
    std::cout<<"_xs0= "<<_er0<<std::endl;
    _xs1 = gamHXSection(s,Egamcut,EgamH,nmc);
    double xs1_err = _er1;
 
   //--- sigularity point x from 0 to 0.0001GeV
    double xb= 2*Esig/_cms;
    double sig_xs = SoftPhoton_xs(s,xb)*(myxsection->getXsection(mx));
    _xs0 += sig_xs;
 
    //for information output
    double bornXS   = myxsection->getXsection(mx);  // std::cout<<"BornXS= "<<bornXS<<std::endl;
    double bornXS_er=myxsection->getErr(mx);
    double obsvXS   = _xs0 + _xs1; //gamHXSection(mth ,mx);
    double obsvXS_er= _er0 + xs1_err;
    double corr = obsvXS/bornXS; 
    double corr_er =corr*sqrt(bornXS_er*bornXS_er/bornXS/bornXS + obsvXS_er*obsvXS_er/obsvXS/obsvXS);
 
 
    if(bornXS==0){std::cout<<"EvtConExc::init : The Born cross section at this point is zero!"<<std::endl;abort();}
    if(obsvXS==0){std::cout<<"EvtConExc::init : The observed cross section at this point is zero!"<<std::endl;abort();}
 
    //_xs0 += bornXS;//  events with very soft photon (Egam<0.025) are taken as the born process 
    //_er0  = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
 
    std::cout<<""<<std::endl;
    std::cout<<"========= Generation using cross section (Egamcut=0.025 GeV) =============="<<std::endl;
    std::cout<<" sqrt(s)= "<<mx<< " GeV "<<std::endl;
    if(_mode>=0 || _mode ==-2 ){
    std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<"+/-"<<_er0<<" in Unit "<<_unit.c_str()<<std::endl;
    std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"+/-"<<xs1_err<<" in Unit "<<_unit.c_str()<<std::endl;
    std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
    std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<","<<std::endl;
    std::cout<<"which is calcualted with these quantities:"<<std::endl;
    std::cout<<"the observed cross section is "<<obsvXS<<"+/-"<< obsvXS_er<<_unit.c_str()<<std::endl;
    std::cout<<"and the Born cross section (s) is  "<<bornXS<<"+/-"<<bornXS_er<<_unit.c_str()<<std::endl;
    std::cout<<"Within the range "<<myxsection->getXlw()<<" ~"<<myxsection->getXup()<<" GeV, "<<myxsection->getMsg().c_str()<<std::endl;
    std::cout<<"==========================================================================="<<std::endl;
    }else if(_mode==-1){
    std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<" *Br_ee"<<" in Unit "<<_unit.c_str()<<std::endl;
    std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"*Br_ee in Unit "<<_unit.c_str()<<std::endl;
    std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
    std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<std::endl;
    std::cout<<"==========================================================================="<<std::endl;
    }
    std::cout<<std::endl;
    std::cout<<std::endl;
    
    findMaxXS(p);
    //--debugging
    //std::cout<<"maxXS= "<<maxXS<<std::endl;
    
    if(_xs0 == 0 && _xs1==0){std::cout<<"EvtConExc::ini() has zero cross section"<<std::endl;abort();}
 
    init_Br_ee();
    double mthrld = EvtPDL::getMeanMass(daugs[0]); //threshold mass of hadrons
    for(int jj=1;jj<_ndaugs;jj++){
      mthrld += EvtPDL::getMeanMass(daugs[jj]);
    }
 
    for(int ii=0;ii<6;ii++){
      double mR = EvtPDL::getMeanMass(ResId[ii]);
      if(mx<mR || mR<mthrld) continue;
      double narRxs=Ros_xs(mx,BR_ee[ii],ResId[ii]);
      std::cout<<"The corss section for gamma "<<EvtPDL::name(ResId[ii]).c_str()<<" is: "<< narRxs<<"*Br (nb)."<<std::endl;
    }
    std::cout<<"==========================================================================="<<std::endl;
    std::cout<<std::endl;
    std::cout<<std::endl;
 
    //--- for debugging
    if(mydbg){
      int nd = myxsection->getYY().size();
      double  xx[10000],yy[10000],xer[10000],yer[10000];
      for(int i=0;i<nd;i++){
    xx[i]=myxsection->getXX()[i];
    yy[i]=myxsection->getYY()[i];
    yer[i]=myxsection->getEr()[i];
    xer[i]=0.;
    //std::cout<<"yy "<<yy[i]<<std::endl;
      }
      myth=new TH1F("myth","Exp.data",200,xx[0],xx[nd]);
      for(int i=0;i<nd;i++){
    myth->Fill(xx[i],yy[i]);
      }
      myth->SetError(yer);
      myth->Write();
    } //if(mydbg)_block
    //-------------------------
 
}

init_Br_ee()

void EvtConExc::init_Br_ee

(

)

Definition at line 1141 of file EvtConExc.cc.

                          {
  double psipp_ee=9.6E-06;
  double psip_ee =7.73E-03;
  double jsi_ee  =5.94*0.01;
  double phi_ee  =2.954E-04;
  double omega_ee=7.28E-05;
  double rho_ee = 4.72E-05;
  EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
  EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
  EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
  EvtId phiId=EvtPDL::getId(std::string("phi"));
  EvtId omegaId=EvtPDL::getId(std::string("omega"));
  EvtId rhoId=EvtPDL::getId(std::string("rho0"));
  BR_ee.clear();
  ResId.clear();
 
  BR_ee.push_back(rho_ee);  
  BR_ee.push_back(omega_ee);
  BR_ee.push_back(phi_ee);
  BR_ee.push_back(jsi_ee);
  BR_ee.push_back(psip_ee);
  BR_ee.push_back(psipp_ee);
 
  ResId.push_back(rhoId);
  ResId.push_back(omegaId);
  ResId.push_back(phiId);
  ResId.push_back(jsiId);
  ResId.push_back(pspId);
  ResId.push_back(psppId);
  
}

Referenced by init().

init_mode()

void EvtConExc::init_mode

(

int

mode

)

Definition at line 337 of file EvtConExc.cc.

                                 {
  if(mode ==0 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("p+"));
    daugs[1]=EvtPDL::getId(std::string("anti-p-"));
  }else  if(mode ==1 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("n0"));
    daugs[1]=EvtPDL::getId(std::string("anti-n0"));
  }else  if(mode ==2 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("Lambda0"));
    daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
  }else  if(mode ==3 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("Sigma0"));
    daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
  }else  if(mode ==4 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("Lambda0"));
    daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
  }else  if(mode ==5 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("Sigma0"));
    daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
  }  else  if(mode ==6 ){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
  }else  if(mode ==7 ){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
  }else  if(mode ==8 ){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
  }else  if(mode ==9 ){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K_S0"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else  if(mode ==10 ){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K_S0"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
  }else  if(mode ==11 ){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("eta"));
  }else  if(mode ==12 ){
    _ndaugs =4;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
  }else  if(mode ==13 ){
    _ndaugs =4;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
    daugs[3]=EvtPDL::getId(std::string("pi0"));
  }else  if(mode ==14 ){
    _ndaugs =4;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
  }else  if(mode ==15 ){
    _ndaugs =4;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
    daugs[3]=EvtPDL::getId(std::string("pi0"));
  }else  if(mode ==16 ){
    _ndaugs =4;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("K+"));
    daugs[3]=EvtPDL::getId(std::string("K-"));
  }else  if(mode ==17 ){ 
    _ndaugs =5;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
    daugs[4]=EvtPDL::getId(std::string("pi0"));
  }else  if(mode ==18 ){
    _ndaugs =5;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
    daugs[4]=EvtPDL::getId(std::string("eta"));
  }else  if(mode ==19 ){
    _ndaugs =5;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
    daugs[4]=EvtPDL::getId(std::string("pi0"));
  }else  if(mode ==20 ){
    _ndaugs =5;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
    daugs[4]=EvtPDL::getId(std::string("eta"));
  }else  if(mode ==21 ){
    _ndaugs =6;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
    daugs[4]=EvtPDL::getId(std::string("pi+"));
    daugs[5]=EvtPDL::getId(std::string("pi-"));
  }else  if(mode ==22 ){
    _ndaugs =6;
    daugs[0]=EvtPDL::getId(std::string("pi+"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
    daugs[3]=EvtPDL::getId(std::string("pi-"));
    daugs[4]=EvtPDL::getId(std::string("pi0"));
    daugs[5]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 23){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("eta"));
    daugs[1]=EvtPDL::getId(std::string("phi"));
  }else if(mode == 24){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("phi"));
    daugs[1]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 25){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K*-"));
  }else if(mode == 26){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("K-"));
    daugs[1]=EvtPDL::getId(std::string("K*+"));
  }else if(mode == 27){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("K_S0"));
    daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
  }else if(mode == 28){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K*0"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 29){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K*0"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
  }else if(mode == 30){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K*+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 31){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K*-"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 32){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K_2*0"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 33){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K_2*0"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
  }else if(mode == 34){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
    daugs[2]=EvtPDL::getId(std::string("rho0"));
  }else if(mode == 35){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("phi"));
    daugs[1]=EvtPDL::getId(std::string("pi-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
  }else if(mode == 36){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("phi"));
    daugs[1]=EvtPDL::getId(std::string("f_0"));
  }else if(mode == 37){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("eta"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 38){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("omega"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 39){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("omega"));
    daugs[1]=EvtPDL::getId(std::string("f_0"));
  }else if(mode == 40){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("eta'"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 41){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("f_1'"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 42){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("omega'"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("K-"));
  }else if(mode == 43){
    _ndaugs =4;
    daugs[0]=EvtPDL::getId(std::string("omega'"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
    daugs[3]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 44){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("Sigma-"));
    daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
  }else if(mode == 45){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("K+"));
    daugs[1]=EvtPDL::getId(std::string("K-"));
  }else if(mode == 46){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("K_S0"));
    daugs[1]=EvtPDL::getId(std::string("K_L0"));
  }else if(mode == 50){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("p+"));
    daugs[1]=EvtPDL::getId(std::string("anti-p-"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 51){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("p+"));
    daugs[1]=EvtPDL::getId(std::string("anti-p-"));
    daugs[2]=EvtPDL::getId(std::string("eta"));
  }else if(mode == 68){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D0"));
    daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
 
  }else if(mode == 69){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("anti-D0"));
    daugs[1]=EvtPDL::getId(std::string("D*0"));
 
  }else if(mode == 70){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D0"));
    daugs[1]=EvtPDL::getId(std::string("anti-D0"));
 
}else if(mode == 71){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D+"));
    daugs[1]=EvtPDL::getId(std::string("D-"));
  }else if(mode == 72){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D+"));
    daugs[1]=EvtPDL::getId(std::string("D*-"));
 
  }else if(mode == 73){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D-"));
    daugs[1]=EvtPDL::getId(std::string("D*+"));
 
  }else if(mode == 74){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D*+"));
    daugs[1]=EvtPDL::getId(std::string("D*-"));
 
  }else if(mode == 75){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("D0"));
    daugs[1]=EvtPDL::getId(std::string("D-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
  }else if(mode == 76){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("anti-D0"));
    daugs[1]=EvtPDL::getId(std::string("D+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 77){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("D0"));
    daugs[1]=EvtPDL::getId(std::string("D*-"));
    daugs[2]=EvtPDL::getId(std::string("pi+"));
  }else if(mode == 78){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("anti-D0"));
    daugs[1]=EvtPDL::getId(std::string("D*+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 79){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
    daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
  }else if(mode == 80){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("eta"));
    daugs[1]=EvtPDL::getId(std::string("J/psi"));
  }else if(mode == 81){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("h_c"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 82){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("h_c"));
    daugs[1]=EvtPDL::getId(std::string("pi0"));
    daugs[2]=EvtPDL::getId(std::string("pi0"));
  }else if(mode == 83){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("J/psi"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("K-"));
  }else if(mode == 84){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("J/psi"));
    daugs[1]=EvtPDL::getId(std::string("K_S0"));
    daugs[2]=EvtPDL::getId(std::string("K_S0"));
  }else if(mode == 90){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("J/psi"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 91){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
    daugs[1]=EvtPDL::getId(std::string("pi+"));
    daugs[2]=EvtPDL::getId(std::string("pi-"));
  }else if(mode == 92){
    _ndaugs =3;
    daugs[0]=EvtPDL::getId(std::string("J/psi"));
    daugs[1]=EvtPDL::getId(std::string("K+"));
    daugs[2]=EvtPDL::getId(std::string("K-"));
  }else if(mode == 93){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D_s+"));
    daugs[1]=EvtPDL::getId(std::string("D_s-"));
  }else if(mode == 94){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D_s*+"));
    daugs[1]=EvtPDL::getId(std::string("D_s-"));
  }else if(mode == 95){
    _ndaugs =2;
    daugs[0]=EvtPDL::getId(std::string("D_s*-"));
    daugs[1]=EvtPDL::getId(std::string("D_s+"));
  }else if(mode == -1){
    _ndaugs = nson;
    for(int i=0;i<nson;i++){ daugs[i]=son[i];  }
    std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
  }else if(mode == -2){
    _ndaugs = nson;
    for(int i=0;i<nson;i++){  daugs[i]=son[i];  }
  }else {
    std::cout<<"Bad mode_index number, please refer to the manual."<<std::endl;
    ::abort();
  }
  // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
 
}

Referenced by init().

initProbMax()

void EvtConExc::initProbMax ( )

virtual

Reimplemented from EvtDecayBase.

Definition at line 711 of file EvtConExc.cc.

                           {
 
  noProbMax();
 
}

Li2()

double EvtConExc::Li2

(

double

x

)

Definition at line 1265 of file EvtConExc.cc.

                             {
  double li2= -x +x*x/4. - x*x*x/9;
  return li2;
}

Referenced by SoftPhoton_xs().

meson_sampling()

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

Definition at line 1025 of file EvtConExc.cc.

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

Referenced by decay().

Rad1()

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

Definition at line 1036 of file EvtConExc.cc.

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

Referenced by Rad1difXs().

Rad1difXs()

double EvtConExc::Rad1difXs

(

EvtParticle *

p

)

Definition at line 1113 of file EvtConExc.cc.

                                         {// // first order xsection
  // where p -> hdarons(0 ~ _ndaugs-1) + gamma
  double summass  = p->getDaug(0)->getP4().mass();
  for(int i=1;i<_ndaugs;i++){
    summass += p->getDaug(i)->getP4().mass();
  }
 
  double cms = p->getP4().mass();
  double mrg = cms - summass;
  double pm= EvtRandom::Flat(0.,1);
  double mhs = pm*mrg + summass;
 
  double s = cms * cms;
  double x = 1 - mhs*mhs/s;
  double wsx = Rad1(s,x);
 
  //  std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<<  mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
 
  double xs = myxsection->getXsection(mhs);
  if(xs>XS_max){XS_max = xs;}
  double difxs = 2*mhs/cms/cms * wsx *xs;
 
  differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
  differ *= mrg; //Jacobi factor for variable m_hds
  difxs *= mrg;
  return difxs;
}

Referenced by gamHXSection().

Rad2()

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

Definition at line 1048 of file EvtConExc.cc.

                                        {
  //radiator correction upto second order, see Int. J. of Moder.Phys. A
  // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
  double alpha = 1./137.;
  double pi=3.1415926;
  double me = 0.5*0.001; //e mass
  double xi2 = 1.64493407;
  double xi3=1.2020569;
  double sme = sqrt(s)/me;
  double L = 2*log(sme);
  double beta = 2*alpha/pi*(L-1);
  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.;
  double ap = alpha/pi;
  double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
  double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
  double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
  wsx = wsx + beta*beta/8. *wsx2;
  return wsx;
}

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

Rad2difXs() [1/2]

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

Definition at line 1101 of file EvtConExc.cc.

                                              {// leading second order xsection
  double wsx = Rad2(s,x);
  double mhs = sqrt(s*(1-x));
  double xs = myxsection->getXsection(mhs);
  //std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
  if(xs>XS_max){XS_max = xs;}
  double difxs =  wsx *xs;
  differ2 =  wsx *(myxsection->getErr(mhs));
  return difxs;
}

Rad2difXs() [2/2]

double EvtConExc::Rad2difXs

(

EvtParticle *

p

)

Definition at line 1070 of file EvtConExc.cc.

                                         {// leading second order xsection
  // where p -> hdarons(0 ~ _ndaugs-1) + gamma
  double summass = p->getDaug(0)->getP4().mass();
  EvtVector4R v4p=p->getDaug(0)->getP4();
  for(int i=1;i<_ndaugs;i++){
    summass += p->getDaug(i)->getP4().mass();
    v4p +=  p->getDaug(i)->getP4();
  }
 
  double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
  double cms = p->getP4().mass();
  double mrg = cms - summass;
  double pm= EvtRandom::Flat(0.,1);
  double mhs = pm*mrg + summass;
 
  double s = cms * cms;
  double x = 2*Egam/cms;
  //double mhs = sqrt( s*(1-x) );
  double wsx = Rad2(s,x);
 
  //  std::cout<<"random : "<<pm<<std::endl;
 
  double xs = myxsection->getXsection(mhs);
  if(xs>XS_max){XS_max = xs;}
  double difxs = 2*mhs/cms/cms * wsx *xs;
  differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
  differ *= mrg; //Jacobi factor for variable m_hds
  difxs *= mrg;
  return difxs;
}

Referenced by gamHXSection().

Ros_xs()

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

Definition at line 1173 of file EvtConExc.cc.

                                                        {// in unit nb
  double pi=3.1415926;
  double s= mx*mx;
  double width = EvtPDL::getWidth(pid);
  double mass = EvtPDL::getMeanMass(pid);
  double xv = 1-mass*mass/s;
  double sigma = 12*pi*pi*bree*width/mass/s;
  sigma *= Rad2(s, xv);
  double nbar = 4E05;  // ! GeV-2 = 4*10^5 nbar
    return sigma*nbar;
}

Referenced by init().

SoftPhoton_xs()

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

Definition at line 1242 of file EvtConExc.cc.

                                                 {
  double alpha = 1./137.;
  double pi=3.1415926;
  double me = 0.5*0.001; //e mass
  double xi2 = 1.64493407;
  double xi3=1.2020569;
  double sme = sqrt(s)/me;
  double L = 2*log(sme);
  double beta = 2*alpha/pi*(L-1);
  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.;
  double ap = alpha/pi;
  double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
 
  double beta2 = beta*beta;
  double b2 = b*b;
 
  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 - 
     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) + 
         16*pow(beta,2)*Li2(b))/32. ;
  return xs;
 
}

Referenced by init().

xs_sampling()

bool EvtConExc::xs_sampling

(

double

xs

)

Definition at line 1006 of file EvtConExc.cc.

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

Referenced by decay().

Member Data Documentation

_cms

double EvtConExc::_cms

static

Definition at line 125 of file EvtConExc.hh.

Referenced by init(), Rad2difXs(), Rad2difXs2(), Rad2difXs_er(), and Rad2difXs_er2().

getMode

int EvtConExc::getMode

static

Definition at line 128 of file EvtConExc.hh.

Referenced by init().

myxsection

EvtXsection * EvtConExc::myxsection

static

Definition at line 124 of file EvtConExc.hh.

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

XS_max

double EvtConExc::XS_max

static

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

• EvtConExc.hh
• EvtConExc.cc