BesSim.cc

Go to the documentation of this file.

#include "BesSim.hh"
#include "BesActionInitializer.hh"
#include "BesDetectorConstruction.hh"
#include "BesEventAction.hh"
#include "BesRunAction.hh"
#include "BesSteppingAction.hh"
#include "BesTrackingAction.hh"
#include "GenSim/BesPrimaryGeneratorAction.hh"
#include "MdcSim/BesMdcGeoParameter.hh"
#include "PhySim/BesPhysicsList.hh"
#include "PhySim/BgsPhysicsList.hh" //caogf
#include "SimUtil/ReadBoostRoot.hh"
#include "TruSim/BesSensitiveManager.hh"
// #include "BesVisManager.hh"
// #include "G4VisManager.hh"
#include "FTFP_BERT.hh"
#include "FTFP_BERT_ATL.hh"
#include "FTFP_BERT_HP.hh"
#include "FTFP_BERT_TRV.hh"
#include "FTFQGSP_BERT.hh"
#include "FTF_BIC.hh"
#include "QBBC.hh"
#include "QGSP_BERT.hh"
#include "QGSP_BERT_HP.hh"
#include "QGSP_BIC.hh"
#include "QGSP_BIC_AllHP.hh"
#include "QGSP_BIC_HP.hh"
#include "QGSP_FTFP_BERT.hh"
#include "QGS_BIC.hh"
 
#include "G4Run.hh"
#include "G4RunManager.hh"
 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/NTuple.h"
#include "GaudiKernel/SmartIF.h"
// #include "GaudiKernel/Property.h"
// #include "G4Svc/G4Svc.h"
#include "BesRndmGenSvc/IBesRndmGenSvc.h"
#include "G4Svc/IG4Svc.h"
 
#include "G4Writer/G4GDMLWriter.h"
// caogf for random seed
#include "CLHEP/Random/Random.h"
// sunss for detector version
#include "DetVerSvc/IDetVerSvc.h"
 
// static const AlgFactory<BesSim>    xxBAFactory;
// const IAlgFactory& BesSimFactory = xxBAFactory;
 
/////////////////////////////////////////////////////////////////////////////
DECLARE_COMPONENT( BesSim )
 
BesSim::BesSim( const std::string& name, ISvcLocator* pSvcLocator )
    : Algorithm( name, pSvcLocator ) {
  declareProperty( "Mdc", m_mdc = 2 );
  declareProperty( "Tof", m_tof = 2 );
  declareProperty( "Emc", m_emc = 2 );
  declareProperty( "Muc", m_muc = 2 );
  declareProperty( "PipeSCM", m_pipeSCM = 2 );
  declareProperty( "Field", m_field = 1 );
  declareProperty( "PhysicsList", m_physicsList = 0 );
  declareProperty( "FormatAR", m_formatAR = 0 );
  declareProperty( "Tuning", m_tuning = 0 );
  declareProperty( "TuningFile", m_tuningFile );
  declareProperty( "HitOut", m_hitOut = 0 );
  declareProperty( "TDSFlag", m_TDSFlag = false );
  declareProperty( "AsciiFlag", m_asciiFlag = "00000000" );
  declareProperty( "AsciiFile", m_asciiFile = "boost.dat" );
  declareProperty( "RootFlag", m_rootFlag = false );
  declareProperty( "RootFile", m_rootFile = "Hit.root" );
  declareProperty( "CutValue", newcut = 0.7 ); // unite is mm
  declareProperty( "LogLevel", m_logLevel = 5000 );
  declareProperty( "BesGenAction", m_besGenAction = false );
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesSim::initialize() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "initialize()" << endmsg;
  m_nEvent = 0;
  m_nRun   = 1;
 
  // caogf for random seed
  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;
  }
 
  StatusCode status;
  IG4Svc*    tmpSvc;
  status = service( "G4Svc", tmpSvc );
  if ( status.isSuccess() )
  {
    log << MSG::INFO << "got the G4Svc" << endmsg;
    // m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
    m_G4Svc = tmpSvc;
  }
  else
  {
    log << MSG::ERROR << "could not get the G4Svc" << endmsg;
    return StatusCode::FAILURE;
  }
 
  IDetVerSvc* detVerSvc;
  StatusCode  sc_det = service( "DetVerSvc", detVerSvc );
  if ( sc_det.isFailure() )
  {
    log << MSG::ERROR << "can't retrieve DetVerSvc instance" << endmsg;
    return sc_det;
  }
 
  int phase = detVerSvc->phase();
  log << MSG::INFO << "** ~~~ZZZ~~~ ** : retrieved DetectorStage = " << phase << endmsg;
 
  // phase 1: Endcap TOF Scintillator; phase 2: 4 Scintillator==> 2 MRPC; phase 3: Endcap TOF
  // MRPC
  if ( phase > 2 ) { m_tof = 4; }
 
  // book Ntuple
  if ( m_G4Svc->MdcRootFlag() ) bookMdcRootFile();
 
  if ( m_G4Svc->TofRootFlag() ) bookTofRootFile();
 
  if ( m_G4Svc->EmcRootFlag() ) bookEmcRootFile();
 
  if ( m_G4Svc->MucRootFlag() ) bookMucRootFile();
 
  log << MSG::INFO << "Four sub-detector construction flags: " << m_mdc << " " << m_tof << " "
      << m_emc << " " << m_muc << endmsg;
 
  ReadBoostRoot* readBoost =
      new ReadBoostRoot( m_mdc, m_tof, m_emc, m_muc, m_field, m_formatAR, m_tuning, m_hitOut );
  BesMdcGeoParameter* mdcGeoService;
  mdcGeoService = new BesMdcGeoParameter;
 
  log << MSG::INFO << "Instantiating BESIII Detector" << endmsg;
 
  BesSensitiveManager* sensitiveManager = new BesSensitiveManager;
  sensitiveManager->SetLogLevel( m_G4Svc->LogLevel() );
 
  BesDetectorConstruction* detectorConstruction = new BesDetectorConstruction();
  detectorConstruction->SetPipeSCM( m_pipeSCM );
  m_G4Svc->SetUserInitialization( detectorConstruction );
 
  switch ( m_physicsList )
  {
  case 1: m_G4Svc->SetUserInitialization( new QBBC ); break;
  case 2: m_G4Svc->SetUserInitialization( new QGS_BIC ); break;
  case 3: m_G4Svc->SetUserInitialization( new QGSP_BERT ); break;
  case 4: m_G4Svc->SetUserInitialization( new QGSP_BERT_HP ); break;
  case 5: m_G4Svc->SetUserInitialization( new QGSP_BIC_AllHP ); break;
  case 6: m_G4Svc->SetUserInitialization( new QGSP_BIC ); break;
  case 7: m_G4Svc->SetUserInitialization( new FTF_BIC ); break;
  case 8: m_G4Svc->SetUserInitialization( new FTFP_BERT_ATL ); break;
  case 9: m_G4Svc->SetUserInitialization( new FTFP_BERT ); break;
  case 10: m_G4Svc->SetUserInitialization( new FTFP_BERT_HP ); break;
  case 11: m_G4Svc->SetUserInitialization( new FTFP_BERT_TRV ); break;
  case 12: m_G4Svc->SetUserInitialization( new FTFQGSP_BERT ); break;
  default: m_G4Svc->SetUserInitialization( new FTFP_BERT ); break;
  }
 
  BesActionInitializer* actionInit = new BesActionInitializer();
  actionInit->SetBesGenActionFlag( m_besGenAction );
  actionInit->SetTDSFlag( m_TDSFlag );
  actionInit->SetAsciiFlag( m_asciiFlag );
  actionInit->SetAsciiFile( m_asciiFile );
  actionInit->SetRootFlag( m_rootFlag );
  actionInit->SetRootFile( m_rootFile );
  actionInit->SetTuningFile( m_tuningFile );
 
  m_G4Svc->SetUserInitialization( actionInit );
 
  // m_G4Svc->SetUserAction(new BesPrimaryGeneratorAction);
  // BesRunAction* runAction = new BesRunAction;
  // runAction->SetTDSFlag(m_TDSFlag);
  // runAction->SetAsciiFlag(m_asciiFlag);
  // runAction->SetAsciiFile(m_asciiFile);
  // runAction->SetRootFlag(m_rootFlag);
  // runAction->SetRootFile(m_rootFile);
  // vector<string>::iterator pd;
  // for (pd = m_tuningFile.begin(); pd != m_tuningFile.end(); pd++){
  //   std::cout << "tuning File: " << *pd << std::endl;
  // }
 
  // runAction->SetTuningFile(m_tuningFile);
  // m_G4Svc->SetUserAction(runAction);
 
  // BesEventAction* evtAction = new BesEventAction(runAction);
  // m_G4Svc->SetUserAction(evtAction);
 
  // m_G4Svc->SetUserAction(new BesSteppingAction);
 
  // m_G4Svc->SetUserAction(new BesTrackingAction(runAction) );
 
  // #ifdef G4VIS_USE
  // visualization manager
  //   G4VisManager* visManager = new BesVisManager;
  //   visManager->Initialize();
  // #endif
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
// StatusCode BesSim::beginRun() {
StatusCode BesSim::beginRun() {
  debug() << "In start()" << endmsg;
  static int flag = -1;
  if ( flag < 0 )
  {
    MsgStream log( msgSvc(), name() );
    log << MSG::DEBUG << "\t start initializing G4 " << endmsg;
    // new visMgr,  StartUISession(),  runMgr->Initialize()
    m_G4Svc->G4Init();
 
    // to do: runMgr->RunInitialization()
    m_G4Svc->RunInitialize();
 
    // G4VPhysicalVolume* g4wv =G4TransportationManager::GetTransportationManager()->
    // GetNavigatorForTracking()->GetWorldVolume(); //Apply this for storing TOF GDML
    // G4GDMLWriter
    // g4writerBes("/afs/ihep.ac.cn/bes3/offline/sw/packages/BesGDML/2.8.0/GDMLSchema/gdml.xsd",
    // "Bes.gdml",2); G4GDMLWriter
    // g4writerMdc("/afs/ihep.ac.cn/bes3/offline/sw/packages/BesGDML/2.8.0/GDMLSchema/gdml.xsd",
    // "Mdc.gdml",2); G4GDMLWriter
    // g4writerTof("/afs/ihep.ac.cn/bes3/offline/sw/packages/BesGDML/2.8.0/GDMLSchema/gdml.xsd",
    // "Tof.gdml",2);//Apply this for storing TOF GDML G4GDMLWriter
    // g4writerEmc("/afs/ihep.ac.cn/bes3/offline/sw/packages/BesGDML/2.8.0/GDMLSchema/gdml.xsd",
    // "Emc.gdml",2); G4GDMLWriter
    // g4writerMuc("/afs/ihep.ac.cn/bes3/offline/sw/packages/BesGDML/2.8.0/GDMLSchema/gdml.xsd",
    // "Muc.gdml",2); try
    //{
    // g4writerBes.DumpGeometryInfo(g4wv);
    // g4writerTof.DumpGeometryInfo(g4wv->GetLogicalVolume()->GetDaughter(0)); //Changed from 0
    // to 3//Apply this for storing TOF GDML
 
    //}
    // catch(std::logic_error &lerr)
    //{
    //  std::cout << "Caught an exception: "
    //    << lerr.what () << std::endl;
    //}
 
    // For random seed added by caogf. Note the position of the code, otherwise it is not
    // available.
    CLHEP::HepRandomEngine* engine = p_BesRndmGenSvc->GetEngine( "SIM" );
    HepRandom::setTheEngine( engine );
    HepRandom::showEngineStatus();
  }
  return StatusCode::SUCCESS;
}
 
StatusCode BesSim::bookMdcRootFile() {
  MsgStream      log( msgSvc(), name() );
  NTuplePtr      nt( ntupleSvc(), "FILE801/n1" );
  NTuple::Tuple* tupleMdc;
  if ( nt ) tupleMdc = nt;
  else
  {
    tupleMdc = ntupleSvc()->book( "FILE801/n1", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleMdc( tupleMdc );
  }
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesSim::bookTofRootFile() {
  MsgStream      log( msgSvc(), name() );
  NTuplePtr      nt1( ntupleSvc(), "FILE802/n1" );
  NTuple::Tuple* tupleTof1;
  if ( nt1 ) tupleTof1 = nt1;
  else
  {
    tupleTof1 = ntupleSvc()->book( "FILE802/n1", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleTof1( tupleTof1 );
  }
  NTuplePtr      nt2( ntupleSvc(), "FILE802/n2" );
  NTuple::Tuple* tupleTof2;
  if ( nt2 ) tupleTof2 = nt2;
  else
  {
    tupleTof2 = ntupleSvc()->book( "FILE802/n2", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleTof2( tupleTof2 );
  }
  NTuplePtr      nt3( ntupleSvc(), "FILE802/n3" );
  NTuple::Tuple* tupleTof3;
  if ( nt3 ) tupleTof3 = nt3;
  else
  {
    tupleTof3 = ntupleSvc()->book( "FILE802/n3", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleTof3( tupleTof3 );
  }
  /*
  NTuplePtr nt4(ntupleSvc(), "FILE802/n4");
  NTuple::Tuple* tupleTof4;
  if(nt4) tupleTof4 = nt4;
  else {
    tupleTof4 = ntupleSvc()->book("FILE802/n4",CLID_ColumnWiseTuple,"BesSim");
    m_G4Svc->SetTupleTof4(tupleTof4);
  }
  */
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesSim::bookEmcRootFile() {
  NTuplePtr      nt1( ntupleSvc(), "FILE803/n1" );
  NTuple::Tuple* tupleEmc1;
  if ( nt1 ) tupleEmc1 = nt1;
  else
  {
    tupleEmc1 = ntupleSvc()->book( "FILE803/n1", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleEmc1( tupleEmc1 );
  }
 
  NTuplePtr      nt2( ntupleSvc(), "FILE803/n2" );
  NTuple::Tuple* tupleEmc2;
  if ( nt2 ) tupleEmc2 = nt2;
  else
  {
    tupleEmc2 = ntupleSvc()->book( "FILE803/n2", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleEmc2( tupleEmc2 );
  }
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesSim::bookMucRootFile() {
  MsgStream      log( msgSvc(), name() );
  NTuplePtr      nt( ntupleSvc(), "FILE804/n1" );
  NTuple::Tuple* tupleMuc;
  if ( nt ) tupleMuc = nt;
  else
  {
    tupleMuc = ntupleSvc()->book( "FILE804/n1", CLID_ColumnWiseTuple, "BesSim" );
    m_G4Svc->SetTupleMuc( tupleMuc );
  }
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesSim::execute() {
  if ( !m_beginRun )
  {
    StatusCode sc = beginRun();
    if ( sc.isFailure() )
    {
      error() << "beginRun failed" << endmsg;
      return StatusCode::FAILURE;
    }
    m_beginRun = true;
  }
 
  MsgStream log( msgSvc(), name() );
 
  // log << MSG::INFO << "event ID:"<<m_nEvent<<endmsg;
  if ( m_nEvent % m_logLevel == 0 )
    std::cout << " BesSim Begin of Event: " << m_nEvent << std::endl;
 
  m_G4Svc->SimulateEvents( m_nEvent );
  m_nEvent++;
 
  return StatusCode::SUCCESS;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 
StatusCode BesSim::finalize() {
 
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "finalize()" << endmsg;
  log << MSG::DEBUG << "\t terminating the current G4 run" << endmsg;
  m_G4Svc->RunTerminate();
 
  std::cout << "BesSim::finalize(), total events in this run: " << m_nEvent << std::endl;
 
  if ( m_G4Svc->runMgr() != 0 ) delete m_G4Svc->runMgr();
 
#ifdef G4VIS_USE
  if ( m_G4Svc->visMgr() != 0 ) delete m_G4Svc->visMgr();
#endif
 
  return StatusCode::SUCCESS;
}