/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Private Attributes
Sim15 Class Reference

#include <Sim15.h>

Collaboration diagram for Sim15:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 Sim15 (const std::string &name, ISvcLocator *pSvcLocator)
virtual ~Sim15 ()
virtual StatusCode initialize ()
virtual StatusCode execute ()
virtual StatusCode finalize ()

Private Attributes

std::string m_TopStageName
 Fifteen begin.
int m_TopStageNum
IStagem_TopStage
TimeStamp m_StartTime
 Fifteen end.
TimeStamp m_CurrTime
double m_TimeRange
IEventProcessor * m_EvtLoopMgr
bool m_monitor
 create a root tree
ValidationTreem_valiTree
int m_counter
IStageDataManagerm_sdmgr

Detailed Description

Definition at line 26 of file Sim15.h.


Constructor & Destructor Documentation

Sim15::Sim15 ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Fifteen begin

Fifteen end

Definition at line 65 of file Sim15.cc.

                                                           :
  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 ( ) [virtual]

Definition at line 79 of file Sim15.cc.

{
}

Member Function Documentation

StatusCode Sim15::initialize ( ) [virtual]

Convert TimeRange to [second]

Definition at line 83 of file Sim15.cc.

{
  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 ( ) [virtual]

reset tree

fill sim hit SimHitHeader

typedef std::map<short int,DayaBay::SimHitCollection*> hc_map

loop over every SimHitCollection

loop over every SimHit

check all kinematic information if available

fill the tree

Definition at line 138 of file Sim15.cc.

{
  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 ( ) [virtual]

Definition at line 487 of file Sim15.cc.

{
  if(m_monitor) {
    m_valiTree->close();
    delete m_valiTree;
  }

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

Member Data Documentation

std::string Sim15::m_TopStageName [private]

Fifteen begin.

Definition at line 40 of file Sim15.h.

int Sim15::m_TopStageNum [private]

Definition at line 41 of file Sim15.h.

Definition at line 42 of file Sim15.h.

Fifteen end.

Definition at line 45 of file Sim15.h.

Definition at line 46 of file Sim15.h.

double Sim15::m_TimeRange [private]

Definition at line 47 of file Sim15.h.

IEventProcessor* Sim15::m_EvtLoopMgr [private]

Definition at line 48 of file Sim15.h.

bool Sim15::m_monitor [private]

create a root tree

Definition at line 51 of file Sim15.h.

Definition at line 52 of file Sim15.h.

int Sim15::m_counter [private]

Definition at line 53 of file Sim15.h.

Definition at line 56 of file Sim15.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:07:33 for Stage by doxygen 1.7.4