In This Package:

ReconDataHistogram.cc

Go to the documentation of this file.

// General stuff

#include <sstream>
#include <string>
#include <stdio.h>

// Root stuff

#include "TGraph.h"
#include "TH2F.h"
#include "TH1F.h"
#include "TH1I.h"
#include <TProfile.h>
#include "TMath.h"
// Local Alg stuff

#include "ReconDataHistogram.h"

// Recon Header readout stuff
#include "StatisticsSvc/IStatisticsSvc.h"
#include "Event/RecHeader.h"
#include "Event/ReadoutHeader.h"
#include "Conventions/Reconstruction.h"
#include "Conventions/Calibration.h"

using namespace std;
using namespace DayaBay;
using namespace DybDaq;

ReconDataHistogram::ReconDataHistogram(const std::string& name, 
                                       ISvcLocator* pSvcLocator)
  : GaudiAlgorithm(name,pSvcLocator),
    m_statsSvc(0),
    m_eventCount(0),
    m_firstTriggerTime(0)
{
}

ReconDataHistogram::~ReconDataHistogram()
{
}

StatusCode ReconDataHistogram::initialize()
{
  // Initialize the necessary services
  StatusCode sc = service("StatisticsSvc",m_statsSvc,true);
  if(sc.isFailure()){
    error() << "Failed to get StatisticsSvc" << endreq;
    return sc;
  }
  return sc;
}

StatusCode ReconDataHistogram::execute()
{
  verbose() << "executing: " << m_eventCount << endreq;
  // Add the current event into histograms
  DayaBay::RecHeader* recHeader = 
    get<DayaBay::RecHeader>("/Event/Rec/AdSimple");
  if(!recHeader){
    error() << "Failed to get reconstructed header:" << "/Event/Rec/AdSimple"
            << endreq;
    return StatusCode::FAILURE;
  }

  const DayaBay::RecTrigger* recTrigger = &(recHeader->recTrigger());
  if(!recTrigger){
    error() << "Failed to get reconstructed trigger: " << "Event/Rec/AdSimple"
            << endreq;
    return StatusCode::FAILURE;
  }
  
  if(!recTrigger->detector().isAD()){
    // Not an AD, continue
    return StatusCode::SUCCESS;
  }
 
  const DayaBay::Detector& detector = recTrigger->detector();
  info() << "Detector = " << detector << endreq;

  if(!m_firstTriggerTime){
    m_firstTriggerTime = new TimeStamp(recTrigger->triggerTime());
  }
  
  TH2F* h_energyVsR = dynamic_cast<TH2F*>(getOrMakeHist(1, detector, kEnergyVsR));
  TH2F* h_energyVsZ = dynamic_cast<TH2F*>(getOrMakeHist(1, detector, kEnergyVsZ));
  TH1F* h_rateVsR = dynamic_cast<TH1F*>(getOrMakeHist(1, detector, kRateVsR));
  TH1F* h_rateVsZ = dynamic_cast<TH1F*>(getOrMakeHist(1, detector, kRateVsZ));
  
  if(recTrigger->positionStatus()==ReconStatus::kGood){
    double x = recTrigger->position().x() / 1000;
    double y = recTrigger->position().y() / 1000;
    double z = recTrigger->position().z() / 1000;
    double r = TMath::Sqrt(x*x+y*y);
    h_rateVsR->Fill(r);
    h_rateVsZ->Fill(z);
    if(recTrigger->energyStatus()==ReconStatus::kGood){
      h_energyVsR->Fill(r, recTrigger->energy());
      h_energyVsZ->Fill(z,recTrigger->energy());
    }
  }

  m_lastTriggerTime = recTrigger->triggerTime();
  
  ++m_eventCount;
  if (0 == (m_eventCount % 1000)) {
    info() << " Processed " << m_eventCount << " events" << endreq;
  }

  return StatusCode::SUCCESS;
}

StatusCode ReconDataHistogram::finalize()
{
  info() << "Finalizing " << endreq;

  for(unsigned int Idx=0; Idx<m_scale.size(); Idx++){
    TH1* hist = m_scale[Idx];
    hist->Sumw2();
    double duration = m_lastTriggerTime.GetSeconds() - m_firstTriggerTime->GetSeconds();
    hist->Scale(1/duration);
  }
  m_scale.clear();

  if( m_statsSvc ) {
    m_statsSvc->release();
  }

  if(m_firstTriggerTime){
    delete m_firstTriggerTime;
    m_firstTriggerTime = 0;
  }

  return StatusCode::SUCCESS;
}

TH1* ReconDataHistogram::getOrMakeHist(int run, const DayaBay::Detector& detector,
                                       int histIndex)
{
  debug() << "Detector=" << detector << "," << detector.fullPackedData() << endreq;

  int siteInt = 0;
  switch(detector.site()){
  case Site::kUnknown:
    siteInt = 0; break;
  case Site::kDayaBay:
    siteInt = 1; break;
  case Site::kLingAo:
    siteInt = 2; break;
  case Site::kFar:
    siteInt = 3; break;
  case Site::kSAB:
    siteInt = 4; break;
  default:
    error() << "Unknown site: " << detector.detName() << endreq;
    return 0;
  }

//  int detectorInt = detector.fullPackedData();
  int detectorInt = int(detector.detectorId());
  std::map<int,TH1**>::iterator histIter = m_hist.find(detectorInt);
  TH1** histograms = 0;
  if( histIter == m_hist.end() ){
    info() << "Create " << MAXRECHISTS << " histograms for " << detector << endreq;
    // Initialize histogram 
    histograms = new TH1*[MAXRECHISTS];
    for(int i=0; i<MAXRECHISTS; i++) histograms[i] = 0;
    m_hist[detectorInt] = histograms;
  }else{
    // Found detector
    histograms = histIter->second;
  }

  int iHist = siteInt * NDETECTORS * NRECHISTOGRAMS + detectorInt * NRECHISTOGRAMS + histIndex;
  //cout<<"siteInt:  "<<siteInt<<"  detectorInt:  "<<detectorInt<<"  iHist:  "<<iHist<<endreq;
  info() << "Histogram index = " << iHist << endreq;
  TH1* hist = histograms[siteInt * NDETECTORS * NRECHISTOGRAMS + detectorInt * NRECHISTOGRAMS + histIndex];
  if(!hist){
    // Make this histogram
    std::string histName;
    if(detector.detectorId()){
       
      // Make detector histIndex
      switch(histIndex){
        case kEnergyVsR:
          hist = new TH2F("h_sum_energyVsR","h_sum_energyVsR",250, 0, 2.5, 110, -1, 10);
          hist->GetXaxis()->SetTitle("R_{recon}/m");
          hist->GetYaxis()->SetTitle("E_{recon}/MeV");
          break;
        case kEnergyVsZ:
          hist = new TH2F("h_sum_energyVsZ","h_sum_energyVsZ",500, -2.5, 2.5, 110, -1, 10);
          hist->GetXaxis()->SetTitle("Z_{recon}/m");
          hist->GetYaxis()->SetTitle("E_{recon}/MeV");
          break;
        case kRateVsR:
          hist = new TH1F("h_sum_rateVsR","h_sum_rateVsR", 250, 0, 2.5);
          hist->GetXaxis()->SetTitle("R_{recon}/m");
          hist->GetYaxis()->SetTitle("EventRate/bin(1cm)");
          m_scale.push_back(hist);
          break;
        case kRateVsZ:
          hist = new TH1F("h_sum_rateVsZ","h_sum_rateVsZ",500, -2.5, 2.5);
          hist->GetXaxis()->SetTitle("Z_{recon}/m");
          hist->GetYaxis()->SetTitle("EventRate/bin(1cm)");
          m_scale.push_back(hist);
          break;
        default:
          error() << "Unknown histogram for event: " << histIndex << endreq;
          return 0;
      }
      
    }
    info() << "Making histogram: " << hist->GetName() << endreq;
    m_statsSvc->put( getPath(run, detector, hist->GetName()), hist );
    histograms[ siteInt * NDETECTORS * NRECHISTOGRAMS + detectorInt * NRECHISTOGRAMS + histIndex] = hist;
  }

  return hist;
}

std::string ReconDataHistogram::getPath(int run, const Detector &detector,
                                        const char* histName)
{
  // Construct histogram path in statistics service
  std::stringstream path;
  path << "/file1/";
  path << Site::AsString(detector.site()) << "/" << DetectorId::AsString(detector.detectorId());
  path << "/Recon";
  path << "/" << histName;
  debug() << "Histogram path = " << path << endreq;
  return path.str();
}

---