/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes
PmtCalibLeadingEdge Class Reference

#include <PmtCalibLeadingEdge.h>

Inheritance diagram for PmtCalibLeadingEdge:
Inheritance graph
[legend]
Collaboration diagram for PmtCalibLeadingEdge:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 PmtCalibLeadingEdge (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~PmtCalibLeadingEdge ()
virtual StatusCode initialize ()
virtual StatusCode finalize ()
virtual StatusCode process (const DayaBay::ReadoutHeader &readout)
 Implementation of IPmtCalibParamTool ////////////////////////////.
virtual StatusCode calibrate ()
 calibrate() is called after processing many readouts.

Static Public Member Functions

static const InterfaceID & interfaceID ()
 Retrieve interface ID.

Private Member Functions

bool hasStats (const DayaBay::Detector &detector)
 Check if the statistics for this detector have been initialized.
StatusCode prepareStats (const Context &context)
 Get the current statistics for a detector, or create if needed.
std::string getPath (const DayaBay::Detector &detector)
std::string getPath (const DayaBay::Detector &detector, int ring)
std::string getPath (const DayaBay::FeeChannelId &channelId)
bool GoodTdc (int Tdc)

Private Attributes

std::string m_cableSvcName
std::string m_calibSvcName
std::string m_floatFeePedesSvcName
std::string m_filepath
std::string m_pmtGeomSvcName
IPmtGeomInfoSvcm_pmtGeomSvc
ICableSvcm_cableSvc
ICalibDataSvcm_calibSvc
IStatisticsSvcm_statsSvc
IFloatingFeePedestalSvcm_floatFeePedesSvc
bool m_useFloatFeePedes
std::vector< DayaBay::Detectorm_processedDetectors
ofstream m_textFile
string m_textFileName
int m_tdclow
int m_tdchigh

Detailed Description

Definition at line 74 of file PmtCalibLeadingEdge.h.


Constructor & Destructor Documentation

PmtCalibLeadingEdge::PmtCalibLeadingEdge ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 29 of file PmtCalibLeadingEdge.cc.

  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_statsSvc(0)
{
  
  declareInterface< IPmtCalibParamTool >(this) ;   
  declareProperty("CableSvcName",m_cableSvcName="CableSvc",
                  "Name of service to map between detector, hardware, and electronic IDs");
  declareProperty("CalibSvcName",m_calibSvcName="CalibDataSvc",
                  "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("TDCCutLow",m_tdclow=950,"low tdc cut for good hits");
  declareProperty("TDCCutHigh",m_tdchigh=1150,"high tdc cut for good hits");
  declareProperty("TextFile", m_textFileName="Sum.txt", "Summary of fitting result in text file");
        declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc", 
      "Name of Pmt Geometry Information Service");

}
PmtCalibLeadingEdge::~PmtCalibLeadingEdge ( ) [virtual]

Definition at line 56 of file PmtCalibLeadingEdge.cc.

{}

Member Function Documentation

StatusCode PmtCalibLeadingEdge::initialize ( ) [virtual]

Definition at line 58 of file PmtCalibLeadingEdge.cc.

{
  info() << "initialize()" << endreq;
  StatusCode sc = this->GaudiTool::initialize();
  if( sc != StatusCode::SUCCESS ) return sc;
  
  //Get PMT geometry service
  m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);
   if(!m_pmtGeomSvc ){
    error() << "Failed to access PMT geometry svc: " << m_pmtGeomSvcName << endreq;
    return StatusCode::FAILURE;
  }
  // 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 ( ) [virtual]

Definition at line 125 of file PmtCalibLeadingEdge.cc.

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

  m_textFile.close();

  return StatusCode::SUCCESS;
}
StatusCode PmtCalibLeadingEdge::process ( const DayaBay::ReadoutHeader readout) [virtual]

Implementation of IPmtCalibParamTool ////////////////////////////.

process(const DayaBay::ReadoutPmtCrate&) process a single readout, extracting the information needed for calibration

Restrict to fine range gain measurement

Implements IPmtCalibParamTool.

Definition at line 137 of file PmtCalibLeadingEdge.cc.

                                                                                {

  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();
  //DetectorId::DetectorId_t detId = detector.detectorId();
  //cout<<detId<<endl; 
 // 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;
  }
    // Calib trigger only
  DayaBay::Trigger::TriggerType_t trigType = daqCrate->triggerType();
  if( (trigType & DayaBay::Trigger::kCalib) != DayaBay::Trigger::kCalib )  
    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::DetectorSensor& sensDetId 
      = m_cableSvc->sensor(chanId, svcMode);
    if (sensDetId.bogus()) {
      warning() << "Got bogus PMT: \"" << sensDetId
                << "\" using channelId: \"" << chanId
                << "\" and context \"" << context.AsString() << "\""
                << endreq;
        continue;
    }

    // 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);
                if(adcClock<3||adcClock>8) continue;   // Bad cycle
                int residue = (tdc-15)%16;
          if(residue<4) continue; // remove 40MHz Electronic noise
      // 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;
      
        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();
                const DayaBay::DetectorSensor& sensDetId 
      = m_cableSvc->sensor(chanId, svcMode);
        if (sensDetId.bogus()) {
      warning() << "Got bogus PMT: \"" << sensDetId
                << "\" using channelId: \"" << chanId
                << "\" and context \"" << context.AsString() << "\""
                << endreq;
        continue;
    }
      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 ( ) [virtual]

calibrate() is called after processing many readouts.

This method is responsible for calculating the calibration parameters.

Adc

Adc

end of channel

end of WP PMTs

only deal with AD and waterpool yuzy

Implements IPmtCalibParamTool.

Definition at line 313 of file PmtCalibLeadingEdge.cc.

                                         {
  // 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()) { //for now this code can only tolerate ADs (jpochoa)    // the following is Ok for AD until meeting WP

    // 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)
                                                +"/occupancyChannel");
    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) 
                                              + "/ADoccupancy");
    if(!occupancyH2D) return StatusCode::FAILURE;
    TH2F* gainH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/ADgain2D");
    if(!gainH2D) return StatusCode::FAILURE;
    TH1F* gainH1D = m_statsSvc->getTH1F( this->getPath(detector) 
                                              + "/gain1D");
    if(!gainH1D) return StatusCode::FAILURE;
   TH2F* Chi2NDF = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/ADChi2NDF");
    if(!Chi2NDF) return StatusCode::FAILURE;
    TH1F* Chi2NDF_1D = m_statsSvc->getTH1F( this->getPath(detector) 
                                              + "/Chi2_1D");
    if(!Chi2NDF_1D) 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;
      const DayaBay::DetectorSensor& sensDetId 
      = m_cableSvc->sensor(chanId, svcMode);
    if (sensDetId.bogus()) {
      warning() << "Got bogus PMT: \"" << sensDetId
                << "\" using channelId: \"" << chanId
                << "\" and context \"" << context.AsString() << "\""
                << endreq;
        continue;
    }
      // 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) );
        }
        meanOccupancyH->SetBinContent(chanId.board()*16+chanId.connector(),occupancy);
        meanOccupancyH->SetBinError(chanId.board()*16+chanId.connector(),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()<=18) {
          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()<=18) {
          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);
        //cout<<adcMean<<' '<<adcSigma<<endl;
        TF1 adcFit("adcFit","gaus");
        //adcFit.SetParameter( 0, adcArea );
        adcFit.SetParameter( 1, adcMean );
        adcFit.SetParameter( 2, adcSigma);
        adcH->Fit( &adcFit,"Q","",
                   10, 50);     

        adcMean = adcFit.GetParameter(1);
        adcMeanUncert = adcFit.GetParError(1);
        adcSigma = adcFit.GetParameter(2);
        adcSigmaUncert = adcFit.GetParError(2);
        
        adcFit.SetParameter( 1, adcMean );
   adcFit.SetParameter( 2, adcSigma);
        adcH->Fit( &adcFit,"Q","",
                   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;
        double Chi2 = adcFit.GetChisquare();
        double ndf = adcFit.GetNDF();
        double fitquality = Chi2/ndf;
        
        

        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);
   Chi2NDF->SetBinContent(Chi2NDF->GetXaxis()->FindBin(column),
                               Chi2NDF->GetYaxis()->FindBin(ring),
                               fitquality);
        if(chanId.board()<=18) {
          expAdcMeanH->SetBinContent(
                                     chanId.board()*16+chanId.connector(),adcMean);
         expAdcMeanH->SetBinError(chanId.board()*16+chanId.connector(),adcMeanUncert);
          expAdcSigmaH->SetBinContent(chanId.board()*16+chanId.connector(),adcSigma);
        expAdcSigmaH->SetBinError(chanId.board()*16+chanId.connector(),adcSigmaUncert);
          Chi2NDF_1D->SetBinContent(chanId.board()*16+chanId.connector(), fitquality);
        }

        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
   } //end of AD
   
   
   
   
   
   if(detector.isWaterShield()){  //deal with WP PMTs, yuzy
       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)
                                                +"/occupancyChannel");
    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;
   TH1F* gainH1D = m_statsSvc->getTH1F( this->getPath(detector) 
                                              + "/gain1D");
    if(!gainH1D) return StatusCode::FAILURE;
   TH2F* WPIoccupancyH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/WPIoccupancy");
    if(!WPIoccupancyH2D) return StatusCode::FAILURE;
    TH2F* WPOinwardoccupancyH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/WPOinwardoccupancy");
         TH2F* WPOoutwardoccupancyH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/WPOoutwardoccupancy");
    TH2F* WPIgainH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/WPIgain2D");
        TH2F* WPOinwardgainH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/WPOinwardgain2D");
        TH2F* WPOoutwardgainH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/WPOoutwardgain2D");

        TH1F* Chi2NDF_1D = m_statsSvc->getTH1F( this->getPath(detector) 
                                              + "/Chi2_1D");
    if(!Chi2NDF_1D) return StatusCode::FAILURE;
        TH2F* WPIChi2 = m_statsSvc->getTH2F( this->getPath(detector) + "/WPIChi2NDF");
        TH2F* WPOinwardChi2 = m_statsSvc->getTH2F( this->getPath(detector) + "/WPOinwardChi2NDF");
        TH2F* WPOoutwardChi2 = m_statsSvc->getTH2F( this->getPath(detector) + "/WPOoutwardChi2NDF");
       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);
        // 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;
      DayaBay::PoolPmtSensor pmtId = m_cableSvc->poolPmtSensor(chanId, svcMode);
  const DayaBay::DetectorSensor& sensDetId 
      = m_cableSvc->sensor(chanId, svcMode);
    if (sensDetId.bogus()) {
      warning() << "Got bogus PMT: \"" << sensDetId
                << "\" using channelId: \"" << chanId
                << "\" and context \"" << context.AsString() << "\""
                << endreq;
        continue;
    }
    double minValidOcc = 0.01;
        int wallnum = pmtId.wallNumber();
        int wallspot = pmtId.wallSpot();
        bool infacing = pmtId.inwardFacing();
        int site1;
        if(site ==1) site1 = 0x01;
        if(site ==2)site1 = 0x02;
        if(site==4)site1 = 0x04;
       int pmtid = (site1<<24) | (detId <<16) | (infacing<<12) | (wallnum<<8) | wallspot;
       //cout<<pmtid<<endl;
       CLHEP::Hep3Vector pos =  m_pmtGeomSvc->get(pmtid)->localPosition();
       double localx = pos.x()/1000.;
       double localy = pos.y()/1000.;
       double localz=pos.z()/1000.;
       //if(detId ==5 )cout<<"iws "<<chanId.board()<<' '<<chanId.connector()<<' '<<wallnum<<' '<<wallspot<<' '<<infacing<<' '<<localx<<' '<<localy<<' '<<localz<<endl;
      //if(detId ==6 )cout<<"ows "<<chanId.board()<<' '<<chanId.connector()<<' '<<wallnum<<' '<<wallspot<<' '<<infacing<<' '<<localx<<' '<<localy<<' '<<localz<<endl;
       if(detId==5) localz += 4.2334;
       if(detId==6)localz += 4.6537;
       if(wallnum==1) {localx += localz;}
       if(wallnum==2) {localx+=0.71*localz; localy+=localz;}
       if(wallnum==3) {localy+=localz;}
       if(wallnum==4) {localx-=0.71*localz; localy+=localz;}
       if(wallnum==5) {localx -= localz;}
       if(wallnum==6) {localx-=0.71*localz; localy-=localz;}
       if(wallnum==7) {localy-=localz;}
       if(wallnum==8) {localx+=0.71*localz; localy-=localz;}
       
       //ofstream out1("position.txt", ofstream::out || ofstream::app);
      //if(detId ==5 )cout<<"iws "<<chanId.board()<<' '<<chanId.connector()<<' '<<wallnum<<' '<<wallspot<<' '<<infacing<<' '<<localx<<' '<<localy<<endl;
      //if(detId ==6 )cout<<"ows "<<chanId.board()<<' '<<chanId.connector()<<' '<<wallnum<<' '<<wallspot<<' '<<infacing<<' '<<localx<<' '<<localy<<endl;
       //out1.close();
     // m_textFile<<-1<<" "<<chanId.board()<<" "<<chanId.connector()<<" ";
        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;
      
           // 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) );
            
        }
        meanOccupancyH->SetBinContent(chanId.board()*16+chanId.connector(),occupancy);
        meanOccupancyH->SetBinError(chanId.board()*16+chanId.connector(),occupancyUncert);
        if(detId==5)
                WPIoccupancyH2D->SetBinContent(WPIoccupancyH2D->GetXaxis()->FindBin(localx),
                                    WPIoccupancyH2D->GetYaxis()->FindBin(localy),
                                    occupancy);
        if(detId==6&&infacing==1)
                WPOinwardoccupancyH2D->SetBinContent(WPOinwardoccupancyH2D->GetXaxis()->FindBin(localx),
                                    WPOinwardoccupancyH2D->GetYaxis()->FindBin(localy),
                                    occupancy);
    if(detId==6&&infacing==0)
                WPOoutwardoccupancyH2D->SetBinContent(WPOoutwardoccupancyH2D->GetXaxis()->FindBin(localx),
                                    WPOoutwardoccupancyH2D->GetYaxis()->FindBin(localy),
                                    occupancy);
        if(occupancy < minValidOcc) occupancyOk = false;
      }
      
       // PreAdc
      double preAdcArea;
      double preAdcMean;
      double preAdcSigm;
      {
        preAdcArea = preAdcH->GetEntries();
        preAdcMean = preAdcH->GetMean(1);
        preAdcSigm = preAdcH->GetRMS(1);
        //if(detId==5) WPIgainH2D->Fill(localx,localy,preAdcSigm);
        //if(detId==6&&infacing==1) WPOinwardgainH2D->Fill(localx,localy,preAdcSigm);
        //if(detId==6&&infacing==0) WPOoutwardgainH2D->Fill(localx,localy,preAdcSigm);
        TF1 preAdcFit("preAdcFit","gaus");
        preAdcFit.SetParameter( 0, preAdcArea );
        preAdcFit.SetParameter( 1, preAdcMean );
        preAdcFit.SetParameter( 2, preAdcSigm );
        preAdcH->Fit( &preAdcFit, "Q", "", 
                      preAdcMean-2*preAdcSigm, preAdcMean+2*preAdcSigm );
        preAdcMean = preAdcFit.GetParameter(1);
        preAdcSigm = preAdcFit.GetParameter(2);
        //
        if(chanId.board()<=25) {
          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, "Q", "",
                      adcRawMean-1.0*adcRawSigm, adcRawMean+1.0*adcRawSigm );
        adcRawMean = adcRawFit.GetParameter(1);
        adcRawSigm = adcRawFit.GetParameter(2);
        
        if(chanId.board()<=25) {
          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
        TF1 adcFit("adcFit","gaus");

        adcArea = adcH->GetEntries();
        adcMean = adcH->GetMean(1);
        adcSigma= adcH->GetRMS(1);
        //cout<<adcMean<<' '<<adcSigma<<endl;
        //TF1 adcFit("adcFit","gaus");
        //adcFit.SetParameter( 0, adcArea );
        adcFit.SetParameter( 1, adcMean );
        adcFit.SetParameter( 2, adcSigma);
        adcH->Fit( &adcFit,"Q","",
                   10,50);      

        adcMean = adcFit.GetParameter(1);
        adcMeanUncert = adcFit.GetParError(1);
        adcSigma = adcFit.GetParameter(2);
        adcSigmaUncert = adcFit.GetParError(2);
        
        adcFit.SetParameter( 1, adcMean );
   adcFit.SetParameter( 2, adcSigma);
        adcH->Fit( &adcFit,"Q","",
                   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;
        double Chi2 = adcFit.GetChisquare();
        double ndf = adcFit.GetNDF();
        double fitquality = Chi2/ndf;
        //gainOk=true;


    
        gainH1D->Fill(adcMean);
        if(chanId.board()<=25) {
          expAdcMeanH->SetBinContent(
                                     chanId.board()*16+chanId.connector(),adcMean);
         expAdcMeanH->SetBinError(
                                     chanId.board()*16+chanId.connector(),adcMeanUncert);
          expAdcSigmaH->SetBinContent(
                                      chanId.board()*16+chanId.connector(),adcSigma);
        expAdcSigmaH->SetBinError(
                                      chanId.board()*16+chanId.connector(),adcSigmaUncert);
          Chi2NDF_1D->SetBinContent(chanId.board()*16+chanId.connector(), fitquality);
        }
        m_textFile<<adcMean<<" "<<adcSigma<<" ";
        if(adcMean<0) adcMean =1; // a trick to hold 2D histogram. only valid to the 2-D histograms.
        if(detId==5) {WPIgainH2D->Fill(localx,localy,adcMean); WPIChi2->Fill(localx,localy,fitquality);}
        if(detId==6&&infacing==1){ WPOinwardgainH2D->Fill(localx,localy,adcMean);WPOinwardChi2->Fill(localx,localy,fitquality);}
        if(detId==6&&infacing==0) {WPOoutwardgainH2D->Fill(localx,localy,adcMean);WPOoutwardChi2->Fill(localx,localy,fitquality);}
      }

      // 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, "Q", "",
                        quasiAdcMean-1.0*quasiAdcSigm, quasiAdcMean+1.0*quasiAdcSigm );
        quasiAdcMean = quasiAdcFit.GetParameter(1);
        quasiAdcSigm = quasiAdcFit.GetParameter(2);

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

        m_textFile<<quasiAdcMean<<" "<<quasiAdcSigm<<endl;
      }
      } 
    }
    if(!(detector.isAD() || detector.isWaterShield())) continue; 
  } // Loop over detectors
  return StatusCode::SUCCESS;
}
bool PmtCalibLeadingEdge::hasStats ( const DayaBay::Detector detector) [private]

Check if the statistics for this detector have been initialized.

Definition at line 1063 of file PmtCalibLeadingEdge.cc.

                                                                 {
  // Check if statistics have been initialized for this detector
  return (std::find(m_processedDetectors.begin(), 
                    m_processedDetectors.end(), 
                    detector) 
          != m_processedDetectors.end());
}
StatusCode PmtCalibLeadingEdge::prepareStats ( const Context context) [private]

Get the current statistics for a detector, or create if needed.

Definition at line 15 of file PmtCalibLeadingEdgePrepare.cc.

                                                                  {
  // 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;
    }
  }

  // NChannel 
  {
    std::ostringstream title, path;
    std::string name = "NChannel";
    title << "NChannel " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* nchannel = new TH1F(name.c_str(),title.str().c_str(),200,0,200);
    nchannel->GetXaxis()->SetTitle("NChannel");
    nchannel->GetYaxis()->SetTitle("Entries");
    if( m_statsSvc->put(path.str(),nchannel).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete nchannel;
      return StatusCode::FAILURE;
    }
  }
  // Raw Adc, mean
  {
    std::ostringstream title, path;
    std::string name = "adcRawMean";
    title << "Raw ADC, mean in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* adcRawMean = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    adcRawMean->GetXaxis()->SetTitle("board*16+connector");
    adcRawMean->GetYaxis()->SetTitle("Raw ADC mean");
    if( m_statsSvc->put(path.str(),adcRawMean).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete adcRawMean;
      return StatusCode::FAILURE;
    }
  }
  // Raw Adc, sigma
  {
    std::ostringstream title, path;
    std::string name = "adcRawSigma";
    title << "Raw ADC, sigma in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* adcRawSigma = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    adcRawSigma->GetXaxis()->SetTitle("board*16+connector");
    adcRawSigma->GetYaxis()->SetTitle("Raw ADC sigma");
    if( m_statsSvc->put(path.str(),adcRawSigma).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete adcRawSigma;
      return StatusCode::FAILURE;
    }
  }
  // Adc - preAdc, mean
  {
    std::ostringstream title, path;
    std::string name = "quasiAdcMean";
    title << "ADC - preAdc, mean in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* quasiAdcMean = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    quasiAdcMean->GetXaxis()->SetTitle("board*16+connector");
    quasiAdcMean->GetYaxis()->SetTitle("QuasiAdc mean");
    if( m_statsSvc->put(path.str(),quasiAdcMean).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete quasiAdcMean;
      return StatusCode::FAILURE;
    }
  }
  // Adc - preAdc, sigma
  {
    std::ostringstream title, path;
    std::string name = "quasiAdcSigma";
    title << "ADC - preAdc, sigma in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* quasiAdcSigma = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    quasiAdcSigma->GetXaxis()->SetTitle("board*16+connector");
    quasiAdcSigma->GetYaxis()->SetTitle("QuasiAdc sigma");
    if( m_statsSvc->put(path.str(),quasiAdcSigma).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete quasiAdcSigma;
      return StatusCode::FAILURE;
    }
  }
  // Adc - Ave preAdc, mean
  {
    std::ostringstream title, path;
    std::string name = "expAdcMean";
    title << "ADC - ave preAdc, mean in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* expAdcMean = new TH1F(name.c_str(),title.str().c_str(),400,0,400); // to hold water PMT fee board 
    expAdcMean->GetXaxis()->SetTitle("board*16+connector");
    expAdcMean->GetYaxis()->SetTitle("ExpAdc mean");
    if( m_statsSvc->put(path.str(),expAdcMean).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete expAdcMean;
      return StatusCode::FAILURE;
    }
  }
  // Adc - Ave preAdc, sigma
  {
    std::ostringstream title, path;
    std::string name = "expAdcSigma";
    title << "ADC - ave preAdc, sigma in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* expAdcSigma = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    expAdcSigma->GetXaxis()->SetTitle("board*16+connector");
    expAdcSigma->GetYaxis()->SetTitle("ExpAdc sigma");
    if( m_statsSvc->put(path.str(),expAdcSigma).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete expAdcSigma;
      return StatusCode::FAILURE;
    }
  }
    // fitting chi2/ndf
  {
    std::ostringstream title, path;
    std::string name = "Chi2_1D";
    title << "Chi2/NDF of fitting " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* Chi2_1D = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    Chi2_1D->GetXaxis()->SetTitle("board*16+connector");
    Chi2_1D->GetYaxis()->SetTitle("Chi2/NDF");
    if( m_statsSvc->put(path.str(),Chi2_1D).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete Chi2_1D;
      return StatusCode::FAILURE;
    }
  }
  // preAdc, mean
  {
    std::ostringstream title, path;
    std::string name = "preAdcMean";
    title << "preAdc, mean in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* preAdcMean = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    preAdcMean->GetXaxis()->SetTitle("board*16+connector");
    preAdcMean->GetYaxis()->SetTitle("preAdc mean");
    if( m_statsSvc->put(path.str(),preAdcMean).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete preAdcMean;
      return StatusCode::FAILURE;
    }
  }
  // preAdc, sigma
  {
    std::ostringstream title, path;
    std::string name = "preAdcSigma";
    title << "preAdc, sigma in " << detector.detName();
    path << this->getPath(detector) << "/" << name;
    TH1F* preAdcSigma = new TH1F(name.c_str(),title.str().c_str(),400,0,400);
    preAdcSigma->GetXaxis()->SetTitle("board*16+connector");
    preAdcSigma->GetYaxis()->SetTitle("preAdc sigma");
    if( m_statsSvc->put(path.str(),preAdcSigma).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete preAdcSigma;
      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(),
                            400,0,400);
    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;
    }
  }
    // AD Occupancy for each channel (jpochoa)
    {
      std::ostringstream title, path;
      std::string name = "ADoccupancy";
      title << "AD Occupancy per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* ADoccupancy = new TH2F(name.c_str(),title.str().c_str(),49,0.25,24.75,19,-0.75,8.75);
      ADoccupancy->GetXaxis()->SetTitle("Column");
      ADoccupancy->GetYaxis()->SetTitle("Ring");
      if( m_statsSvc->put(path.str(),ADoccupancy).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete ADoccupancy;
        return StatusCode::FAILURE;
      }
    }
    // WPI Occupancy for each channel (Zeyuan)
    {
      std::ostringstream title, path;
      std::string name = "WPIoccupancy";
      title << "WPI Occupancy per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      info()<<"WPI 2D::::::::::::::::"<< this->getPath(detector) << "/" << name<<endl;
      TH2F* WPIoccupancy = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPIoccupancy->GetXaxis()->SetTitle("PMT local X: m");
      WPIoccupancy->GetYaxis()->SetTitle("PMT local Y: m");
      if( m_statsSvc->put(path.str(),WPIoccupancy).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPIoccupancy;
        return StatusCode::FAILURE;
      }
    }
    // WPO Occupancy for each channel (InwardSeeing PMTs)(Zeyuan)
    {
      std::ostringstream title, path;
      std::string name = "WPOinwardoccupancy";
      title << "WPOinward Occupancy per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      info()<<"WPOinward 2D::::::::::::::::"<< this->getPath(detector) << "/" << name<<endl;
      TH2F* WPOinwardoccupancy = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPOinwardoccupancy->GetXaxis()->SetTitle("PMT local X: m");
      WPOinwardoccupancy->GetYaxis()->SetTitle("PMT local Y: m");
      if( m_statsSvc->put(path.str(),WPOinwardoccupancy).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPOinwardoccupancy;
        return StatusCode::FAILURE;
      }
    }
        // WPO Occupancy for each channel (OutwardSeeing PMTs)(Zeyuan)
    {
      std::ostringstream title, path;
      std::string name = "WPOoutwardoccupancy";
      title << "WPOoutward Occupancy per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      info()<<"WPOoutward 2D::::::::::::::::"<< this->getPath(detector) << "/" << name<<endl;
      TH2F* WPOoutwardoccupancy = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPOoutwardoccupancy->GetXaxis()->SetTitle("PMT local X: m");
      WPOoutwardoccupancy->GetYaxis()->SetTitle("PMT local Y: m");
      if( m_statsSvc->put(path.str(),WPOoutwardoccupancy).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPOoutwardoccupancy;
        return StatusCode::FAILURE;
      }
    }
  // AD Gain for each channel (jpochoa)
  {
      std::ostringstream title, path;
      std::string name = "ADgain2D";
      title << "AD Gain per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* ADgain2D = new TH2F(name.c_str(),title.str().c_str(),49,0.25,24.75,19,-0.75,8.75);
      ADgain2D->GetXaxis()->SetTitle("Column");
      ADgain2D->GetYaxis()->SetTitle("Ring");
      if( m_statsSvc->put(path.str(),ADgain2D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete ADgain2D;
        return StatusCode::FAILURE;
      }
    }
      // WPI Gain for each channel (Zeyuan)
  {
      std::ostringstream title, path;
      std::string name = "WPIgain2D";
      title << "WPI Gain per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* WPIgain2D = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPIgain2D->GetXaxis()->SetTitle("PMT local X: m");
      WPIgain2D->GetYaxis()->SetTitle("PMT local Y: m");
      if( m_statsSvc->put(path.str(),WPIgain2D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPIgain2D;
        return StatusCode::FAILURE;
      }
    }
    // WPO Gain for each channel Inwarding seeing PMTs(Zeyuan)
  {
      std::ostringstream title, path;
      std::string name = "WPOinwardgain2D";
      title << "WPOinward Gain per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* WPOinwardgain2D = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPOinwardgain2D->GetXaxis()->SetTitle("PMT local X: m");
      WPOinwardgain2D->GetYaxis()->SetTitle("PMT local Y: m");
      if( m_statsSvc->put(path.str(),WPOinwardgain2D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPOinwardgain2D;
        return StatusCode::FAILURE;
      }
    }
    // WPO Gain for each channel Outwarding seeing PMTs(Zeyuan)
  {
      std::ostringstream title, path;
      std::string name = "WPOoutwardgain2D";
      title << "WPOoutward Gain per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* WPOoutwardgain2D = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPOoutwardgain2D->GetXaxis()->SetTitle("PMT local X: m");
      WPOoutwardgain2D->GetYaxis()->SetTitle("PMT local Y: m");
      if( m_statsSvc->put(path.str(),WPOoutwardgain2D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPOoutwardgain2D;
        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(),400,-200,200);
      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;
      }
    }
    
    // AD Chi2/ndf for each channel (yuzy)
  {
      std::ostringstream title, path;
      std::string name = "ADChi2NDF";
      title << "AD Chi2/NDF per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* ADChi2NDF = new TH2F(name.c_str(),title.str().c_str(),49,0.25,24.75,19,-0.75,8.75);
      ADChi2NDF->GetXaxis()->SetTitle("Column");
      ADChi2NDF->GetYaxis()->SetTitle("Ring");
      if( m_statsSvc->put(path.str(),ADChi2NDF).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete ADChi2NDF;
        return StatusCode::FAILURE;
      }
    }
        // WPI Chi2/ndf for each channel (yuzy)
  {
      std::ostringstream title, path;
      std::string name = "WPIChi2NDF";
      title << "WPI Chi2/NDF per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* WPIChi2NDF = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPIChi2NDF->GetXaxis()->SetTitle("X :m");
      WPIChi2NDF->GetYaxis()->SetTitle("Y: m");
      if( m_statsSvc->put(path.str(),WPIChi2NDF).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPIChi2NDF;
        return StatusCode::FAILURE;
      }
    }
   // WPO infacing Chi2/ndf for each channel (yuzy)
  {
      std::ostringstream title, path;
      std::string name = "WPOinwardChi2NDF";
      title << "WPOinward Chi2/NDF per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* WPOinwardChi2NDF = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPOinwardChi2NDF->GetXaxis()->SetTitle("X:m");
      WPOinwardChi2NDF->GetYaxis()->SetTitle("Y:m");
      if( m_statsSvc->put(path.str(),WPOinwardChi2NDF).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPOinwardChi2NDF;
        return StatusCode::FAILURE;
      }
    }
       // WPO outfacing Chi2/ndf for each channel (yuzy)
  {
      std::ostringstream title, path;
      std::string name = "WPOoutwardChi2NDF";
      title << "WPOoutward Chi2/NDF per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* WPOoutwardChi2NDF = new TH2F(name.c_str(),title.str().c_str(),101,-20,20,101,-20,20);
      WPOoutwardChi2NDF->GetXaxis()->SetTitle("X:m");
      WPOoutwardChi2NDF->GetYaxis()->SetTitle("Y:m");
      if( m_statsSvc->put(path.str(),WPOoutwardChi2NDF).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete WPOoutwardChi2NDF;
        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;
    }
  }
  // TDC Mean spectrum
  {
    std::ostringstream title, path;
    std::string name = "tdcMean";
    title << "Mean TDC for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* tdcMean = new TH1F(name.c_str(),title.str().c_str(),
                               2000,0,2000);
    tdcMean->GetXaxis()->SetTitle("TDC value");
    tdcMean->GetYaxis()->SetTitle("Entries");
    if( m_statsSvc->put(path.str(),tdcMean).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete tdcMean;
      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 Mean by ring
    {
      std::ostringstream title, path;
      std::string name = "adcMean";
      title << "ADC Mean for PMTs in " << detector.detName() 
            << " ring " << ring; 
      path << this->getPath(detector, ring) << "/" << name;
      TH1F* adcMean = new TH1F(name.c_str(),title.str().c_str(),
                                 24,1,25);
      adcMean->GetXaxis()->SetTitle("Column");
      adcMean->GetYaxis()->SetTitle("Mean ADC");
      if( m_statsSvc->put(path.str(),adcMean).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcMean;
        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(),
                              2000,0,2000);
      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, mean corrected
    {
      std::string name = "tdcByMean";
      std::ostringstream title, path;
      title << "Mean-corrected TDC values for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* tdcByMean = new TH1F(name.c_str(),title.str().c_str(),
                                   600,-300,300);
      tdcByMean->GetXaxis()->SetTitle("TDC value");
      tdcByMean->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),tdcByMean).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete tdcByMean;
        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;
      }
    }
    
    // Quasi ADC spectrum, baseline subtracted ( -preAdc)
    {
      std::ostringstream title, path;
      std::string name = "quasiAdc";
      title << "ADC - preAdc " << detector.detName()
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* quasiAdc = new TH1F(name.c_str(),title.str().c_str(),
                           1200,-200,1000);
      quasiAdc->GetXaxis()->SetTitle("QuasiADC value");
      quasiAdc->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),quasiAdc).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete quasiAdc;
        return StatusCode::FAILURE;
      }
    }    

    // ADC spectrum, average pedestal substracted
    {
      std::ostringstream title, path;
      std::string name = "adc";
      title << "ADC - Average preAdc " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* adc = new TH1F(name.c_str(),title.str().c_str(),
                           1200,-200,1000);
      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;
      }
    }

    // PreADC spectrum, baseline subtracted
    {
      std::ostringstream title, path;
      std::string name = "preAdc";
      title << "PreADC values for " << detector.detName()
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* preAdc = new TH1F(name.c_str(),title.str().c_str(),
                           1200,-200,1000);
      preAdc->GetXaxis()->SetTitle("PreADC value");
      preAdc->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),preAdc).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete preAdc;
        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 PmtCalibLeadingEdge::getPath ( const DayaBay::Detector detector) [private]

Definition at line 1071 of file PmtCalibLeadingEdge.cc.

                                                                       {
  return m_filepath + "/" + detector.detName();
}
std::string PmtCalibLeadingEdge::getPath ( const DayaBay::Detector detector,
int  ring 
) [private]

Definition at line 1075 of file PmtCalibLeadingEdge.cc.

                                                  {
  std::ostringstream path;
  path << m_filepath << "/" << detector.detName() << "/ring_" << ring; 
  return path.str();
}
std::string PmtCalibLeadingEdge::getPath ( const DayaBay::FeeChannelId channelId) [private]

Definition at line 1082 of file PmtCalibLeadingEdge.cc.

                                                                       {
  std::ostringstream path;
  path << m_filepath << "/" << chanId.detName() 
       << "/chan_" << chanId.board() << "_" << chanId.connector(); 
  return path.str();
}
bool PmtCalibLeadingEdge::GoodTdc ( int  Tdc) [inline, private]

Definition at line 111 of file PmtCalibLeadingEdge.h.

    {
      if(Tdc>m_tdclow&&Tdc<m_tdchigh) return true;
      else return false;
    };
const InterfaceID & IPmtCalibParamTool::interfaceID ( ) [static, inherited]

Retrieve interface ID.

Definition at line 8 of file IPmtCalibParamTool.cc.

{ 
    return IID_IPmtCalibParamTool; 
}

Member Data Documentation

std::string PmtCalibLeadingEdge::m_cableSvcName [private]

Definition at line 115 of file PmtCalibLeadingEdge.h.

std::string PmtCalibLeadingEdge::m_calibSvcName [private]

Definition at line 123 of file PmtCalibLeadingEdge.h.

Definition at line 126 of file PmtCalibLeadingEdge.h.

std::string PmtCalibLeadingEdge::m_filepath [private]

Definition at line 129 of file PmtCalibLeadingEdge.h.

Definition at line 131 of file PmtCalibLeadingEdge.h.

Definition at line 134 of file PmtCalibLeadingEdge.h.

Definition at line 137 of file PmtCalibLeadingEdge.h.

Definition at line 140 of file PmtCalibLeadingEdge.h.

Definition at line 143 of file PmtCalibLeadingEdge.h.

Definition at line 146 of file PmtCalibLeadingEdge.h.

Definition at line 149 of file PmtCalibLeadingEdge.h.

Definition at line 152 of file PmtCalibLeadingEdge.h.

ofstream PmtCalibLeadingEdge::m_textFile [private]

Definition at line 155 of file PmtCalibLeadingEdge.h.

Definition at line 156 of file PmtCalibLeadingEdge.h.

Definition at line 158 of file PmtCalibLeadingEdge.h.

Definition at line 159 of file PmtCalibLeadingEdge.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:04:28 for CalibParam by doxygen 1.7.4