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

#include <PmtCalibLeadingEdgeWithCuts.h>

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

List of all members.

Public Member Functions

 PmtCalibLeadingEdgeWithCuts (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~PmtCalibLeadingEdgeWithCuts ()
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)
virtual StatusCode calibrateRaw ()
 used for firstly calibration of PMT before calibration with cuts calibrateRaw() is called after processing many readouts.

Private Attributes

std::string m_cableSvcName
std::string m_calibSvcName
std::string m_filepath
ICableSvcm_cableSvc
ICalibDataSvcm_calibSvc
IStatisticsSvcm_statsSvc
std::vector< DayaBay::Detectorm_processedDetectors
int adc_buffer [max_PMT_number][max_event_number][max_FEE_hit]
int tdc_buffer [max_PMT_number][max_event_number][max_FEE_hit]
int adc_clock_buffer [max_PMT_number][max_event_number][max_FEE_hit]
double adc_baseline_buffer [max_PMT_number][max_event_number]
double m_ADC_Cut_Min
double m_ADC_Cut_Max
std::map< int, double > m_ADC_Cut_Min_Channel
std::map< int, double > m_ADC_Cut_Max_Channel
double m_TDC_Cut_Min
double m_TDC_Cut_Max
std::map< int, double > m_TDC_Cut_Min_Channel
std::map< int, double > m_TDC_Cut_Max_Channel

Detailed Description

Definition at line 80 of file PmtCalibLeadingEdgeWithCuts.h.


Constructor & Destructor Documentation

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

Definition at line 13 of file PmtCalibLeadingEdgeWithCuts.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("FilePath",m_filepath="/file0/PmtCalibLeadingEdgeWithCuts",
                  "File path of registered histograms.");

  m_ADC_Cut_Min=0.25; // in ratio to single PE
  m_ADC_Cut_Max=-1; // in ratio to single PE

  m_TDC_Cut_Min=15; //in ns
  m_TDC_Cut_Max=15; //in ns

  for(int i=0;i<max_PMT_number;i++)
  for(int j=0;j<max_event_number;j++)
  {
    adc_baseline_buffer[i][j]=0;
    for(int k=0;k<max_FEE_hit;k++)
    {
      adc_buffer[i][j][k]=0;
      tdc_buffer[i][j][k]=0;
      adc_clock_buffer[i][j][k]=0;
    }
  }
  
}
PmtCalibLeadingEdgeWithCuts::~PmtCalibLeadingEdgeWithCuts ( ) [virtual]

Definition at line 48 of file PmtCalibLeadingEdgeWithCuts.cc.

{}

Member Function Documentation

StatusCode PmtCalibLeadingEdgeWithCuts::initialize ( ) [virtual]

Definition at line 50 of file PmtCalibLeadingEdgeWithCuts.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 the histogram service
  m_statsSvc = svc<IStatisticsSvc>("StatisticsSvc",true);
  if(!m_statsSvc) {
    error() << " No StatisticsSvc available." << endreq;
    return StatusCode::FAILURE;
  }

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

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

Definition at line 97 of file PmtCalibLeadingEdgeWithCuts.cc.

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

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

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

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

Implements IPmtCalibParamTool.

Definition at line 107 of file PmtCalibLeadingEdgeWithCuts.cc.

{
  /*
  DayaBay::Detector detector = readout.detector();
  // Only process AD detectors
  if(detector.detectorId() != DetectorId::kAD1
     && detector.detectorId() != DetectorId::kAD2
     && detector.detectorId() != DetectorId::kAD3
     && detector.detectorId() != DetectorId::kAD4){
    debug() << "process(): Do not process detector: " << detector.detName() 
            << endreq;
    return StatusCode::SUCCESS;
  }

  const DayaBay::ReadoutHeader* header = readout.header();
  if(!header){
    error() << "process(): Readout has no header!" << endreq;
    return StatusCode::FAILURE;
  }
  Context context = header->context();
  if( !hasStats(detector) ){
    // New detector; initialize all the detector statistics
    StatusCode sc = prepareStats(context);
    if( sc != StatusCode::SUCCESS ) return sc;
  }
  ServiceMode svcMode(context, 0);

  // Retrieve the parameters and histograms
  TObject* nReadoutsObj = 
    m_statsSvc->get(this->getPath(detector) + "/nReadouts");
  TParameter<int>* nReadoutsPar= dynamic_cast<TParameter<int>*>(nReadoutsObj);
  if(!nReadoutsPar) return StatusCode::FAILURE;
  TH1F* adcSumH = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSum");
  if(!adcSumH) return StatusCode::FAILURE;
  TH1F* tdcMedianH = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedian");
  if(!tdcMedianH) return StatusCode::FAILURE;

  // Add the current readout to the calibration statistics
  std::vector<int> readoutTdc;
  double adcSum = 0;
  int nReadouts = nReadoutsPar->GetVal();
  // Loop over channels
  DayaBay::ReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter,
    chanEnd = readout.channelReadout().end();
  for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
      chanIter++){
    const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
    const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
    const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId, 
                                                                     svcMode);

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

      //info()<<"test 7 chan_num="<<chan_num<<", ring="<<ring<<",column="<<column<<endreq;


    // Get Parameters and Histograms
    TObject* nHitsObj = m_statsSvc->get(this->getPath(chanId) + "/nHits");
    TParameter<int>* nHitsPar = dynamic_cast<TParameter<int>*>(nHitsObj);
    if(!nHitsPar) return StatusCode::FAILURE;
    TH1F* tdcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcRaw");
    if(!tdcRawH) return StatusCode::FAILURE;
    TH1F* adcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRaw");
    if(!adcRawH) return StatusCode::FAILURE;
    TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
    if(!adcH) return StatusCode::FAILURE;
    TH1F* adcByClockH = m_statsSvc->getTH1F( this->getPath(chanId) 
                                             + "/adcByClock");
    if(!adcByClockH) return StatusCode::FAILURE;

    std::vector<int>::const_iterator tdcIter, tdcEnd = pmtChan.tdc().end();
    // Loop over tdc values
    int tdc_hit_count=0;
    for(tdcIter = pmtChan.tdc().begin(); tdcIter!=tdcEnd; tdcIter++){
      int tdc = *tdcIter;
      readoutTdc.push_back(tdc);
      tdcRawH->Fill(tdc);
      tdc_buffer[chan_num][nReadouts%max_event_number][tdc_hit_count%max_FEE_hit]=tdc;
      tdc_hit_count++;      
    }
    // Loop over adc values
    double adcBaseline = feeCalib->m_adcBaselineLow;
    adc_baseline_buffer[chan_num][nReadouts%max_event_number]=adcBaseline;
    int adc_hit_count=0;
    for(int i=0; i<pmtChan.size(); i++) {
      int adcClock = pmtChan.adcCycle(i);
      int adc = pmtChan.adc(i);
      adcRawH->Fill(adc);
      adcH->Fill(adc-adcBaseline);
      adcByClockH->Fill(adcClock,adc-adcBaseline);
      adcSum += adc-adcBaseline;
      adc_buffer[chan_num][nReadouts%max_event_number][adc_hit_count%max_FEE_hit]=adc;
     //info()<<"test 4 chan_num="<<chan_num<<",nReadouts= "<<nReadouts<<", adc_hit_count="<<adc_hit_count<<",adc="<<adc_buffer[chan_num][nReadouts%max_event_number][adc_hit_count%max_FEE_hit]<<endreq;
      adc_clock_buffer[chan_num][nReadouts%max_event_number][adc_hit_count%max_FEE_hit]=adcClock;
      adc_hit_count++;
    }
    // Increment number of hits on this channel
    nHitsPar->SetVal( nHitsPar->GetVal() + 1 );
  }


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

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




        //info() << "test first data read over" << endreq;

  // used for preliminary calibrate the PMT and get preliminary results
  // the condition for this is the nReadouts+1==max_event_number
  // if the total nReadouts+1 of whole analysis is less than max_event_number, the  calibrateRaw() will invoke with finalize()
  // wangzhm@ihep.ac.cn 2009/11/29
  if(nReadouts+1==max_event_number) 
  if(calibrateRaw()==StatusCode::FAILURE)
  {
        warning() << "calibrateRaw() error" 
                    << endreq;
        return StatusCode::FAILURE;
  }

 //info() << "test first calibration over" << endreq;


// read the data again with the cuts from firstly calibration
// wangzhm@ihep.ac.cn 2009/11/28


  // the condition for this is the nReadouts>max_event_number
  // if the total nReadouts+1 of whole analysis is less than max_event_number, the  calibrateRaw() will invoke with finalize()
  // wangzhm@ihep.ac.cn 2009/11/29
  if(nReadouts+1>max_event_number)
  {
        info()<<"test 1 nReadouts="<<nReadouts<<endreq;

  TH1F* adcSumHWithCuts = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSumWithCuts");
  if(!adcSumHWithCuts) return StatusCode::FAILURE;
  TH1F* tdcMedianHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedianWithCuts");
  if(!tdcMedianHWithCuts) return StatusCode::FAILURE;

  // Add the current readout to the calibration statistics
  std::vector<int> readoutTdcWithCuts;
  double adcSumWithCuts = 0;
  // Loop over channels
  for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
      chanIter++){
    const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
    const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
    const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId, 
                                                                     svcMode);

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

//info()<<"process ring="<<ring<<"      column="<<column<<endreq;


    // Get Parameters and Histograms
    TObject* nHitsObjWithCuts = m_statsSvc->get(this->getPath(chanId) + "/nHitsWithCuts");
    TParameter<int>* nHitsParWithCuts = dynamic_cast<TParameter<int>*>(nHitsObjWithCuts);
    if(!nHitsParWithCuts) return StatusCode::FAILURE;
    TH1F* tdcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcWithCuts");
    if(!tdcHWithCuts) return StatusCode::FAILURE;
    TH1F* adcRawHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRawWithCuts");
    if(!adcRawHWithCuts) return StatusCode::FAILURE;
    TH1F* adcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcWithCuts");
    if(!adcHWithCuts) return StatusCode::FAILURE;
    TH1F* adcByClockHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) 
                                             + "/adcByClockWithCuts");
    if(!adcByClockHWithCuts) return StatusCode::FAILURE;

    bool hit_ok_flag=0;
    double adcBaseline = feeCalib->m_adcBaselineLow;
    // Loop over tdc values
    for(int i=0; i<pmtChan.size(); i++) {
      int tdc = pmtChan.tdc(i);

      int adcClock = pmtChan.adcCycle(i);
      int adc = pmtChan.adc(i);

      if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
      if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
      {
        readoutTdcWithCuts.push_back(tdc);
        tdcHWithCuts->Fill(tdc);

        adcRawHWithCuts->Fill(adc);
        adcHWithCuts->Fill(adc-adcBaseline);
        adcByClockHWithCuts->Fill(adcClock,adc-adcBaseline);
        adcSumWithCuts += adc-adcBaseline;
        hit_ok_flag=1;
      }

    }
    // Increment number of hits on this channel
    if(hit_ok_flag)
      nHitsParWithCuts->SetVal( nHitsParWithCuts->GetVal() + 1 );
  }

  if(readoutTdcWithCuts.size() > 0){
    // Find the median TDC
    std::sort(readoutTdcWithCuts.begin(), readoutTdcWithCuts.end());
    int medianIdx = readoutTdcWithCuts.size()/2;
    int medianTdc = readoutTdcWithCuts[medianIdx];
    for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
        chanIter++){
      const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
      const DayaBay::FeeChannelId& chanId = pmtChan.channelId();

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


      TH1F* tdcByMedianHWithCuts = 
        m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMedianWithCuts");
      if(!tdcByMedianHWithCuts) return StatusCode::FAILURE;
      // Loop over tdc values
    for(int i=0; i<pmtChan.size(); i++) {
      int tdc = pmtChan.tdc(i);
      int adc = pmtChan.adc(i);

      if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
      if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
        tdcByMedianHWithCuts->Fill(tdc-medianTdc);
      }
    }
    tdcMedianHWithCuts->Fill(medianTdc);
  }

  adcSumHWithCuts->Fill(adcSumWithCuts);


  }//over for if(nReadouts>=max_event_number)
//info() << "test second read data over" << endreq;


  */
  return StatusCode::SUCCESS;  
}
StatusCode PmtCalibLeadingEdgeWithCuts::calibrate ( ) [virtual]

calibrate() is called after processing many readouts.

This method is responsible for calculating the calibration parameters.

Implements IPmtCalibParamTool.

Definition at line 845 of file PmtCalibLeadingEdgeWithCuts.cc.

                                                 {
  // Loop over detectors in summary data, and calculate parameters

  // before the final calibrate of PMT
  // we need to fit the raw data without any cuts
  // and if the total nreadouts number less than max_event_number,
  // the data with cuts need to be fill firstly which will realized in the same program
  // 2009/11/29
  if(calibrateRaw()==StatusCode::FAILURE)
  {
        warning() << "calibrateRaw() error" 
                    << endreq;
        return StatusCode::FAILURE;
  }


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

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

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

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


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

    TH1F* tdcOffsetChannelHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/tdcOffsetChannelWithCuts");
    if(!tdcOffsetChannelHWithCuts) return StatusCode::FAILURE;
    TH1F* occupancyChannelHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/occupancyChannelWithCuts");
    if(!occupancyChannelHWithCuts) return StatusCode::FAILURE;


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

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

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

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

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

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

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

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

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

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

      // Get Parameters and Histograms
      TObject* nHitsObjWithCuts = m_statsSvc->get(this->getPath(chanId) + "/nHitsWithCuts");
      TParameter<int>* nHitsParWithCuts = dynamic_cast<TParameter<int>*>(nHitsObjWithCuts);
      if(!nHitsParWithCuts) return StatusCode::FAILURE;
      TH1F* tdcByMedianHWithCuts= m_statsSvc->getTH1F( this->getPath(chanId)
                                               +"/tdcByMedianWithCuts");
      if(!tdcByMedianHWithCuts) return StatusCode::FAILURE;
      TH1F* adcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcWithCuts");
      if(!adcHWithCuts) return StatusCode::FAILURE;
      TH1F* occupancyHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/occupancyWithCuts");
      if(!occupancyHWithCuts) return StatusCode::FAILURE;
      TH1F* tdcOffsetHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/tdcOffsetWithCuts");
      if(!tdcOffsetHWithCuts) return StatusCode::FAILURE;
      TH1F* adcMedianHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/adcMedianWithCuts");
      if(!adcMedianHWithCuts) return StatusCode::FAILURE;
      TH1F* adcSigmaHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                             +"/adcSigmaWithCuts");
      if(!adcSigmaHWithCuts) return StatusCode::FAILURE;

      // Channel Occupancy
      double occupancy = 0;
      double occupancyUncert = 0;
      {
        int nHits = nHitsParWithCuts->GetVal();
        if(nReadouts > 0){
          occupancy = nHits / ((double) nReadouts);
          occupancyUncert = 
            sqrt( occupancy*(1-occupancy)/((double) nReadouts) );
        }
        occupancyHWithCuts->SetBinContent(occupancyHWithCuts->FindBin(column),occupancy);
        occupancyHWithCuts->SetBinError(occupancyHWithCuts->FindBin(column),occupancyUncert);
        if(occupancy < minValidOcc) occupancyOk = false;

        occupancyChannelHWithCuts->SetBinContent(ring*24+column,occupancy);
        occupancyChannelHWithCuts->SetBinError(ring*24+column,occupancyUncert);


      }
      // TDC offset
      double tdcOffset = 0;
      double tdcOffsetUncert = 0;
      {
        // Use a fit to the leading edge to determine the time offset
        int maxBin = tdcByMedianHWithCuts->GetMaximumBin();
        double fitMin = tdcByMedianHWithCuts->GetBinCenter(maxBin - 2);
        double fitMax = tdcByMedianHWithCuts->GetBinCenter(maxBin + 10);
        tdcByMedianHWithCuts->Fit("gaus","","",fitMin, fitMax);
        TF1* fitFunc = tdcByMedianHWithCuts->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());
          
          tdcOffsetHWithCuts->SetBinContent(tdcOffsetHWithCuts->FindBin(column),tdcOffset);
          tdcOffsetHWithCuts->SetBinError(tdcOffsetHWithCuts->FindBin(column),tdcOffsetUncert);

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


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

      // ADC median and sigma
      double adcMean = 0;
      double adcMeanUncert = 0;
      double adcSigma = 0;
      double adcSigmaUncert = 0;
      {
        // Find ADC SPE peak, width with simple gaussian
        TF1 adcFit("adcFit","gaus");
        double adcPeak = adcHWithCuts->GetBinCenter(adcHWithCuts->GetMaximumBin());
        double adcSig = 5.0;
        adcHWithCuts->Fit(&adcFit,"","",adcPeak - adcSig, adcPeak + adcSig);
        adcMean = adcFit.GetParameter(1);
        if(fabs(adcMean-adcPeak) > adcSig){
          warning() << "Issues with fit on ring, column: " << ring << " / " 
                    << column << " mean/peak: " << adcMean << " / " << adcPeak 
                    << endreq;
          gainOk = false;
        }
        adcMeanUncert = adcFit.GetParError(1);
        adcSigma = adcFit.GetParameter(2);
        adcSigmaUncert = adcFit.GetParError(2);
        adcMedianHWithCuts->SetBinContent(adcMedianHWithCuts->FindBin(column),adcMean);
        adcMedianHWithCuts->SetBinError(adcMedianHWithCuts->FindBin(column),adcMeanUncert);
        adcSigmaHWithCuts->SetBinContent(adcSigmaHWithCuts->FindBin(column),adcSigma);
        adcSigmaHWithCuts->SetBinError(adcSigmaHWithCuts->FindBin(column),adcSigmaUncert);

        meanAdcHWithCuts->SetBinContent(ring,adcMean);
        meanAdcHWithCuts->SetBinError(ring, adcMeanUncert);

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

      }

      // 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 ){
        meanOccupancyHWithCuts->SetBinContent(ring,
                                      meanOccRing[ring] / meanOccWeight[ring]);
        meanOccupancyHWithCuts->SetBinError(ring, 
                                    sqrt( 1.0 / meanOccWeight[ring] ));
      }
      if( meanTdcWeight[ring] > 0 ){
        meanTdcOffsetHWithCuts->SetBinContent(ring,
                                      meanTdcRing[ring] / meanTdcWeight[ring]);
        meanTdcOffsetHWithCuts->SetBinError(ring, 
                                    sqrt( 1.0 / meanTdcWeight[ring] ));
      }
    } // Loop over rings
 
  } // Loop over detectors

info() << "test second calibrate over" << endreq;

  return StatusCode::SUCCESS;

}
bool PmtCalibLeadingEdgeWithCuts::hasStats ( const DayaBay::Detector detector) [private]

Check if the statistics for this detector have been initialized.

Definition at line 1108 of file PmtCalibLeadingEdgeWithCuts.cc.

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

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

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

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


  // Gain vs channel spectrum
  {
    std::ostringstream title, path;
    std::string name = "gainChannel";
    title << "gainChannel for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* gainChannel = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    gainChannel->GetXaxis()->SetTitle("Channel");
    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;
    }
  }

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

  // tdc offset vs channel spectrum
  {
    std::ostringstream title, path;
    std::string name = "tdcOffsetChannel";
    title << "tdcOffsetChannel for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* tdcOffsetChannel = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    tdcOffsetChannel->GetXaxis()->SetTitle("Channel");
    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;
    }
  }

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


  // occupancy vs channel spectrum
  {
    std::ostringstream title, path;
    std::string name = "occupancyChannel";
    title << "occupancyChannel for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* occupancyChannel = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    occupancyChannel->GetXaxis()->SetTitle("Channel");
    occupancyChannel->GetYaxis()->SetTitle("occupancy");
    // 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;
    }
  }

  // occupancy vs channel with cuts
  {
    std::ostringstream title, path;
    std::string name = "occupancyChannelWithCuts";
    title << "occupancyChannel for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* occupancyChannelWithCuts = new TH1F(name.c_str(),title.str().c_str(),
                            300,0,300);
    occupancyChannelWithCuts->GetXaxis()->SetTitle("channel");
    occupancyChannelWithCuts->GetYaxis()->SetTitle("occupancy");
    // DEBUG
    info() << "name= " << occupancyChannelWithCuts->GetName() 
           << " at path = \"" << path.str() << "\"" << endreq;
    if( m_statsSvc->put(path.str(),occupancyChannelWithCuts).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete occupancyChannelWithCuts;
      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(),
                               800,0,800);
    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;
    }
  }

  // TDC Median spectrum with cuts
  {
    std::ostringstream title, path;
    std::string name = "tdcMedianWithCuts";
    title << "Median TDC WithCuts for readouts in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* tdcMedianWithCuts = new TH1F(name.c_str(),title.str().c_str(),
                               800,0,800);
    tdcMedianWithCuts->GetXaxis()->SetTitle("TDC value");
    tdcMedianWithCuts->GetYaxis()->SetTitle("Entries");
    if( m_statsSvc->put(path.str(),tdcMedianWithCuts).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete tdcMedianWithCuts;
      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 TDC Offset WithCuts by AD ring
  {
    std::ostringstream title, path;
    std::string name = "meanTdcOffsetWithCuts";
    title << "Mean TDC offset by ring in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* meanTdcOffsetWithCuts = new TH1F(name.c_str(),title.str().c_str(),
                                   8,1,9);
    meanTdcOffsetWithCuts->GetXaxis()->SetTitle("AD Ring");
    meanTdcOffsetWithCuts->GetYaxis()->SetTitle("Mean TDC Offset With Cuts");
    if( m_statsSvc->put(path.str(),meanTdcOffsetWithCuts).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete meanTdcOffsetWithCuts;
      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;
    }
  }

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

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

  // Mean Adc With Cuts by AD ring 
  {
    std::ostringstream title, path;
    std::string name = "meanAdcWithCuts";
    title << "Mean Adc With Cuts by ring in " << detector.detName(); 
    path << this->getPath(detector) << "/" << name;
    TH1F* meanAdcWithCuts = new TH1F(name.c_str(),title.str().c_str(),
                                   8,1,9);
    meanAdcWithCuts->GetXaxis()->SetTitle("AD Ring");
    meanAdcWithCuts->GetYaxis()->SetTitle("Mean Adc With Cuts");
    if( m_statsSvc->put(path.str(),meanAdcWithCuts).isFailure() ) {
      error() << "prepareStats(): Could not register " << path << endreq;
      delete meanAdcWithCuts;
      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;
      }
    }
    // TDC Raw offset by ring
    {
      std::ostringstream title, path;
      std::string name = "tdcRawOffset";
      title << "Time Raw Offset for PMTs in " << detector.detName() 
            << " ring " << ring; 
      path << this->getPath(detector, ring) << "/" << name;
      TH1F* tdcRawOffset = new TH1F(name.c_str(),title.str().c_str(),
                           24,1,25);
      tdcRawOffset->GetXaxis()->SetTitle("Column");
      tdcRawOffset->GetYaxis()->SetTitle("Time Offset [tdc]");
      if( m_statsSvc->put(path.str(),tdcRawOffset).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete tdcRawOffset;
        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;
      }
    }
  }

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

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


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

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


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

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

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

    // ADC spectrum, baseline subtracted with cuts
    {
      std::ostringstream title, path;
      std::string name = "adcWithCuts";
      title << "Baseline-subtracted ADC values With Cuts for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* adcWithCuts = new TH1F(name.c_str(),title.str().c_str(),
                           5096,-1000,4096);
      adcWithCuts->GetXaxis()->SetTitle("ADC value");
      adcWithCuts->GetYaxis()->SetTitle("Entries");
      if( m_statsSvc->put(path.str(),adcWithCuts).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcWithCuts;
        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;
      }
    }
    // ADC spectrum by Clock Cycle with cuts
    {
      std::ostringstream title, path;
      std::string name = "adcByClockWithCuts";
      title << "ADC values by clock cycle for " << detector.detName() 
            << " channel " << chanId.board() << "_"
            << chanId.connector();
      path << this->getPath(chanId) << "/" << name;
      TH1F* adcByClockWithCuts = new TH1F(name.c_str(),title.str().c_str(),
                                  20,0,20);
      adcByClockWithCuts->GetXaxis()->SetTitle("ADC Clock Cycle");
      adcByClockWithCuts->GetYaxis()->SetTitle("ADC");
      if( m_statsSvc->put(path.str(),adcByClockWithCuts).isFailure() ) {
        error() << "prepareStats(): Could not register " << path << endreq;
        delete adcByClockWithCuts;
        return StatusCode::FAILURE;
      }
    }


  }  

  m_processedDetectors.push_back(detector);

  return StatusCode::SUCCESS;
}
std::string PmtCalibLeadingEdgeWithCuts::getPath ( const DayaBay::Detector detector) [private]

Definition at line 1952 of file PmtCalibLeadingEdgeWithCuts.cc.

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

Definition at line 1956 of file PmtCalibLeadingEdgeWithCuts.cc.

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

Definition at line 1963 of file PmtCalibLeadingEdgeWithCuts.cc.

                                                                               {
  std::ostringstream path;
  path << m_filepath << "/" << chanId.detName() 
       << "/chan_" << chanId.board() << "_" << chanId.connector(); 
  return path.str();
}
StatusCode PmtCalibLeadingEdgeWithCuts::calibrateRaw ( ) [private, virtual]

used for firstly calibration of PMT before calibration with cuts calibrateRaw() is called after processing many readouts.

This method is responsible for calculating the calibration parameters. wangzhm@ihep.ac.cn 2009/11/28

Definition at line 391 of file PmtCalibLeadingEdgeWithCuts.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;

    // 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* meanAdcH = m_statsSvc->getTH1F( this->getPath(detector)
                                                +"/meanAdc");
    if(!meanAdcH) return StatusCode::FAILURE;

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


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

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

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

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

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

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

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

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

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

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

      // Get Parameters and Histograms
      TObject* nHitsObj = m_statsSvc->get(this->getPath(chanId) + "/nHits");
      TParameter<int>* nHitsPar = dynamic_cast<TParameter<int>*>(nHitsObj);
      if(!nHitsPar) return StatusCode::FAILURE;
      TH1F* tdcByMedianH= m_statsSvc->getTH1F( this->getPath(chanId)
                                               +"/tdcByMedian");
      if(!tdcByMedianH) return StatusCode::FAILURE;
      TH1F* tdcRawH= m_statsSvc->getTH1F( this->getPath(chanId)
                                               +"/tdcRaw");
      if(!tdcRawH) return StatusCode::FAILURE;
      TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
      if(!adcH) return StatusCode::FAILURE;
      TH1F* occupancyH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/occupancy");
      if(!occupancyH) return StatusCode::FAILURE;
      TH1F* tdcOffsetH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/tdcOffset");
      if(!tdcOffsetH) return StatusCode::FAILURE;
      TH1F* tdcRawOffsetH = m_statsSvc->getTH1F( this->getPath(detector,ring)
                                              +"/tdcRawOffset");
      if(!tdcRawOffsetH) 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;

      // Channel Occupancy
      double occupancy = 0;
      double occupancyUncert = 0;
      {
        int nHits = nHitsPar->GetVal();
        if(nReadouts > 0){
          occupancy = nHits / ((double) nReadouts);
          occupancyUncert = 
            sqrt( occupancy*(1-occupancy)/((double) nReadouts) );
        }
        occupancyH->SetBinContent(occupancyH->FindBin(column),occupancy);
        occupancyH->SetBinError(occupancyH->FindBin(column),occupancyUncert);
        if(occupancy < minValidOcc) occupancyOk = false;

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

      }

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

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

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

      }

     // TDC offset Raw
      tdcOffset = 0;
      tdcOffsetUncert = 0;
      {
        // Use a fit to the leading edge to determine the time offset
        int maxBin = tdcRawH->GetMaximumBin();
        double fitMin = tdcRawH->GetBinCenter(maxBin - 2);
        double fitMax = tdcRawH->GetBinCenter(maxBin + 10);
        tdcRawH->Fit("gaus","","",fitMin, fitMax);
        TF1* fitFunc = tdcRawH->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());
          
          tdcRawOffsetH->SetBinContent(tdcOffsetH->FindBin(column),tdcOffset);
          tdcRawOffsetH->SetBinError(tdcOffsetH->FindBin(column),tdcOffsetUncert);
        }else{
          warning() << "PMT  " << pmtId << " has no TDC data for fitting." 
                    << endreq;
        }
          //set for cuts of TDC
        if(m_TDC_Cut_Min>=0)
        {
          m_TDC_Cut_Min_Channel[channel_number]=tdcOffset-m_TDC_Cut_Min;
        }
        else
        {
          m_TDC_Cut_Min_Channel[channel_number]=0;
        }
        if(m_TDC_Cut_Max>=0)
        {
          m_TDC_Cut_Max_Channel[channel_number]=tdcOffset+m_TDC_Cut_Max;
        }
        else
        {
          m_TDC_Cut_Max_Channel[channel_number]=10000;
        }

          //info() << "Issues with fit on ring, column: " << ring << " / " << column << " m_TDC_Cut_Min_Channel/m_TDC_Cut_Max_Channel: " << m_TDC_Cut_Min_Channel[channel_number] << " / " << m_TDC_Cut_Max_Channel[channel_number] << endreq;

      }



      // ADC median and sigma
      double adcMean = 0;
      double adcMeanUncert = 0;
      double adcSigma = 0;
      double adcSigmaUncert = 0;
        const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId, svcMode);
        double adcBaseline = feeCalib->m_adcBaselineLow;
      {
        // Find ADC SPE peak, width with simple gaussian
        TF1 adcFit("adcFit","gaus");
        double adcPeak = adcH->GetBinCenter(adcH->GetMaximumBin());
        double adcSig = 5.0;
        adcH->Fit(&adcFit,"","",adcPeak - adcSig, adcPeak + adcSig);
        adcMean = adcFit.GetParameter(1);
        if(fabs(adcMean-adcPeak) > adcSig){
          warning() << "Issues with fit on ring, column: " << ring << " / " 
                    << column << " mean/peak: " << adcMean << " / " << adcPeak 
                    << endreq;
          gainOk = false;
        }
          //set for cuts of ADC
        if(m_ADC_Cut_Min>=0)
        {
          m_ADC_Cut_Min_Channel[channel_number]=(adcMean-m_ADC_Cut_Min*(adcMean-adcBaseline));
        }
        else
        {
          m_ADC_Cut_Min_Channel[channel_number]=0;
        }
        if(m_ADC_Cut_Max>=0)
        {
          m_ADC_Cut_Max_Channel[channel_number]=(adcMean+m_ADC_Cut_Max*(adcMean-adcBaseline));
        }
        else
        {
          m_ADC_Cut_Max_Channel[channel_number]=4096;
        }

        adcMeanUncert = adcFit.GetParError(1);
        adcSigma = adcFit.GetParameter(2);
        adcSigmaUncert = adcFit.GetParError(2);
        adcMedianH->SetBinContent(adcMedianH->FindBin(column),adcMean);
        adcMedianH->SetBinError(adcMedianH->FindBin(column),adcMeanUncert);
        adcSigmaH->SetBinContent(adcSigmaH->FindBin(column),adcSigma);
        adcSigmaH->SetBinError(adcSigmaH->FindBin(column),adcSigmaUncert);

        meanAdcH->SetBinContent(ring,adcMean);
        meanAdcH->SetBinError(ring, adcMeanUncert);

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



      }

      // 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

//fill the data after cuts if the total events number less than max_event_number
  if(nReadouts<=max_event_number)
  for(int event_count=0;event_count<nReadouts;event_count++)
  {

        //info()<<"test 2 event_count="<<event_count<<endreq;

        TH1F* adcSumHWithCuts = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSumWithCuts");
        if(!adcSumHWithCuts) return StatusCode::FAILURE;
        TH1F* tdcMedianHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedianWithCuts");
        if(!tdcMedianHWithCuts) return StatusCode::FAILURE;

        // Add the current readout to the calibration statistics
        std::vector<int> readoutTdcWithCuts;
        double adcSumWithCuts = 0;
    // Initialize statistics for each channel
        for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
        {
                DayaBay::FeeChannelId chanId = *chanIter;

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

                // Get Parameters and Histograms
                TObject* nHitsObjWithCuts = m_statsSvc->get(this->getPath(chanId) + "/nHitsWithCuts");
                TParameter<int>* nHitsParWithCuts = dynamic_cast<TParameter<int>*>(nHitsObjWithCuts);
                if(!nHitsParWithCuts) return StatusCode::FAILURE;
                TH1F* tdcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcWithCuts");
                if(!tdcHWithCuts) return StatusCode::FAILURE;
                TH1F* adcRawHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRawWithCuts");
                if(!adcRawHWithCuts) return StatusCode::FAILURE;
                TH1F* adcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcWithCuts");
                if(!adcHWithCuts) return StatusCode::FAILURE;
                TH1F* adcByClockHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) 
                                             + "/adcByClockWithCuts");
                if(!adcByClockHWithCuts) return StatusCode::FAILURE;

                // Loop over tdc values
                int chan_num=ring*24+column;
                double baseline=adc_baseline_buffer[chan_num%max_PMT_number][event_count];
                for(int k=0;k<max_FEE_hit;k++)
                {
                        int tdc=tdc_buffer[chan_num%max_PMT_number][event_count][k];
                        int adc=adc_buffer[chan_num%max_PMT_number][event_count][k];
                        int adcClock=adc_clock_buffer[chan_num%max_PMT_number][event_count][k];
                        bool hit_ok_flag=0;
                        //if(tdc>0)
                        //{
                        //info()<<"tdc="<<tdc<<",m_TDC_Cut_Min_Channel="<<m_TDC_Cut_Min_Channel[channel_number]<<",m_TDC_Cut_Max_Channel="<<m_TDC_Cut_Max_Channel[channel_number]<<endreq;
                        //info()<<"adc="<<adc<<",m_ADC_Cut_Min_Channel="<<m_ADC_Cut_Min_Channel[channel_number]<<",m_ADC_Cut_Max_Channel="<<m_ADC_Cut_Max_Channel[channel_number]<<endreq;
                        //}
                                if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
                                if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
                                {
                                        //info()<<"test 5 chan_num="<<chan_num<<",event_count="<<event_count<<",k="<<k<<",adc="<<adc<<endreq;
                                        readoutTdcWithCuts.push_back(tdc);
                                        tdcHWithCuts->Fill(tdc);

                                        adcRawHWithCuts->Fill(adc);
                                        adcHWithCuts->Fill(adc-baseline);
                                        adcByClockHWithCuts->Fill(adcClock,adc-baseline);
                                        adcSumWithCuts += adc-baseline;
                                        hit_ok_flag=1;

                                        //info()<<"test 8 adcSumWithCuts="<<adcSumWithCuts<<endreq;
                                }


                                // Increment number of hits on this channel
                                if(hit_ok_flag)
                                                nHitsParWithCuts->SetVal( nHitsParWithCuts->GetVal() + 1 );
                } //over for  for(int k=0;k<max_FEE_hit;k++)

                //info()<<"test 9 adcSumWithCuts="<<adcSumWithCuts<<endreq;


         }// over for for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)

                if(readoutTdcWithCuts.size() > 0)
                {
                                // Find the median TDC
                                std::sort(readoutTdcWithCuts.begin(), readoutTdcWithCuts.end());
                                int medianIdx = readoutTdcWithCuts.size()/2;
                                int medianTdc = readoutTdcWithCuts[medianIdx];
                                for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
                                {
                                        DayaBay::FeeChannelId chanId = *chanIter;

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


                                        TH1F* tdcByMedianHWithCuts = 
                                        m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMedianWithCuts");
                                        if(!tdcByMedianHWithCuts) return StatusCode::FAILURE;

                                        int chan_num=ring*24+column;
                                        for(int k=0;k<max_FEE_hit;k++)
                                        {
                                                int tdc=tdc_buffer[chan_num%max_PMT_number][event_count][k];
                                                int adc=adc_buffer[chan_num%max_PMT_number][event_count][k];

                                                if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
                                                if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
                                                tdcByMedianHWithCuts->Fill(tdc-medianTdc);
                                        }
                                } //over for for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
                                tdcMedianHWithCuts->Fill(medianTdc);
                } //over for if(readoutTdcWithCuts.size() > 0)

        //info()<<"test 3 event_count="<<event_count<<" adcSumWithCuts="<<adcSumWithCuts<<endreq;
        adcSumHWithCuts->Fill(adcSumWithCuts);

  } //over if(nReadouts<=max_event_number) and for(int event_count=0;event_count<nReadouts;event_count++)


 
  } // Loop over detectors


  return StatusCode::SUCCESS;
}
const InterfaceID & IPmtCalibParamTool::interfaceID ( ) [static, inherited]

Retrieve interface ID.

Definition at line 8 of file IPmtCalibParamTool.cc.

{ 
    return IID_IPmtCalibParamTool; 
}

Member Data Documentation

Definition at line 128 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 131 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 134 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 137 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 140 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 143 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 146 of file PmtCalibLeadingEdgeWithCuts.h.

int PmtCalibLeadingEdgeWithCuts::adc_buffer[max_PMT_number][max_event_number][max_FEE_hit] [private]

Definition at line 156 of file PmtCalibLeadingEdgeWithCuts.h.

int PmtCalibLeadingEdgeWithCuts::tdc_buffer[max_PMT_number][max_event_number][max_FEE_hit] [private]

Definition at line 157 of file PmtCalibLeadingEdgeWithCuts.h.

int PmtCalibLeadingEdgeWithCuts::adc_clock_buffer[max_PMT_number][max_event_number][max_FEE_hit] [private]

Definition at line 158 of file PmtCalibLeadingEdgeWithCuts.h.

double PmtCalibLeadingEdgeWithCuts::adc_baseline_buffer[max_PMT_number][max_event_number] [private]

Definition at line 159 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 166 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 167 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 171 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 172 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 179 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 180 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 183 of file PmtCalibLeadingEdgeWithCuts.h.

Definition at line 184 of file PmtCalibLeadingEdgeWithCuts.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