In This Package:

AdPlotterAlg.cc

Go to the documentation of this file.

#include "AdPlotterAlg.h"

#include "Context/Context.h"
#include "Conventions/Detectors.h"
#include "DataSvc/ICableSvc.h"
#include "DetHelpers/ICoordSysSvc.h"
#include "DetDesc/IGeometryInfo.h"
#include "Event/ReadoutHeader.h"
#include "Event/Readout.h"
#include "Event/ReadoutPmtCrate.h"
#include "Event/ReadoutPmtChannel.h"
#include "Event/SimReadoutHeader.h"
#include "Event/ElecHeader.h"
#include "Event/SimHeader.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/GenHeader.h"
#include "HepMC/GenEvent.h"

#include "GaudiKernel/ITHistSvc.h"
#include "GaudiKernel/SystemOfUnits.h"

#include "TH2F.h"
#include "TH1F.h"
#include "TStyle.h"

using namespace std;
using namespace DayaBay;

struct AdSummary {
  TimeStamp triggerTime;
  double adcSum;
};

class AdExamplePlots {
  DayaBay::Detector detector;
  MsgStream& log;
  bool firstReadout;
  TimeStamp lastTriggerTime;
  vector<AdSummary> adHistory;
public:
  
  // public histograms
  // Generated Vertices
  TH2F* genVertexZRho;
  // Simulated Hits
  TH1F* nPhotoElectrons;
  // Readouts
  TH1F* hitPmts;
  TH1F* dtReadout;
  TH1F* dtReadoutShort;
  TH1F* adcSpectrum;
  TH1F* tdcSpectrum;
  TH1F* chargeSum;
  TH2F* hitPmtVsChargeSum;
  
  // generate coincidence figures
  bool findCoincidence;
  TH1F* dtNeutronCapture;
  TH1F* promptEnergy;
  TH1F* delayedEnergy;
  TH2F* delayedVsPromptEnergy;
  
  AdExamplePlots(DayaBay::Detector det, MsgStream& msg) 
    : detector(det), 
      log(msg),
      firstReadout(true),
      findCoincidence(false)
  { }
  
  StatusCode book(ITHistSvc* hsvc, string filepath) {
    // Initialize and register all the histograms
    gStyle->SetHistLineWidth(2);
    gStyle->SetHistLineColor(4);
    {
      string name = detector.detName() + "_HitPmts";
      string title = "Number of Hit Pmts for " + detector.detName();
      hitPmts = new TH1F(name.c_str(),title.c_str(),200,0,200);
      hitPmts->GetXaxis()->SetTitle("Number of Hit PMTs");
      hitPmts->GetYaxis()->SetTitle("Readouts");
      name = filepath+"/"+detector.detName()+"/hitPmts";
      if (hsvc->regHist(name,hitPmts).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete hitPmts; hitPmts = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_DtReadout";
      string title = "Time between readouts for " + detector.detName();
      dtReadout = new TH1F(name.c_str(),title.c_str(),200,0,0.02);
      dtReadout->GetXaxis()->SetTitle("#Deltat [s]");
      dtReadout->GetYaxis()->SetTitle("Readouts / bin");
      name = filepath+"/"+detector.detName()+"/dtReadout";
      if (hsvc->regHist(name,dtReadout).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete dtReadout; dtReadout = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_DtReadoutShort";
      string title = "Short time between readouts for " + detector.detName();
      dtReadoutShort = new TH1F(name.c_str(),title.c_str(),200,0,100.e-6);
      dtReadoutShort->GetXaxis()->SetTitle("#Deltat [s]");
      dtReadoutShort->GetYaxis()->SetTitle("Readouts / bin");
      name = filepath+"/"+detector.detName()+"/dtReadoutShort";
      if (hsvc->regHist(name,dtReadoutShort).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete dtReadoutShort; dtReadoutShort = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_AdcSpectrum";
      string title = "ADC spectrum for all channels in " + detector.detName();
      adcSpectrum = new TH1F(name.c_str(),title.c_str(),4096,0,4096);
      adcSpectrum->GetXaxis()->SetTitle("ADC value");
      adcSpectrum->GetYaxis()->SetTitle("Entries");
      name = filepath+"/"+detector.detName()+"/adcSpectrum";
      if (hsvc->regHist(name,adcSpectrum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete adcSpectrum; adcSpectrum = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_TdcSpectrum";
      string title = "TDC spectrum for all channels in " + detector.detName();
      tdcSpectrum = new TH1F(name.c_str(),title.c_str(),300,0,300);
      tdcSpectrum->GetXaxis()->SetTitle("TDC value");
      tdcSpectrum->GetYaxis()->SetTitle("Entries");
      name = filepath+"/"+detector.detName()+"/tdcSpectrum";
      if (hsvc->regHist(name,tdcSpectrum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete tdcSpectrum; tdcSpectrum = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_ChargeSum";
      string title = "ADC Charge sum for " + detector.detName();
      chargeSum = new TH1F(name.c_str(),title.c_str(),500,0,20000);
      chargeSum->GetXaxis()->SetTitle("Sum(ADC - Baseline)");
      chargeSum->GetYaxis()->SetTitle("Readouts");
      name = filepath+"/"+detector.detName()+"/chargeSum";
      if (hsvc->regHist(name,chargeSum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete chargeSum; chargeSum = 0;
        return StatusCode::FAILURE;
      }
    }
    {
      string name = detector.detName() + "_HitPmtVsChargeSum";
      string title = "Number of Hit PMTs versus the Charge Sum for" 
        + detector.detName();
      hitPmtVsChargeSum = new TH2F(name.c_str(),title.c_str(),
                                   200,0,20000,
                                   200,0,200);
      hitPmtVsChargeSum->GetXaxis()->SetTitle("Sum(ADC - Baseline)");
      hitPmtVsChargeSum->GetYaxis()->SetTitle("Number of Hit PMTs");
      name = filepath+"/"+detector.detName()+"/hitPmtVsChargeSum";
      if (hsvc->regHist(name,hitPmtVsChargeSum).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete hitPmtVsChargeSum; hitPmtVsChargeSum = 0;
        return StatusCode::FAILURE;
      }
    }
    // Simulated Hits
    {
      string name = detector.detName() + "_NPhotoElectrons";
      string title = "Number of simulated photoelectrons for " 
        + detector.detName();
      nPhotoElectrons = new TH1F(name.c_str(),title.c_str(),1000,0,1000);
      nPhotoElectrons->GetXaxis()->SetTitle("Number of Photoelectrons");
      nPhotoElectrons->GetYaxis()->SetTitle("Readouts");
      name = filepath+"/"+detector.detName()+"/nPhotoElectrons";
      if (hsvc->regHist(name,nPhotoElectrons).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete nPhotoElectrons; nPhotoElectrons = 0;
        return StatusCode::FAILURE;
      }
    }
    // Generated Vertices
    {
      string name = detector.detName() + "_GenVertexZRho";
      string title = "Position of generated vertices for " + detector.detName();
      genVertexZRho = new TH2F(name.c_str(),title.c_str(),
                               200,0,3,
                               200,-5,5);
      genVertexZRho->GetXaxis()->SetTitle("Particle (#rho / 3 m)^{2}");
      genVertexZRho->GetYaxis()->SetTitle("Particle z [m]");
      name = filepath+"/"+detector.detName()+"/genVertexZRho";
      if (hsvc->regHist(name,genVertexZRho).isFailure()) {
        log << MSG::ERROR << "Could not register " << name << endreq;
        delete genVertexZRho; genVertexZRho = 0;
        return StatusCode::FAILURE;
      }
    }
    if(findCoincidence){
      // add coincidence histograms
      {
        // IBD Neutron Capture Time
        string name = detector.detName() + "_DtNeutronCapture";
        string title = "IBD Neutron capture spectrum for " + detector.detName();
        dtNeutronCapture = new TH1F(name.c_str(),title.c_str(),200,0,100.e-6);
        dtNeutronCapture->GetXaxis()->SetTitle("#Deltat [s]");
        dtNeutronCapture->GetYaxis()->SetTitle("Readouts");
        name = filepath+"/"+detector.detName()+"/dtNeutronCapture";
        if (hsvc->regHist(name,dtNeutronCapture).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete dtNeutronCapture; dtNeutronCapture = 0;
          return StatusCode::FAILURE;
        }       
      }
      {
        // IBD Prompt Energy
        string name = detector.detName() + "_PromptEnergy";
        string title = "IBD Prompt Energy for " + detector.detName();
        promptEnergy = new TH1F(name.c_str(),title.c_str(),200,0,12);
        promptEnergy->GetXaxis()->SetTitle("Prompt Energy [MeV]");
        promptEnergy->GetYaxis()->SetTitle("Readouts");
        name = filepath+"/"+detector.detName()+"/promptEnergy";
        if (hsvc->regHist(name,promptEnergy).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete promptEnergy; promptEnergy = 0;
          return StatusCode::FAILURE;
        }       
      }
      {
        // IBD Delayed Energy
        string name = detector.detName() + "_DelayedEnergy";
        string title = "IBD Delayed Energy for " + detector.detName();
        delayedEnergy = new TH1F(name.c_str(),title.c_str(),200,0,12);
        delayedEnergy->GetXaxis()->SetTitle("Delayed Energy [MeV]");
        delayedEnergy->GetYaxis()->SetTitle("Readouts");
        name = filepath+"/"+detector.detName()+"/delayedEnergy";
        if (hsvc->regHist(name,delayedEnergy).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete delayedEnergy; delayedEnergy = 0;
          return StatusCode::FAILURE;
        }       
      }
      {
        // IBD Delayed vs. Prompt Energy
        string name = detector.detName() + "_DelayedVsPromptEnergy";
        string title = "IBD Delayed vs. Prompt Energy for " + detector.detName();
        delayedVsPromptEnergy = new TH2F(name.c_str(),title.c_str(),
                                         200,0,12,
                                         200,0,12);
        delayedVsPromptEnergy->GetXaxis()->SetTitle("Prompt Energy [MeV]");
        delayedVsPromptEnergy->GetYaxis()->SetTitle("Delayed Energy [MeV]");
        name = filepath+"/"+detector.detName()+"/delayedVsPromptEnergy";
        if (hsvc->regHist(name,delayedVsPromptEnergy).isFailure()) {
          log << MSG::ERROR << "Could not register " << name << endreq;
          delete delayedVsPromptEnergy; delayedVsPromptEnergy = 0;
          return StatusCode::FAILURE;
        }       
      }
    }
    return StatusCode::SUCCESS;
  }

  StatusCode fill(ICableSvc* /*cable*/, const ServiceMode &/*sm*/,
                  ICoordSysSvc* coord, const ReadoutHeader* roh) 
  {
    
    const DayaBay::ReadoutPmtCrate* ro = 
      dynamic_cast<const DayaBay::ReadoutPmtCrate*>(roh->readout());
    if (!ro) {
      log << MSG::DEBUG << "Got a non-PmtCrate readout, skipping..." << endreq;
      return StatusCode::SUCCESS;
    }
  
    double adcChargeSum = 0;
    int nHitPmts = 0;
    const ReadoutPmtCrate::PmtChannelReadouts& 
      channelMap = ro->channelReadout();
    ReadoutPmtCrate::PmtChannelReadouts::const_iterator
      it, done = channelMap.end();
    for (it=channelMap.begin(); it != done; ++it) {
      // Ignore calibration PMTs
      if(it->first.board() > 12) continue;
      nHitPmts++;
    
      // ADC data
      double baseline = 500.;
      float charge = 0;
      const vector<int> &adc = it->second.adc();
      vector<int>::const_iterator cit, cdone=adc.end();
      for (cit=adc.begin(); cit != cdone; ++cit) {
        adcSpectrum->Fill(*cit);
        charge += *cit - baseline;
      }
      adcChargeSum += charge;      

      // TDC data
      const vector<int> &tdc = it->second.tdc();
      vector<int>::const_iterator tdcIter, tdcDone = tdc.end();
      for (tdcIter=tdc.begin(); tdcIter != tdcDone; ++tdcIter) {
        tdcSpectrum->Fill(*tdcIter);
      }

    }
    hitPmts->Fill(nHitPmts);
    chargeSum->Fill(adcChargeSum);
    hitPmtVsChargeSum->Fill(adcChargeSum, nHitPmts);
  
    // Time between readouts
    if( !firstReadout ){
      TimeStamp dt = ro->triggerTime() - lastTriggerTime;
      dtReadout->Fill( dt.GetSeconds() );
      dtReadoutShort->Fill( dt.GetSeconds() );
    }else{
      firstReadout = false;
    }
    lastTriggerTime = ro->triggerTime();

    // If Simulation headers are available, fill additional figures
    const SimReadoutHeader* simReadoutHeader = 0;
    const ElecHeader* elecHeader = 0;
    vector<const SimHeader*> simHeaders;
    vector<const GenHeader*> genHeaders;
    vector<const IHeader*> inputHeaders = roh->inputHeaders();
    vector<const IHeader*>::const_iterator headerIter, 
      headerDone = inputHeaders.end();
    log << MSG::INFO << "Processing Input Headers: " << inputHeaders.size() 
        << endreq;
    for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
         ++headerIter) {
      simReadoutHeader = dynamic_cast<const SimReadoutHeader*>(*headerIter);
      if(simReadoutHeader){
        log << MSG::INFO << " Found SimReadoutHeader" << endreq;
        break;
      }
    }
    if(simReadoutHeader){
      inputHeaders = simReadoutHeader->inputHeaders();
      headerDone = inputHeaders.end();
      for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
           ++headerIter) {
        elecHeader = dynamic_cast<const ElecHeader*>(*headerIter);
        if(elecHeader){
          log << MSG::INFO << " Found ElecHeader" << endreq;
          break;
        }
      }
    }
    if(elecHeader){
      inputHeaders = elecHeader->inputHeaders();
      headerDone = inputHeaders.end();
      for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
         ++headerIter) {
        const SimHeader* simHeader 
          = dynamic_cast<const SimHeader*>(*headerIter);
        if(simHeader){
          log << MSG::INFO << " Found SimHeader" << endreq;
          simHeaders.push_back(simHeader);
        }
      }
    }
    if(simHeaders.size() > 0){
      for(unsigned int simHeaderIdx=0; simHeaderIdx<simHeaders.size(); 
          simHeaderIdx++){
        inputHeaders = (simHeaders[simHeaderIdx])->inputHeaders();
        headerDone = inputHeaders.end();
        for (headerIter=inputHeaders.begin(); headerIter != headerDone; 
             ++headerIter) {
          const GenHeader* genHeader 
            = dynamic_cast<const GenHeader*>(*headerIter);
          if(genHeader){
            log << MSG::INFO << " Found GenHeader" << endreq;
            genHeaders.push_back(genHeader);
          }
        }
      }
    }
    
    // Fill additional histograms
    // Loop over simulated hits
    for(unsigned int simHeaderIdx=0; simHeaderIdx<simHeaders.size(); 
        simHeaderIdx++){
      log << MSG::INFO << " simHeaderIdx = " << simHeaderIdx << endreq;
      const SimHeader* simHeader = simHeaders[simHeaderIdx];
      const SimHitHeader* shh = simHeader->hits();
      if(!shh){
        log << MSG::INFO << "Can't find hit header." << endreq;
        continue;
      }
      short int detId = detector.siteDetPackedData();
      SimHitHeader::hc_map::const_iterator hcmapIter 
        = shh->hitCollection().find(detId);
      if(hcmapIter == shh->hitCollection().end()){
        log << MSG::INFO << "Can't find hit collection for detector " 
            << detId << endreq;
        log << MSG::INFO << "Have hit collections for " 
            << shh->hitCollection().size() << " detectors: "; 
        SimHitHeader::hc_map::const_iterator hcmapDone 
          = shh->hitCollection().end();
        for(hcmapIter = shh->hitCollection().begin(); hcmapIter != hcmapDone;
            hcmapIter++){
          log << MSG::INFO << hcmapIter->first << "\t";
        }
        continue;
      }
      const SimHitCollection* hc = hcmapIter->second;
      if(!hc) continue;
      int nPeHits = 0;
      const SimHitCollection::hit_container& hits = hc->collection();
      SimHitCollection::hit_container::const_iterator hitIter, 
        hitEnd = hits.end();
      for(hitIter = hits.begin(); hitIter != hitEnd; hitIter++){
        const SimHit* hit = *hitIter;
        if(!hit) continue;
        // Ignore calibration pmts
        if(AdPmtSensor(hit->sensDetId()).ring()<1) continue;
        nPeHits++;
      }
      log << MSG::INFO << " Npe = " << nPeHits << endreq;
      nPhotoElectrons->Fill(nPeHits);
    }
    
    for(unsigned int genHeaderIdx=0; genHeaderIdx<genHeaders.size(); 
        genHeaderIdx++){
      const GenHeader* genHeader = genHeaders[genHeaderIdx];
      const HepMC::GenEvent* event = genHeader->event();
      if(!event) continue;
      HepMC::GenEvent::vertex_const_iterator vtxIter, 
        vtxDone = event->vertices_end();
      for(vtxIter = event->vertices_begin(); vtxIter != vtxDone; vtxIter++){
        const HepMC::GenVertex* vtx = *vtxIter;
        if(!vtx) continue;
        double x = vtx->position().x();
        double y = vtx->position().y();
        double z = vtx->position().z();
        log << MSG::INFO << "\n" << vtx->position().t() << endreq;
        log << MSG::INFO << " Vertex time: " << vtx->position().t() << endreq;
        log << MSG::INFO << "\n" << vtx->position().t() << endreq;
        Gaudi::XYZPoint globalPoint(x,y,z);
        IDetectorElement* coordDE = coord->coordSysDE(globalPoint);
        if(!coordDE){
          log << MSG::ERROR << "Failed to find detector element for point: "
              << x << "\t" << y << "\t" << z 
              << endreq;
          continue;
        }
        Gaudi::XYZPoint localPoint = coordDE->geometry()->toLocal(globalPoint);
        double rho_3meters = localPoint.Rho()/(3.*Gaudi::Units::meter);
        double rhoSq = rho_3meters * rho_3meters;
        double z_meters = localPoint.Z()/Gaudi::Units::meter;
        genVertexZRho->Fill(rhoSq, z_meters);
        if( rhoSq < 0 || rhoSq > 3. || z_meters>5. || z_meters<-5. ){
          log << MSG::VERBOSE << "vtx position: " << localPoint.x() << "\t" 
              << localPoint.y() << "\t" << localPoint.z() 
              << endreq;
        }
        // Example loop over particles in vertex
        /*
        HepMC::GenVertex::particles_out_const_iterator pIter, 
          pDone = vtx->particles_out_const_end();
        for(pIter = vtx->particles_out_const_begin(); pIter != pDone; pIter++){
          const HepMC::GenParticle* particle = *pIter;
          if(!particle) continue;
          log << MSG::VERBOSE << "particle energy: " 
              << particle->momentum().e() << endreq;    
        }
        */
        // End example
      }
    }

    if(findCoincidence){
      double dtCut_s = 100.0e-6; 
      double epCutLow = 0.7;
      double epCutHigh = 8.0;
      double edCutLow = 6.0;
      double edCutHigh = 10.0;
      double chargeScale = 14000./8.;
      double ed = adcChargeSum / chargeScale;

      vector<AdSummary> newHistory;
      vector<AdSummary>::iterator histIter, histEnd = adHistory.end();
      for(histIter = adHistory.begin(); histIter != histEnd; histIter++){
        double dt = ro->triggerTime() - histIter->triggerTime;
        if(dt < dtCut_s){
          AdSummary promptSummary = *histIter;
          newHistory.push_back(promptSummary);
          double ep = promptSummary.adcSum / chargeScale;
          dtNeutronCapture->Fill(dt);
          if( ed > edCutLow && ed < edCutHigh ) promptEnergy->Fill(ep);
          if( ep > epCutLow && ep < epCutHigh ) delayedEnergy->Fill(ed);
          delayedVsPromptEnergy->Fill(ep,ed);
        }
      }
      AdSummary currentSummary;
      currentSummary.triggerTime = ro->triggerTime();
      currentSummary.adcSum = adcChargeSum;
      newHistory.push_back(currentSummary);
      adHistory = newHistory;
    }

    return StatusCode::SUCCESS;
  }
  
};


AdPlotterAlg::AdPlotterAlg(const std::string& name, ISvcLocator* pSvcLocator)
  : GaudiAlgorithm(name,pSvcLocator)
  , m_hsvc(0), m_cable(0), m_coord(0)
{
  declareProperty("Location",m_location=DayaBay::ReadoutHeaderLocation::Default,
                  "Location in the TES to ReadoutHeaders");
  declareProperty("FilePath",m_filepath="/file1/adplotter/",
                  "File path of registered histograms.");
  declareProperty("FindCoincidence",m_findCoincidence=false,
                  "Generate coincidence figures");
}

AdPlotterAlg::~AdPlotterAlg()
{
}

StatusCode AdPlotterAlg::initialize()
{
  this->GaudiAlgorithm::initialize();
  
  // Get the histogram service
  if ( service("THistSvc", m_hsvc).isFailure()) {
    error() << " No THistSvc available." << endreq;
    return StatusCode::FAILURE;
  } 
  // Get the cable service
  if ( service("StaticCableSvc", m_cable).isFailure()) {
    error() << " No CableService available." << endreq;
    return StatusCode::FAILURE;
  } 

  // Get the coord sys service
  if ( service("CoordSysSvc", m_coord).isFailure()) {
    error() << " No CoordSysSvc available." << endreq;
    return StatusCode::FAILURE;
  } 
  
  return StatusCode::SUCCESS;
}

StatusCode AdPlotterAlg::execute()
{
  DayaBay::ReadoutHeader* roh = get<DayaBay::ReadoutHeader>(m_location);
  
  const Context &ctx = roh->context();
  DayaBay::Detector det(ctx.GetSite(),ctx.GetDetId());
  AdExamplePlots* plots = m_plots[det];
  if (!plots) {
    plots = new AdExamplePlots(det,info());
    if(m_findCoincidence){
      info() << "Adding IBD coincidence figure generation" << endreq;
      plots->findCoincidence = m_findCoincidence;
    }else{
      info() << "Do not add IBD coincidence figure generation" << endreq;
    }
    m_plots[det] = plots;
    plots->book(m_hsvc, m_filepath);
  }
  
  return plots->fill(m_cable,ServiceMode(ctx,0),m_coord,roh);
}

StatusCode AdPlotterAlg::finalize()
{
  return this->GaudiAlgorithm::finalize();
}

---