EvtDecay Class Reference

#include <EvtDecay.h>

Inheritance diagram for EvtDecay:

Public Member Functions
	EvtDecay (const string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
double	ampsLenu (vector< double > Vexp, std::vector< double > vpars)
double	ampsLbenu (vector< double > Vexp, std::vector< double > vpars)
double	HV (double i, double j, vector< double > Vexp, std::vector< double > vpars)
double	HA (double i, double j, vector< double > Vexp, std::vector< double > vpars)
double	getObsXsection (double mhds, int mode)
double	energySpread (double mu, double sigma)

Public Attributes
IDataInfoSvc *	tmpInfoSvc
DataInfoSvc *	dataInfoSvc

Detailed Description

Definition at line 65 of file EvtDecay.h.

Constructor & Destructor Documentation

EvtDecay()

EvtDecay::EvtDecay	(	const string &	name,
		ISvcLocator *	pSvcLocator )

Definition at line 88 of file EvtDecay.cxx.

                                                              :Algorithm( name, pSvcLocator )
{
 
 
 
  // these can be used to specify alternative locations or filenames
  // for the EvtGen particle and channel definition files.
 
  declareProperty("DecayDecDir", m_DecayDec             = "");
  declareProperty("PdtTableDir", m_PdtTable             = "");
  declareProperty("userDecayTableName", userDecFileName = "NONE");
  declareProperty("DecayTopology", m_DecayTop           = "");    // output decay topology to a file specified by user    
  declareProperty("DecayRecord", m_DecayRec             = "");    // output decay topology of one particle specified by user to a file                   
  declareProperty("ParentParticle", m_ParentPart        = "NONE");// Mothor particle name in pdt.table for only running BesEvtGen
  declareProperty("Multiplicity", m_Ncharge             = false); // output ncharge number of an event
  declareProperty("NutupleFile",m_NtupleFile            = false); // output Ntuple file 
  declareProperty("mDIY",_mDIY= false);                           // mDIY model 
  declareProperty("FDPdata",m_SB3File                   = "");    // Fit Dalitz Plot (FDP) data file name (*.root) 
  declareProperty("FDPHisT",m_SB3HT                     = "");    // histogram title of Dalitz plot in  data file  
  declareProperty("FDPpart",m_FDPparticle               ="");     //to assign the FDP parent particle name (string)
  declareProperty("TruthFile",m_truthFile               =""); 
  declareProperty("TruthPart",m_truthPart               =""); 
  declareProperty("Psi3SopenCharm",m_Psi4040OpenCharm=false);
  declareProperty("Psi2openCharm", m_Psi2openCharm=false);
  declareProperty("SetMthrForConExc",m_SetMthr=0.0);
  declareProperty("statDecays",m_statDecays=false);
  declareProperty("TagLundCharmModel", m_tagLundModel=false);
  m_mystring.clear();
  declareProperty("FileForTrackGen", m_mystring);
  //for polarized charmonium production
  m_polarization.clear();//= diag{1+Pe, 0, 1-Pe} with electron polarization Pe
  declareProperty("polarization", m_polarization);
  declareProperty("eBeamPolarization", m_eBeamPolarization=-1);
  m_cluster0.clear();m_wind0.clear();
  m_cluster1.clear();m_wind1.clear();
  m_cluster2.clear();m_wind2.clear();
  declareProperty("cluster0",m_cluster0);
  declareProperty("cluster1",m_cluster1);
  declareProperty("cluster2",m_cluster2);
  declareProperty("masswin0",m_wind0);
  declareProperty("masswin1",m_wind1);
  declareProperty("masswin2",m_wind2);
  declareProperty("OutputP4",m_outputp4="");
  declareProperty("ReadMeasuredEcms", m_RdMeasuredEcms = false);
  declareProperty("beamEnergySpread", m_beamEnergySpread = 0);
  m_ReadTruth.clear();
  declareProperty("ReadTruthAtCM",_truthAtCM=false);
  declareProperty("ReadTruth",m_ReadTruth); 
  //ReadTruth={{"ParentName"},{"i0","i1","i2"},{"j0","j1","j2","j3"}}, where the first part. is Parent->getDaug(i0)->getDaug(i1)->getDaug(i2),
  //and second particle is Parent ->getDaug(j0)->getDaug(j1)->getDaug(j2)->getDaug(j3);
  declareProperty("RvalueTag",_RvalueTag=false);
  declareProperty("PrintFSFor",m_printfs="");
}

Member Function Documentation

ampsLbenu()

double EvtDecay::ampsLbenu	(	vector< double >	Vexp,
		std::vector< double >	vpars )

ampsLenu()

double EvtDecay::ampsLenu	(	vector< double >	Vexp,
		std::vector< double >	vpars )

energySpread()

double EvtDecay::energySpread	(	double	mu,
		double	sigma )

Definition at line 1616 of file EvtDecay.cxx.

                                                   {
  //mu: mean value in Gaussian
  //sigma: variance in Gaussian 
 rloop:
  int n=12;
  double ri=0;
  for(int i=0;i<n;i++){
    double pm= EvtRandom::Flat(0.,1);
    ri += pm;
  }
  double eta=sqrt(n*12.0)*(ri/12-0.5);
  double xsig= sigma*eta+mu;
  double mx0=EvtConExc::mlow;
  double mx1=EvtConExc::mup;
  if(xsig<mx0 || xsig>mx1) goto rloop;
  return xsig;
}

execute()

StatusCode EvtDecay::execute

(

)

Definition at line 343 of file EvtDecay.cxx.

{ 
  //new Added
  EvtParticle::_NextLevelDauNum=0; //////////add for Model:FromParent
 
  MsgStream log(messageService(), name());
  // log << MSG::INFO << "EvtDecay executing" << endreq;
  //------------------
  string   key = "GEN_EVENT";
  // retrieve event from Transient Store (Storegate)
  SmartDataPtr<McGenEventCol> McEvtColl(eventSvc(), "/Event/Gen");  
 
  m_numberEvent += 1;  
  writeFlag = 0;
  //std::cout<<"EvtNumber= "<<m_numberEvent<<std::endl;
  if ( McEvtColl == 0 ) 
  {
    //std::cout<<"McEvtColl == 0 INFO: Using the callBesEvtGen!!"<<std::endl;
    HepMC::GenEvent* evt=new GenEvent();  
    evt->set_event_number(m_numberEvent);
 
    //read Ecms from the database
    SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
    int runNo=eventHeader->runNumber();
    int event=eventHeader->eventNumber();
    bool newRunFlag=0;
    if(runNo != 0 && runNo != m_runNo){m_runNo=runNo;newRunFlag = true;}else{newRunFlag=false;}
    if(m_RdMeasuredEcms&& newRunFlag){// using cms energy of beam read from database
      runNo = std::abs(runNo);
      ReadME theME(runNo);
      if(theME.isRunNoValid()){
    dbEcms=theME.getEcms();
    std::cout<<"Read Ecms= "<<dbEcms<<std::endl;
    if(dbEcms!=0){
      std::cout<<"INFO: Read the MeasuredEcms"<<std::endl;
    }
    else{
      std::cout<<"ERROR: Can not read the MeasuredEcms, try to turn off the ReadEcmsDB"<<std::endl;
      return StatusCode::FAILURE;
    }
      }
    }
    //end of read Ecms 
 
    callBesEvtGen( evt );
    if(!(m_DecayTop==""))  evt->print(outfile);  
 
    //create a Transient store 
    McGenEventCol *mcColl = new McGenEventCol;
    McGenEvent* mcEvent = new McGenEvent(evt);
    mcColl->push_back(mcEvent);
    StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
    if(sc != SUCCESS)  return StatusCode::FAILURE;    
 
  } else  //  (McEvtColl != 0) 
 {
 
  //std::cout<<"McEvtColl == 1 INFO: Using the callEvtGen!!"<<std::endl;
  McGenEventCol::iterator mcItr;
  for( mcItr = McEvtColl->begin(); mcItr != McEvtColl->end(); mcItr++ )   
  {
     HepMC::GenEvent* hEvt = (*mcItr)->getGenEvt();
    // MeVToGeV( hEvt );
 
     callEvtGen( hEvt );
 
     if(!(m_DecayTop==""))  hEvt->print(outfile); 
    // GeVToMeV( hEvt );
     //   if(!(m_DecayRec=="")) outfile2<<std::endl;
  if(!(m_DecayRec=="")) std::cout<<std::endl;
  };
 }
  if( m_NtupleFile ){ m_tuple->write();}
 
  if(_mDIY){
    std::cout<<"helangs ";
    for(int i=0;i<m_vangs.size();i++) std::cout<<m_vangs[i]<<" ";
    std::cout<<std::endl;
   }
  return StatusCode::SUCCESS; 
}

finalize()

StatusCode EvtDecay::finalize

(

)

Definition at line 1008 of file EvtDecay.cxx.

{
 
  if(EvtCalHelAmp::nevt>0){
    double H2=EvtCalHelAmp::_H2;
    double H1=EvtCalHelAmp::_H1;
    double H0=EvtCalHelAmp::_H0;
    double H2err=EvtCalHelAmp::_H2err;
    double H1err=EvtCalHelAmp::_H1err;
    double H0err=EvtCalHelAmp::_H0err;
    int nd = EvtCalHelAmp::nevt;
    std::cout<<"Calculated helicity amplitude sqaured are:"<<std::endl;
    std::cout<<"|H_{2}|^2=|H_{-2}|^2= "<<H2/nd<<" +/- "<<sqrt(H2err/nd)<<endl;
    std::cout<<"|H_{1}|^2=|H_{-1}|^2= "<<H1/nd<<" +/- "<<sqrt(H1err/nd)<<endl;
    //std::cout<<"|H_{0}|^2= "<<H0/nd           <<" +/- "<<sqrt(H0err/nd)<<endl;
  }
  MsgStream log(messageService(), name());
  truth.close();
  log << MSG::INFO << "EvtDecay finalized" << endreq;
   fstream Forfile;
   Forfile.open("fort.16",ios::in);
   char delcmd[300];
   strcpy(delcmd,"rm -f ");
   strcat(delcmd,"fort.16");
   // if(Forfile)system(delcmd);
 
   fstream Forfile2;
   Forfile2.open("fort.22",ios::in);
   strcpy(delcmd,"rm -f ");
   strcat(delcmd,"fort.22");
   // if(Forfile2)system(delcmd);
 
   //   To get the branching fraction of the specified mode in Lundcharm Model
   /*    
   EvtLundCharm lundcharmEvt;
   int nevt=lundcharmEvt.getTotalEvt();
   std::cout<<"The total number of the specified mode is  "<<nevt<<std::endl;
   */ 
   int totalEvt=0;
   if(!(m_SB3File=="" || m_SB3HT=="")){
   for(int i=0;i<500;i++){totalEvt=totalEvt+br[i];}
   for(int i=0;i<500;i++){
     double bi=1.0*br[i]/totalEvt;
     std::cout<<"Branching fraction of channel "<<i<<" is: "<<bi<<std::endl;
     if(br[i]==0) break;
   }
   }
 
   if(m_statDecays){
     int totalEvt=0;
     for(int i=0;i<=totalChannels;i++){totalEvt=totalEvt+br[i];}
     std::cout<<"==================Statistical first chain decay ==============================="<<std::endl;
     std::cout<<"i-th channel      Events Generated         Branching fraction generated        "<<std::endl;
     for(int i=0;i<=totalChannels;i++){
      std::cout<<"|  "<<i<<"              "<<br[i]<<"                        "<<br[i]*1.0/totalEvt<<std::endl;      
 
     }
     std::cout<<"--------------------------------------------------------------------------------"<<std::endl;
     std::cout<<"  sum:            "<<totalEvt<<"                     "<<"1"<<std::endl;
     std::cout<<"================== End of statistical first chain decay ========================"<<std::endl;
   }
 
   return StatusCode::SUCCESS;
}

getObsXsection()

double EvtDecay::getObsXsection	(	double	mhds,
		int	mode )

HA()

double EvtDecay::HA	(	double	i,
		double	j,
		vector< double >	Vexp,
		std::vector< double >	vpars )

HV()

double EvtDecay::HV	(	double	i,
		double	j,
		vector< double >	Vexp,
		std::vector< double >	vpars )

initialize()

StatusCode EvtDecay::initialize

(

)

Definition at line 143 of file EvtDecay.cxx.

                               {
 
  MsgStream log(messageService(), name());
  //system("cat $BESEVTGENROOT/share/phokhara_9.1.fferr>phokhara_9.1.fferr");
  //system("cat $BESEVTGENROOT/share/phokhara_9.1.ffwarn>phokhara_9.1.ffwarn");
  system("cat $BESEVTGENROOT/share/phokhara_10.0.fferr>phokhara_10.0.fferr");    //modified by luojiashun 2022.2.21
  system("cat $BESEVTGENROOT/share/phokhara_10.0.ffwarn>phokhara_10.0.ffwarn");  //modified by luojiashun 2022.2.21
 
  if(m_truthFile!=""){truth.open(m_truthFile.c_str());}
  log << MSG::INFO << "EvtDecay initialize" << endreq;
  static const bool CREATEIFNOTTHERE(true);
  StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
  if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
  {
    log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
    return RndmStatus;
  }
 
  EvtGlobalSet::SV.clear();
  for(int i=0;i<m_mystring.size();i++){EvtGlobalSet::SV.push_back(m_mystring[i]);}
 
  EvtGlobalSet::iVV.clear();
  EvtGlobalSet::dVV.clear();
  std::vector<std::vector<int> >myivv;
  std::vector<std::vector<double> >mydvv;
 
  EvtConExc::SetMthr=m_SetMthr;
 
  myivv.push_back(m_cluster0);
  myivv.push_back(m_cluster1);
  myivv.push_back(m_cluster2);
 
  mydvv.push_back(m_wind0);
  mydvv.push_back(m_wind1);
  mydvv.push_back(m_wind2);
 
  for(int i=0;i<myivv.size();i++){
    std::vector<int> theivv;
    for(int j=0;j<myivv[i].size();j++){
      theivv.push_back(myivv[i][j]);
    }
    if(theivv.size()) EvtGlobalSet::iVV.push_back(theivv);
  }
 
  for(int i=0;i<mydvv.size();i++){
    std::vector<double> thedvv;
    for(int j=0;j<mydvv[i].size();j++){
      thedvv.push_back(mydvv[i][j]);
    }
    if(thedvv.size()) EvtGlobalSet::dVV.push_back(thedvv);
  }
 
 
  if(m_eBeamPolarization>1) {std::cout<<"The e-beam polaziation should be in 0 to 1"<<std::endl;abort();}
  m_numberEvent=0; 
  first = true; 
  m_ampscalflag = true;
  // Save the random number seeds in the event
  CLHEP::HepRandomEngine* engine  = p_BesRndmGenSvc->GetEngine("EVTGEN");
  const long    s  = engine->getSeed();
  p_BesRndmGenSvc->setGenseed(s+1);
 
  m_seeds.clear();
  m_seeds.push_back(s);
  totalChannels = 0;
 
 // open an interface to EvtGen
 
  if ( m_DecayDec == "") {  //use default DECAY.DEC and pdt.table if not defined by user
   m_DecayDec=getenv("BESEVTGENROOT");
   m_DecayDec +="/share/DECAY.DEC";
   }
 
  if ( m_PdtTable == "") {
   m_PdtTable =getenv("BESEVTGENROOT");
   m_PdtTable +="/share/pdt.table";
   }  
  
  char catcmd[300];  //output user decay cards to log file
  std::cout<<"===================== user decay card ========================"<<std::endl;
  strcpy(catcmd, "cat ");
  strcat(catcmd, userDecFileName.c_str());
  system(catcmd);
  
   IDataInfoSvc *tmpInfoSvc;
   DataInfoSvc* dataInfoSvc;
  
  // write decay cards to the root file: pingrg@2009-09-09
  StatusCode status;
  status = service("DataInfoSvc",tmpInfoSvc);
  if (status.isSuccess()) {
    log << MSG::INFO << "got the DataInfoSvc" << endreq;
    dataInfoSvc=dynamic_cast<DataInfoSvc *>(tmpInfoSvc);
    dataInfoSvc->setDecayCard(userDecFileName);
  } else { 
    log << MSG::ERROR << "could not get the DataInfoSvc" << endreq;
    return StatusCode::FAILURE;
  }
 
  
 
  m_RandomEngine  = new EvtBesRandom(engine);
  m_Gen = new EvtGen( m_DecayDec.c_str(), m_PdtTable.c_str(), m_RandomEngine );
      
  if(userDecFileName =="NONE")  log << MSG::INFO << "EvtDecay User did not define his Decay table EvtGen will use standart" << endreq;   
  if(!(userDecFileName =="NONE")) m_Gen->readUDecay(userDecFileName.c_str());
 
 if(!(m_DecayTop==" ")) {
   outfile.open(m_DecayTop.c_str()); 
  }
 
 //for output Ntuple file, pingrg-2009-09-07
 
 
  if(m_NtupleFile) {
       NTuplePtr nt1(ntupleSvc(),"MYROOTFILE/Trk");
       if(nt1) {m_tuple=nt1;}
       else {
          m_tuple = ntupleSvc()->book ("MYROOTFILE/Trk", CLID_ColumnWiseTuple, "Generator-trk-Ntuple");
          if(m_tuple){
          status = m_tuple->addItem ("TotNumTrk", TotNumTrk, 0,100);
          status = m_tuple->addIndexedItem ("Trk_index", TotNumTrk, m_Trk_index);    
          status = m_tuple->addIndexedItem ("m_px_trk",  TotNumTrk, m_px_trk); 
          status = m_tuple->addIndexedItem ("m_py_trk",  TotNumTrk, m_py_trk); 
          status = m_tuple->addIndexedItem ("m_pz_trk",  TotNumTrk, m_pz_trk); 
          status = m_tuple->addIndexedItem ("m_en_trk",  TotNumTrk, m_en_trk); 
          status = m_tuple->addIndexedItem ("FST",       TotNumTrk, m_fst); 
 
          status = m_tuple->addItem ("nchr",     m_nchr); 
          status = m_tuple->addItem ("nchr_e",   m_nchr_e); 
          status = m_tuple->addItem ("nchr_mu",  m_nchr_mu); 
          status = m_tuple->addItem ("nchr_pi",  m_nchr_pi); 
          status = m_tuple->addItem ("nchr_k",   m_nchr_k); 
          status = m_tuple->addItem ("nchr_p",   m_nchr_p); 
          status = m_tuple->addItem ("N_gamma",  m_gamma);
          status = m_tuple->addItem ("N_gammaFSR",  m_gammaFSR); 
          status = m_tuple->addItem ("Flag1",  m_flag1); 
      } else {   
               log << MSG::ERROR << "    Cannot book N-tuple:"<< long(m_tuple) << endmsg;  
               return StatusCode::FAILURE;
      }    
       }
 
       NTuplePtr nt2(ntupleSvc(),"MYROOTFILE/mass");
       if(nt2) {mass_tuple=nt2;}
       else {
          mass_tuple = ntupleSvc()->book ("MYROOTFILE/mass", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
          if(mass_tuple){
          status = mass_tuple->addItem ("m12", m_m12);
          status = mass_tuple->addItem ("m13", m_m13);
          status = mass_tuple->addItem ("m23", m_m23);
          status = mass_tuple->addItem ("m1", m_m1);
          status = mass_tuple->addItem ("m2", m_m2);
          status = mass_tuple->addItem ("m3", m_m3);
          status = mass_tuple->addItem ("costheta1", m_cos1);
          status = mass_tuple->addItem ("costheta2", m_cos2);
          status = mass_tuple->addItem ("costheta3", m_cos3);
          status = mass_tuple->addItem ("ichannel", m_ich);
      } else {   
               log << MSG::ERROR << "    Cannot book N-tuple:"<< long(m_tuple) << endmsg;  
               return StatusCode::FAILURE;
      }    
       }
 
       NTuplePtr nt3(ntupleSvc(),"MYROOTFILE/massGen");
       if(nt3) {massgen_tuple=nt3;}
       else {
          massgen_tuple = ntupleSvc()->book ("MYROOTFILE/massGen", CLID_ColumnWiseTuple, "Generator-mass-Ntuple");
          if(massgen_tuple){
          status = massgen_tuple->addItem ("m12", _m12);
          status = massgen_tuple->addItem ("m13", _m13);
          status = massgen_tuple->addItem ("m23", _m23);
          status = massgen_tuple->addItem ("m1", _m1);
          status = massgen_tuple->addItem ("m2", _m2);
          status = massgen_tuple->addItem ("m3", _m3);
          status = massgen_tuple->addItem ("costheta1", _cos1);
          status = massgen_tuple->addItem ("costheta2", _cos2);
          status = massgen_tuple->addItem ("costheta3", _cos3);
          status = massgen_tuple->addItem ("ichannel", _ich);
      } else {   
               log << MSG::ERROR << "    Cannot book N-tuple:"<< long(m_tuple) << endmsg;  
               return StatusCode::FAILURE;
      }    
       }
 
 
  }
  for(int i=0;i<500;i++){br[i]=0;}
  if(!(m_SB3File=="" && m_SB3HT=="")) {
    const char *filename, *histitle;
    filename=(char*)m_SB3File.c_str();
    histitle=(char*)m_SB3HT.c_str();
    SuperBody3decay.setFile(filename);
    SuperBody3decay.setHTitle(histitle);
    SuperBody3decay.init();
  }
 
  return StatusCode::SUCCESS;
}

Member Data Documentation

dataInfoSvc

DataInfoSvc* EvtDecay::dataInfoSvc

Definition at line 75 of file EvtDecay.h.

Referenced by initialize().

tmpInfoSvc

IDataInfoSvc* EvtDecay::tmpInfoSvc

Definition at line 74 of file EvtDecay.h.

Referenced by initialize().

---

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

• EvtDecay.h
• EvtDecay.cxx
• UsermDIY.cc