In This Package:

PmtCalibLeadingEdge.cc

Go to the documentation of this file.

#include "PmtCalibLeadingEdge.h"

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

PmtCalibLeadingEdge::PmtCalibLeadingEdge(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=false,
                  "Use FloatFeePedes or not? In case of false, it will use CalibSvc.");
  declareProperty("FilePath",m_filepath="/file0/pmtCalibLeadingEdge",
                  "File path of registered histograms.");

  declareProperty("TextFile", m_textFileName="Sum.txt", "Summary of fitting result in text file");

}

PmtCalibLeadingEdge::~PmtCalibLeadingEdge(){}

StatusCode PmtCalibLeadingEdge::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;
  }

  if(m_useFloatFeePedes){
    // 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);
  }

  m_textFile.open( m_textFileName.c_str(), ios_base::out );

  m_textFile
    <<"run/I:board/I:channel/I:preAdc/D:preAdcSigma/D:rawAdc/D:rawAdcSigma/D:expAdc/D:expAdcSigma/D:quasiAdc/D:quasiAdcSigma/D"<<endl;

  return StatusCode::SUCCESS;
}

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

  m_textFile.close();

  return StatusCode::SUCCESS;
}


StatusCode PmtCalibLeadingEdge::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* tdcMeanH = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMean");
  if(!tdcMeanH) return StatusCode::FAILURE;
  TH1F* NChannelH = m_statsSvc->getTH1F( this->getPath(detector) + "/NChannel");
  if(!NChannelH) return StatusCode::FAILURE;
  TH2F* occupancyH = m_statsSvc->getTH2F(this->getPath(detector) + "/occupancy");
  if(!occupancyH) return StatusCode::FAILURE;

  // Add the current readout to the calibration statistics
  std::vector<int> readoutTdc;
  double adcSum = 0;
  // Loop over channels
  int NChannel = 0;

  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* preAdcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/preAdc");
    if(!preAdcH) return StatusCode::FAILURE;
    TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");    
    if(!adcH) return StatusCode::FAILURE;
    TH1F* quasiAdcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/quasiAdc");
    if(!quasiAdcH) return StatusCode::FAILURE;
    TH1F* adcByClockH = m_statsSvc->getTH1F( this->getPath(chanId) 
                                             + "/adcByClock");
    if(!adcByClockH) return StatusCode::FAILURE;
    
    // Loop over each hit
    bool FindFirstHit = false;
    for(unsigned int i=0; i<channel.hitCount(); i++) {
      int adcClock = channel.peakCycle(i);
      int adc = channel.adc(i);
      float preAdc = channel.preAdcAvg(i);
      int tdc = channel.tdc(i);

      // Only hit with good timing can join the fit and only the first one.
      if( FindFirstHit ) break;
      if( !GoodTdc(tdc) ) continue;
      FindFirstHit = true;

      // Loop over tdc values
      readoutTdc.push_back(tdc);
      tdcRawH->Fill(tdc);
    
      // Loop over adc values
      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 enough
        adcBaseline = channel.preAdcAvg(i);  // Use preAdc is much better
      }

      if( adc>0 && channel.isHighGainAdc(i) ) {
        adcRawH->Fill(adc);
        adcH->Fill(adc-adcBaseline);
        quasiAdcH->Fill(adc-preAdc);
        adcByClockH->Fill(adcClock,adc-adcBaseline);
        preAdcH->Fill(preAdc);
      }
      adcSum += adc-adcBaseline;
    }
    // Increment number of hits on this channel
    if( FindFirstHit ) {
      nHitsPar->SetVal( nHitsPar->GetVal() + 1 );
     }
  }
  NChannel = readoutTdc.size();
  if( NChannel<=0 ) {
    info()<<" Too few hits for this event. Skipped. "<<endreq;
    return StatusCode::SUCCESS;
  }

  if(readoutTdc.size() > 0){
    // Find the mean TDC
    double meanTdc = accumulate(readoutTdc.begin(), readoutTdc.end(),0);
    meanTdc = meanTdc/NChannel;

    for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
      const DayaBay::DaqPmtChannel& pmtChan = *(*channelIter);
      const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
      TH1F* tdcByMeanH = 
        m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMean");
      if(!tdcByMeanH) return StatusCode::FAILURE;
      // Loop over tdc values
      bool FirstHit = false;
      for(unsigned int i=0; i<pmtChan.hitCount(); i++) {
        int tdc = pmtChan.tdc(i);
        // Only hit with good timing can join the fit and only the first one.
        if( FirstHit ) break;
        if( !GoodTdc(tdc) ) continue;
        FirstHit = true;
        tdcByMeanH->Fill(tdc-meanTdc);
      }
    }
    tdcMeanH->Fill(meanTdc);
  }

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

StatusCode PmtCalibLeadingEdge::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;

    if(!detector.isAD()) continue;//for now this code can only tolerate ADs (jpochoa)    

    // 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* adcRawMeanH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/adcRawMean");
    if(!adcRawMeanH) return StatusCode::FAILURE;
    TH1F* adcRawSigmaH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/adcRawSigma");
    if(!adcRawSigmaH) return StatusCode::FAILURE;
    //
    TH1F* quasiAdcMeanH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/quasiAdcMean");
    if(!quasiAdcMeanH) return StatusCode::FAILURE;
    TH1F* quasiAdcSigmaH = m_statsSvc->getTH1F( this->getPath(detector)
                                                 +"/quasiAdcSigma");
    if(!quasiAdcSigmaH) return StatusCode::FAILURE;
    //
    TH1F* expAdcMeanH = m_statsSvc->getTH1F( this->getPath(detector)
                                                  +"/expAdcMean");
    if(!expAdcMeanH) return StatusCode::FAILURE;
    TH1F* expAdcSigmaH = m_statsSvc->getTH1F( this->getPath(detector)
                                                   +"/expAdcSigma");
    if(!expAdcSigmaH) return StatusCode::FAILURE;
    //
    TH1F* preAdcMeanH = m_statsSvc->getTH1F( this->getPath(detector)
                                             +"/preAdcMean");
    if(!preAdcMeanH) return StatusCode::FAILURE;
    TH1F* preAdcSigmaH = m_statsSvc->getTH1F( this->getPath(detector)
                                              +"/preAdcSigma");
    if(!preAdcSigmaH) 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 (not from data, but just from cableSvc)
    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();

      m_textFile<<-1<<" "<<chanId.board()<<" "<<chanId.connector()<<" ";

      // 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* tdcByMeanH= m_statsSvc->getTH1F( this->getPath(chanId)
                                               +"/tdcByMean");
      if(!tdcByMeanH) return StatusCode::FAILURE;
      TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
      if(!adcH) return StatusCode::FAILURE;
      TH1F* adcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRaw");
      if(!adcRawH) return StatusCode::FAILURE;
      TH1F* preAdcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/preAdc");
      if(!preAdcH) return StatusCode::FAILURE;
      TH1F* quasiAdcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/quasiAdc");
      if(!quasiAdcH) 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* adcMeanH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/adcMean");
      if(!adcMeanH) return StatusCode::FAILURE;
      TH1F* adcSigmaH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                             +"/adcSigma");
      if(!adcSigmaH) return StatusCode::FAILURE;

      // Channel Occupancy
      double occupancy = 0;
      double occupancyUncert = 0;
      {
        int nHits = nHitsPar->GetVal();
        if(nReadouts > 0){
          occupancy = nHits / ((double) nReadouts);
          occupancyUncert = 
            sqrt( occupancy*(1-occupancy)/((double) nReadouts) );
        }
        occupancyH->SetBinContent(occupancyH->FindBin(column),occupancy);
        occupancyH->SetBinError(occupancyH->FindBin(column),occupancyUncert);
        occupancyH2D->SetBinContent(occupancyH2D->GetXaxis()->FindBin(column),
                                    occupancyH2D->GetYaxis()->FindBin(ring),
                                    occupancy);
        occupancyH2D->SetBinError(occupancyH2D->GetXaxis()->FindBin(column),
                                    occupancyH2D->GetYaxis()->FindBin(ring),
                                    occupancyUncert);
        
        if(occupancy < minValidOcc) occupancyOk = false;
      }
      // TDC offset
      double tdcOffset = 0;
      double tdcOffsetUncert = 0;
      {
        // Use a fit to the leading edge to determine the time offset
        int maxbin = tdcByMeanH->GetMaximumBin();
        double fitMin = tdcByMeanH->GetBinCenter(maxbin) - 3; //tdcByMeanH->GetRMS(1);
        double fitMax = tdcByMeanH->GetBinCenter(maxbin) + 3; //tdcByMeanH->GetRMS(1);
        tdcByMeanH->Fit("gaus","","",fitMin, fitMax);
        TF1* fitFunc = tdcByMeanH->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);
        }else{
          warning() << "PMT  " << pmtId << " has no TDC data for fitting." 
                    << endreq;
        }
      }

      // PreAdc
      double preAdcArea;
      double preAdcMean;
      double preAdcSigm;
      {
        preAdcArea = preAdcH->GetEntries();
        preAdcMean = preAdcH->GetMean(1);
        preAdcSigm = preAdcH->GetRMS(1);
        TF1 preAdcFit("preAdcFit","gaus");
        preAdcFit.SetParameter( 0, preAdcArea );
        preAdcFit.SetParameter( 1, preAdcMean );
        preAdcFit.SetParameter( 2, preAdcSigm );
        preAdcH->Fit( &preAdcFit, "", "", 
                      preAdcMean-2*preAdcSigm, preAdcMean+2*preAdcSigm );
        preAdcMean = preAdcFit.GetParameter(1);
        preAdcSigm = preAdcFit.GetParameter(2);
        //
        if(chanId.board()<=16) {
          preAdcMeanH->SetBinContent(
                                     chanId.board()*16+chanId.connector(),preAdcMean);
          preAdcSigmaH->SetBinContent(
                                      chanId.board()*16+chanId.connector(),preAdcSigm);
        }
        
        m_textFile<<preAdcMean<<" "<<preAdcSigm<<" ";
      }

      // Raw Adc
      double adcRawArea;
      double adcRawMean;
      double adcRawSigm;
      {
        adcRawArea = adcRawH->GetEntries();
        adcRawMean = adcRawH->GetMean(1);
        adcRawSigm = adcRawH->GetRMS(1);
        TF1 adcRawFit("adcRawFit","gaus");
        adcRawFit.SetParameter( 0, adcRawArea );
        adcRawFit.SetParameter( 1, adcRawMean );
        adcRawFit.SetParameter( 2, adcRawSigm );
        adcRawH->Fit( &adcRawFit, "", "",
                      adcRawMean-1.0*adcRawSigm, adcRawMean+1.0*adcRawSigm );
        adcRawMean = adcRawFit.GetParameter(1);
        adcRawSigm = adcRawFit.GetParameter(2);
        
        if(chanId.board()<=16) {
          adcRawMeanH->SetBinContent(
                                     chanId.board()*16+chanId.connector(),adcRawMean);
          adcRawSigmaH->SetBinContent(
                                      chanId.board()*16+chanId.connector(),adcRawSigm);
        }

        m_textFile<<adcRawMean<<" "<<adcRawSigm<<" ";
      }

      // ADC (Adc -Ave Ped)
      double adcArea = 0;
      double adcMean = 0;
      double adcMeanUncert = 0;
      double adcSigma = 0;
      double adcSigmaUncert = 0;
      {
        // Find ADC SPE peak, width with simple gaussian
        adcArea = adcH->GetEntries();
        adcMean = adcH->GetMean(1);
        adcSigma= adcH->GetRMS(1);
        TF1 adcFit("adcFit","gaus");
        adcFit.SetParameter( 0, adcArea );
        adcFit.SetParameter( 1, adcMean );
        adcFit.SetParameter( 2, adcSigma);
        adcH->Fit( &adcFit,"","",
                   adcMean-1.0*adcSigma, adcMean+1.0*adcSigma); 

        adcMean = adcFit.GetParameter(1);
        adcMeanUncert = adcFit.GetParError(1);
        adcSigma = adcFit.GetParameter(2);
        adcSigmaUncert = adcFit.GetParError(2);
        gainOk=true;

        adcMeanH->SetBinContent(adcMeanH->FindBin(column),adcMean);
        adcMeanH->SetBinError(adcMeanH->FindBin(column),adcMeanUncert);
        adcSigmaH->SetBinContent(adcSigmaH->FindBin(column),adcSigma);
        adcSigmaH->SetBinError(adcSigmaH->FindBin(column),adcSigmaUncert);
        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);

        if(chanId.board()<=16) {
          expAdcMeanH->SetBinContent(
                                     chanId.board()*16+chanId.connector(),adcMean);
          expAdcSigmaH->SetBinContent(
                                      chanId.board()*16+chanId.connector(),adcSigma);
        }

        m_textFile<<adcMean<<" "<<adcSigma<<" ";
      }

      // Quasi Adc
      double quasiAdcArea;
      double quasiAdcMean;
      double quasiAdcSigm;
      {
        quasiAdcArea = quasiAdcH->GetEntries();
        quasiAdcMean = quasiAdcH->GetMean(1);
        quasiAdcSigm = quasiAdcH->GetRMS(1);
        TF1 quasiAdcFit("quasiAdcFit","gaus");
        quasiAdcFit.SetParameter( 0, quasiAdcArea );
        quasiAdcFit.SetParameter( 1, quasiAdcMean );
        quasiAdcFit.SetParameter( 2, quasiAdcSigm );
        quasiAdcH->Fit( &quasiAdcFit, "", "",
                        quasiAdcMean-1.0*quasiAdcSigm, quasiAdcMean+1.0*quasiAdcSigm );
        quasiAdcMean = quasiAdcFit.GetParameter(1);
        quasiAdcSigm = quasiAdcFit.GetParameter(2);

        if(chanId.board()<=16) {
          quasiAdcMeanH->SetBinContent(
                                       chanId.board()*16+chanId.connector(),quasiAdcMean);
          quasiAdcSigmaH->SetBinContent(
                                        chanId.board()*16+chanId.connector(),quasiAdcSigm);
        }

        m_textFile<<quasiAdcMean<<" "<<quasiAdcSigm<<endl;
      }


      // 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 PmtCalibLeadingEdge::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());
}

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

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

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

---