In This Package:

DsPullEvent.cc

Go to the documentation of this file.

#include "DsPullEvent.h"

#include "G4DataHelpers/G4DhHit.h"

#include "Event/GenHeader.h"

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

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

DsPullEvent::DsPullEvent(const std::string& name, ISvcLocator* pSvcLocator)
    : DybAlgorithm<DayaBay::SimHeader>(name,pSvcLocator)
    , m_giga(0)
    , m_historyKeeper(0)
      //, m_converter(0)
{
    declareProperty("GenLocation",m_genLocation=DayaBay::GenHeaderLocation::Default,
                    "Location in the TES where the input GenHeader is to be found.");
}

DsPullEvent::~DsPullEvent()
{
}

StatusCode DsPullEvent::initialize()
{
    m_giga = svc<IGiGaSvc>("GiGa",true);
    m_historyKeeper = svc<IHistoryKeeper>("HistoryKeeper",true);

    m_capinfo = tool<INeutronCaptureInfo>("G4DhNeutronCaptureInfoTool");

    return StatusCode::SUCCESS;
}

StatusCode DsPullEvent::execute()
{
    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;
    }

    // reset Capture
    G4DhNeutronCapture capture;
    m_capinfo->addCapture(capture);

    int nverts = g4event->GetNumberOfPrimaryVertex();
    if( nverts == 0 ) {
        warning() << "The g4event has zero primary vertices!" << endreq;
        return StatusCode::SUCCESS;
    }


    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);

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

    // introduce the headers to each other
    DayaBay::SimHitHeader* hit_header = new DayaBay::SimHitHeader(header);
    header->setHits(hit_header);

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

    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();
          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);
    }

    // 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;
    context.SetTimeStamp(gen_header->timeStamp());
    header->setContext(context);

    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());
    }

    return StatusCode::SUCCESS;
}

StatusCode DsPullEvent::finalize()
{
    return StatusCode::SUCCESS;
}

---