In This Package:

Sim15.cc

Go to the documentation of this file.

#include "Sim15.h"
#include "Stage/HeaderStageData.h"
#include "Stage/IStageDataManager.h"

// for GenHeader
#include "Event/GenHeader.h"

// for GenEvent
#include "HepMC/GenEvent.h"

// detector element
#include "DetDesc/IGeometryInfo.h"
#include "DetDesc/DetectorElement.h"
#include "DetDesc/SolidBase.h"

// for SimHit
#include "Event/SimHeader.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/SimPmtHit.h"
#include "Event/SimHit.h"

// for ElecEvent
#include "Event/ElecHeader.h"

// SimTrig
#include "Event/SimTrigHeader.h"

// SimReadout
#include "Event/SimReadoutHeader.h"

// for Readout
#include "Event/ReadoutHeader.h"
#include "Event/ReadoutPmtCrate.h"
#include "Event/ReadoutRpcCrate.h"

// for history info
//#include "Event/SimParticleHistoryHeader.h"
#include "Event/SimParticleHistory.h"
#include "Event/SimUnobservableStatisticsHeader.h"
#include "Event/SimStatistic.h"

// for Detector ID etc.
#include "Conventions/Detectors.h"

#include "GaudiKernel/ITHistSvc.h"

#include <Event/HeaderObject.h>

#include "CLHEP/Units/SystemOfUnits.h"

using namespace std;

Sim15::Sim15(const std::string& name, ISvcLocator* pSvcLocator):
  GaudiAlgorithm(name, pSvcLocator)
  , m_TopStage(0)
  , m_sdmgr(0)
{
  declareProperty("TopStage", m_TopStageName,"Name of top stage");  
  
  declareProperty("TimeRange", m_TimeRange = 365*24*60*60*CLHEP::second, "Time range to simulate, default 1 year.");
  
  declareProperty("Monitor", m_monitor=false,"Generate a monitor.root");
}

Sim15::~Sim15()
{
}

StatusCode Sim15::initialize()
{
  StatusCode sc = this->GaudiAlgorithm::initialize();
  if (sc.isFailure()) return sc;
  
  sc= toolSvc()->retrieveTool("Stage",m_TopStageName,m_TopStage);
  if( sc.isFailure() ) {
    error() << "Error retrieving the tool"<< m_TopStageName << endreq;
  }  
  debug() << "Got top Stage \"" << m_TopStageName << "\" at " << (void*) m_TopStage << endreq;
  // should probably make this configurable one day
  sc = service("StageDataManager",m_sdmgr);
  if (sc.isFailure()) return sc;
  debug() << "Got StageDataManager at " << (void*) m_sdmgr << endreq;

  // init root tree
  if(m_monitor) {
    m_valiTree=new ValidationTree;
    m_valiTree->create();
  }

  if(m_TopStageName=="SingleLoader") {
    // ReadoutHeader
    m_TopStageNum = 5;
  } else if(m_TopStageName=="TrigRead") {
    // SimReadoutHeader
    m_TopStageNum = 4;
  } else if(m_TopStageName=="Electronic") {
    // ElecHeader
    m_TopStageNum = 3;
  } else if (m_TopStageName=="Detector") {
    // SimHeader
    m_TopStageNum = 2;
  } else if (m_TopStageName=="Kinematic") {
    // GenHeader
    m_TopStageNum = 1;
  }

  m_counter = 0;

  sc = service("EventLoopMgr",m_EvtLoopMgr,false);
  if( sc.isFailure() ) {
    error() << "Failed to retrive EventLoopMgr service"<< endreq;
  }
  info()<<"Got EventLoopMgr service"<<endreq;

  info()<<m_TimeRange/CLHEP::second<<" seconds detector simulation is required."<<endreq;

  return StatusCode::SUCCESS;
}

StatusCode Sim15::execute()
{
  IStageData*  pIStageData = 0;
  StatusCode sc = m_TopStage->nextElement(pIStageData); // get a new data
  if (sc.isFailure()) return sc;
  debug() << "consuming stage data at " << pIStageData->time() << endreq;
  StatusCode sc2 = m_sdmgr->consumeData(*pIStageData);
  if (sc2.isFailure()) return sc2;


  m_counter++;
  if( m_counter==1 ) m_StartTime = pIStageData->time();
  m_CurrTime = pIStageData->time();

  info()<<"New element #"<<m_counter<<" at time "<<pIStageData->time()<<endreq;

  // if no monitor.root is required, good to go now.
  if(m_monitor) {
    typedef HeaderStageData<DayaBay::GenHeader>        GenData;
    typedef HeaderStageData<DayaBay::SimHeader>        SimData;
    typedef HeaderStageData<DayaBay::ElecHeader>       ElecData;
    typedef HeaderStageData<DayaBay::SimReadoutHeader> SimReadoutData;
    typedef HeaderStageData<DayaBay::ReadoutHeader>    ReadoutData;

    GenData        *pGenData=0;
    SimData        *pSimData=0;
    ElecData       *pElecData=0;
    SimReadoutData *pSimReadoutData=0;
    ReadoutData    *pReadoutData=0;
    
    const DayaBay::GenHeader        *pGenHeader=0;
    const DayaBay::SimHeader        *pSimHeader=0;
    const DayaBay::ElecHeader       *pElecHeader=0;
    const DayaBay::SimTrigHeader    *pSimTrigHeader=0;
    const DayaBay::SimReadoutHeader *pSimReadoutHeader=0;
    const DayaBay::ReadoutHeader    *pReadoutHeader=0;
    
    const DayaBay::IHeader* firstHeader=0;
    const DayaBay::HeaderObject* pHeader=0;
    
    // get the header object pointer of the top stage
    if(m_TopStageNum==5) {
      // ReadoutHeader
      pReadoutData=dynamic_cast<ReadoutData*>(pIStageData);
      pReadoutHeader=&(pReadoutData->header());
      pHeader=pReadoutHeader;
    } else if(m_TopStageNum==4) { 
      // SimReadoutHeader
      pSimReadoutData=dynamic_cast<SimReadoutData*>(pIStageData);
      pSimReadoutHeader=&(pSimReadoutData->header());
      pHeader=pSimReadoutHeader;
    } else if(m_TopStageNum==3) {
      // ElecHeader
      pElecData=dynamic_cast<ElecData*>(pIStageData);
      pElecHeader=&(pElecData->header());
      pHeader=pElecHeader;
    } else if (m_TopStageNum==2) {
      // SimHeader
      pSimData=dynamic_cast<SimData*>(pIStageData);
      pSimHeader=&(pSimData->header());
      pHeader=pSimHeader;
    } else if (m_TopStageNum==1) {
      // GenHeader
      pGenData=dynamic_cast<GenData*>(pIStageData);
      pGenHeader=&(pGenData->header()); 
      pHeader=pGenHeader;
    }
    
    // monitor, do the rest
    double m_e=5.109989e-1;  // MeV
    double m_alpha=3727.379; // MeV
    
    m_valiTree->reset();
    
    // get the header object pointer of it lower stages
    // only the first one, for low event rate, it is enough
    if(pReadoutHeader!=0) {
      // SimReadoutHeader
      if(pReadoutHeader->inputHeaders().size() != 0) {
        firstHeader=*((pReadoutHeader->inputHeaders()).begin());
        pSimReadoutHeader= dynamic_cast<const DayaBay::SimReadoutHeader*>(firstHeader);
      }
    }
    if(pSimReadoutHeader!=0) {
      debug() << "Number of input SimTrigHeaders:" <<
        pSimReadoutHeader->inputHeaders().size() << endreq;
      // SimTrigHeader
      if(pSimReadoutHeader->inputHeaders().size() != 0) {
        firstHeader=*((pSimReadoutHeader->inputHeaders()).begin());
        pSimTrigHeader= dynamic_cast<const DayaBay::SimTrigHeader*>(firstHeader);
      }
    }
    if(pSimTrigHeader!=0) {
      // ElecHeader
      if(pSimTrigHeader->inputHeaders().size() != 0) {
        firstHeader=*((pSimTrigHeader->inputHeaders()).begin());
        pElecHeader= dynamic_cast<const DayaBay::ElecHeader*>(firstHeader);
      }
    }
    if(pElecHeader!=0) {
      // SimHeader
      if(pElecHeader->inputHeaders().size() != 0) {
        firstHeader=*((pElecHeader->inputHeaders()).begin());
        pSimHeader= dynamic_cast<const DayaBay::SimHeader*>(firstHeader);
      }
    }
    if(pSimHeader!=0) {
      // GenHeader
      if(pSimHeader->inputHeaders().size() != 0) {
        firstHeader=*((pSimHeader->inputHeaders()).begin());
        pGenHeader= dynamic_cast<const DayaBay::GenHeader*>(firstHeader);
      }
    }
    
    
    // event information
    m_valiTree->time=pIStageData->time().GetSeconds();
    
    // check ReadoutHeader information
    if(pReadoutHeader!=0) {
      m_valiTree->execSing=pReadoutHeader->execNumber();
      
      // adc tdc
      const DayaBay::Readout* pReadout = pReadoutHeader->readout();
      DayaBay::Detector det = pReadout->detector();
      
      if(det.detectorId() == DetectorId::kAD1 ||  // AD information
         det.detectorId() == DetectorId::kAD2 ||
         det.detectorId() == DetectorId::kAD3 ||
         det.detectorId() == DetectorId::kAD4) {
        
        // convert it to PMT Crate Readout
        const DayaBay::ReadoutPmtCrate* pmtcrate = 0;
        pmtcrate = dynamic_cast<const DayaBay::ReadoutPmtCrate*>(pReadout);
        
        // channnel map
        const DayaBay::ReadoutPmtCrate::PmtChannelReadouts& channelMap
          = pmtcrate->channelReadout();
        
        DayaBay::ReadoutPmtCrate::PmtChannelReadouts::const_iterator cmci;
        
        for(cmci=channelMap.begin();cmci!=channelMap.end();++cmci) {
          //const DayaBay::FeeChannelId&      channelId = cmci->first;
          const DayaBay::ReadoutPmtChannel& aChannel  = cmci->second;
          
          //modificaton for new readout data model. Please fix me if not right.
          m_valiTree->adc_sum=aChannel.sumAdc();
          //const vector<int>&          tdc=aChannel.tdc();
          //const std::map<int,int>&    adc=aChannel.adc();
          //int tdc_peak,adc_peak;
          
          //for(unsigned long ii=0;ii<tdc.size();++ii) {
          //  tdc[ii];
          //  tdc_peak=tdc[ii];
          //  adc_peak=adc.find(tdc_peak)->second;
          //  m_valiTree->adc_sum+=adc_peak;
          //}
          //verbose()<<channelId<<" tdc "<<tdc_peak<<" adc "<<adc_peak<<endreq;
        }
      }  // end of AD crate
    } // end of pReadoutHeader
    
    // check SimReadoutHeader information
    if(pSimReadoutHeader!=0) {
      m_valiTree->execTrig=pSimReadoutHeader->execNumber();
      
      const DayaBay::SimReadoutHeader::SimReadoutContainer 
        simReadout=pSimReadoutHeader->readouts();
      debug()<<"No. of SimReadout: "<<simReadout.size()<<endreq;
      debug()<<pSimReadoutHeader->execNumber()<<endreq;
    } // end of pSimReadoutHeader
    
    // check electrontic simulation information if available
    if(pElecHeader!=0) {
      m_valiTree->execElec=pElecHeader->execNumber();
    } // end of pElecHeader
    
    // check trigger header information if available
    if(pSimTrigHeader!=0) {
      debug()<<pSimTrigHeader->commandHeader()->collections().size();
    }
    
    // check detector simulation information if available
    if(pSimHeader!=0) {
      m_valiTree->execDets=pSimHeader->execNumber();
      const DayaBay::SimHitHeader* simhitheader = pSimHeader->hits();
      if(simhitheader)  {
        const DayaBay::SimHitHeader::hc_map& hcmap = simhitheader->hitCollection();
        DayaBay::SimHitHeader::hc_map::const_iterator it;
        
        m_valiTree->SimHitsEntries=0;
        int index(0);
      
        debug()<<"size of hit collection "<< hcmap.size()<<endreq;
        for (it=hcmap.begin(); it != hcmap.end(); ++it) {
          
          DayaBay::Detector det(it->first);
          debug() << "Got hit collection from " << det.detName()<<" id= "
                  << det.siteDetPackedData()<<endreq;
          
          const std::vector<DayaBay::SimHit*>& hitvec=it->second->collection();
          std::vector<DayaBay::SimHit*>::const_iterator it_hvec;
          
          debug() <<"size of hit vector "<< hitvec.size()<<endreq;
          for (it_hvec=hitvec.begin(); it_hvec!=hitvec.end(); ++it_hvec) {
            
            index=m_valiTree->SimHitsEntries;
            if(index<MaxSimHitEntries){
              
              m_valiTree->hitTime[index]  = (*it_hvec)->hitTime();
              m_valiTree->hitx[index]     = (*it_hvec)->localPos().x();
              m_valiTree->hity[index]     = (*it_hvec)->localPos().y();
              m_valiTree->hitz[index]     = (*it_hvec)->localPos().z();
              m_valiTree->sensID[index]   = (*it_hvec)->sensDetId();
              m_valiTree->weight[index]   = (*it_hvec)->weight();
              
              // optical photon's ancestor's pdg
              if((*it_hvec)->ancestor().track()) {
                m_valiTree->ancestorPdg[index] = (*it_hvec)->ancestor().track()->particle();
              }
              else {
                m_valiTree->ancestorPdg[index] = 0;
              }
              //optical photon's wavelength
              const DayaBay::SimPmtHit* pmthit = static_cast<const DayaBay::SimPmtHit*>(*it_hvec);
              m_valiTree->wavelength[index]   = (*pmthit).wavelength();
              m_valiTree->polx[index]   = (*pmthit).pol().x();
              m_valiTree->poly[index]   = (*pmthit).pol().y();
              m_valiTree->polz[index]   = (*pmthit).pol().z();
              m_valiTree->px[index]   = (*pmthit).dir().x();
              m_valiTree->py[index]   = (*pmthit).dir().y();
              m_valiTree->pz[index]   = (*pmthit).dir().z();
              
              m_valiTree->SimHitsEntries = index+1;
              
            }
            else {
              warning()<< " Reached the max hit limit!!!!!" <<endl;
            }
          }
        }
      }
    } // end of pSimHeader
    
    //  if(0==1) {
    if(pGenHeader!=0) {
      m_valiTree->execKine=pGenHeader->execNumber();
      debug()<<"Generator name: "<<pGenHeader->generatorName()<<endreq;
      debug()<<"execNum: "<<pGenHeader->execNumber()<<endreq;
      
      const HepMC::GenEvent* event=pGenHeader->event();
      //  event->print();
      
      HepMC::GenEvent::particle_const_iterator p;
      int count=0;
      for ( p = event->particles_begin(); p!=event->particles_end(); ++p ) {      
        (*p)->print();
        if((*p)->production_vertex()) {   // only outgoing particle has a production vertex
          if(  (*p)->status() == 1  ||       // Normal track
               abs((*p)->status()) > 1000 ) {     // MuonProphet muon
            count++;
            
            if(count==1) {
              
              m_valiTree->trk1_pdgId=(*p)->pdg_id();
              
              m_valiTree->trk1_px=(*p)->momentum().px();
              m_valiTree->trk1_py=(*p)->momentum().py();
              m_valiTree->trk1_pz=(*p)->momentum().pz();
              m_valiTree->trk1_e =(*p)->momentum().e();
              
              m_valiTree->trk1_x=(*p)->production_vertex()->position().x();
              m_valiTree->trk1_y=(*p)->production_vertex()->position().y();
              m_valiTree->trk1_z=(*p)->production_vertex()->position().z();
              m_valiTree->trk1_t=(*p)->production_vertex()->position().t();
              
              debug()<<"trk1_pdgId: "<<m_valiTree->trk1_pdgId<<endreq;
              
              if(m_valiTree->trk1_pdgId==22) { // gamma
                m_valiTree->trk1_ke=m_valiTree->trk1_e;
                m_valiTree->trk1_ve=m_valiTree->trk1_e;
              }else if (m_valiTree->trk1_pdgId==11) { // e-
                m_valiTree->trk1_ke=m_valiTree->trk1_e-m_e;
                m_valiTree->trk1_ve=m_valiTree->trk1_e-m_e;
              }else if (m_valiTree->trk1_pdgId==-11) { // e+
                m_valiTree->trk1_ke=m_valiTree->trk1_e-m_e;
                m_valiTree->trk1_ve=m_valiTree->trk1_e+m_e;
              }else if (m_valiTree->trk1_pdgId==1000020040) { // alpha
                m_valiTree->trk1_ke=m_valiTree->trk1_e-m_alpha;
                m_valiTree->trk1_ve=m_valiTree->trk1_e-m_alpha;
              }
            }  // outgoing particle 1
            if(count==2) {
              m_valiTree->trk2_pdgId=(*p)->pdg_id();
              
              m_valiTree->trk2_px=(*p)->momentum().px();
              m_valiTree->trk2_py=(*p)->momentum().py();
              m_valiTree->trk2_pz=(*p)->momentum().pz();
              m_valiTree->trk2_e =(*p)->momentum().e();
              
              m_valiTree->trk2_x=(*p)->production_vertex()->position().x();
              m_valiTree->trk2_y=(*p)->production_vertex()->position().y();
              m_valiTree->trk2_z=(*p)->production_vertex()->position().z();
              m_valiTree->trk2_t=(*p)->production_vertex()->position().t();
              if(m_valiTree->trk2_pdgId==22) { // gamma
                m_valiTree->trk2_ke=m_valiTree->trk2_e;
                m_valiTree->trk2_ve=m_valiTree->trk2_e;
              }else if (m_valiTree->trk2_pdgId==11) { // e-
                m_valiTree->trk2_ke=m_valiTree->trk2_e-m_e;
                m_valiTree->trk2_ve=m_valiTree->trk2_e-m_e;
              }else if (m_valiTree->trk2_pdgId==-11) { // e+
                m_valiTree->trk2_ke=m_valiTree->trk2_e-m_e;
                m_valiTree->trk2_ve=m_valiTree->trk2_e+m_e;
              }else if (m_valiTree->trk2_pdgId==1000020040) { // alpha
                m_valiTree->trk2_ke=m_valiTree->trk2_e-m_alpha;
                m_valiTree->trk2_ve=m_valiTree->trk2_e-m_alpha;
              }
            }  // outgoing particle 2
          }
        } // end of outgoing particles
      }
    } // end of pGenHeader

    m_valiTree->n = m_counter;
    m_valiTree->fill();
  }
    
  delete pIStageData;
  pIStageData=0;

  if( (m_CurrTime-m_StartTime).GetSeconds()> (m_TimeRange/CLHEP::second) ) {
    info()<<m_TimeRange<<" seconds detector simulation finished. Issue scheduled stop."<<endreq;
    m_EvtLoopMgr->stopRun();
  }

  return StatusCode::SUCCESS;
}

StatusCode Sim15::finalize()
{
  if(m_monitor) {
    m_valiTree->close();
    delete m_valiTree;
  }

  return this->GaudiAlgorithm::finalize();
}

---