In This Package:

AdReadoutFigs.cc

Go to the documentation of this file.

#include "AdReadoutFigs.h"

#include "StatisticsSvc/IStatisticsSvc.h"
#include "Event/ReadoutHeader.h"
#include "Event/DaqPmtCrate.h"
#include "Event/DaqPmtChannel.h"
#include "Conventions/Electronics.h"

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

#include <sstream>


char timeFormat[] = "#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. :(

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

AdReadoutFigs::~AdReadoutFigs()
{
}

StatusCode AdReadoutFigs::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 AdReadoutFigs::execute()
{
  // 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* daqCrate = readoutHeader->daqCrate();
  if(!daqCrate){
    error() << "Failed to get DAQ crate from header" << endreq;
    return StatusCode::FAILURE;
  }
  
  // Convert to PMT crate readout
  const DayaBay::DaqPmtCrate* pmtCrate
      = daqCrate->asPmtCrate();
  if(!pmtCrate){
    // Not an AD, continue
    return StatusCode::SUCCESS;
  }

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

  // Find histograms
  int runNumber = pmtCrate->runNumber();
  const DayaBay::Detector& detector = pmtCrate->detector();
  TH1F* triggerRate = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, TRIGGERRATE));
  TH1F* dtTrigger = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, DTTRIGGER));
  TH1F* logDtTrigger = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                               detector,
                                                               0, 0, 
                                                               LOGDTTRIGGER));
  TH1F* nChannels = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, 
                                                            NCHANNELS));
  TH2F* nChannelsVsTime =  dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, 
                                                            NCHANNELSVSTIME));
  TH2F* channelHits = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              0, 0, 
                                                              CHANNELHITS));
  TH2F* channelDadcFine = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, 
                                                            CHANNELDADCFINE));
  TH2F* channelDadcCoarse = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            0, 0, 
                                                            CHANNELDADCCOARSE));
  TH2F* tdcVsBoard = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                             detector,
                                                             0, 0, 
                                                             TDCVSBOARD));

  if(pmtCrate->triggerTime().GetSeconds()-triggerRate->GetXaxis()->GetXmax() > 1800)
  {
    // Catch large gaps in data, and fail
    error() << "An invalid gap of " 
            << pmtCrate->triggerTime().GetSeconds()-triggerRate->GetXaxis()->GetXmax() 
            << " seconds has been found in the data.  Bailing..." << endreq;
    return  StatusCode::FAILURE;
  }

  while(pmtCrate->triggerTime().GetSeconds() > 
        triggerRate->GetXaxis()->GetXmax()){
    this->extendRange(triggerRate, 60);
  }
  while(pmtCrate->triggerTime().GetSeconds() > 
        nChannelsVsTime->GetXaxis()->GetXmax()){
    this->extendRange(nChannelsVsTime, 60);
  }

  const DayaBay::DaqPmtCrate::PmtChannelPtrList& channels
    = pmtCrate->channelReadouts();
    
  DayaBay::DaqPmtCrate::PmtChannelPtrList::const_iterator channelIter, 
    channelEnd = channels.end();

  triggerRate->Fill(pmtCrate->triggerTime().GetSeconds());
  nChannelsVsTime->Fill(pmtCrate->triggerTime().GetSeconds(), channels.size() );

  if(m_lastTriggerTime.find(pmtCrate->detector()) != m_lastTriggerTime.end()){
    TimeStamp dtTriggerTime = pmtCrate->triggerTime();
    dtTriggerTime.Subtract(m_lastTriggerTime[pmtCrate->detector()]);
    dtTrigger->Fill(dtTriggerTime.GetSeconds());
    logDtTrigger->Fill(TMath::Log10(dtTriggerTime.GetSeconds()));
  }

  nChannels->Fill( channels.size() );

  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
    const DayaBay::DaqPmtChannel& channel = *(*channelIter);

    int board = channel.channelId().board();
    int connector = channel.channelId().connector();
    TH1F* tdc = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                        detector,
                                                        board, connector, 
                                                        TDC));
    TH1F* adcFine = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                            detector,
                                                            board, connector, 
                                                            ADCFINE));
    TH1F* adcCoarse = dynamic_cast<TH1F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              board, connector, 
                                                              ADCCOARSE));
    TH2F* dadcVsTimeFine = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              board, 
                                                              connector, 
                                                              DADCVSTIMEFINE));
    TH2F* dadcVsPedFine = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              board, 
                                                              connector, 
                                                              DADCVSPEDFINE));
    TH2F* dadcFineVsTdc = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              board, 
                                                              connector, 
                                                              DADCFINEVSTDC));
    TH2F* pedFineVsTdc = dynamic_cast<TH2F*>(this->getOrMakeHist(runNumber,
                                                              detector,
                                                              board, 
                                                              connector, 
                                                              PEDFINEVSTDC));
    while(pmtCrate->triggerTime().GetSeconds() > 
          dadcVsTimeFine->GetXaxis()->GetXmax()){
      this->extendRange(dadcVsTimeFine, 60);
    }

    // Channel Hit Map
    channelHits->Fill(board, connector);

    // First hit histograms
    if(channel.hitCount()>0){
        if(channel.isHighGainAdc(0)){
        channelDadcFine->Fill(board, connector,
                              channel.deltaAdc(0));
        channelDadcFine->SetNormFactor( channelDadcFine->GetNormFactor()+1 );
      }else{
        channelDadcCoarse->Fill(board, connector,
                                channel.deltaAdc(0));
        channelDadcCoarse->SetNormFactor(channelDadcCoarse->GetNormFactor()+1);
      }
    }
    // All hits
    for(unsigned int hitIdx=0; hitIdx<channel.hitCount(); hitIdx++){
      double epsilon = 0.01;
      tdcVsBoard->Fill(board + connector/double(NCONNECTORS) + epsilon,
                       channel.tdc(hitIdx));
      tdc->Fill(channel.tdc(hitIdx));
      float deltaAdc = channel.deltaAdc(hitIdx);
      if(channel.isHighGainAdc(hitIdx)){
        adcFine->Fill(channel.adc(hitIdx));
        if(dadcVsTimeFine->GetYaxis()->GetXmin()<deltaAdc
           && dadcVsTimeFine->GetYaxis()->GetXmax()>deltaAdc){
          dadcVsTimeFine->Fill(pmtCrate->triggerTime().GetSeconds(), deltaAdc);
        }
        dadcVsPedFine->Fill(channel.preAdcAvg(hitIdx), deltaAdc);
        dadcFineVsTdc->Fill(channel.tdc(hitIdx), deltaAdc);
        pedFineVsTdc->Fill(channel.tdc(hitIdx), channel.preAdcAvg(hitIdx));
      }else{
        adcCoarse->Fill(channel.adc(hitIdx));
      }
    }
  }

  channelHits->SetNormFactor( channelHits->GetNormFactor()+1 );

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

  return StatusCode::SUCCESS;
}

StatusCode AdReadoutFigs::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* AdReadoutFigs::getOrMakeHist(int run, const DayaBay::Detector& detector,
                                  int board, int connector, 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*[MAXHISTS];
    for(int i=0; i<MAXHISTS; i++) histograms[i] = 0;
    m_shortCuts[run] = histograms;
  }else{
    // Found run
    histograms = histIter->second;
  }
  
  TH1* hist = 
    histograms[siteInt * NDETECTORS * NBOARDS * NCONNECTORS * NHISTOGRAMS
               + detectorInt * NBOARDS * NCONNECTORS * NHISTOGRAMS
               + board * NCONNECTORS * NHISTOGRAMS
               + connector * NHISTOGRAMS
               + histogram];
  if(!hist){
    // Make this histogram
    std::string histName;
    if(detector.detectorId()){
      if(board){
        // Make channel histograms
        switch(histogram){
        case TDC:
          // TDC Values
          histName = "tdc";
          hist = new TH1F(histName.c_str(),
                          "TDC",
                          1000,300,1300);
          hist->GetXaxis()->SetTitle("TDC");
          hist->GetYaxis()->SetTitle("Entries");
          break;
        case ADCFINE:
          // ADC Values (Fine Range)
          histName = "adcFine";
          hist = new TH1F(histName.c_str(),
                          "ADC (Fine Range)",
                          4096,0,4096);
          hist->GetXaxis()->SetTitle("ADC");
          hist->GetYaxis()->SetTitle("Entries");
          break;
        case ADCCOARSE:
          // ADC Values (Coarse Range)
          histName = "adcCoarse";
          hist = new TH1F(histName.c_str(),
                          "ADC (Coarse Range)",
                          4096,0,4096);
          hist->GetXaxis()->SetTitle("ADC");
          hist->GetYaxis()->SetTitle("Entries");
          break;
        case DADCVSTIMEFINE:
          // DeltaADC Values Vs Time (Fine Range)
          histName = "dadcVsTimeFine";
          hist = new TH2F(histName.c_str(),
                          "#DeltaADC vs Time (Fine Range)",
                          60,m_firstTriggerTime->GetSec(),
                          m_firstTriggerTime->GetSec()+60,
                          200,-100,100);
          hist->GetXaxis()->SetTitle("Run Time");
          hist->GetYaxis()->SetTitle("#DeltaADC");
          hist->GetXaxis()->SetTimeDisplay(1);
          hist->GetXaxis()->SetTimeFormat(timeFormat);
          hist->GetXaxis()->SetTimeOffset(0,"gmt");
          hist->GetXaxis()->SetNdivisions(505);
          break;
        case DADCVSPEDFINE:
          // DeltaADC Values Vs. PreADC (Fine Range)
          histName = "dadcVsPedFine";
          hist = new TH2F(histName.c_str(),
                          "#DeltaADC vs. PreADC (Fine Range)",
                          200,0,200,200,-100,100);
          hist->GetXaxis()->SetTitle("PreADC");
          hist->GetYaxis()->SetTitle("#DeltaADC");
          break;
        case DADCFINEVSTDC:
          // DeltaADC Values Vs. TDC (Fine Range)
          histName = "dadcFineVsTdc";
          hist = new TH2F(histName.c_str(),
                          "#DeltaADC vs. TDC (Fine Range)",
                          100,300,1300,200,-100,100);
          hist->GetXaxis()->SetTitle("TDC");
          hist->GetYaxis()->SetTitle("#DeltaADC");
          break;
        case PEDFINEVSTDC:
          // PreADC Values Vs. TDC (Fine Range)
          histName = "pedFineVsTdc";
          hist = new TH2F(histName.c_str(),
                          "PreADC vs. TDC (Fine Range)",
                          100,300,1300,200,0,200);
          hist->GetXaxis()->SetTitle("TDC");
          hist->GetYaxis()->SetTitle("PreADC");
          break;
        default:
          error() << "Unknown Histogram: " << histogram << endreq;
          return 0;
        }
      }else{
        // Make detector histogram
        switch(histogram){
        case TRIGGERRATE:
          // Trigger rate
          histName = "triggerRate";
          hist = new TH1F(histName.c_str(),"Trigger Rate",
                          60,m_firstTriggerTime->GetSec(),
                          m_firstTriggerTime->GetSec()+60);
          hist->GetXaxis()->SetTitle("Run Time");
          hist->GetYaxis()->SetTitle("Rate [Hz]");
          hist->GetXaxis()->SetTimeDisplay(1);
          hist->GetXaxis()->SetTimeFormat(timeFormat);
          hist->GetXaxis()->SetTimeOffset(0,"gmt");
          hist->GetXaxis()->SetNdivisions(505);
          break;
        case DTTRIGGER:
          // Time Between triggers
          histName = "dtTrigger";
          hist = new TH1F(histName.c_str(),"Time Between Triggers",
                          1000,0,0.015);
          hist->GetXaxis()->SetTitle("#Deltat [s]");
          break;
        case LOGDTTRIGGER:
          // Logarithm of Time Between triggers
          histName = "logDtTrigger";
          hist = new TH1F(histName.c_str(),"Time Between Triggers",
                          900,-8,1);
          hist->GetXaxis()->SetTitle("log_{10}(#Deltat / 1s)");
          break;
        case NCHANNELS:
          // Number of channels in event
          histName = "nChannels";
          hist = new TH1F(histName.c_str(),
                          "Number of Channels in Event",200,0,200);
          hist->GetXaxis()->SetTitle("Number of Channels");
          break;
        case NCHANNELSVSTIME:
          // Number of channels in event vs. time
          histName = "nChannelsVsTime";
          hist = new TH2F(histName.c_str(),
                          "Number of Channels in Event vs. Time",
                          60,m_firstTriggerTime->GetSec(),
                          m_firstTriggerTime->GetSec()+60,200,0,200);
          hist->GetXaxis()->SetTitle("Run Time");
          hist->GetYaxis()->SetTitle("Number of Channels");
          hist->GetXaxis()->SetTimeDisplay(1);
          hist->GetXaxis()->SetTimeFormat(timeFormat);
          hist->GetXaxis()->SetTimeOffset(0,"gmt");
          hist->GetXaxis()->SetNdivisions(505);
          break;
        case CHANNELHITS:
          // Channel hit rates
          histName = "channelHits";
          hist = new TH2F(histName.c_str(),"Occupancy of FEE Channels",
                          20,0.5,20.5,16,0.5,16.5);
          hist->GetXaxis()->SetTitle("FEE Board");
          hist->GetYaxis()->SetTitle("FEE Connector");
          m_normalize.push_back(hist);
          hist->SetNormFactor(0);
          hist->SetStats(0);
          break;
        case CHANNELDADCFINE:
          // Channel dADC (Fine ADC)
          histName = "channelDadcFine";
          hist = new TH2F(histName.c_str(),
                          "Mean Fine #DeltaADC on FEE Channels",
                          20,0.5,20.5,16,0.5,16.5);
          hist->GetXaxis()->SetTitle("FEE Board");
          hist->GetYaxis()->SetTitle("FEE Connector");
          m_normalize.push_back(hist);
          hist->SetNormFactor(0);
          hist->SetStats(0);
          break;
        case CHANNELDADCCOARSE:
          // Channel charge (Coarse ADC)
          histName = "channelDadcCoarse";
          hist = new TH2F(histName.c_str(),
                          "Mean Coarse #DeltaADC on FEE Channels",
                          20,0.5,20.5,16,0.5,16.5);
          hist->GetXaxis()->SetTitle("FEE Board");
          hist->GetYaxis()->SetTitle("FEE Connector");
          m_normalize.push_back(hist);
          hist->SetNormFactor(0);
          hist->SetStats(0);
          break;
        case TDCVSBOARD:
          // TDC values vs board
          histName = "tdcVsBoard";
          hist = new TH2F(histName.c_str(),
                          "TDC vs. Board",
                          NBOARDS*NCONNECTORS,0,NBOARDS,1000,300,1300);
          hist->GetXaxis()->SetTitle("FEE Board");
          hist->GetYaxis()->SetTitle("TDC");
          break;
        default:
          error() << "Unknown Histogram: " << histogram << endreq;
          return 0;
        }
      }
    }
    debug() << "Making histogram: " << histName << endreq;
    m_statsSvc->put( this->getPath(run, detector.detName().c_str(), board, 
                                   connector, histName.c_str()), 
                     hist );
    histograms[siteInt * NDETECTORS * NBOARDS * NCONNECTORS * NHISTOGRAMS
               + detectorInt * NBOARDS * NCONNECTORS * NHISTOGRAMS
               + board * NCONNECTORS * NHISTOGRAMS
               + connector * NHISTOGRAMS
               + histogram] = hist;
  }

  return hist;
}

std::string AdReadoutFigs::getPath(int run, const char* detectorName,
                                 int board, int connector, 
                                 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(board){
      path << "/channel_board" << std::setfill('0') << std::setw(2) 
           << board << "_connector" << std::setfill('0') << std::setw(2)
           << connector;
    }
  }
  path << "/" << histName;
  return path.str();
}

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

---