In This Package:

AdReconFigs.cc

Go to the documentation of this file.

#include "AdReconFigs.h"

#include "StatisticsSvc/IStatisticsSvc.h"
#include "Event/RecHeader.h"
#include "Event/ReadoutHeader.h"
#include "Conventions/Reconstruction.h"
#include "Conventions/Calibration.h"

#include "TH2.h"
#include "TMath.h"
#include "TList.h"

#include <sstream>

AdReconFigs::AdReconFigs(const std::string& name, 
                         ISvcLocator* pSvcLocator)
  : GaudiAlgorithm(name,pSvcLocator),
    m_statsSvc(0),
    m_firstTriggerTime(0)
{
  declareProperty("ReconLocation", m_reconLocation="/Event/Rec/AdSimple",
                  "Generate histograms for this reconstruction TES path");
}

AdReconFigs::~AdReconFigs()
{
}

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

StatusCode AdReconFigs::execute()
{
  // Add the current event into histograms
  DayaBay::RecHeader* recHeader = 
    get<DayaBay::RecHeader>(m_reconLocation);
  if(!recHeader){
    error() << "Failed to get reconstructed header: " << m_reconLocation
            << endreq;
    return StatusCode::FAILURE;
  }

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

  // Wow... Look at the trouble we have to go through to get the run number
  int runNumber = 0;
  {
    // Add the current event into histograms
    DayaBay::ReadoutHeader* readoutHeader = 
      get<DayaBay::ReadoutHeader>("/Event/Readout/ReadoutHeader");
    if(!readoutHeader){
      error() << "Failed to get readout header." << endreq;
      return StatusCode::FAILURE;
    }
    
    const DayaBay::DaqCrate* ro = readoutHeader->daqCrate();
    if(!ro){
      error() << "Failed to get readout from header" << endreq;
      return StatusCode::FAILURE;
    }
    
    runNumber = ro->runNumber();
  }

  if(!m_firstTriggerTime){
    m_firstTriggerTime = new TimeStamp(recTrigger->triggerTime());
  }

  // Find histograms
  const DayaBay::Detector& detector = recTrigger->detector();
  TH1F* adEnergy = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                           detector,
                                                           ADENERGY));
  TH2F* adZVsRho = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                           detector,
                                                           ADZVSRHO));
  TH2F* adYVsX = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                         detector,
                                                         ADYVSX));
  TH2F* adZVsXcal = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            ADZVSXCAL));

  if(recTrigger->energyStatus()==ReconStatus::kGood){
    adEnergy->Fill(recTrigger->energy());
  }
  if(recTrigger->positionStatus()==ReconStatus::kGood){
    // FIXME: Use CLHEP units
    double x = recTrigger->position().x() / 1000.;
    double y = recTrigger->position().y() / 1000.;
    double z = recTrigger->position().z() / 1000.;
    double rhoSq = x*x+y*y;
    double xCal = x*TMath::Cos(DayaBay::Calibration::AD_AxisB_phi)
      + y*TMath::Sin(DayaBay::Calibration::AD_AxisB_phi);
    adZVsRho->Fill( rhoSq/4., z );
    adYVsX->Fill(x, y);
    adZVsXcal->Fill(xCal, z);
  }

  m_lastTriggerTime[recTrigger->detector()] = recTrigger->triggerTime();

  return StatusCode::SUCCESS;
}

StatusCode AdReconFigs::finalize()
{
  // Handle Normalized histograms
  for(unsigned int normIdx=0; normIdx<m_normalize.size(); normIdx++){
    TH1* hist = m_normalize[normIdx];
    hist->Sumw2();
    if(hist->GetNormFactor()>0){
      double scale = 1./hist->GetNormFactor();
      hist->SetNormFactor(1);
      hist->Scale( scale );
    }
    hist->SetBit(TH1::kIsAverage);
  }
  m_normalize.clear();
  // Clean-up histogram shortcuts
  std::map<int,TH1**>::iterator histIter, histEnd = m_shortCuts.end();
  for(histIter = m_shortCuts.begin(); histIter!=histEnd; histIter++){
    delete [] (histIter->second);
    histIter->second = 0;
  }
  if( m_statsSvc ) m_statsSvc->release();
  if(m_firstTriggerTime){
    delete m_firstTriggerTime;
    m_firstTriggerTime = 0;
  }
  return StatusCode::SUCCESS;
}


TH1* AdReconFigs::getOrMakeHist(int run, const DayaBay::Detector& detector,
                                int histogram)
{
  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 = int(detector.detectorId());
  std::map<int,TH1**>::iterator histIter = m_shortCuts.find(run);
  TH1** histograms = 0;
  if( histIter == m_shortCuts.end() ){
    // Initialize histogram shortcuts
    histograms = new TH1*[MAXRECHISTS];
    for(int i=0; i<MAXRECHISTS; i++) histograms[i] = 0;
    m_shortCuts[run] = histograms;
  }else{
    // Found run
    histograms = histIter->second;
  }
  
  TH1* hist = 
    histograms[siteInt * NDETECTORS * NRECHISTOGRAMS
               + detectorInt * NRECHISTOGRAMS
               + histogram];
  if(!hist){
    // Make this histogram
    std::string histName;
    std::string histTitle;
    if(detector.detectorId()){
      // Make detector histogram
      std::string reconStyle = "Unknown"; 
      for(int strIdx=m_reconLocation.size()-1; strIdx>=0; strIdx--){
        if(m_reconLocation[strIdx]=='/'){
          reconStyle=std::string(m_reconLocation,strIdx+1);
          break;
        }
      }
      switch(histogram){
      case ADENERGY:
        // AD Energy 
        histName = std::string("adEnergy_") + reconStyle;
        histTitle = std::string("AD Energy (") + reconStyle
          + std::string(")");
        hist = new TH1F(histName.c_str(),histTitle.c_str(),
                        1100,-1,10);
        hist->GetXaxis()->SetTitle("Energy [MeV]");
        hist->GetYaxis()->SetTitle("Events / keV");
        break;
      case ADZVSRHO:
        // AD Z vs Rho 
        histName = std::string("adZVsRho_") + reconStyle;
        histTitle = std::string("AD Reconstructed Position (") + reconStyle
          + std::string(")");
        hist = new TH2F(histName.c_str(),histTitle.c_str(),
                        400,0,1.6,400,-2,2);
        hist->GetXaxis()->SetTitle("( #rho / 2 m )^{2}");
        hist->GetYaxis()->SetTitle("z [m]");
        break;
      case ADYVSX:
        // AD Y vs X 
        histName = std::string("adYVsX_") + reconStyle;
        histTitle = std::string("AD Reconstructed Position (") + reconStyle
          + std::string(")");
        hist = new TH2F(histName.c_str(),histTitle.c_str(),
                        400,-3,3,400,-3,3);
        hist->GetXaxis()->SetTitle("x [m]");
        hist->GetYaxis()->SetTitle("y [m]");
        break;
      case ADZVSXCAL:
        // AD Z vs X_calibration  
        histName = std::string("adZVsXcal_") + reconStyle;
        histTitle = std::string("AD Reconstructed Position (") + reconStyle
          + std::string(")");
        hist = new TH2F(histName.c_str(),histTitle.c_str(),
                        400,-3,3,400,-2,2);
        hist->GetXaxis()->SetTitle("x_{calib} [m]");
        hist->GetYaxis()->SetTitle("z [m]");
        break;
      default:
        error() << "Unknown Histogram: " << histogram << endreq;
        return 0;
      }
    }
    debug() << "Making histogram: " << histName << endreq;
    m_statsSvc->put( this->getPath(run, detector.detName().c_str(), 
                                   histName.c_str()), 
                     hist );
    histograms[siteInt * NDETECTORS * NRECHISTOGRAMS
               + detectorInt * NRECHISTOGRAMS
               + histogram] = hist;
  }

  return hist;
}

std::string AdReconFigs::getPath(int run, const char* detectorName,
                                 const char* histName)
{
  // Construct histogram path in statistics service
  std::stringstream path;
  path << "/file1/diagnostics/run_" << std::setfill('0') << std::setw(7) 
       << run;
  if(detectorName){
    path << "/detector_" << detectorName;
  }
  path << "/" << histName;
  return path.str();
}

void AdReconFigs::extendRange(TH1* hist, int xBinsAdd, int yBinsAdd/*=0*/)
{
  // Automatically extend the histogram range, and replace
  //  Currently works for TH1F and TH2F
  info() << "Extending histogram range: " << hist->GetName() << endreq;
  TH1* extension = 0;
  if(hist->IsA() == TH2F::Class()){
    extension = new TH2F("","",
                         xBinsAdd + hist->GetNbinsX(),
                         hist->GetXaxis()->GetXmin(),
                         hist->GetXaxis()->GetXmax() 
                         + xBinsAdd*hist->GetXaxis()->GetBinWidth(0),
                         yBinsAdd + hist->GetNbinsY(),
                         hist->GetYaxis()->GetXmin(),
                         hist->GetYaxis()->GetXmax() 
                         + yBinsAdd*hist->GetYaxis()->GetBinWidth(0));
  }else if(hist->IsA() == TH1F::Class()){
    extension = new TH1F("","",
                         xBinsAdd + hist->GetNbinsX(),
                         hist->GetXaxis()->GetXmin(),
                         hist->GetXaxis()->GetXmax() 
                         + xBinsAdd*hist->GetXaxis()->GetBinWidth(0));
  }else{
    error() << "Cannot extend range of type: " << hist->ClassName() << endreq;
    return;
  }
  TList list;
  list.Add(extension);
  hist->Merge(&list);
  delete extension;
  return;
}

---