00001 //-------------------------------------------------------------------------- 00002 // pingrg-2014-10-16 00003 // Model REXC : used to decay vgam to the final state as in ConExc model 00004 // Name: REXC: R-scan exclusive decay model 00005 // Decay cards: 00006 // Decay vgam 00007 // 1 REXC; 00008 // Enddecay 00009 //------------------------------------------------------------------------ 00010 00011 #include "EvtGenBase/EvtPatches.hh" 00012 #include <stdlib.h> 00013 #include "EvtGenBase/EvtParticle.hh" 00014 #include "EvtGenBase/EvtGenKine.hh" 00015 #include "EvtGenBase/EvtPDL.hh" 00016 #include "EvtGenBase/EvtReport.hh" 00017 #include "EvtGenModels/EvtRexc.hh" 00018 #include "EvtGenBase/EvtRandom.hh" 00019 #include <string> 00020 00021 EvtRexc::~EvtRexc() {} 00022 00023 void EvtRexc::getName(std::string& model_name){ 00024 00025 model_name="REXC"; //R-scan exclusive decay model 00026 00027 } 00028 00029 EvtDecayBase* EvtRexc::clone(){ 00030 00031 return new EvtRexc; 00032 00033 } 00034 00035 00036 void EvtRexc::init(){ 00037 00038 // check that there are 0 arguments 00039 checkNArg(0); 00040 myconexc = new EvtConExc(); 00041 } 00042 00043 void EvtRexc::initProbMax(){ 00044 00045 noProbMax(); 00046 00047 } 00048 00049 void EvtRexc::decay( EvtParticle *p ){ 00050 double mhds = p->mass(); 00051 int mymode = EvtConExc::conexcmode; 00052 myconexc->init_mode(mymode); 00053 //std::cout<<"EvtRexc:: A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging 00054 int _ndaugs = myconexc->getNdaugs(); 00055 EvtId *daugs=myconexc->getDaugId(); 00056 //debugging 00057 /* 00058 std::cout<<"Ndaugs= "<<_ndaugs<<std::endl; 00059 for(int i=0;i<_ndaugs;i++){ 00060 std::cout<<i<<" "<<EvtPDL::getStdHep(daugs[i])<<std::endl; 00061 } 00062 */ 00063 Loop: 00064 double totmass=0; 00065 p->makeDaughters(_ndaugs,daugs); 00066 for(int i=0;i< _ndaugs;i++){ 00067 EvtParticle* di=p->getDaug(i); 00068 double xmi=EvtPDL::getMass(di->getId()); 00069 di->setMass(xmi); 00070 totmass += xmi; 00071 } 00072 if(totmass > p->mass()) goto Loop; 00073 00074 double weight = p->initializePhaseSpace( _ndaugs , daugs); 00075 // std::cout<<"weight= "<<weight<<std::endl; 00076 if( (2.5< mhds && mhds<3.092 || mhds>3.1012) && !angularSampling(p)) goto Loop; 00077 return ; 00078 } 00079 00080 bool EvtRexc::angularSampling(EvtParticle* par){ 00081 bool tagp,tagk; 00082 tagk=0; 00083 tagp=0; 00084 int nds = par->getNDaug(); 00085 for(int i=0;i<par->getNDaug();i++){ 00086 EvtId di=par->getDaug(i)->getId(); 00087 EvtVector4R p4i = par->getDaug(i)->getP4Lab(); 00088 int pdgcode =EvtPDL::getStdHep( di ); 00089 double alpha=1; 00090 double angmax = 1+alpha; 00091 double costheta = cos(p4i.theta()); 00092 double ang=1+alpha*costheta*costheta; 00093 double xratio = ang/angmax; 00094 double xi=EvtRandom::Flat(0.,1); 00095 //std::cout<<"pdgcode "<<pdgcode<<std::endl; 00096 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl; 00097 if(xi>xratio) return false; 00098 }//loop over duaghters 00099 return true; 00100 } 00101 00102 double EvtRexc::baryonAng(double mx){ 00103 double mp=0.938; 00104 double u = 0.938/mx; 00105 u = u*u; 00106 double u2 = (1+u)*(1+u); 00107 double uu = u*(1+6*u); 00108 double alpha = (u2-uu)/(u2+uu); 00109 return alpha; 00110 } 00111