In This Package:

PmtCalibFullModel.cc

Go to the documentation of this file.

#include "PmtCalibFullModel.h"

#include "Event/ReadoutPmtCrate.h"
#include "StatisticsSvc/IStatisticsSvc.h"
#include "DataSvc/ICableSvc.h"
#include "DataSvc/ICalibDataSvc.h"
#include "Context/ServiceMode.h"
#include "TH1F.h"
#include "TH2.h"
#include "TF1.h"
#include "TParameter.h"
#include <math.h>

//#include "pmt_full_model_fit.h"
#include "NIMmodel.h"

PmtCalibFullModel::PmtCalibFullModel(const std::string& type,
                                         const std::string& name, 
                                         const IInterface* parent)
  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_statsSvc(0)
{
  declareInterface< IPmtCalibParamTool >(this) ;   
  declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc",
                  "Name of service to map between detector, hardware, and electronic IDs");
  declareProperty("CalibSvcName",m_calibSvcName="StaticCalibDataSvc",
                  "Name of service to access FEE channel baselines");
  declareProperty("FloatingFeePedestalSvcName", m_floatFeePedesSvcName="FloatingFeePedestalSvc",
                  "Name of service to access floating FEE channel baselines");
  declareProperty("UseFloatFeePedes", m_useFloatFeePedes=true,
                  "Use FloatFeePedes or not? In case of false, it will use CalibSvc.");
  declareProperty("FilePath",m_filepath="/file0/PmtCalibFullModel",
                  "File path of registered histograms.");
}

PmtCalibFullModel::~PmtCalibFullModel(){}

StatusCode PmtCalibFullModel::initialize()
{
  info() << "initialize()" << endreq;
  StatusCode sc = this->GaudiTool::initialize();
  if( sc != StatusCode::SUCCESS ) return sc;

  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
  if(!m_cableSvc){
    error() << "Failed to access cable svc: " << m_cableSvcName << endreq;
    return StatusCode::FAILURE;
  }

  // Get Calibration Service
  m_calibSvc = svc<ICalibDataSvc>(m_calibSvcName,true);
  if(!m_calibSvc){
    error() << "Failed to access calibration svc: " << m_calibSvcName << endreq;
    return StatusCode::FAILURE;
  }

  // Get floating fee pedestal service
  m_floatFeePedesSvc = svc<IFloatingFeePedestalSvc>(m_floatFeePedesSvcName,true);
  if(!m_floatFeePedesSvc){
    error() << "Failed to access FloatingFeePedestalSvc: " << m_floatFeePedesSvcName << endreq;
    return StatusCode::FAILURE;
  }

  // Get the histogram service
  m_statsSvc = svc<IStatisticsSvc>("StatisticsSvc",true);
  if(!m_statsSvc) {
    error() << " No StatisticsSvc available." << endreq;
    return StatusCode::FAILURE;
  }

  // Initialize data from input file.
  // Processed detectors are identified by directory name.
  std::vector<std::string> detDirs = m_statsSvc->getSubFolders(m_filepath);
  std::vector<std::string>::iterator dirIter, dirEnd = detDirs.end();
  for(dirIter=detDirs.begin(); dirIter != dirEnd; dirIter++){
    info() << "Processing input directory: " << m_filepath << "/" 
           << *dirIter << endreq;
    DayaBay::Detector detector( 
                 DayaBay::Detector::siteDetPackedFromString(*dirIter) );
    if( detector.site() == Site::kUnknown 
        || detector.detectorId() == DetectorId::kUnknown ){
      warning() << "Input directory: " << m_filepath << "/" 
                << *dirIter << " not processed." << endreq;
    }
    m_processedDetectors.push_back(detector);
  }

  return StatusCode::SUCCESS;
}

StatusCode PmtCalibFullModel::finalize()
{
  StatusCode sc = this->GaudiTool::finalize();
  if( sc != StatusCode::SUCCESS ) return sc;

  return StatusCode::SUCCESS;
}


StatusCode PmtCalibFullModel::process(const DayaBay::ReadoutHeader& readoutHeader){
  
  const DayaBay::DaqCrate* daqCrate = readoutHeader.daqCrate();
  if(!daqCrate){
    error() << "Failed to get DAQ crate from header" << endreq;
    return StatusCode::FAILURE;
  }
  
  const DayaBay::Detector& detector = daqCrate->detector();
  
  // Convert to PMT crate readout
  const DayaBay::DaqPmtCrate* pmtCrate
    = daqCrate->asPmtCrate();
  if(!pmtCrate){
    // Not an AD, continue
    debug() << "process(): Do not process detector: " << detector.detName() 
            << endreq;
    return StatusCode::SUCCESS;
  }
  
  Context context = readoutHeader.context();
  if( !hasStats(detector) ){
    // New detector; initialize all the detector statistics
    StatusCode sc = prepareStats(context);
    if( sc != StatusCode::SUCCESS ) return sc;
  }
  ServiceMode svcMode(context, 0);
  
  // Retrieve the parameters and histograms
  TObject* nReadoutsObj = 
    m_statsSvc->get(this->getPath(detector) + "/nReadouts");
  TParameter<int>* nReadoutsPar= dynamic_cast<TParameter<int>*>(nReadoutsObj);
  if(!nReadoutsPar) return StatusCode::FAILURE;
  TH1F* adcSumH = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSum");
  if(!adcSumH) return StatusCode::FAILURE;
  TH1F* tdcMedianH = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedian");
  if(!tdcMedianH) return StatusCode::FAILURE;
  
  // Add the current readout to the calibration statistics
  std::vector<int> readoutTdc;
  double adcSum = 0;
  // Loop over channels
  const DayaBay::DaqPmtCrate::PmtChannelPtrList& channels
    = pmtCrate->channelReadouts();
  
  DayaBay::DaqPmtCrate::PmtChannelPtrList::const_iterator channelIter, 
    channelEnd = channels.end();
  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
    const DayaBay::DaqPmtChannel& channel = *(*channelIter);
    const DayaBay::FeeChannelId& chanId = channel.channelId();
    //const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId,svcMode);
    
    
    // Get Parameters and Histograms
    TObject* nHitsObj = m_statsSvc->get(this->getPath(chanId) + "/nHits");
    TParameter<int>* nHitsPar = dynamic_cast<TParameter<int>*>(nHitsObj);
    if(!nHitsPar) return StatusCode::FAILURE;
    TH1F* tdcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcRaw");
    if(!tdcRawH) return StatusCode::FAILURE;
    TH1F* adcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRaw");
    if(!adcRawH) return StatusCode::FAILURE;
    TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
    if(!adcH) return StatusCode::FAILURE;
    TH1F* adcByClockH = m_statsSvc->getTH1F( this->getPath(chanId) 
                                             + "/adcByClock");
    if(!adcByClockH) return StatusCode::FAILURE;
    
    //std::vector<int>::const_iterator tdcIter, tdcEnd = pmtChan.tdc().end();
    // Loop over tdc values
    for(unsigned int i=0; i<channel.hitCount(); i++) {
      int tdc=channel.tdc(i);
      readoutTdc.push_back(tdc);
      tdcRawH->Fill(tdc);
    }
    
    // Loop over adc values
    //double adcBaseline = feeCalib->m_adcBaselineLow;
    for(unsigned int i=0; i<channel.hitCount(); i++) {
      int adcClock = channel.peakCycle(i);
      int adc = channel.adc(i);

      double adcBaseline;
      if( m_useFloatFeePedes ) {
        if(channel.isHighGainAdc(i)){
          adcBaseline = m_floatFeePedesSvc->pedestal( chanId, DayaBay::FeeGain::kHigh);
        }else{
          adcBaseline = m_floatFeePedesSvc->pedestal( chanId, DayaBay::FeeGain::kLow);
        }
        
        // In case floatFeePedesSvc didn't get enough statistics yet, use preAdc instead.
        if( adcBaseline<0 ) adcBaseline = channel.preAdcAvg(i);
      } else {
        // adcBaseline = feeCalib->m_adcBaselineLow;  // The pedestal run result is not good enoughreating CalibParamConfDbMerge.lock and Starting CalibParamConfDbMerge
        adcBaseline = channel.preAdcAvg(i);  // Use preAdc is much better
      }
  
      // Restrict to fine range gain measurement
      if( adc>0 && channel.isHighGainAdc(i) ) {
        adcRawH->Fill(adc);
        adcH->Fill(adc-adcBaseline);
        adcByClockH->Fill(adcClock,adc-adcBaseline);
        adcSum += adc-adcBaseline;
      }
    }
    // Increment number of hits on this channel
    nHitsPar->SetVal( nHitsPar->GetVal() + 1 );

  }
  if(readoutTdc.size() > 0){
    // Find the median TDC
    std::sort(readoutTdc.begin(), readoutTdc.end());
    int medianIdx = readoutTdc.size()/2;
    int medianTdc = readoutTdc[medianIdx];
    for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
      const DayaBay::DaqPmtChannel& channel = *(*channelIter);
      const DayaBay::FeeChannelId& chanId = channel.channelId();
      TH1F* tdcByMedianH = 
        m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMedian");
      if(!tdcByMedianH) return StatusCode::FAILURE;
      //std::vector<int>::const_iterator tdcIter, tdcEnd = pmtChan.tdc().end();
      // Loop over tdc values
      for(unsigned int i=0; i<channel.hitCount(); i++) {
        int tdc = channel.tdc(i);
        tdcByMedianH->Fill(tdc-medianTdc);
      }
    }
    tdcMedianH->Fill(medianTdc);
  }

  adcSumH->Fill(adcSum);
  // Increment number of readouts from this detector
  nReadoutsPar->SetVal( nReadoutsPar->GetVal() + 1 );
  return StatusCode::SUCCESS;  
}

StatusCode PmtCalibFullModel::calibrate(){
  // Loop over detectors in summary data, and calculate parameters

  std::vector<DayaBay::Detector>::iterator detIter, 
    detEnd = m_processedDetectors.end();
  for(detIter = m_processedDetectors.begin(); detIter != detEnd; detIter++){
    DayaBay::Detector detector = *detIter;

    // Retrieve the data description and histograms
    TObject* nReadoutsObj = 
      m_statsSvc->get(this->getPath(detector) + "/nReadouts");
    TParameter<int>* nReadoutsPar= dynamic_cast<TParameter<int>*>(nReadoutsObj);
    if(!nReadoutsPar) return StatusCode::FAILURE;

    TH1F* meanOccupancyH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/meanOccupancy");
    if(!meanOccupancyH) return StatusCode::FAILURE;

    TH1F* meanTdcOffsetH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/meanTdcOffset");
    if(!meanTdcOffsetH) return StatusCode::FAILURE;

    TH1F* gainChannelH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/gainChannel");
    if(!gainChannelH) return StatusCode::FAILURE;

    TH1F* tdcOffsetChannelH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/tdcOffsetChannel");
    if(!tdcOffsetChannelH) return StatusCode::FAILURE;

    TH1F* occupancyChannelH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/occupancyChannel");
    if(!occupancyChannelH) return StatusCode::FAILURE;

    TH2F* occupancyH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/occupancy");
    if(!occupancyH2D) return StatusCode::FAILURE;
    TH2F* gainH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/gain2D");
    if(!gainH2D) return StatusCode::FAILURE;
    TH1F* gainH1D = m_statsSvc->getTH1F( this->getPath(detector) 
                                              + "/gain1D");
    if(!gainH1D) return StatusCode::FAILURE;

    TObject* simFlagObj = m_statsSvc->get(this->getPath(detector) +"/simFlag");
    TParameter<int>* simFlagPar = dynamic_cast<TParameter<int>*>(simFlagObj);
    if(!simFlagPar) return StatusCode::FAILURE;

    TObject* timeSecObj = m_statsSvc->get(this->getPath(detector) +"/timeSec");
    TParameter<int>* timeSecPar = dynamic_cast<TParameter<int>*>(timeSecObj);
    if(!timeSecPar) return StatusCode::FAILURE;
    
    TObject* timeNanoSecObj = 
      m_statsSvc->get(this->getPath(detector) +"/timeNanoSec");
    TParameter<int>* timeNanoSecPar = 
      dynamic_cast<TParameter<int>*>(timeNanoSecObj);
    if(!timeNanoSecPar) return StatusCode::FAILURE;

    Site::Site_t site = detector.site();
    DetectorId::DetectorId_t detId = detector.detectorId();
    SimFlag::SimFlag_t simFlag = (SimFlag::SimFlag_t)(simFlagPar->GetVal());
    time_t timeSec = (time_t)(timeSecPar->GetVal());
    int timeNanoSec = timeNanoSecPar->GetVal();
    int nReadouts = nReadoutsPar->GetVal();

    // Context, ServiceMode
    Context context(site,simFlag,TimeStamp(timeSec,timeNanoSec),detId);
    ServiceMode svcMode(context, 0);

    // Prepare data for ring-to-ring comparison
    int nRings = 8;
    std::map<int,double> meanOccRing;
    std::map<int,double> meanOccWeight;
    std::map<int,double> meanTdcRing;
    std::map<int,double> meanTdcWeight;
    for(int ring = 1; ring <= nRings; ring++){
      meanOccRing[ring] = 0.0;
      meanOccWeight[ring] = 0.0;
      meanTdcRing[ring] = 0.0;
      meanTdcWeight[ring] = 0.0;
    }
    double minValidOcc = 0.01;  // Cut on PMTs with bad/low occupancy
    std::vector<DayaBay::FeeChannelId> badChannels;

    // Get list of all FEE channels
    const std::vector<DayaBay::FeeChannelId>& channelList 
      = m_cableSvc->feeChannelIds( svcMode );
    std::vector<DayaBay::FeeChannelId>::const_iterator chanIter, 
      chanEnd = channelList.end();
    // Initialize statistics for each channel
    for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++){
      DayaBay::FeeChannelId chanId = *chanIter;

      // Analyze the channel data
      // Determine ring and column
      DayaBay::AdPmtSensor pmtId = 
        m_cableSvc->adPmtSensor(chanId, svcMode);
      int ring = pmtId.ring();
      int column = pmtId.column();

      // Channel status
      bool occupancyOk = true;
      bool gainOk = true;

      // Get Parameters and Histograms
      TObject* nHitsObj = m_statsSvc->get(this->getPath(chanId) + "/nHits");
      TParameter<int>* nHitsPar = dynamic_cast<TParameter<int>*>(nHitsObj);
      if(!nHitsPar) return StatusCode::FAILURE;
      TH1F* tdcByMedianH= m_statsSvc->getTH1F( this->getPath(chanId)
                                               +"/tdcByMedian");
      if(!tdcByMedianH) return StatusCode::FAILURE;
      TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
      if(!adcH) return StatusCode::FAILURE;
      TH1F* occupancyH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/occupancy");
      if(!occupancyH) return StatusCode::FAILURE;
      TH1F* tdcOffsetH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/tdcOffset");
      if(!tdcOffsetH) return StatusCode::FAILURE;
      TH1F* adcMedianH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/adcMedian");
      if(!adcMedianH) return StatusCode::FAILURE;
      TH1F* adcSigmaH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                             +"/adcSigma");
      if(!adcSigmaH) return StatusCode::FAILURE;


      // ADC median and sigma
      double adcMean = 0;
      double adcMeanUncert = 0;
      double adcSigma = 0;
      double adcSigmaUncert = 0;

      double occupancytosave=0;
      double occupancyerrtosave=0;
      {
        {
          int nHits = nHitsPar->GetVal();
          if(nReadouts > 0){
            occupancytosave = nHits / ((double) nReadouts);
            occupancyerrtosave = sqrt(occupancytosave*(1-occupancytosave)*1./nReadouts);
          }
        }

        //Fit code modified by jpochoa on 03/18/11
        TF1 *nimmodelfunc = new TF1("nimmodelfunc",NIMmodel,5,300,8);
        nimmodelfunc->SetParNames(  "N","Q_{0}","#sigma_{0}","Q_{1}","#sigma_{1}","w", "#alpha","#mu");
        double AreaGuess, Q1MeanGuess, Q1SigmaGuess;
        float q0=0;//for a pedestal substracted distribution
        float sigma0=0;//for a pedestal substracted distribution
        AreaGuess = adcH->GetEntries();
        Q1MeanGuess = adcH->GetMean() - q0;
        Q1SigmaGuess = adcH->GetRMS();
        nimmodelfunc->SetParameters(AreaGuess, 90,1.9,Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03,    0.1);
        //fix four parameters
        nimmodelfunc->FixParameter(1,q0);
        nimmodelfunc->FixParameter(2,sigma0);
        nimmodelfunc->FixParameter(5,0.);
        nimmodelfunc->FixParameter(6,0.);
       
        adcH->Fit(nimmodelfunc,"R");
        if(nimmodelfunc!=NULL){
          
          adcMean=nimmodelfunc->GetParameter(3);
          adcMeanUncert =nimmodelfunc->GetParError(3);
          adcSigma = nimmodelfunc->GetParameter(4);
          adcSigmaUncert = nimmodelfunc->GetParError(3);
                  
          if(fabs(adcH->GetMean()-nimmodelfunc->GetParameter(3)) > adcH->GetRMS()){
            warning() << "Issues with fit on ring, column: " << ring << " / " 
                      << column << " mean/peak: " << adcMean << " / " << adcH->GetMean() 
                      << endreq;
            gainOk = false;
          }//print warning if issues
          
          
        } else {

          adcMean=0;
          adcMeanUncert=0;
          adcSigma = 0;
          adcSigmaUncert = 0;

          warning() << "PMT full model fit error on ring/column: " << ring << "," << column << endreq;
          gainOk = false;
        }
               
        adcMedianH->SetBinContent(adcMedianH->FindBin(column),adcMean);
        adcMedianH->SetBinError(adcMedianH->FindBin(column),adcMeanUncert);
        adcSigmaH->SetBinContent(adcSigmaH->FindBin(column),adcSigma);
        adcSigmaH->SetBinError(adcSigmaH->FindBin(column),adcSigmaUncert);

        gainChannelH->SetBinContent(ring*24+column,adcMean);
        gainChannelH->SetBinError(ring*24+column,adcMeanUncert);

        gainH2D->SetBinContent(gainH2D->GetXaxis()->FindBin(column),
                               gainH2D->GetYaxis()->FindBin(ring),
                               adcMean);
        gainH2D->SetBinError(gainH2D->GetXaxis()->FindBin(column),
                               gainH2D->GetYaxis()->FindBin(ring),
                               adcMeanUncert);
        gainH1D->Fill(adcMean);

      }


      // Channel Occupancy
      double occupancy = occupancytosave;
      double occupancyUncert = occupancyerrtosave;
      {
        occupancyH->SetBinContent(occupancyH->FindBin(column),occupancy);
        occupancyH->SetBinError(occupancyH->FindBin(column),occupancyUncert);
        if(occupancy < minValidOcc) occupancyOk = false;

        occupancyChannelH->SetBinContent(ring*24+column,occupancy);
        occupancyChannelH->SetBinError(ring*24+column,occupancyUncert);
        occupancyH2D->SetBinContent(occupancyH2D->GetXaxis()->FindBin(column),
                                    occupancyH2D->GetYaxis()->FindBin(ring),
                                    occupancy);
        occupancyH2D->SetBinError(occupancyH2D->GetXaxis()->FindBin(column),
                                    occupancyH2D->GetYaxis()->FindBin(ring),
                                    occupancyUncert);

      }
      // TDC offset
      double tdcOffset = 0;
      double tdcOffsetUncert = 0;
      {
        // Use a fit to the leading edge to determine the time offset
        int maxBin = tdcByMedianH->GetMaximumBin();
        double fitMin = tdcByMedianH->GetBinCenter(maxBin - 2);
        double fitMax = tdcByMedianH->GetBinCenter(maxBin + 10);
        tdcByMedianH->Fit("gaus","","",fitMin, fitMax);
        TF1* fitFunc = tdcByMedianH->GetFunction("gaus");
        if( fitFunc ){
          tdcOffset = fitFunc->GetParameter(1);
          tdcOffsetUncert = fitFunc->GetParError(1);
          
          if( fitFunc->GetNDF() < 3 ) 
            // Catch bad fits
            tdcOffsetUncert = 0;
          else
            // Scale uncertainty by the quality of the fit
            tdcOffsetUncert *= sqrt(fitFunc->GetChisquare() / fitFunc->GetNDF());
          
          tdcOffsetH->SetBinContent(tdcOffsetH->FindBin(column),tdcOffset);
          tdcOffsetH->SetBinError(tdcOffsetH->FindBin(column),tdcOffsetUncert);

          tdcOffsetChannelH->SetBinContent(ring*24+column,tdcOffset);
          tdcOffsetChannelH->SetBinError(ring*24+column,tdcOffsetUncert);

        }else{
          warning() << "PMT  " << pmtId << " has no TDC data for fitting." 
                    << endreq;
        }
      }



      // Record ring averages, ignoring bad channels
      if( occupancyOk && gainOk ){
        if( occupancyUncert > 0 ){
          meanOccRing[ring] += occupancy / (occupancyUncert*occupancyUncert); 
          meanOccWeight[ring] += 1.0 / (occupancyUncert*occupancyUncert); 
        }
        if( tdcOffsetUncert > 0 ){
          meanTdcRing[ring] += tdcOffset / (tdcOffsetUncert*tdcOffsetUncert); 
          meanTdcWeight[ring] += 1.0 / (tdcOffsetUncert*tdcOffsetUncert); 
        }
      }

    }

    // Record mean time offset and mean occupancy by ring
    for(int ring = 1; ring <= nRings; ring++){
      if( meanOccWeight[ring] > 0 ){
        meanOccupancyH->SetBinContent(ring,
                                      meanOccRing[ring] / meanOccWeight[ring]);
        meanOccupancyH->SetBinError(ring, 
                                    sqrt( 1.0 / meanOccWeight[ring] ));
      }
      if( meanTdcWeight[ring] > 0 ){
        meanTdcOffsetH->SetBinContent(ring,
                                      meanTdcRing[ring] / meanTdcWeight[ring]);
        meanTdcOffsetH->SetBinError(ring, 
                                    sqrt( 1.0 / meanTdcWeight[ring] ));
      }
    } // Loop over rings
 
  } // Loop over detectors
 
  return StatusCode::SUCCESS;
}

bool PmtCalibFullModel::hasStats(const DayaBay::Detector& detector){
  // Check if statistics have been initialized for this detector
  return (std::find(m_processedDetectors.begin(), 
                    m_processedDetectors.end(), 
                    detector) 
          != m_processedDetectors.end());
}

StatusCode PmtCalibFullModel::prepareStats(const Context& context){
  // Create the histograms and parameters for this detector's statistics
  DayaBay::Detector detector(context.GetSite(), context.GetDetId());

  // Site
  {
    std::string name = "site";
    std::ostringstream path;
    path << this->getPath(detector) << "/" << name;
    TParameter<int>* par = new TParameter<int>(name.c_str(), 
                                               context.GetSite());
    if( m_statsSvc->put(path.str(),par).isFailure() ) {
      error() << "prepareStats(): Could not register " << name 
              << " at " << path << endreq;
      delete par;
      par = 0;
      return StatusCode::FAILURE;
    }
  }

  // Detector ID
  {
    std::string name = "detectorId";
    std::ostringstream path;
    path << this->getPath(detector) << "/" << name;
    TParameter<int>* par = new TParameter<int>(name.c_str(), 
                                               context.GetDetId());
    if( m_statsSvc->put(path.str(),par).isFailure() ) {
      error() << "prepareStats(): Could not register " << name 
              << " at " << path << endreq;
      delete par;
      par = 0;
      return StatusCode::FAILURE;
    }
  }

  // SimFlag
  {
    std::string name = "simFlag";
    std::ostringstream path;
    path << this->getPath(detector) << "/" << name;
    TParameter<int>* par = new TParameter<int>(name.c_str(), 
                                               context.GetSimFlag());
    if( m_statsSvc->put(path.str(),par).isFailure() ) {
      error() << "prepareStats(): Could not register " << name 
              << " at " << path << endreq;
      delete par;
      par = 0;
      return StatusCode::FAILURE;
    }
  }

  // Timestamp: seconds
  {
    std::string name = "timeSec";
    std::ostringstream path;
    path << this->getPath(detector) << "/" << name;
    TParameter<int>* par = new TParameter<int>(name.c_str(), 
                                               context.GetTimeStamp().GetSec());
    if( m_statsSvc->put(path.str(),par).isFailure() ) {
      error() << "prepareStats(): Could not register " << name 
              << " at " << path << endreq;
      delete par;
      par = 0;
      return StatusCode::FAILURE;
    }
  }

  // Timestamp: nanoseconds
  {
    std::string name = "timeNanoSec";
    std::ostringstream path;
    path << this->getPath(detector) << "/" << name;
    TParameter<int>* par = 
      new TParameter<int>(name.c_str(), context.GetTimeStamp().GetNanoSec());
    if( m_statsSvc->put(path.str(),par).isFailure() ) {
      error() << "prepareStats(): Could not register " << name 
              << " at " << path << endreq;
      delete par;
      par = 0;
      return StatusCode::FAILURE;
    }
  }

  // Number of Readouts
  {
    std::string name = "nReadouts";
    std::ostringstream path;
    path << this->getPath(detector) << "/" << name;
    TParameter<int>* par = new TParameter<int>(name.c_str(), 0);
    if( m_statsSvc->put(path.str(),par).isFailure() ) {
      error() << "prepareStats(): Could not register " << name 
              << " at " << path << endreq;
      delete par;
      par = 0;
      return StatusCode::FAILURE;
    }
  }

  // ADC Sum spectrum
  {
    std::ostringstream title, path;
    std::string name = "adcSum";
    title << "ADC Sum for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* adcSum = new TH1F(name.c_str(),title.str().c_str(),
                            2000,0,20000);
    adcSum->GetXaxis()->SetTitle("ADC Sum value");
    adcSum->GetYaxis()->SetTitle("Entries");
    // DEBUG
    info() << "name= " << adcSum->GetName() 
           << " at path = \"" << path.str() << "\"" << endreq;
    if( m_statsSvc->put(path.str(),adcSum).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete adcSum;
      return StatusCode::FAILURE;
    }
  }

  // channel Gain
  {
    std::ostringstream title, path;
    std::string name = "gainChannel";
    title << "gainChannel " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* gainChannel = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    gainChannel->GetXaxis()->SetTitle("channel ID");
    gainChannel->GetYaxis()->SetTitle("gain in ADC bin");
    // DEBUG
    info() << "name= " << gainChannel->GetName() 
           << " at path = \"" << path.str() << "\"" << endreq;
    if( m_statsSvc->put(path.str(),gainChannel).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete gainChannel;
      return StatusCode::FAILURE;
    }
  }

  // channel tdc offset
  {
    std::ostringstream title, path;
    std::string name = "tdcOffsetChannel";
    title << "tdcOffsetChannel " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* tdcOffsetChannel = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    tdcOffsetChannel->GetXaxis()->SetTitle("channel ID");
    tdcOffsetChannel->GetYaxis()->SetTitle("tdc offset in TDC bin");
    // DEBUG
    info() << "name= " << tdcOffsetChannel->GetName() 
           << " at path = \"" << path.str() << "\"" << endreq;
    if( m_statsSvc->put(path.str(),tdcOffsetChannel).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete tdcOffsetChannel;
      return StatusCode::FAILURE;
    }
  }

  // channel occupancy
  {
    std::ostringstream title, path;
    std::string name = "occupancyChannel";
    title << "occupancyChannel " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* occupancyChannel = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    occupancyChannel->GetXaxis()->SetTitle("channel ID");
    occupancyChannel->GetYaxis()->SetTitle("occupancy in mean pe");
    // DEBUG
    info() << "name= " << occupancyChannel->GetName() 
           << " at path = \"" << path.str() << "\"" << endreq;
    if( m_statsSvc->put(path.str(),occupancyChannel).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete occupancyChannel;
      return StatusCode::FAILURE;
    }
  }


  // TDC Median spectrum
  {
    std::ostringstream title, path;
    std::string name = "tdcMedian";
    title << "Median TDC for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* tdcMedian = new TH1F(name.c_str(),title.str().c_str(),
                               300,0,300);
    tdcMedian->GetXaxis()->SetTitle("TDC value");
    tdcMedian->GetYaxis()->SetTitle("Entries");
    if( m_statsSvc->put(path.str(),tdcMedian).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete tdcMedian;
      return StatusCode::FAILURE;
    }
  }
  // Occupancy for each channel (jpochoa)
    {
      std::ostringstream title, path;
      std::string name = "occupancy";
      title << "Occupancy per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* occupancy = new TH2F(name.c_str(),title.str().c_str(),49,0.25,24.75,19,-0.75,8.75);
      occupancy->GetXaxis()->SetTitle("Column");
      occupancy->GetYaxis()->SetTitle("Ring");
      if( m_statsSvc->put(path.str(),occupancy).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete occupancy;
        return StatusCode::FAILURE;
      }
    }
  // Gain for each channel (jpochoa)
  {
      std::ostringstream title, path;
      std::string name = "gain2D";
      title << "Gain per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* gain2D = new TH2F(name.c_str(),title.str().c_str(),49,0.25,24.75,19,-0.75,8.75);
      gain2D->GetXaxis()->SetTitle("Column");
      gain2D->GetYaxis()->SetTitle("Ring");
      if( m_statsSvc->put(path.str(),gain2D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete gain2D;
        return StatusCode::FAILURE;
      }
    }
    // Gain for all channels (jpochoa)
  {
      std::ostringstream title, path;
      std::string name = "gain1D";
      title << "Gain for all channels " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH1F* gain1D = new TH1F(name.c_str(),title.str().c_str(),50,20,60);
      gain1D->GetXaxis()->SetTitle("ADC/PE");
      gain1D->GetYaxis()->SetTitle("Number of entries");
      if( m_statsSvc->put(path.str(),gain1D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete gain1D;
        return StatusCode::FAILURE;
      }
    }
    

  // Mean TDC Offset by AD ring
  {
    std::ostringstream title, path;
    std::string name = "meanTdcOffset";
    title << "Mean TDC offset by ring in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* meanTdcOffset = new TH1F(name.c_str(),title.str().c_str(),
                                   8,1,9);
    meanTdcOffset->GetXaxis()->SetTitle("AD Ring");
    meanTdcOffset->GetYaxis()->SetTitle("Mean TDC Offset");
    if( m_statsSvc->put(path.str(),meanTdcOffset).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete meanTdcOffset;
      return StatusCode::FAILURE;
    }
  }
  // Mean Occupancy by AD ring
  {
    std::ostringstream title, path;
    std::string name = "meanOccupancy";
    title << "Mean occupancy by ring in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* meanOccupancy = new TH1F(name.c_str(),title.str().c_str(),
                                   8,1,9);
    meanOccupancy->GetXaxis()->SetTitle("AD Ring");
    meanOccupancy->GetYaxis()->SetTitle("Mean Occupancy");
    if( m_statsSvc->put(path.str(),meanOccupancy).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete meanOccupancy;
      return StatusCode::FAILURE;
    }
  }

  // Make sure summary histograms for each ring exist in detector statistics
  for(int ring = 0; ring <= 8; ring++){
    // Occupancy by ring
    {
      std::ostringstream title, path;
      std::string name = "occupancy";
      title << "Occupancy for PMTs in " << detector.detName() 
            << " ring " << ring; 
      path << this->getPath(detector, ring) << "/" << name;
      TH1F* occupancy = new TH1F(name.c_str(),title.str().c_str(),
                                 24,1,25);
      occupancy->GetXaxis()->SetTitle("Column");
      occupancy->GetYaxis()->SetTitle("Occupancy");
      if( m_statsSvc->put(path.str(),occupancy).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete occupancy;
        return StatusCode::FAILURE;
      }
    }
    // TDC offset by ring
    {
      std::ostringstream title, path;
      std::string name = "tdcOffset";
      title << "Time Offset for PMTs in " << detector.detName() 
            << " ring " << ring; 
      path << this->getPath(detector, ring) << "/" << name;
      TH1F* tdcOffset = new TH1F(name.c_str(),title.str().c_str(),
                           24,1,25);
      tdcOffset->GetXaxis()->SetTitle("Column");
      tdcOffset->GetYaxis()->SetTitle("Time Offset [tdc]");
      if( m_statsSvc->put(path.str(),tdcOffset).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete tdcOffset;
        return StatusCode::FAILURE;
      }
    }
    // ADC Median by ring
    {
      std::ostringstream title, path;
      std::string name = "adcMedian";
      title << "ADC Median for PMTs in " << detector.detName() 
            << " ring " << ring; 
      path << this->getPath(detector, ring) << "/" << name;
      TH1F* adcMedian = new TH1F(name.c_str(),title.str().c_str(),
                                 24,1,25);
      adcMedian->GetXaxis()->SetTitle("Column");
      adcMedian->GetYaxis()->SetTitle("Median ADC");
      if( m_statsSvc->put(path.str(),adcMedian).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcMedian;
        return StatusCode::FAILURE;
      }
    }
    // ADC Sigma by ring
    {
      std::ostringstream title, path;
      std::string name = "adcSigma";
      title << "ADC Sigma for PMTs in " << detector.detName() 
            << " ring " << ring;
      path << this->getPath(detector, ring) << "/" << name;
      TH1F* adcSigma = new TH1F(name.c_str(),title.str().c_str(),
                                24,1,25);
      adcSigma->GetXaxis()->SetTitle("Column");
      adcSigma->GetYaxis()->SetTitle("ADC Sigma");
      if( m_statsSvc->put(path.str(),adcSigma).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcSigma;
        return StatusCode::FAILURE;
      }
    }
  }

  ServiceMode svcMode(context, 0);

  // Get list of all FEE channels
  const std::vector<DayaBay::FeeChannelId>& channelList 
    = m_cableSvc->feeChannelIds( svcMode );
  std::vector<DayaBay::FeeChannelId>::const_iterator chanIter, 
    chanEnd = channelList.end();
  // Initialize statistics for each channel
  for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++){
    DayaBay::FeeChannelId chanId = *chanIter;

    // Board Number
    {
      std::string name = "board";
      std::ostringstream path;
      path << this->getPath(chanId) << "/" << name;
      TParameter<int>* par = new TParameter<int>(name.c_str(), chanId.board());
      if( m_statsSvc->put(path.str(),par).isFailure() ) {
        error() << "prepareStats(): Could not register " << name 
                << " at " << path << endreq;
        delete par;
        par = 0;
        return StatusCode::FAILURE;
      }
    }

    // Connector Number
    {
      std::string name = "connector";
      std::ostringstream path;
      path << this->getPath(chanId) << "/" << name;
      TParameter<int>* par = new TParameter<int>(name.c_str(), 
                                                 chanId.connector());
      if( m_statsSvc->put(path.str(),par).isFailure() ) {
        error() << "prepareStats(): Could not register " << name 
                << " at " << path << endreq;
        delete par;
        par = 0;
        return StatusCode::FAILURE;
      }
    }

    // Number of Hits
    {
      std::string name = "nHits";
      std::ostringstream path;
      path << this->getPath(chanId) << "/" << name;
      TParameter<int>* par = new TParameter<int>(name.c_str(), 0);
      if( m_statsSvc->put(path.str(),par).isFailure() ) {
        error() << "prepareStats(): Could not register " << name 
                << " at " << path << endreq;
        delete par;
        par = 0;
        return StatusCode::FAILURE;
      }
    }

    // Raw TDC spectrum
    {
      std::string name = "tdcRaw";
      std::ostringstream title, path;
      title << "Raw TDC values for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* tdcRaw = new TH1F(name.c_str(),title.str().c_str(),
                              300,0,300);
      tdcRaw->GetXaxis()->SetTitle("TDC value");
      tdcRaw->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),tdcRaw).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete tdcRaw;
        return StatusCode::FAILURE;
      }
    }

    // TDC spectrum, median corrected
    {
      std::string name = "tdcByMedian";
      std::ostringstream title, path;
      title << "Median-corrected TDC values for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* tdcByMedian = new TH1F(name.c_str(),title.str().c_str(),
                                   600,-300,300);
      tdcByMedian->GetXaxis()->SetTitle("TDC value");
      tdcByMedian->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),tdcByMedian).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete tdcByMedian;
        return StatusCode::FAILURE;
      }
    }

    // Raw ADC spectrum
    {
      std::ostringstream title, path;
      std::string name = "adcRaw";
      title << "ADC values for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* adcRaw = new TH1F(name.c_str(),title.str().c_str(),
                              4096,0,4096);
      adcRaw->GetXaxis()->SetTitle("ADC value");
      adcRaw->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),adcRaw).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcRaw;
        return StatusCode::FAILURE;
      }
    }
    
    // ADC spectrum, baseline subtracted
    {
      std::ostringstream title, path;
      std::string name = "adc";
      title << "Baseline-subtracted ADC values for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* adc = new TH1F(name.c_str(),title.str().c_str(),
                           5096,-1000,4096);
      adc->GetXaxis()->SetTitle("ADC value");
      adc->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),adc).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adc;
        return StatusCode::FAILURE;
      }
    }

    // ADC spectrum by Clock Cycle
    {
      std::ostringstream title, path;
      std::string name = "adcByClock";
      title << "ADC values by clock cycle for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* adcByClock = new TH1F(name.c_str(),title.str().c_str(),
                                  20,0,20);
      adcByClock->GetXaxis()->SetTitle("ADC Clock Cycle");
      adcByClock->GetYaxis()->SetTitle("ADC");
      if( m_statsSvc->put(path.str(),adcByClock).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcByClock;
        return StatusCode::FAILURE;
      }
    }
  }  

  m_processedDetectors.push_back(detector);

  return StatusCode::SUCCESS;
}

std::string PmtCalibFullModel::getPath( const DayaBay::Detector& detector ){
  return m_filepath + "/" + detector.detName();
}

std::string PmtCalibFullModel::getPath(const DayaBay::Detector& detector,
                                         int ring){
  std::ostringstream path;
  path << m_filepath << "/" << detector.detName() << "/ring_" << ring; 
  return path.str();
}

std::string PmtCalibFullModel::getPath(const DayaBay::FeeChannelId& chanId){
  std::ostringstream path;
  path << m_filepath << "/" << chanId.detName() 
       << "/chan_" << chanId.board() << "_" << chanId.connector(); 
  return path.str();
}

---