/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
PmtCalibFullModel Class Reference

#include <PmtCalibFullModel.h>

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

List of all members.

Public Member Functions

 PmtCalibFullModel (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~PmtCalibFullModel ()
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
ICableSvcm_cableSvc
ICalibDataSvcm_calibSvc
std::string m_pmtGeomSvcName
 PMT geometry.
IPmtGeomInfoSvcm_pmtGeomSvc
IStatisticsSvcm_statsSvc
IFloatingFeePedestalSvcm_floatFeePedesSvc
bool m_useFloatFeePedes
int m_tdclow
int m_tdchigh
std::vector< DayaBay::Detectorm_processedDetectors

Detailed Description

Definition at line 73 of file PmtCalibFullModel.h.


Constructor & Destructor Documentation

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

Definition at line 31 of file PmtCalibFullModel.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=true,
                  "Use FloatFeePedes or not? In case of false, it will use CalibSvc.");
  declareProperty("FilePath",m_filepath="/file0/PmtCalibFullModel",
                  "File path of registered histograms.");
                  declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc", 
      "Name of Pmt Geometry Information Service");
  declareProperty("TDCCutLow",m_tdclow=950,"low tdc cut for good hits");
  declareProperty("TDCCutHigh",m_tdchigh=1150,"high tdc cut for good hits");
}
PmtCalibFullModel::~PmtCalibFullModel ( ) [virtual]

Definition at line 55 of file PmtCalibFullModel.cc.

{}

Member Function Documentation

StatusCode PmtCalibFullModel::initialize ( ) [virtual]

Definition at line 57 of file PmtCalibFullModel.cc.

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

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

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

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

  // Get the histogram service
  m_statsSvc = svc<IStatisticsSvc>("StatisticsSvc",true);
  if(!m_statsSvc) {
    error() << " No StatisticsSvc available." << endreq;
    return StatusCode::FAILURE;
  }
  //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;
  }
  // Initialize data from input file.
  // Processed detectors are identified by directory name.
  std::vector<std::string> detDirs = m_statsSvc->getSubFolders(m_filepath);
  std::vector<std::string>::iterator dirIter, dirEnd = detDirs.end();
  for(dirIter=detDirs.begin(); dirIter != dirEnd; dirIter++){
    info() << "Processing input directory: " << m_filepath << "/" 
           << *dirIter << endreq;
    DayaBay::Detector detector( 
                 DayaBay::Detector::siteDetPackedFromString(*dirIter) );
    if( detector.site() == Site::kUnknown 
        || detector.detectorId() == DetectorId::kUnknown ){
      warning() << "Input directory: " << m_filepath << "/" 
                << *dirIter << " not processed." << endreq;
    }
    m_processedDetectors.push_back(detector);
  }

  return StatusCode::SUCCESS;
}
StatusCode PmtCalibFullModel::finalize ( ) [virtual]

Definition at line 116 of file PmtCalibFullModel.cc.

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

  return StatusCode::SUCCESS;
}
StatusCode PmtCalibFullModel::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

quasiAdcH->Fill(adc-preAdc);

Implements IPmtCalibParamTool.

Definition at line 126 of file PmtCalibFullModel.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();
  // Calib trigger only
  DayaBay::Trigger::TriggerType_t trigType = daqCrate->triggerType();
  if( (trigType & DayaBay::Trigger::kCalib) != DayaBay::Trigger::kCalib )  
    return StatusCode::SUCCESS;

  // 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;
  
  // Add the current readout to the calibration statistics
  std::vector<int> readoutTdc;
  double adcSum = 0;
  // Loop over channels
  const DayaBay::DaqPmtCrate::PmtChannelPtrList& channels
    = pmtCrate->channelReadouts();
  
  DayaBay::DaqPmtCrate::PmtChannelPtrList::const_iterator channelIter, 
    channelEnd = channels.end();
  for(channelIter = channels.begin(); channelIter!=channelEnd; channelIter++) { 
    const DayaBay::DaqPmtChannel& channel = *(*channelIter);
    const DayaBay::FeeChannelId& chanId = channel.channelId();
    const DayaBay::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* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
    if(!adcH) return StatusCode::FAILURE;
    TH1F* pedH = m_statsSvc->getTH1F( this->getPath(chanId) + "/ped");
    if(!pedH) return StatusCode::FAILURE;
    TH1F* adcByClockH = m_statsSvc->getTH1F( this->getPath(chanId) 
                                             + "/adcByClock");
    if(!adcByClockH) return StatusCode::FAILURE;

    // Loop over adc values
        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
      tdcRawH->Fill(tdc);
      // Loop over adc values

      if( adc>0 && channel.isHighGainAdc(i) ) {
        adcRawH->Fill(adc);
        adcH->Fill(adc-preAdc);
        pedH->Fill(preAdc);
        adcByClockH->Fill(adcClock,adc-preAdc);
      }
      adcSum += adc-preAdc;
    }
    // Increment number of hits on this channel
    if( FindFirstHit ) {
      nHitsPar->SetVal( nHitsPar->GetVal() + 1 );
     }
  }
  
  adcSumH->Fill(channels.size());
  // Increment number of readouts from this detector
  if(channels.size()>10)nReadoutsPar->SetVal( nReadoutsPar->GetVal() + 1 );
  return StatusCode::SUCCESS;  
}
StatusCode PmtCalibFullModel::calibrate ( ) [virtual]

calibrate() is called after processing many readouts.

This method is responsible for calculating the calibration parameters.

end of WP PMTs

Implements IPmtCalibParamTool.

Definition at line 247 of file PmtCalibFullModel.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()) {
    // 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* gainChannelH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/gainChannel");
    if(!gainChannelH) return StatusCode::FAILURE;
    TH1F* gainChannelHSigma = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/gainChannelSigma");
    if(!gainChannelHSigma) return StatusCode::FAILURE;

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

    TH2F* occupancyH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/ADoccupancy");
    if(!occupancyH2D) return StatusCode::FAILURE;
    TH2F* gainH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/ADgain2D");
    if(!gainH2D) return StatusCode::FAILURE;
    TH2F* pedH2D = m_statsSvc->getTH2F( this->getPath(detector) 
                                              + "/ADped2D");
    if(!pedH2D) 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
    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();

      // 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* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
      if(!adcH) return StatusCode::FAILURE;
      TH1F* pedH = m_statsSvc->getTH1F( this->getPath(chanId) + "/ped");
      if(!pedH) return StatusCode::FAILURE;
      TH1F* occupancyH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/occupancy");
      if(!occupancyH) return StatusCode::FAILURE;
      TH1F* adcMedianH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/adcMedian");
      if(!adcMedianH) return StatusCode::FAILURE;
      TH1F* adcSigmaH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                             +"/adcSigma");
      if(!adcSigmaH) return StatusCode::FAILURE;


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

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

        //Fit code modified by jpochoa on 03/18/11
        TF1 *nimmodelfunc = new TF1("nimmodelfunc",NIMmodel,5,100,8);
        
        nimmodelfunc->SetParNames(  "N","Q_{0}","#sigma_{0}","Q_{1}","#sigma_{1}","w", "#alpha","#mu");
        double AreaGuess, Q1MeanGuess, Q1SigmaGuess;
        float q0=0;//for a pedestal substracted distribution
        float sigma0=0;//for a pedestal substracted distribution
        AreaGuess = adcH->GetEntries();
        Q1MeanGuess = adcH->GetMean() - q0;
        Q1SigmaGuess = adcH->GetRMS();
        nimmodelfunc->SetParameters(AreaGuess, 90,1.9,Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03,    0.05);
        //fix four parameters
        nimmodelfunc->FixParameter(1,q0);
        nimmodelfunc->FixParameter(2,sigma0);
        nimmodelfunc->FixParameter(5,0.);
        nimmodelfunc->FixParameter(6,0.);
       
        adcH->Fit(nimmodelfunc,"RQ");
        if(nimmodelfunc!=NULL){
          
          adcMean=nimmodelfunc->GetParameter(3);
          adcMeanUncert =nimmodelfunc->GetParError(3);
          adcSigma = nimmodelfunc->GetParameter(4);
          adcSigmaUncert = nimmodelfunc->GetParError(4);
          double Chi2 = nimmodelfunc->GetChisquare();
          double ndf = nimmodelfunc->GetNDF();
          double fitquality = Chi2/ndf;
          if(adcMean<10 || adcMean>50 || adcMeanUncert>1 || adcSigma>15 || adcSigmaUncert>1 || fitquality>3){
                nimmodelfunc->SetParameters(AreaGuess, 90,1.9,Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03,    0.02);
                nimmodelfunc->FixParameter(1,0);
                nimmodelfunc->FixParameter(2,0);
                nimmodelfunc->FixParameter(5,0.);
                nimmodelfunc->FixParameter(6,0.);
                adcH->Fit(nimmodelfunc,"QR");
                if(adcMean<10 || adcMean>50 || adcMeanUncert>1 || adcSigma>15 || adcSigmaUncert>1 || fitquality>3){
                        nimmodelfunc->SetParameters(AreaGuess, 90,1.9,Q1MeanGuess, Q1SigmaGuess, 0.01, 0.03,    0.02);
                        nimmodelfunc->FixParameter(1,0);
                        nimmodelfunc->FixParameter(2,0);
                        nimmodelfunc->FixParameter(5,0.);
                        nimmodelfunc->FixParameter(6,0.);
                        adcH->Fit(nimmodelfunc,"RQ","",10,100);
                }
                adcMean=nimmodelfunc->GetParameter(3);
                adcMeanUncert =nimmodelfunc->GetParError(3);
                adcSigma = nimmodelfunc->GetParameter(4);
                adcSigmaUncert = nimmodelfunc->GetParError(4);
          }//print warning if issues
          gainOk = true;
          
        } else {

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

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

        gainChannelH->SetBinContent(chanId.board()*16+chanId.connector(),adcMean);
        gainChannelH->SetBinError(chanId.board()*16+chanId.connector(),adcMeanUncert);
        gainChannelHSigma->SetBinContent(chanId.board()*16+chanId.connector(),adcSigma);
        gainChannelHSigma->SetBinError(chanId.board()*16+chanId.connector(),adcSigmaUncert);
        Chi2NDF_1D->SetBinContent(chanId.board()*16+chanId.connector(), fitquality);
        gainH2D->SetBinContent(gainH2D->GetXaxis()->FindBin(column),
                               gainH2D->GetYaxis()->FindBin(ring),
                               adcMean);
        gainH2D->SetBinError(gainH2D->GetXaxis()->FindBin(column),
                               gainH2D->GetYaxis()->FindBin(ring),
                               adcMeanUncert);
        //gainH2D->SetBinContent(gainH2D->GetXaxis()->FindBin(column),
        //                     gainH2D->GetYaxis()->FindBin(ring),
        //                     adcH->GetMean());
        //gainH2D->SetBinError(gainH2D->GetXaxis()->FindBin(column),
        //                     gainH2D->GetYaxis()->FindBin(ring),
        //                     adcH->GetRMS());
        pedH2D->SetBinContent(pedH2D->GetXaxis()->FindBin(column),
                               pedH2D->GetYaxis()->FindBin(ring),
                               pedH->GetMean());
        pedH2D->SetBinError(pedH2D->GetXaxis()->FindBin(column),
                               pedH2D->GetYaxis()->FindBin(ring),
                               pedH->GetRMS());
        gainH1D->Fill(adcMean);
        Chi2NDF->SetBinContent(Chi2NDF->GetXaxis()->FindBin(column),
                               Chi2NDF->GetYaxis()->FindBin(ring),
                               fitquality);
                              

      }


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

        occupancyChannelH->SetBinContent(chanId.board()*16+chanId.connector(),occupancy);
        occupancyChannelH->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);

      }
   } //end of channels    
 } 
 
 
 
   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)
                                                +"/meanOccupancy");
    if(!meanOccupancyH) return StatusCode::FAILURE;

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

    TH1F* occupancyChannelH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/occupancyChannel");
    if(!occupancyChannelH) 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");
    TH1F* gainH1D = m_statsSvc->getTH1F( this->getPath(detector) 
                                              + "/gain1D");
        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");
    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);
        // 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();
      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) 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;}
       
       // 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* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
      if(!adcH) return StatusCode::FAILURE;

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

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

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

          warning() << "PMT full model fit error on : Board:channel " << chanId.board() << " . " 
                      <<chanId.connector() << endreq;
          gainOk = false;
        }
        double Chi2 = nimmodelfunc->GetChisquare();
        double ndf = nimmodelfunc->GetNDF();
        double fitquality = Chi2/ndf;
        Chi2NDF_1D->SetBinContent(chanId.board()*16+chanId.connector(), fitquality);
   gainH1D->Fill(adcMean);
        if(chanId.board()<=25) {
         gainChannelH->SetBinContent(chanId.board()*16+chanId.connector(),adcMean);
        gainChannelH->SetBinError(chanId.board()*16+chanId.connector(),adcMeanUncert);
        gainChannelSigmaH->SetBinContent(chanId.board()*16+chanId.connector(),adcSigma);
        gainChannelSigmaH->SetBinError(chanId.board()*16+chanId.connector(),adcSigmaUncert);
        }
        //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);}
        
            // Channel Occupancy
      double occupancy = occupancytosave;
      double occupancyUncert = occupancyerrtosave;
      {
        if(occupancy < minValidOcc) occupancyOk = false;
        occupancyChannelH->SetBinContent(chanId.board()*16+chanId.connector(),occupancy);
        occupancyChannelH->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);
      }
      }//end of channel
   }
  } // Loop over detectors
 
  return StatusCode::SUCCESS;
}
bool PmtCalibFullModel::hasStats ( const DayaBay::Detector detector) [private]

Check if the statistics for this detector have been initialized.

Definition at line 748 of file PmtCalibFullModel.cc.

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

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

Definition at line 756 of file PmtCalibFullModel.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;
    }
  }

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

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

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

  // channel occupancy
  {
    std::ostringstream title, path;
    std::string name = "occupancyChannel";
    title << "occupancyChannel " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* occupancyChannel = new TH1F(name.c_str(),title.str().c_str(),
                            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;
    }
  }


  // TDC Median spectrum
  {
    std::ostringstream title, path;
    std::string name = "tdcMedian";
    title << "Median TDC for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* tdcMedian = new TH1F(name.c_str(),title.str().c_str(),
                               400,0,400);
    tdcMedian->GetXaxis()->SetTitle("TDC value");
    tdcMedian->GetYaxis()->SetTitle("Entries");
    if( m_statsSvc->put(path.str(),tdcMedian).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete tdcMedian;
      return StatusCode::FAILURE;
    }
  }
// 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<<endreq;
      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<<endreq;
      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<<endreq;
      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;
      }
    }
  
  // AD Gain for each channel (jpochoa)
  {
      std::ostringstream title, path;
      std::string name = "ADped2D";
      title << "AD Gain per channel " << detector.detName();
      path << this->getPath(detector) << "/" << name;
      TH2F* ADped2D = new TH2F(name.c_str(),title.str().c_str(),49,0.25,24.75,19,-0.75,8.75);
      ADped2D->GetXaxis()->SetTitle("Column");
      ADped2D->GetYaxis()->SetTitle("Ring");
      if( m_statsSvc->put(path.str(),ADped2D).isFailure() ) {
          error() << "prepareStats(): Could not register " << path << endreq;
          delete ADped2D;
          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;
      }
    }
        // fitting chi2/ndf 1D 
  {
    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;
    }
  }
      // 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;
      }
    }
    // 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(),500,-100,100);
      gain1D->GetXaxis()->SetTitle("ADC/PE");
      gain1D->GetYaxis()->SetTitle("Number of entries");
      if( m_statsSvc->put(path.str(),gain1D).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete gain1D;
        return StatusCode::FAILURE;
      }
    }
    

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

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

  ServiceMode svcMode(context, 0);

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

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

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

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

    // Raw TDC spectrum
    {
      std::string name = "tdcRaw";
      std::ostringstream title, path;
      title << "Raw TDC values for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* tdcRaw = new TH1F(name.c_str(),title.str().c_str(),
                              400,0,400);
      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;
      }
    }


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

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

  m_processedDetectors.push_back(detector);
 
 
  return StatusCode::SUCCESS;
}
std::string PmtCalibFullModel::getPath ( const DayaBay::Detector detector) [private]

Definition at line 1467 of file PmtCalibFullModel.cc.

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

Definition at line 1471 of file PmtCalibFullModel.cc.

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

Definition at line 1478 of file PmtCalibFullModel.cc.

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

Definition at line 108 of file PmtCalibFullModel.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 PmtCalibFullModel::m_cableSvcName [private]

Definition at line 112 of file PmtCalibFullModel.h.

std::string PmtCalibFullModel::m_calibSvcName [private]

Definition at line 119 of file PmtCalibFullModel.h.

Definition at line 122 of file PmtCalibFullModel.h.

std::string PmtCalibFullModel::m_filepath [private]

Definition at line 125 of file PmtCalibFullModel.h.

Definition at line 128 of file PmtCalibFullModel.h.

Definition at line 131 of file PmtCalibFullModel.h.

std::string PmtCalibFullModel::m_pmtGeomSvcName [private]

PMT geometry.

Definition at line 133 of file PmtCalibFullModel.h.

Definition at line 136 of file PmtCalibFullModel.h.

Definition at line 139 of file PmtCalibFullModel.h.

Definition at line 142 of file PmtCalibFullModel.h.

Definition at line 145 of file PmtCalibFullModel.h.

Definition at line 146 of file PmtCalibFullModel.h.

Definition at line 147 of file PmtCalibFullModel.h.

Definition at line 150 of file PmtCalibFullModel.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