In This Package:

DetSimProc.cc

Go to the documentation of this file.

#include "DetSimProc.h"

#include "G4DataHelpers/G4DhHit.h"

#include "Event/GenHeader.h"

#include "GiGa/IGiGaSvc.h"
#include "G4DataHelpers/IHistoryKeeper.h"

#include "G4Event.hh"
#include "G4PrimaryVertex.hh"

#include "G4DataHelpers/IHepMCtoG4.h"

#include "Event/GenHeader.h"

#include "MuonProphet/MpMuonFate.h"

DetSimProc::DetSimProc(const std::string& name, ISvcLocator* pSvcLocator)
    : StageProcessor<DayaBay::SimHeader>(name,pSvcLocator)
    , m_giga(0)
    , m_converter(0)
    , m_historyKeeper(0)
    , m_genPrune(0)
{
  DS_constructor();

  // Auxiliary
  declareProperty("GenPruneName",m_genPruneName = "GenPruneTool", "Name of GenHeader pruner");

}

DetSimProc::~DetSimProc()
{
  DS_destructor();
}

StatusCode DetSimProc::initialize()
{
  StatusCode sc = this->StageProcessor<DayaBay::SimHeader>::initialize();
  if (sc.isFailure()) return sc;

  debug()<<"DetSimProc"<<endreq;

  sc = DS_initialize();
  if (sc.isFailure()) return sc;

  // Although initialized to 0, however, no special meaning
  m_CurrentTime=0;
  m_SimDataList.clear();

  m_Start=true;

  // Auxiliary
  m_genPrune = tool<IGenPruneTool>(m_genPruneName);
  if( !m_genPrune ) {
    error() << "Failed to retrive GenPruneTool" << endreq;
    return StatusCode::FAILURE;
  }

  return StatusCode::SUCCESS;
}

StatusCode DetSimProc::execute()
{
  debug()<<"DetSimProc"<<endreq;

  SimData::DataList::iterator psimdatalist;

  FFTimeStamp tm(0,0);
  FFTimeStamp dt(0,0); 
  FFTimeStamp localtime(0,0);

  if(lowerStage()) {
    debug() << "Consumer&Producer" << endreq;
    if(m_Start||m_CurrentTime<=thisStage()->currentTime()) {

      if(m_Start) {
        // nothing to do here
      } else {
        psimdatalist=m_SimDataList.begin();
        localtime=psimdatalist->first;
      }

      //debug()<<"A: local time= "<<localtime<<endreq;
      //debug()<<"A: lower stage time= "<<m_LowerStage->currentTime()<<endreq;
      while(m_Start||localtime>=lowerStage()->currentTime()) {
        m_Start=false;
        //debug() << "In while loop, pulling data from lower Stage ......" << endreq;

        IStageData* pIStageData=0;
        StatusCode sc = lowerStage()->nextElement(pIStageData);  // pulling out
        if (sc.isFailure()) {
          error() << "Failed to pull kinematics" << endreq;
          return sc;
        }
        
        GnrtrData* pGnrtrData=0;
        try{
          pGnrtrData = dynamic_cast<GnrtrData*>(pIStageData);
        }
        catch(...) {
          error() << "Failed to get GnrtrData pointer" <<endreq;
        }

        //>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        // delete the used one
        tm=pGnrtrData->time();
        //debug() << "data pulled out from lower stage has time: "<<tm<<endreq;

        preExecute();
        SimData*   pSimData=0;   
        DS_execute(pGnrtrData,pSimData);
        this->AppendInputHeader(&pGnrtrData->header());
        //m_genPrune->prune( &pGnrtrData->header() );
        //
        postExecute();
        delete pGnrtrData; pGnrtrData=0;

        //debug() << "new data generated at time: "<< pSimData->time()<< endreq;
        //<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<

        // push back into a local store
        m_SimDataList.insert(SimData::DataList::value_type(pSimData->time(),pSimData));

        psimdatalist=m_SimDataList.begin();
        localtime=psimdatalist->first;
        //debug()<<"B: local time= "<<localtime<<endreq;
        //debug()<<"B: lower stage time= "<<m_LowerStage->currentTime()<<endreq;        
      }

      m_CurrentTime=localtime;
      //debug() << "m_CurrentTime= "<< m_CurrentTime << endreq;

      psimdatalist=m_SimDataList.begin();
      debug() << "to grep: new data pushed out at time " << psimdatalist->first << endreq;
      thisStage()->pushElement(psimdatalist->second);
      this->registerData(*(psimdatalist->second));
      m_SimDataList.erase(psimdatalist);

    }
  }
  return StatusCode::SUCCESS;  
}

StatusCode DetSimProc::finalize()
{
  DS_finalize();

  return StatusCode::SUCCESS;
}

void DetSimProc::DS_constructor()
{
  declareProperty("Converter",m_converterName="HepMCtoG4",
                  "Name of tool to convert HepMCEvents to G4PrimaryVertex");

  declareProperty("GenLocation",m_genLocation=DayaBay::GenHeaderLocation::Default,
                  "Location in the TES where the input GenHeader is to be found.");
}

void DetSimProc::DS_destructor()
{}

StatusCode DetSimProc::DS_initialize()
{

  m_giga = svc<IGiGaSvc>("GiGa",true);
  try {
    m_converter = tool<IHepMCtoG4>(m_converterName);
  }
  catch (const GaudiException& gex) {
    fatal() << "Failed to get converter tool \""
            << m_converterName << "\"" << endreq;
    return StatusCode::FAILURE;
  }

  m_historyKeeper = svc<IHistoryKeeper>("HistoryKeeper",true);

  return StatusCode::SUCCESS;
}

StatusCode DetSimProc::DS_execute(GnrtrData* p_input, SimData*& p_output)
{

//DayaBay::GenHeader* gen_header = get<DayaBay::GenHeader>(m_location);
  DayaBay::GenHeader* gen_header = &p_input->header();

  //debug()<<"genheader, timeStamp: "<<gen_header->timeStamp()<<endreq;
  //debug()<<"genheader, earliest: "<<gen_header->earliest()<<endreq;

  HepMC::GenEvent* genevt = gen_header->event();
  std::vector<G4PrimaryVertex*> g4verts;
  m_converter->convert(*genevt,g4verts);

  debug() << "Got " << g4verts.size() << " primary vertices" << endreq;

  for (size_t ind=0; ind < g4verts.size(); ++ind) {

    *m_giga << g4verts[ind];

  }

  debug() << "Pushed " << gen_header->generatorName()
         << " to giga" << endreq;

  DayaBay::SimHeader* header = MakeHeaderObject();

  // Just pass through GenHeader's timestamp.  This also causes
  // GenHeader to be registered as input, something that would
  // normally just happen if DsPushKine and DsPullEvent were the
  // same algorithm.
  //  DayaBay::GenHeader* gen_header = getTES<DayaBay::GenHeader>(m_genLocation);
  header->setTimeStamp(gen_header->timeStamp());

  // Primary event vertices.
  const G4Event* g4event = 0;
  m_giga->retrieveEvent(g4event);
  if (!g4event) {
    error() << "No G4Event!" << endreq;
    return StatusCode::FAILURE;
  }

  int nverts = g4event->GetNumberOfPrimaryVertex();
  if( nverts == 0 ) {
    warning() << "The g4event has zero primary vertices!" << endreq;
  }
  debug() << "Pulled event with " << nverts
         << " primary vertices, event id:" << g4event->GetEventID() << endreq;

  G4PrimaryVertex* g4vtx = g4event->GetPrimaryVertex(0);
  while (g4vtx) {
    debug() << "\n\tat (" << g4vtx->GetX0() << "," << g4vtx->GetY0() << "," << g4vtx->GetZ0() << ")";
    g4vtx = g4vtx->GetNext();
    break;
  }
  debug() << endreq;

  // particle histories.
  // Do this first so we can use it below.
  DayaBay::SimParticleHistory* history =0;
  m_historyKeeper->ClaimCurrentHistory(history); // This takes ownership from the Keeper.
  header->setParticleHistory(history);

  // Unobservable Statistics
  DayaBay::SimUnobservableStatisticsHeader* unobs =0;
  m_historyKeeper->ClaimCurrentUnobservable(unobs); // This takes ownership from the Keeper.
  header->setUnobservableStatistics(unobs);

  // introduce the headers to each other
  DayaBay::SimHitHeader* hit_header = 0;

  double earliestTime = 0;
  double latestTime = 0;
  Context context;
  context.SetSimFlag(SimFlag::kMC);
  bool firstDetector = true;
  int hitcount=0;  // deal with no hits situation

  // Hit collections.
  G4HCofThisEvent* hcs = g4event->GetHCofThisEvent();
  if (!hcs) {
    debug() << "No HitCollections in this event" << endreq;
  }  else  {
    int nhc = hcs->GetNumberOfCollections();
    if (!nhc) {
      debug() << "Number of HitCollections is zero" << endreq;
    }  else  {
      debug () << "# HitCollections = " << nhc << endreq;

      hit_header = new DayaBay::SimHitHeader(header);
      for (int ihc=0; ihc<nhc; ++ihc) {
        G4DhHitCollection* g4hc = dynamic_cast<G4DhHitCollection*>(hcs->GetHC(ihc));
        if (!g4hc) {
          error() << "Failed to get hit collection #" << ihc << endreq;
          return StatusCode::FAILURE;
        }
        
        // DetSim produces hit collections even for unsimulated detectors
        size_t nhits = g4hc->GetSize();
        hitcount+=nhits;
        if (!nhits) continue;
        
        bool firstHit = true;
        DayaBay::SimHitCollection::hit_container hits;
        DayaBay::Detector detector;
        DayaBay::SimHitCollection* shc =
          new DayaBay::SimHitCollection(hit_header,detector,hits);
        for (size_t ihit=0; ihit<nhits; ++ihit) {
          DayaBay::SimHit* simhit = (*g4hc)[ihit]->get();

          // A protection on time overflow or underflow.
          // Kill the hit over two years.
          // It is impossible to simulate them with full Geant4 simulation.
          // See Trac#489.
          if( abs(simhit->hitTime()) > 6.3e16 ) {
            warning() << "Eccentric long hit time "<<simhit->hitTime()<<"ns. Skipped."<<endreq;
            hitcount--;
            continue;
          }

          if(history) {
            int trackid = (*g4hc)[ihit]->trackId();
            simhit->setAncestor(history->track(trackid));
          }
          // Set detector and site using first hit
          if(firstHit){
            detector = DayaBay::Detector(simhit->sensDetId()).siteDetPackedData();
            if(firstDetector){
              // Keep track of sites/detectors in this simulation for context
              context.SetSite(detector.site());
              context.SetDetId(detector.detectorId());
            }
            if(context.GetSite() != detector.site()){
              // Simulation contains multiple sites, unset
              context.SetSite(Site::kUnknown);
            }
            if(context.GetDetId() != detector.detectorId()){
              // Simulation contains multiple detectors, unset
              context.SetDetId(DetectorId::kUnknown);
            }
          }
          // Record earliest/latest simhit
          debug() << "Hit Time: " << simhit->hitTime() << endreq;
          if((firstHit && firstDetector) || simhit->hitTime() < earliestTime)
            earliestTime = simhit->hitTime();
          if((firstHit && firstDetector) || simhit->hitTime() > latestTime)
            latestTime = simhit->hitTime();
          firstHit = false;

          simhit->setHc(shc);
          hits.push_back(simhit);
        }
        firstDetector = false;
        
        shc->setDetector(detector);
        shc->setCollection(hits);
        debug() << "Adding " << detector.detName() << " hits ("
                << shc->collection().size() << ") to header." << endreq;
        hit_header->addHitCollection(shc);
      }
    }
  }

  header->setHits(hit_header);
  context.SetTimeStamp(gen_header->timeStamp());
  header->setContext(context);

  // Set header time range based on hit times
  debug() << "Header Time sec, nanosec: " << header->timeStamp().GetSec()
          << "  " << header->timeStamp().GetNanoSec() << endreq;
  debug() << "Earliest Hit: " << earliestTime << endreq;
  debug() << "Latest Hit: " << latestTime << endreq;
  TimeStamp earliestTS(header->timeStamp());
  TimeStamp latestTS(header->timeStamp());
  earliestTS.Add(earliestTime * 1e-9);
  latestTS.Add(latestTime * 1e-9);
  debug() << "Earliest Time sec, nanosec: " << earliestTS.GetSec()
          << "  " << earliestTS.GetNanoSec() << endreq;
  debug() << "Latest Time sec, nanosec: " << latestTS.GetSec()
          << "  " << latestTS.GetNanoSec() << endreq;
  header->setEarliest(earliestTS);
  header->setLatest(latestTS);
  debug() << "Earliest Time sec, nanosec: " << header->earliest().GetSec()
          << "  " << header->earliest().GetNanoSec() << endreq;
  debug() << "Latest Time sec, nanosec: " << header->latest().GetSec()
          << "  " << header->latest().GetNanoSec() << endreq;

  HepMC::GenEvent::particle_const_iterator pci, pci_end = genevt->particles_end();

  for( pci = genevt->particles_begin(); pci != pci_end; ++pci ) {
    if( (*pci)->pdg_id() == 13 || (*pci)->pdg_id() == -13 ) {   // muon                                            
      if( (*pci)->status() != MpMuonFate::kNeedSim ) {   //  geant4-noble particle
        header->setEarliest(gen_header->timeStamp());
      }
    }
  }

  if(hitcount==0) {  // when there is no hit, set time stamp to genheader's earliest() to save causality
    header->setTimeStamp(gen_header->earliest());
    header->setEarliest(gen_header->earliest());
    header->setLatest(gen_header->earliest());
  }

  p_output=new SimData(*header);

  return StatusCode::SUCCESS;

}

StatusCode DetSimProc::DS_finalize()
{
  return StatusCode::SUCCESS;
}


// Local Variables:
// c-basic-offset: 2
// End:

---