In This Package:

AdCalibFigs.cc

Go to the documentation of this file.

#include "AdCalibFigs.h"

#include "StatisticsSvc/IStatisticsSvc.h"
#include "Event/ReadoutHeader.h"
#include "Event/CalibReadoutHeader.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/CalibReadoutPmtChannel.h"
#include "Conventions/Electronics.h"

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

#include <sstream>

char timeFormat_AdCalibFigs[] = "#splitline{%H:%M:%S}{%b %d %Y}";

// Another frustrating feature of ROOT: timezones are a pain.
// It doesn't matter if you set the china local time here,
// it will always be drawn in the local time of the
// machine which draws the histogram. :(

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

AdCalibFigs::~AdCalibFigs()
{
}

StatusCode AdCalibFigs::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 AdCalibFigs::execute()
{
  // Add the current event into histograms
  DayaBay::CalibReadoutHeader* calibHeader = 
    get<DayaBay::CalibReadoutHeader>("/Event/CalibReadout/CalibReadoutHeader");
  if(!calibHeader){
    error() << "Failed to get calibrated readout header." << endreq;
    return StatusCode::FAILURE;
  }

  const DayaBay::CalibReadout* readout = calibHeader->calibReadout();
  if(!readout){
    error() << "Failed to get calibrated readout from header" << endreq;
    return StatusCode::FAILURE;
  }
  
  if(!readout->detector().isAD()){
    // Not an AD, continue
    return StatusCode::SUCCESS;
  }

  // Convert to PMT crate readout
  const DayaBay::CalibReadoutPmtCrate* pmtReadout
    = dynamic_cast<const DayaBay::CalibReadoutPmtCrate*>(readout);
  if(!pmtReadout){
    error() << "Failed to get PMT readout" << endreq;
    return StatusCode::FAILURE;
  }

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

  // Find histograms
  const DayaBay::Detector& detector = readout->detector();
  TH1F* adCharge = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                           detector,
                                                           0, 0, ADCHARGE));
  TH1F* logAdCharge = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              0, 0, 
                                                              LOGADCHARGE));
  TH2F* pmtCharge = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, 
                                                            PMTCHARGE));
  TH2F* adChargeVsTime =  dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                               detector,
                                                               0, 0, 
                                                               ADCHARGEVSTIME));
  TH2F* adChargeVsDtTrigger =  dynamic_cast<TH2F*>(this->getOrMakeHist(
                                                         runNumber,
                                                         detector,
                                                         0, 0, 
                                                         ADCHARGEVSDTTRIGGER));
  TH2F* maxPmtChargeVsAdCharge =  dynamic_cast<TH2F*>(this->getOrMakeHist(
                                                      runNumber,
                                                      detector,
                                                      0, 0, 
                                                      MAXPMTCHARGEVSADCHARGE));
  TH2F* nChannelsVsAdCharge = dynamic_cast<TH2F*>(this->getOrMakeHist(
                                                 runNumber,
                                                 detector,
                                                 0, 0, 
                                                 NCHANNELSVSADCHARGE));
  while(readout->triggerTime().GetSeconds() > 
        adChargeVsTime->GetXaxis()->GetXmax()){
    this->extendRange(adChargeVsTime, 60);
  }

  const DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts& channels
    = pmtReadout->channelReadout();
    
  DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator channelIter,
    channelEnd = channels.end();

  double adSumCharge = 0;
  double maxPmtCharge = 0;
  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
    const DayaBay::CalibReadoutPmtChannel& channel = *channelIter;

    DayaBay::AdPmtSensor pmtId( channel.pmtSensorId().fullPackedData() );
    int column = pmtId.column();
    int ring = pmtId.ring();

    // Channel Hit Map
    double charge = channel.sumCharge();
    pmtCharge->Fill(column, ring, charge);
    adSumCharge += charge;
    if( charge > maxPmtCharge ) maxPmtCharge = charge;
  }

  if(m_lastTriggerTime.find(readout->detector()) != m_lastTriggerTime.end()){
    TimeStamp dtTriggerTime = readout->triggerTime();
    dtTriggerTime.Subtract(m_lastTriggerTime[readout->detector()]);
    if(adSumCharge>0){
      adChargeVsDtTrigger->Fill(dtTriggerTime.GetSeconds()*1e6,
                                TMath::Log10( adSumCharge ) );
    }
  }

  adCharge->Fill( adSumCharge );
  if(adSumCharge>0){
    double logAdSumCharge = TMath::Log10( adSumCharge );
    logAdCharge->Fill( logAdSumCharge );
    if(adChargeVsTime->GetYaxis()->GetXmin()<logAdSumCharge
       && adChargeVsTime->GetYaxis()->GetXmax()>logAdSumCharge){
      adChargeVsTime->Fill(readout->triggerTime().GetSeconds(), logAdSumCharge );
    }
    maxPmtChargeVsAdCharge->Fill(logAdSumCharge, maxPmtCharge/adSumCharge);
    nChannelsVsAdCharge->Fill(logAdSumCharge, channels.size());
  }

  pmtCharge->SetNormFactor( pmtCharge->GetNormFactor()+1 );
  m_lastTriggerTime[readout->detector()] = readout->triggerTime();

  return StatusCode::SUCCESS;
}

StatusCode AdCalibFigs::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* AdCalibFigs::getOrMakeHist(int run, const DayaBay::Detector& detector,
                                int column, int ring, 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*[MAXCALIBHISTS];
    for(int i=0; i<MAXCALIBHISTS; i++) histograms[i] = 0;
    m_shortCuts[run] = histograms;
  }else{
    // Found run
    histograms = histIter->second;
  }
  
  TH1* hist = 
    histograms[siteInt * NDETECTORS * NCOLUMNS * NRINGS * NCALIBHISTOGRAMS
               + detectorInt * NCOLUMNS * NRINGS * NCALIBHISTOGRAMS
               + column * NRINGS * NCALIBHISTOGRAMS
               + ring * NCALIBHISTOGRAMS
               + histogram];
  if(!hist){
    // Make this histogram
    std::string histName;
    if(detector.detectorId()){
      // Make detector histogram
      switch(histogram){
      case ADCHARGE:
        // AD Calibrated Charge
        histName = "adCharge";
        hist = new TH1F(histName.c_str(),"AD Calibrated Charge",
                        1100,-100,1000);
        hist->GetXaxis()->SetTitle("Photoelectrons");
        hist->GetYaxis()->SetTitle("Events");
        break;
      case LOGADCHARGE:
        // Log of AD Calibrated Charge
        histName = "logAdCharge";
        hist = new TH1F(histName.c_str(),"AD Calibrated Charge",
                        800,-1,7);
        hist->GetXaxis()->SetTitle("log_{10}( Photoelectrons )");
        hist->GetYaxis()->SetTitle("Events");
        break;
      case PMTCHARGE:
        // Sum of PMT charge
        histName = "pmtCharge";
        hist = new TH2F(histName.c_str(),"Mean PMT Charge by Column/Ring",
                        49,0.25,24.75,19,-0.75,8.75);
        hist->GetXaxis()->SetTitle("PMT Column");
        hist->GetYaxis()->SetTitle("PMT Ring");
        hist->SetStats(0);
        m_normalize.push_back(hist);
        hist->SetNormFactor(0);
        break;
      case ADCHARGEVSTIME:
        // AD Charge vs. time
        histName = "adChargeVsTime";
        hist = new TH2F(histName.c_str(),
                        "AD Charge vs. Time",
                        60,m_firstTriggerTime->GetSec(),
                        m_firstTriggerTime->GetSec()+60,
                        200,-1,7);
        hist->GetXaxis()->SetTitle("Run Time");
        hist->GetYaxis()->SetTitle("log_{10}( Photoelectrons )");
        hist->GetXaxis()->SetTimeDisplay(1);
        hist->GetXaxis()->SetTimeFormat(timeFormat_AdCalibFigs);
        hist->GetXaxis()->SetTimeOffset(0,"gmt");
        hist->GetXaxis()->SetNdivisions(505);
        break;
      case ADCHARGEVSDTTRIGGER:
        // AD Charge vs DT trigger
        histName = "adChargeVsDtTrigger";
        hist = new TH2F(histName.c_str(),"AD Charge vs time to last event",
                        200,0,100,200,-1,7);
        hist->GetXaxis()->SetTitle("#DeltaT [us]");
        hist->GetYaxis()->SetTitle("log_{10}( Photoelectrons )");
        break;
      case MAXPMTCHARGEVSADCHARGE:
        // Max PMT Charge vs. AD Charge
        histName = "maxPmtChargeVsAdCharge";
        hist = new TH2F(histName.c_str(),"Fraction of Charge in one PMT",
                        200,-1,7,110,0,1.1);
        hist->GetXaxis()->SetTitle("log_{10}( Photoelectrons )");
        hist->GetYaxis()->SetTitle("PMT / Total");
        break;
      case NCHANNELSVSADCHARGE:
        // Number of Hit Channels vs. AD Charge
        histName = "nChannelsVsAdCharge";
        hist = new TH2F(histName.c_str(),
                        "Number of hit channels vs. Total AD Charge",
                        200,-1,7,200,0,200);
        hist->GetXaxis()->SetTitle("log_{10}( Photoelectrons )");
        hist->GetYaxis()->SetTitle("Hit Channels");
        break;
      default:
        error() << "Unknown Histogram: " << histogram << endreq;
        return 0;
      }
    }

    debug() << "Making histogram: " << histName << endreq;
    m_statsSvc->put( this->getPath(run, detector.detName().c_str(), column, 
                                   ring, histName.c_str()), 
                     hist );
    histograms[siteInt * NDETECTORS * NCOLUMNS * NRINGS * NCALIBHISTOGRAMS
               + detectorInt * NCOLUMNS * NRINGS * NCALIBHISTOGRAMS
               + column * NRINGS * NCALIBHISTOGRAMS
               + ring * NCALIBHISTOGRAMS
               + histogram] = hist;
  }

  return hist;
}

std::string AdCalibFigs::getPath(int run, const char* detectorName,
                                 int column, int ring, 
                                 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;
    if(column){
      path << "/channel_column" << std::setfill('0') << std::setw(2) 
           << column << "_ring" << std::setfill('0') << std::setw(2)
           << ring;
    }
  }
  path << "/" << histName;
  return path.str();
}

void AdCalibFigs::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;
}

---