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
IDataInfoSvc *	dataInfoSvc

Detailed Description

Definition at line 63 of file EvtDecay.h.

Constructor & Destructor Documentation

EvtDecay()

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

Definition at line 90 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 = "" );
}

Referenced by EvtDecay().

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 1846 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 375 of file EvtDecay.cxx.

                             {
  // new Added
  EvtParticle::_NextLevelDauNum = 0; //////////add for Model:FromParent
 
  MsgStream log( msgSvc(), name() );
  // log << MSG::INFO << "EvtDecay executing" << endmsg;
  //------------------
  string key = "GEN_EVENT";
  // retrieve event from Transient Store (Storegate)
  SmartDataPtr<McGenEventCol> McEvtColl( eventSvc(), "/Event/Gen" );
 
  m_numberEvent += 1;
  writeFlag = 0;
  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 != StatusCode::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 1105 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( msgSvc(), name() );
  truth.close();
  log << MSG::INFO << "EvtDecay finalized" << endmsg;
  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 153 of file EvtDecay.cxx.

                                {
 
  MsgStream log( msgSvc(), 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" << endmsg;
  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" << endmsg;
    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* dataInfoSvc;
 
  // write decay cards to the root file: pingrg@2009-09-09
  StatusCode status;
  status = service( "DataInfoSvc", dataInfoSvc );
  if ( status.isSuccess() )
  {
    log << MSG::INFO << "got the DataInfoSvc" << endmsg;
    dataInfoSvc->setDecayCard( userDecFileName );
  }
  else
  {
    log << MSG::ERROR << "could not get the DataInfoSvc" << endmsg;
    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"
        << endmsg;
  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

IDataInfoSvc* EvtDecay::dataInfoSvc

Definition at line 72 of file EvtDecay.h.

Referenced by initialize().

tmpInfoSvc

IDataInfoSvc* EvtDecay::tmpInfoSvc

Definition at line 71 of file EvtDecay.h.

---

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

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