CalibParam: /home/dayabay/installs/trunk_2011_04_11_opt/NuWa-trunk/dybgaudi/Calibration/CalibParam/src/components/PmtCalibLeadingEdgeWithCuts.cc Source File

00001 #include "PmtCalibLeadingEdgeWithCuts.h"
00002 
00003 #include "Event/ReadoutHeader.h"
00004 #include "StatisticsSvc/IStatisticsSvc.h"
00005 #include "DataSvc/ICableSvc.h"
00006 #include "DataSvc/ICalibDataSvc.h"
00007 #include "Context/ServiceMode.h"
00008 #include "TH1F.h"
00009 #include "TF1.h"
00010 #include "TParameter.h"
00011 #include <math.h>
00012 
00013 PmtCalibLeadingEdgeWithCuts::PmtCalibLeadingEdgeWithCuts(const std::string& type,
00014                                          const std::string& name, 
00015                                          const IInterface* parent)
00016   : GaudiTool(type,name,parent)
00017   , m_cableSvc(0)
00018   , m_statsSvc(0)
00019 {
00020   declareInterface< IPmtCalibParamTool >(this) ;   
00021   declareProperty("CableSvcName",m_cableSvcName="StaticCableSvc",
00022                   "Name of service to map between detector, hardware, and electronic IDs");
00023   declareProperty("CalibSvcName",m_calibSvcName="StaticCalibDataSvc",
00024                   "Name of service to access FEE channel baselines");
00025   declareProperty("FilePath",m_filepath="/file0/PmtCalibLeadingEdgeWithCuts",
00026                   "File path of registered histograms.");
00027 
00028   m_ADC_Cut_Min=0.25; // in ratio to single PE
00029   m_ADC_Cut_Max=-1; // in ratio to single PE
00030 
00031   m_TDC_Cut_Min=15; //in ns
00032   m_TDC_Cut_Max=15; //in ns
00033 
00034   for(int i=0;i<max_PMT_number;i++)
00035   for(int j=0;j<max_event_number;j++)
00036   {
00037     adc_baseline_buffer[i][j]=0;
00038     for(int k=0;k<max_FEE_hit;k++)
00039     {
00040       adc_buffer[i][j][k]=0;
00041       tdc_buffer[i][j][k]=0;
00042       adc_clock_buffer[i][j][k]=0;
00043     }
00044   }
00045   
00046 }
00047 
00048 PmtCalibLeadingEdgeWithCuts::~PmtCalibLeadingEdgeWithCuts(){}
00049 
00050 StatusCode PmtCalibLeadingEdgeWithCuts::initialize()
00051 {
00052   info() << "initialize()" << endreq;
00053   StatusCode sc = this->GaudiTool::initialize();
00054   if( sc != StatusCode::SUCCESS ) return sc;
00055 
00056   // Get Cable Service
00057   m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
00058   if(!m_cableSvc){
00059     error() << "Failed to access cable svc: " << m_cableSvcName << endreq;
00060     return StatusCode::FAILURE;
00061   }
00062 
00063   // Get Calibration Service
00064   m_calibSvc = svc<ICalibDataSvc>(m_calibSvcName,true);
00065   if(!m_calibSvc){
00066     error() << "Failed to access calibration svc: " << m_calibSvcName << endreq;
00067     return StatusCode::FAILURE;
00068   }
00069 
00070   // Get the histogram service
00071   m_statsSvc = svc<IStatisticsSvc>("StatisticsSvc",true);
00072   if(!m_statsSvc) {
00073     error() << " No StatisticsSvc available." << endreq;
00074     return StatusCode::FAILURE;
00075   }
00076 
00077   // Initialize data from input file.
00078   // Processed detectors are identified by directory name.
00079   std::vector<std::string> detDirs = m_statsSvc->getSubFolders(m_filepath);
00080   std::vector<std::string>::iterator dirIter, dirEnd = detDirs.end();
00081   for(dirIter=detDirs.begin(); dirIter != dirEnd; dirIter++){
00082     info() << "Processing input directory: " << m_filepath << "/" 
00083            << *dirIter << endreq;
00084     DayaBay::Detector detector( 
00085                  DayaBay::Detector::siteDetPackedFromString(*dirIter) );
00086     if( detector.site() == Site::kUnknown 
00087         || detector.detectorId() == DetectorId::kUnknown ){
00088       warning() << "Input directory: " << m_filepath << "/" 
00089                 << *dirIter << " not processed." << endreq;
00090     }
00091     m_processedDetectors.push_back(detector);
00092   }
00093 
00094   return StatusCode::SUCCESS;
00095 }
00096 
00097 StatusCode PmtCalibLeadingEdgeWithCuts::finalize()
00098 {
00099   StatusCode sc = this->GaudiTool::finalize();
00100   if( sc != StatusCode::SUCCESS ) return sc;
00101 
00102   return StatusCode::SUCCESS;
00103 }
00104 
00106 
00107 StatusCode PmtCalibLeadingEdgeWithCuts::process(const DayaBay::ReadoutHeader& readout)
00108 {
00109   /*
00110   DayaBay::Detector detector = readout.detector();
00111   // Only process AD detectors
00112   if(detector.detectorId() != DetectorId::kAD1
00113      && detector.detectorId() != DetectorId::kAD2
00114      && detector.detectorId() != DetectorId::kAD3
00115      && detector.detectorId() != DetectorId::kAD4){
00116     debug() << "process(): Do not process detector: " << detector.detName() 
00117             << endreq;
00118     return StatusCode::SUCCESS;
00119   }
00120 
00121   const DayaBay::ReadoutHeader* header = readout.header();
00122   if(!header){
00123     error() << "process(): Readout has no header!" << endreq;
00124     return StatusCode::FAILURE;
00125   }
00126   Context context = header->context();
00127   if( !hasStats(detector) ){
00128     // New detector; initialize all the detector statistics
00129     StatusCode sc = prepareStats(context);
00130     if( sc != StatusCode::SUCCESS ) return sc;
00131   }
00132   ServiceMode svcMode(context, 0);
00133 
00134   // Retrieve the parameters and histograms
00135   TObject* nReadoutsObj = 
00136     m_statsSvc->get(this->getPath(detector) + "/nReadouts");
00137   TParameter<int>* nReadoutsPar= dynamic_cast<TParameter<int>*>(nReadoutsObj);
00138   if(!nReadoutsPar) return StatusCode::FAILURE;
00139   TH1F* adcSumH = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSum");
00140   if(!adcSumH) return StatusCode::FAILURE;
00141   TH1F* tdcMedianH = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedian");
00142   if(!tdcMedianH) return StatusCode::FAILURE;
00143 
00144   // Add the current readout to the calibration statistics
00145   std::vector<int> readoutTdc;
00146   double adcSum = 0;
00147   int nReadouts = nReadoutsPar->GetVal();
00148   // Loop over channels
00149   DayaBay::ReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter,
00150     chanEnd = readout.channelReadout().end();
00151   for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
00152       chanIter++){
00153     const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
00154     const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
00155     const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId, 
00156                                                                      svcMode);
00157 
00158       // Analyze the channel data
00159       // Determine ring and column
00160       DayaBay::AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chanId, svcMode);
00161       int ring = pmtId.ring();
00162       int column = pmtId.column();
00163       int chan_num=(ring*24+column)%max_PMT_number;
00164 
00165       //info()<<"test 7 chan_num="<<chan_num<<", ring="<<ring<<",column="<<column<<endreq;
00166 
00167 
00168     // Get Parameters and Histograms
00169     TObject* nHitsObj = m_statsSvc->get(this->getPath(chanId) + "/nHits");
00170     TParameter<int>* nHitsPar = dynamic_cast<TParameter<int>*>(nHitsObj);
00171     if(!nHitsPar) return StatusCode::FAILURE;
00172     TH1F* tdcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcRaw");
00173     if(!tdcRawH) return StatusCode::FAILURE;
00174     TH1F* adcRawH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRaw");
00175     if(!adcRawH) return StatusCode::FAILURE;
00176     TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
00177     if(!adcH) return StatusCode::FAILURE;
00178     TH1F* adcByClockH = m_statsSvc->getTH1F( this->getPath(chanId) 
00179                                              + "/adcByClock");
00180     if(!adcByClockH) return StatusCode::FAILURE;
00181 
00182     std::vector<int>::const_iterator tdcIter, tdcEnd = pmtChan.tdc().end();
00183     // Loop over tdc values
00184     int tdc_hit_count=0;
00185     for(tdcIter = pmtChan.tdc().begin(); tdcIter!=tdcEnd; tdcIter++){
00186       int tdc = *tdcIter;
00187       readoutTdc.push_back(tdc);
00188       tdcRawH->Fill(tdc);
00189       tdc_buffer[chan_num][nReadouts%max_event_number][tdc_hit_count%max_FEE_hit]=tdc;
00190       tdc_hit_count++;      
00191     }
00192     // Loop over adc values
00193     double adcBaseline = feeCalib->m_adcBaselineLow;
00194     adc_baseline_buffer[chan_num][nReadouts%max_event_number]=adcBaseline;
00195     int adc_hit_count=0;
00196     for(int i=0; i<pmtChan.size(); i++) {
00197       int adcClock = pmtChan.adcCycle(i);
00198       int adc = pmtChan.adc(i);
00199       adcRawH->Fill(adc);
00200       adcH->Fill(adc-adcBaseline);
00201       adcByClockH->Fill(adcClock,adc-adcBaseline);
00202       adcSum += adc-adcBaseline;
00203       adc_buffer[chan_num][nReadouts%max_event_number][adc_hit_count%max_FEE_hit]=adc;
00204      //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;
00205       adc_clock_buffer[chan_num][nReadouts%max_event_number][adc_hit_count%max_FEE_hit]=adcClock;
00206       adc_hit_count++;
00207     }
00208     // Increment number of hits on this channel
00209     nHitsPar->SetVal( nHitsPar->GetVal() + 1 );
00210   }
00211 
00212 
00213   if(readoutTdc.size() > 0){
00214     // Find the median TDC
00215     std::sort(readoutTdc.begin(), readoutTdc.end());
00216     int medianIdx = readoutTdc.size()/2;
00217     int medianTdc = readoutTdc[medianIdx];
00218     for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
00219         chanIter++){
00220       const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
00221       const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
00222       TH1F* tdcByMedianH = 
00223         m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMedian");
00224       if(!tdcByMedianH) return StatusCode::FAILURE;
00225       std::vector<int>::const_iterator tdcIter, tdcEnd = pmtChan.tdc().end();
00226       // Loop over tdc values
00227       for(tdcIter = pmtChan.tdc().begin(); tdcIter!=tdcEnd; tdcIter++){
00228         int tdc = *tdcIter;
00229         tdcByMedianH->Fill(tdc-medianTdc);
00230       }
00231     }
00232     tdcMedianH->Fill(medianTdc);
00233   }
00234 
00235   adcSumH->Fill(adcSum);
00236   // Increment number of readouts from this detector
00237   nReadoutsPar->SetVal( nReadoutsPar->GetVal() + 1 );
00238 
00239 
00240 
00241 
00242         //info() << "test first data read over" << endreq;
00243 
00244   // used for preliminary calibrate the PMT and get preliminary results
00245   // the condition for this is the nReadouts+1==max_event_number
00246   // if the total nReadouts+1 of whole analysis is less than max_event_number, the  calibrateRaw() will invoke with finalize()
00247   // wangzhm@ihep.ac.cn 2009/11/29
00248   if(nReadouts+1==max_event_number) 
00249   if(calibrateRaw()==StatusCode::FAILURE)
00250   {
00251         warning() << "calibrateRaw() error" 
00252                     << endreq;
00253         return StatusCode::FAILURE;
00254   }
00255 
00256  //info() << "test first calibration over" << endreq;
00257 
00258 
00260 // read the data again with the cuts from firstly calibration
00261 // wangzhm@ihep.ac.cn 2009/11/28
00263 
00264 
00265   // the condition for this is the nReadouts>max_event_number
00266   // if the total nReadouts+1 of whole analysis is less than max_event_number, the  calibrateRaw() will invoke with finalize()
00267   // wangzhm@ihep.ac.cn 2009/11/29
00268   if(nReadouts+1>max_event_number)
00269   {
00270         info()<<"test 1 nReadouts="<<nReadouts<<endreq;
00271 
00272   TH1F* adcSumHWithCuts = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSumWithCuts");
00273   if(!adcSumHWithCuts) return StatusCode::FAILURE;
00274   TH1F* tdcMedianHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedianWithCuts");
00275   if(!tdcMedianHWithCuts) return StatusCode::FAILURE;
00276 
00277   // Add the current readout to the calibration statistics
00278   std::vector<int> readoutTdcWithCuts;
00279   double adcSumWithCuts = 0;
00280   // Loop over channels
00281   for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
00282       chanIter++){
00283     const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
00284     const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
00285     const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId, 
00286                                                                      svcMode);
00287 
00288       // Analyze the channel data
00289       // Determine ring and column
00290       DayaBay::AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chanId, svcMode);
00291       int ring = pmtId.ring();
00292       int column = pmtId.column();
00293         int channel_number=ring*100+column;
00294 
00295 //info()<<"process ring="<<ring<<"      column="<<column<<endreq;
00296 
00297 
00298     // Get Parameters and Histograms
00299     TObject* nHitsObjWithCuts = m_statsSvc->get(this->getPath(chanId) + "/nHitsWithCuts");
00300     TParameter<int>* nHitsParWithCuts = dynamic_cast<TParameter<int>*>(nHitsObjWithCuts);
00301     if(!nHitsParWithCuts) return StatusCode::FAILURE;
00302     TH1F* tdcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcWithCuts");
00303     if(!tdcHWithCuts) return StatusCode::FAILURE;
00304     TH1F* adcRawHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRawWithCuts");
00305     if(!adcRawHWithCuts) return StatusCode::FAILURE;
00306     TH1F* adcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcWithCuts");
00307     if(!adcHWithCuts) return StatusCode::FAILURE;
00308     TH1F* adcByClockHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) 
00309                                              + "/adcByClockWithCuts");
00310     if(!adcByClockHWithCuts) return StatusCode::FAILURE;
00311 
00312     bool hit_ok_flag=0;
00313     double adcBaseline = feeCalib->m_adcBaselineLow;
00314     // Loop over tdc values
00315     for(int i=0; i<pmtChan.size(); i++) {
00316       int tdc = pmtChan.tdc(i);
00317 
00318       int adcClock = pmtChan.adcCycle(i);
00319       int adc = pmtChan.adc(i);
00320 
00321       if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
00322       if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
00323       {
00324         readoutTdcWithCuts.push_back(tdc);
00325         tdcHWithCuts->Fill(tdc);
00326 
00327         adcRawHWithCuts->Fill(adc);
00328         adcHWithCuts->Fill(adc-adcBaseline);
00329         adcByClockHWithCuts->Fill(adcClock,adc-adcBaseline);
00330         adcSumWithCuts += adc-adcBaseline;
00331         hit_ok_flag=1;
00332       }
00333 
00334     }
00335     // Increment number of hits on this channel
00336     if(hit_ok_flag)
00337       nHitsParWithCuts->SetVal( nHitsParWithCuts->GetVal() + 1 );
00338   }
00339 
00340   if(readoutTdcWithCuts.size() > 0){
00341     // Find the median TDC
00342     std::sort(readoutTdcWithCuts.begin(), readoutTdcWithCuts.end());
00343     int medianIdx = readoutTdcWithCuts.size()/2;
00344     int medianTdc = readoutTdcWithCuts[medianIdx];
00345     for(chanIter = readout.channelReadout().begin(); chanIter != chanEnd; 
00346         chanIter++){
00347       const DayaBay::ReadoutPmtChannel& pmtChan = chanIter->second;
00348       const DayaBay::FeeChannelId& chanId = pmtChan.channelId();
00349 
00350       // Analyze the channel data
00351       // Determine ring and column
00352       DayaBay::AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chanId, svcMode);
00353       int ring = pmtId.ring();
00354       int column = pmtId.column();
00355         int channel_number=ring*100+column;
00356 
00357 
00358       TH1F* tdcByMedianHWithCuts = 
00359         m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMedianWithCuts");
00360       if(!tdcByMedianHWithCuts) return StatusCode::FAILURE;
00361       // Loop over tdc values
00362     for(int i=0; i<pmtChan.size(); i++) {
00363       int tdc = pmtChan.tdc(i);
00364       int adc = pmtChan.adc(i);
00365 
00366       if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
00367       if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
00368         tdcByMedianHWithCuts->Fill(tdc-medianTdc);
00369       }
00370     }
00371     tdcMedianHWithCuts->Fill(medianTdc);
00372   }
00373 
00374   adcSumHWithCuts->Fill(adcSumWithCuts);
00375 
00376 
00377   }//over for if(nReadouts>=max_event_number)
00378 //info() << "test second read data over" << endreq;
00379 
00382 
00383   */
00384   return StatusCode::SUCCESS;  
00385 }
00386 
00387 
00389 //wangzhm@ihep.ac.cn 2009/11/28
00391 StatusCode PmtCalibLeadingEdgeWithCuts::calibrateRaw(){
00392   // Loop over detectors in summary data, and calculate parameters
00393 
00394   std::vector<DayaBay::Detector>::iterator detIter, 
00395     detEnd = m_processedDetectors.end();
00396   for(detIter = m_processedDetectors.begin(); detIter != detEnd; detIter++){
00397     DayaBay::Detector detector = *detIter;
00398 
00399     // Retrieve the data description and histograms
00400     TObject* nReadoutsObj = 
00401       m_statsSvc->get(this->getPath(detector) + "/nReadouts");
00402     TParameter<int>* nReadoutsPar= dynamic_cast<TParameter<int>*>(nReadoutsObj);
00403     if(!nReadoutsPar) return StatusCode::FAILURE;
00404 
00405     TH1F* meanOccupancyH = m_statsSvc->getTH1F( this->getPath(detector)
00406                                                 +"/meanOccupancy");
00407     if(!meanOccupancyH) return StatusCode::FAILURE;
00408 
00409     TH1F* meanAdcH = m_statsSvc->getTH1F( this->getPath(detector)
00410                                                 +"/meanAdc");
00411     if(!meanAdcH) return StatusCode::FAILURE;
00412 
00413     TH1F* gainChannelH = m_statsSvc->getTH1F( this->getPath(detector)
00414                                                 +"/gainChannel");
00415     if(!gainChannelH) return StatusCode::FAILURE;
00416     TH1F* tdcOffsetChannelH = m_statsSvc->getTH1F( this->getPath(detector)
00417                                                 +"/tdcOffsetChannel");
00418     if(!tdcOffsetChannelH) return StatusCode::FAILURE;
00419     TH1F* occupancyChannelH = m_statsSvc->getTH1F( this->getPath(detector)
00420                                                 +"/occupancyChannel");
00421     if(!occupancyChannelH) return StatusCode::FAILURE;
00422 
00423 
00424     TH1F* meanTdcOffsetH = m_statsSvc->getTH1F( this->getPath(detector)
00425                                                 +"/meanTdcOffset");
00426     if(!meanTdcOffsetH) return StatusCode::FAILURE;
00427 
00428     TObject* simFlagObj = m_statsSvc->get(this->getPath(detector) +"/simFlag");
00429     TParameter<int>* simFlagPar = dynamic_cast<TParameter<int>*>(simFlagObj);
00430     if(!simFlagPar) return StatusCode::FAILURE;
00431 
00432     TObject* timeSecObj = m_statsSvc->get(this->getPath(detector) +"/timeSec");
00433     TParameter<int>* timeSecPar = dynamic_cast<TParameter<int>*>(timeSecObj);
00434     if(!timeSecPar) return StatusCode::FAILURE;
00435 
00436     TObject* timeNanoSecObj = 
00437       m_statsSvc->get(this->getPath(detector) +"/timeNanoSec");
00438     TParameter<int>* timeNanoSecPar = 
00439       dynamic_cast<TParameter<int>*>(timeNanoSecObj);
00440     if(!timeNanoSecPar) return StatusCode::FAILURE;
00441 
00442     Site::Site_t site = detector.site();
00443     DetectorId::DetectorId_t detId = detector.detectorId();
00444     SimFlag::SimFlag_t simFlag = (SimFlag::SimFlag_t)(simFlagPar->GetVal());
00445     time_t timeSec = (time_t)(timeSecPar->GetVal());
00446     int timeNanoSec = timeNanoSecPar->GetVal();
00447     int nReadouts = nReadoutsPar->GetVal();
00448 
00449     // Context, ServiceMode
00450     Context context(site,simFlag,TimeStamp(timeSec,timeNanoSec),detId);
00451     ServiceMode svcMode(context, 0);
00452 
00453     // Prepare data for ring-to-ring comparison
00454     int nRings = 8;
00455     std::map<int,double> meanOccRing;
00456     std::map<int,double> meanOccWeight;
00457     std::map<int,double> meanTdcRing;
00458     std::map<int,double> meanTdcWeight;
00459     for(int ring = 1; ring <= nRings; ring++){
00460       meanOccRing[ring] = 0.0;
00461       meanOccWeight[ring] = 0.0;
00462       meanTdcRing[ring] = 0.0;
00463       meanTdcWeight[ring] = 0.0;
00464     }
00465     double minValidOcc = 0.01;  // Cut on PMTs with bad/low occupancy
00466     std::vector<DayaBay::FeeChannelId> badChannels;
00467 
00468     // Get list of all FEE channels
00469     const std::vector<DayaBay::FeeChannelId>& channelList 
00470       = m_cableSvc->feeChannelIds( svcMode );
00471     std::vector<DayaBay::FeeChannelId>::const_iterator chanIter, 
00472       chanEnd = channelList.end();
00473     // Initialize statistics for each channel
00474     for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++){
00475       DayaBay::FeeChannelId chanId = *chanIter;
00476 
00477       // Analyze the channel data
00478       // Determine ring and column
00479       DayaBay::AdPmtSensor pmtId = 
00480         m_cableSvc->adPmtSensor(chanId, svcMode);
00481       int ring = pmtId.ring();
00482       int column = pmtId.column();
00483         int channel_number=ring*100+column;
00484 
00485       // Channel status
00486       bool occupancyOk = true;
00487       bool gainOk = true;
00488 
00489       // Get Parameters and Histograms
00490       TObject* nHitsObj = m_statsSvc->get(this->getPath(chanId) + "/nHits");
00491       TParameter<int>* nHitsPar = dynamic_cast<TParameter<int>*>(nHitsObj);
00492       if(!nHitsPar) return StatusCode::FAILURE;
00493       TH1F* tdcByMedianH= m_statsSvc->getTH1F( this->getPath(chanId)
00494                                                +"/tdcByMedian");
00495       if(!tdcByMedianH) return StatusCode::FAILURE;
00496       TH1F* tdcRawH= m_statsSvc->getTH1F( this->getPath(chanId)
00497                                                +"/tdcRaw");
00498       if(!tdcRawH) return StatusCode::FAILURE;
00499       TH1F* adcH = m_statsSvc->getTH1F( this->getPath(chanId) + "/adc");
00500       if(!adcH) return StatusCode::FAILURE;
00501       TH1F* occupancyH = m_statsSvc->getTH1F( this->getPath(detector,ring)
00502                                               +"/occupancy");
00503       if(!occupancyH) return StatusCode::FAILURE;
00504       TH1F* tdcOffsetH = m_statsSvc->getTH1F( this->getPath(detector,ring)
00505                                               +"/tdcOffset");
00506       if(!tdcOffsetH) return StatusCode::FAILURE;
00507       TH1F* tdcRawOffsetH = m_statsSvc->getTH1F( this->getPath(detector,ring)
00508                                               +"/tdcRawOffset");
00509       if(!tdcRawOffsetH) return StatusCode::FAILURE;
00510       TH1F* adcMedianH = m_statsSvc->getTH1F( this->getPath(detector,ring)
00511                                               +"/adcMedian");
00512       if(!adcMedianH) return StatusCode::FAILURE;
00513       TH1F* adcSigmaH = m_statsSvc->getTH1F( this->getPath(detector,ring)
00514                                              +"/adcSigma");
00515       if(!adcSigmaH) return StatusCode::FAILURE;
00516 
00517       // Channel Occupancy
00518       double occupancy = 0;
00519       double occupancyUncert = 0;
00520       {
00521         int nHits = nHitsPar->GetVal();
00522         if(nReadouts > 0){
00523           occupancy = nHits / ((double) nReadouts);
00524           occupancyUncert = 
00525             sqrt( occupancy*(1-occupancy)/((double) nReadouts) );
00526         }
00527         occupancyH->SetBinContent(occupancyH->FindBin(column),occupancy);
00528         occupancyH->SetBinError(occupancyH->FindBin(column),occupancyUncert);
00529         if(occupancy < minValidOcc) occupancyOk = false;
00530 
00531         occupancyChannelH->SetBinContent(ring*24+column,occupancy);
00532         occupancyChannelH->SetBinError(ring*24+column,occupancyUncert);
00533 
00534       }
00535 
00536       // TDC offset by Median
00537       double tdcOffset = 0;
00538       double tdcOffsetUncert = 0;
00539       {
00540         // Use a fit to the leading edge to determine the time offset
00541         int maxBin = tdcByMedianH->GetMaximumBin();
00542         double fitMin = tdcByMedianH->GetBinCenter(maxBin - 2);
00543         double fitMax = tdcByMedianH->GetBinCenter(maxBin + 10);
00544         tdcByMedianH->Fit("gaus","","",fitMin, fitMax);
00545         TF1* fitFunc = tdcByMedianH->GetFunction("gaus");
00546         if( fitFunc ){
00547           tdcOffset = fitFunc->GetParameter(1);
00548           tdcOffsetUncert = fitFunc->GetParError(1);
00549           
00550           if( fitFunc->GetNDF() < 3 ) 
00551             // Catch bad fits
00552             tdcOffsetUncert = 0;
00553           else
00554             // Scale uncertainty by the quality of the fit
00555             tdcOffsetUncert *= sqrt(fitFunc->GetChisquare() / fitFunc->GetNDF());
00556           
00557           tdcOffsetH->SetBinContent(tdcOffsetH->FindBin(column),tdcOffset);
00558           tdcOffsetH->SetBinError(tdcOffsetH->FindBin(column),tdcOffsetUncert);
00559 
00560         tdcOffsetChannelH->SetBinContent(ring*24+column,tdcOffset);
00561         tdcOffsetChannelH->SetBinError(ring*24+column,tdcOffsetUncert);
00562 
00563         }else{
00564           warning() << "PMT  " << pmtId << " has no TDC data for fitting." 
00565                     << endreq;
00566         }
00567 
00568       }
00569 
00570      // TDC offset Raw
00571       tdcOffset = 0;
00572       tdcOffsetUncert = 0;
00573       {
00574         // Use a fit to the leading edge to determine the time offset
00575         int maxBin = tdcRawH->GetMaximumBin();
00576         double fitMin = tdcRawH->GetBinCenter(maxBin - 2);
00577         double fitMax = tdcRawH->GetBinCenter(maxBin + 10);
00578         tdcRawH->Fit("gaus","","",fitMin, fitMax);
00579         TF1* fitFunc = tdcRawH->GetFunction("gaus");
00580         if( fitFunc ){
00581           tdcOffset = fitFunc->GetParameter(1);
00582           tdcOffsetUncert = fitFunc->GetParError(1);
00583           
00584           if( fitFunc->GetNDF() < 3 ) 
00585             // Catch bad fits
00586             tdcOffsetUncert = 0;
00587           else
00588             // Scale uncertainty by the quality of the fit
00589             tdcOffsetUncert *= sqrt(fitFunc->GetChisquare() / fitFunc->GetNDF());
00590           
00591           tdcRawOffsetH->SetBinContent(tdcOffsetH->FindBin(column),tdcOffset);
00592           tdcRawOffsetH->SetBinError(tdcOffsetH->FindBin(column),tdcOffsetUncert);
00593         }else{
00594           warning() << "PMT  " << pmtId << " has no TDC data for fitting." 
00595                     << endreq;
00596         }
00597           //set for cuts of TDC
00598         if(m_TDC_Cut_Min>=0)
00599         {
00600           m_TDC_Cut_Min_Channel[channel_number]=tdcOffset-m_TDC_Cut_Min;
00601         }
00602         else
00603         {
00604           m_TDC_Cut_Min_Channel[channel_number]=0;
00605         }
00606         if(m_TDC_Cut_Max>=0)
00607         {
00608           m_TDC_Cut_Max_Channel[channel_number]=tdcOffset+m_TDC_Cut_Max;
00609         }
00610         else
00611         {
00612           m_TDC_Cut_Max_Channel[channel_number]=10000;
00613         }
00614 
00615           //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;
00616 
00617       }
00618 
00619 
00620 
00621       // ADC median and sigma
00622       double adcMean = 0;
00623       double adcMeanUncert = 0;
00624       double adcSigma = 0;
00625       double adcSigmaUncert = 0;
00626         const DayaBay::FeeCalibData* feeCalib = m_calibSvc->feeCalibData(chanId, svcMode);
00627         double adcBaseline = feeCalib->m_adcBaselineLow;
00628       {
00629         // Find ADC SPE peak, width with simple gaussian
00630         TF1 adcFit("adcFit","gaus");
00631         double adcPeak = adcH->GetBinCenter(adcH->GetMaximumBin());
00632         double adcSig = 5.0;
00633         adcH->Fit(&adcFit,"","",adcPeak - adcSig, adcPeak + adcSig);
00634         adcMean = adcFit.GetParameter(1);
00635         if(fabs(adcMean-adcPeak) > adcSig){
00636           warning() << "Issues with fit on ring, column: " << ring << " / " 
00637                     << column << " mean/peak: " << adcMean << " / " << adcPeak 
00638                     << endreq;
00639           gainOk = false;
00640         }
00641           //set for cuts of ADC
00642         if(m_ADC_Cut_Min>=0)
00643         {
00644           m_ADC_Cut_Min_Channel[channel_number]=(adcMean-m_ADC_Cut_Min*(adcMean-adcBaseline));
00645         }
00646         else
00647         {
00648           m_ADC_Cut_Min_Channel[channel_number]=0;
00649         }
00650         if(m_ADC_Cut_Max>=0)
00651         {
00652           m_ADC_Cut_Max_Channel[channel_number]=(adcMean+m_ADC_Cut_Max*(adcMean-adcBaseline));
00653         }
00654         else
00655         {
00656           m_ADC_Cut_Max_Channel[channel_number]=4096;
00657         }
00658 
00659         adcMeanUncert = adcFit.GetParError(1);
00660         adcSigma = adcFit.GetParameter(2);
00661         adcSigmaUncert = adcFit.GetParError(2);
00662         adcMedianH->SetBinContent(adcMedianH->FindBin(column),adcMean);
00663         adcMedianH->SetBinError(adcMedianH->FindBin(column),adcMeanUncert);
00664         adcSigmaH->SetBinContent(adcSigmaH->FindBin(column),adcSigma);
00665         adcSigmaH->SetBinError(adcSigmaH->FindBin(column),adcSigmaUncert);
00666 
00667         meanAdcH->SetBinContent(ring,adcMean);
00668         meanAdcH->SetBinError(ring, adcMeanUncert);
00669 
00670         gainChannelH->SetBinContent(ring*24+column,adcMean);
00671         gainChannelH->SetBinError(ring*24+column,adcMeanUncert);
00672 
00673 
00674 
00675       }
00676 
00677       // Record ring averages, ignoring bad channels
00678       if( occupancyOk && gainOk ){
00679         if( occupancyUncert > 0 ){
00680           meanOccRing[ring] += occupancy / (occupancyUncert*occupancyUncert); 
00681           meanOccWeight[ring] += 1.0 / (occupancyUncert*occupancyUncert); 
00682         }
00683         if( tdcOffsetUncert > 0 ){
00684           meanTdcRing[ring] += tdcOffset / (tdcOffsetUncert*tdcOffsetUncert); 
00685           meanTdcWeight[ring] += 1.0 / (tdcOffsetUncert*tdcOffsetUncert); 
00686         }
00687       }
00688 
00689     }
00690 
00691     // Record mean time offset and mean occupancy by ring
00692     for(int ring = 1; ring <= nRings; ring++){
00693       if( meanOccWeight[ring] > 0 ){
00694         meanOccupancyH->SetBinContent(ring,
00695                                       meanOccRing[ring] / meanOccWeight[ring]);
00696         meanOccupancyH->SetBinError(ring, 
00697                                     sqrt( 1.0 / meanOccWeight[ring] ));
00698       }
00699       if( meanTdcWeight[ring] > 0 ){
00700         meanTdcOffsetH->SetBinContent(ring,
00701                                       meanTdcRing[ring] / meanTdcWeight[ring]);
00702         meanTdcOffsetH->SetBinError(ring, 
00703                                     sqrt( 1.0 / meanTdcWeight[ring] ));
00704       }
00705     } // Loop over rings
00706 
00708 //fill the data after cuts if the total events number less than max_event_number
00711   if(nReadouts<=max_event_number)
00712   for(int event_count=0;event_count<nReadouts;event_count++)
00713   {
00714 
00715         //info()<<"test 2 event_count="<<event_count<<endreq;
00716 
00717         TH1F* adcSumHWithCuts = m_statsSvc->getTH1F( this->getPath(detector) + "/adcSumWithCuts");
00718         if(!adcSumHWithCuts) return StatusCode::FAILURE;
00719         TH1F* tdcMedianHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)+"/tdcMedianWithCuts");
00720         if(!tdcMedianHWithCuts) return StatusCode::FAILURE;
00721 
00722         // Add the current readout to the calibration statistics
00723         std::vector<int> readoutTdcWithCuts;
00724         double adcSumWithCuts = 0;
00725     // Initialize statistics for each channel
00726         for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00727         {
00728                 DayaBay::FeeChannelId chanId = *chanIter;
00729 
00730                 // Analyze the channel data
00731                 // Determine ring and column
00732                 DayaBay::AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chanId, svcMode);
00733                 int ring = pmtId.ring();
00734                 int column = pmtId.column();
00735                 int channel_number=ring*100+column;
00736 
00737                 // Get Parameters and Histograms
00738                 TObject* nHitsObjWithCuts = m_statsSvc->get(this->getPath(chanId) + "/nHitsWithCuts");
00739                 TParameter<int>* nHitsParWithCuts = dynamic_cast<TParameter<int>*>(nHitsObjWithCuts);
00740                 if(!nHitsParWithCuts) return StatusCode::FAILURE;
00741                 TH1F* tdcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcWithCuts");
00742                 if(!tdcHWithCuts) return StatusCode::FAILURE;
00743                 TH1F* adcRawHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcRawWithCuts");
00744                 if(!adcRawHWithCuts) return StatusCode::FAILURE;
00745                 TH1F* adcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcWithCuts");
00746                 if(!adcHWithCuts) return StatusCode::FAILURE;
00747                 TH1F* adcByClockHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) 
00748                                              + "/adcByClockWithCuts");
00749                 if(!adcByClockHWithCuts) return StatusCode::FAILURE;
00750 
00751                 // Loop over tdc values
00752                 int chan_num=ring*24+column;
00753                 double baseline=adc_baseline_buffer[chan_num%max_PMT_number][event_count];
00754                 for(int k=0;k<max_FEE_hit;k++)
00755                 {
00756                         int tdc=tdc_buffer[chan_num%max_PMT_number][event_count][k];
00757                         int adc=adc_buffer[chan_num%max_PMT_number][event_count][k];
00758                         int adcClock=adc_clock_buffer[chan_num%max_PMT_number][event_count][k];
00759                         bool hit_ok_flag=0;
00760                         //if(tdc>0)
00761                         //{
00762                         //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;
00763                         //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;
00764                         //}
00765                                 if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
00766                                 if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
00767                                 {
00768                                         //info()<<"test 5 chan_num="<<chan_num<<",event_count="<<event_count<<",k="<<k<<",adc="<<adc<<endreq;
00769                                         readoutTdcWithCuts.push_back(tdc);
00770                                         tdcHWithCuts->Fill(tdc);
00771 
00772                                         adcRawHWithCuts->Fill(adc);
00773                                         adcHWithCuts->Fill(adc-baseline);
00774                                         adcByClockHWithCuts->Fill(adcClock,adc-baseline);
00775                                         adcSumWithCuts += adc-baseline;
00776                                         hit_ok_flag=1;
00777 
00778                                         //info()<<"test 8 adcSumWithCuts="<<adcSumWithCuts<<endreq;
00779                                 }
00780 
00781 
00782                                 // Increment number of hits on this channel
00783                                 if(hit_ok_flag)
00784                                                 nHitsParWithCuts->SetVal( nHitsParWithCuts->GetVal() + 1 );
00785                 } //over for  for(int k=0;k<max_FEE_hit;k++)
00786 
00787                 //info()<<"test 9 adcSumWithCuts="<<adcSumWithCuts<<endreq;
00788 
00789 
00790          }// over for for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00791 
00792                 if(readoutTdcWithCuts.size() > 0)
00793                 {
00794                                 // Find the median TDC
00795                                 std::sort(readoutTdcWithCuts.begin(), readoutTdcWithCuts.end());
00796                                 int medianIdx = readoutTdcWithCuts.size()/2;
00797                                 int medianTdc = readoutTdcWithCuts[medianIdx];
00798                                 for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00799                                 {
00800                                         DayaBay::FeeChannelId chanId = *chanIter;
00801 
00802                                         // Analyze the channel data
00803                                         // Determine ring and column
00804                                         DayaBay::AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chanId, svcMode);
00805                                         int ring = pmtId.ring();
00806                                         int column = pmtId.column();
00807                                         int channel_number=ring*100+column;
00808 
00809 
00810                                         TH1F* tdcByMedianHWithCuts = 
00811                                         m_statsSvc->getTH1F( this->getPath(chanId) + "/tdcByMedianWithCuts");
00812                                         if(!tdcByMedianHWithCuts) return StatusCode::FAILURE;
00813 
00814                                         int chan_num=ring*24+column;
00815                                         for(int k=0;k<max_FEE_hit;k++)
00816                                         {
00817                                                 int tdc=tdc_buffer[chan_num%max_PMT_number][event_count][k];
00818                                                 int adc=adc_buffer[chan_num%max_PMT_number][event_count][k];
00819 
00820                                                 if(tdc>=m_TDC_Cut_Min_Channel[channel_number]&&tdc<=m_TDC_Cut_Max_Channel[channel_number])
00821                                                 if(adc>=m_ADC_Cut_Min_Channel[channel_number]&&adc<=m_ADC_Cut_Max_Channel[channel_number])
00822                                                 tdcByMedianHWithCuts->Fill(tdc-medianTdc);
00823                                         }
00824                                 } //over for for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++)
00825                                 tdcMedianHWithCuts->Fill(medianTdc);
00826                 } //over for if(readoutTdcWithCuts.size() > 0)
00827 
00828         //info()<<"test 3 event_count="<<event_count<<" adcSumWithCuts="<<adcSumWithCuts<<endreq;
00829         adcSumHWithCuts->Fill(adcSumWithCuts);
00830 
00831   } //over if(nReadouts<=max_event_number) and for(int event_count=0;event_count<nReadouts;event_count++)
00832 
00835 
00836  
00837   } // Loop over detectors
00838 
00839 
00840   return StatusCode::SUCCESS;
00841 }
00842 
00843 
00844 
00845 StatusCode PmtCalibLeadingEdgeWithCuts::calibrate(){
00846   // Loop over detectors in summary data, and calculate parameters
00847 
00848   // before the final calibrate of PMT
00849   // we need to fit the raw data without any cuts
00850   // and if the total nreadouts number less than max_event_number,
00851   // the data with cuts need to be fill firstly which will realized in the same program
00852   // 2009/11/29
00853   if(calibrateRaw()==StatusCode::FAILURE)
00854   {
00855         warning() << "calibrateRaw() error" 
00856                     << endreq;
00857         return StatusCode::FAILURE;
00858   }
00859 
00860 
00861   std::vector<DayaBay::Detector>::iterator detIter, 
00862     detEnd = m_processedDetectors.end();
00863   for(detIter = m_processedDetectors.begin(); detIter != detEnd; detIter++){
00864     DayaBay::Detector detector = *detIter;
00865 
00866     // Retrieve the data description and histograms
00867     TObject* nReadoutsObj = 
00868       m_statsSvc->get(this->getPath(detector) + "/nReadouts");
00869     TParameter<int>* nReadoutsPar= dynamic_cast<TParameter<int>*>(nReadoutsObj);
00870     if(!nReadoutsPar) return StatusCode::FAILURE;
00871 
00872     TH1F* meanOccupancyHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
00873                                                 +"/meanOccupancyWithCuts");
00874     if(!meanOccupancyHWithCuts) return StatusCode::FAILURE;
00875 
00876     TH1F* meanAdcHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
00877                                                 +"/meanAdcWithCuts");
00878     if(!meanAdcHWithCuts) return StatusCode::FAILURE;
00879 
00880 
00881     TH1F* gainChannelHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
00882                                                 +"/gainChannelWithCuts");
00883     if(!gainChannelHWithCuts) return StatusCode::FAILURE;
00884 
00885     TH1F* tdcOffsetChannelHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
00886                                                 +"/tdcOffsetChannelWithCuts");
00887     if(!tdcOffsetChannelHWithCuts) return StatusCode::FAILURE;
00888     TH1F* occupancyChannelHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
00889                                                 +"/occupancyChannelWithCuts");
00890     if(!occupancyChannelHWithCuts) return StatusCode::FAILURE;
00891 
00892 
00893     TH1F* meanTdcOffsetHWithCuts = m_statsSvc->getTH1F( this->getPath(detector)
00894                                                 +"/meanTdcOffsetWithCuts");
00895     if(!meanTdcOffsetHWithCuts) return StatusCode::FAILURE;
00896 
00897     TObject* simFlagObj = m_statsSvc->get(this->getPath(detector) +"/simFlag");
00898     TParameter<int>* simFlagPar = dynamic_cast<TParameter<int>*>(simFlagObj);
00899     if(!simFlagPar) return StatusCode::FAILURE;
00900 
00901     TObject* timeSecObj = m_statsSvc->get(this->getPath(detector) +"/timeSec");
00902     TParameter<int>* timeSecPar = dynamic_cast<TParameter<int>*>(timeSecObj);
00903     if(!timeSecPar) return StatusCode::FAILURE;
00904 
00905     TObject* timeNanoSecObj = 
00906       m_statsSvc->get(this->getPath(detector) +"/timeNanoSec");
00907     TParameter<int>* timeNanoSecPar = 
00908       dynamic_cast<TParameter<int>*>(timeNanoSecObj);
00909     if(!timeNanoSecPar) return StatusCode::FAILURE;
00910 
00911     Site::Site_t site = detector.site();
00912     DetectorId::DetectorId_t detId = detector.detectorId();
00913     SimFlag::SimFlag_t simFlag = (SimFlag::SimFlag_t)(simFlagPar->GetVal());
00914     time_t timeSec = (time_t)(timeSecPar->GetVal());
00915     int timeNanoSec = timeNanoSecPar->GetVal();
00916     int nReadouts = nReadoutsPar->GetVal();
00917 
00918     // Context, ServiceMode
00919     Context context(site,simFlag,TimeStamp(timeSec,timeNanoSec),detId);
00920     ServiceMode svcMode(context, 0);
00921 
00922     // Prepare data for ring-to-ring comparison
00923     int nRings = 8;
00924     std::map<int,double> meanOccRing;
00925     std::map<int,double> meanOccWeight;
00926     std::map<int,double> meanTdcRing;
00927     std::map<int,double> meanTdcWeight;
00928     for(int ring = 1; ring <= nRings; ring++){
00929       meanOccRing[ring] = 0.0;
00930       meanOccWeight[ring] = 0.0;
00931       meanTdcRing[ring] = 0.0;
00932       meanTdcWeight[ring] = 0.0;
00933     }
00934     double minValidOcc = 0.01;  // Cut on PMTs with bad/low occupancy
00935     std::vector<DayaBay::FeeChannelId> badChannels;
00936 
00937     // Get list of all FEE channels
00938     const std::vector<DayaBay::FeeChannelId>& channelList 
00939       = m_cableSvc->feeChannelIds( svcMode );
00940     std::vector<DayaBay::FeeChannelId>::const_iterator chanIter, 
00941       chanEnd = channelList.end();
00942     // Initialize statistics for each channel
00943     for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++){
00944       DayaBay::FeeChannelId chanId = *chanIter;
00945 
00946       // Analyze the channel data
00947       // Determine ring and column
00948       DayaBay::AdPmtSensor pmtId = 
00949         m_cableSvc->adPmtSensor(chanId, svcMode);
00950       int ring = pmtId.ring();
00951       int column = pmtId.column();
00952 
00953       // Channel status
00954       bool occupancyOk = true;
00955       bool gainOk = true;
00956 
00957       // Get Parameters and Histograms
00958       TObject* nHitsObjWithCuts = m_statsSvc->get(this->getPath(chanId) + "/nHitsWithCuts");
00959       TParameter<int>* nHitsParWithCuts = dynamic_cast<TParameter<int>*>(nHitsObjWithCuts);
00960       if(!nHitsParWithCuts) return StatusCode::FAILURE;
00961       TH1F* tdcByMedianHWithCuts= m_statsSvc->getTH1F( this->getPath(chanId)
00962                                                +"/tdcByMedianWithCuts");
00963       if(!tdcByMedianHWithCuts) return StatusCode::FAILURE;
00964       TH1F* adcHWithCuts = m_statsSvc->getTH1F( this->getPath(chanId) + "/adcWithCuts");
00965       if(!adcHWithCuts) return StatusCode::FAILURE;
00966       TH1F* occupancyHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
00967                                               +"/occupancyWithCuts");
00968       if(!occupancyHWithCuts) return StatusCode::FAILURE;
00969       TH1F* tdcOffsetHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
00970                                               +"/tdcOffsetWithCuts");
00971       if(!tdcOffsetHWithCuts) return StatusCode::FAILURE;
00972       TH1F* adcMedianHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
00973                                               +"/adcMedianWithCuts");
00974       if(!adcMedianHWithCuts) return StatusCode::FAILURE;
00975       TH1F* adcSigmaHWithCuts = m_statsSvc->getTH1F( this->getPath(detector,ring)
00976                                              +"/adcSigmaWithCuts");
00977       if(!adcSigmaHWithCuts) return StatusCode::FAILURE;
00978 
00979       // Channel Occupancy
00980       double occupancy = 0;
00981       double occupancyUncert = 0;
00982       {
00983         int nHits = nHitsParWithCuts->GetVal();
00984         if(nReadouts > 0){
00985           occupancy = nHits / ((double) nReadouts);
00986           occupancyUncert = 
00987             sqrt( occupancy*(1-occupancy)/((double) nReadouts) );
00988         }
00989         occupancyHWithCuts->SetBinContent(occupancyHWithCuts->FindBin(column),occupancy);
00990         occupancyHWithCuts->SetBinError(occupancyHWithCuts->FindBin(column),occupancyUncert);
00991         if(occupancy < minValidOcc) occupancyOk = false;
00992 
00993         occupancyChannelHWithCuts->SetBinContent(ring*24+column,occupancy);
00994         occupancyChannelHWithCuts->SetBinError(ring*24+column,occupancyUncert);
00995 
00996 
00997       }
00998       // TDC offset
00999       double tdcOffset = 0;
01000       double tdcOffsetUncert = 0;
01001       {
01002         // Use a fit to the leading edge to determine the time offset
01003         int maxBin = tdcByMedianHWithCuts->GetMaximumBin();
01004         double fitMin = tdcByMedianHWithCuts->GetBinCenter(maxBin - 2);
01005         double fitMax = tdcByMedianHWithCuts->GetBinCenter(maxBin + 10);
01006         tdcByMedianHWithCuts->Fit("gaus","","",fitMin, fitMax);
01007         TF1* fitFunc = tdcByMedianHWithCuts->GetFunction("gaus");
01008         if( fitFunc ){
01009           tdcOffset = fitFunc->GetParameter(1);
01010           tdcOffsetUncert = fitFunc->GetParError(1);
01011           
01012           if( fitFunc->GetNDF() < 3 ) 
01013             // Catch bad fits
01014             tdcOffsetUncert = 0;
01015           else
01016             // Scale uncertainty by the quality of the fit
01017             tdcOffsetUncert *= sqrt(fitFunc->GetChisquare() / fitFunc->GetNDF());
01018           
01019           tdcOffsetHWithCuts->SetBinContent(tdcOffsetHWithCuts->FindBin(column),tdcOffset);
01020           tdcOffsetHWithCuts->SetBinError(tdcOffsetHWithCuts->FindBin(column),tdcOffsetUncert);
01021 
01022         tdcOffsetChannelHWithCuts->SetBinContent(ring*24+column,tdcOffset);
01023         tdcOffsetChannelHWithCuts->SetBinError(ring*24+column,tdcOffsetUncert);
01024 
01025 
01026         }else{
01027           warning() << "PMT  " << pmtId << " has no TDC data for fitting." 
01028                     << endreq;
01029         }
01030       }
01031 
01032       // ADC median and sigma
01033       double adcMean = 0;
01034       double adcMeanUncert = 0;
01035       double adcSigma = 0;
01036       double adcSigmaUncert = 0;
01037       {
01038         // Find ADC SPE peak, width with simple gaussian
01039         TF1 adcFit("adcFit","gaus");
01040         double adcPeak = adcHWithCuts->GetBinCenter(adcHWithCuts->GetMaximumBin());
01041         double adcSig = 5.0;
01042         adcHWithCuts->Fit(&adcFit,"","",adcPeak - adcSig, adcPeak + adcSig);
01043         adcMean = adcFit.GetParameter(1);
01044         if(fabs(adcMean-adcPeak) > adcSig){
01045           warning() << "Issues with fit on ring, column: " << ring << " / " 
01046                     << column << " mean/peak: " << adcMean << " / " << adcPeak 
01047                     << endreq;
01048           gainOk = false;
01049         }
01050         adcMeanUncert = adcFit.GetParError(1);
01051         adcSigma = adcFit.GetParameter(2);
01052         adcSigmaUncert = adcFit.GetParError(2);
01053         adcMedianHWithCuts->SetBinContent(adcMedianHWithCuts->FindBin(column),adcMean);
01054         adcMedianHWithCuts->SetBinError(adcMedianHWithCuts->FindBin(column),adcMeanUncert);
01055         adcSigmaHWithCuts->SetBinContent(adcSigmaHWithCuts->FindBin(column),adcSigma);
01056         adcSigmaHWithCuts->SetBinError(adcSigmaHWithCuts->FindBin(column),adcSigmaUncert);
01057 
01058         meanAdcHWithCuts->SetBinContent(ring,adcMean);
01059         meanAdcHWithCuts->SetBinError(ring, adcMeanUncert);
01060 
01061         gainChannelHWithCuts->SetBinContent(ring*24+column,adcMean);
01062         gainChannelHWithCuts->SetBinError(ring*24+column,adcMeanUncert);
01063 
01064       }
01065 
01066       // Record ring averages, ignoring bad channels
01067       if( occupancyOk && gainOk ){
01068         if( occupancyUncert > 0 ){
01069           meanOccRing[ring] += occupancy / (occupancyUncert*occupancyUncert); 
01070           meanOccWeight[ring] += 1.0 / (occupancyUncert*occupancyUncert); 
01071         }
01072         if( tdcOffsetUncert > 0 ){
01073           meanTdcRing[ring] += tdcOffset / (tdcOffsetUncert*tdcOffsetUncert); 
01074           meanTdcWeight[ring] += 1.0 / (tdcOffsetUncert*tdcOffsetUncert); 
01075         }
01076       }
01077 
01078     }
01079 
01080     // Record mean time offset and mean occupancy by ring
01081     for(int ring = 1; ring <= nRings; ring++){
01082       if( meanOccWeight[ring] > 0 ){
01083         meanOccupancyHWithCuts->SetBinContent(ring,
01084                                       meanOccRing[ring] / meanOccWeight[ring]);
01085         meanOccupancyHWithCuts->SetBinError(ring, 
01086                                     sqrt( 1.0 / meanOccWeight[ring] ));
01087       }
01088       if( meanTdcWeight[ring] > 0 ){
01089         meanTdcOffsetHWithCuts->SetBinContent(ring,
01090                                       meanTdcRing[ring] / meanTdcWeight[ring]);
01091         meanTdcOffsetHWithCuts->SetBinError(ring, 
01092                                     sqrt( 1.0 / meanTdcWeight[ring] ));
01093       }
01094     } // Loop over rings
01095  
01096   } // Loop over detectors
01097 
01098 info() << "test second calibrate over" << endreq;
01099 
01100   return StatusCode::SUCCESS;
01101 
01102 }
01105 
01106 
01108 bool PmtCalibLeadingEdgeWithCuts::hasStats(const DayaBay::Detector& detector){
01109   // Check if statistics have been initialized for this detector
01110   return (std::find(m_processedDetectors.begin(), 
01111                     m_processedDetectors.end(), 
01112                     detector) 
01113           != m_processedDetectors.end());
01114 }
01115 
01116 StatusCode PmtCalibLeadingEdgeWithCuts::prepareStats(const Context& context){
01117   // Create the histograms and parameters for this detector's statistics
01118   DayaBay::Detector detector(context.GetSite(), context.GetDetId());
01119 
01120   // Site
01121   {
01122     std::string name = "site";
01123     std::ostringstream path;
01124     path << this->getPath(detector) << "/" << name;
01125     TParameter<int>* par = new TParameter<int>(name.c_str(), 
01126                                                context.GetSite());
01127     if( m_statsSvc->put(path.str(),par).isFailure() ) {
01128       error() << "prepareStats(): Could not register " << name 
01129               << " at " << path << endreq;
01130       delete par;
01131       par = 0;
01132       return StatusCode::FAILURE;
01133     }
01134   }
01135 
01136   // Detector ID
01137   {
01138     std::string name = "detectorId";
01139     std::ostringstream path;
01140     path << this->getPath(detector) << "/" << name;
01141     TParameter<int>* par = new TParameter<int>(name.c_str(), 
01142                                                context.GetDetId());
01143     if( m_statsSvc->put(path.str(),par).isFailure() ) {
01144       error() << "prepareStats(): Could not register " << name 
01145               << " at " << path << endreq;
01146       delete par;
01147       par = 0;
01148       return StatusCode::FAILURE;
01149     }
01150   }
01151 
01152   // SimFlag
01153   {
01154     std::string name = "simFlag";
01155     std::ostringstream path;
01156     path << this->getPath(detector) << "/" << name;
01157     TParameter<int>* par = new TParameter<int>(name.c_str(), 
01158                                                context.GetSimFlag());
01159     if( m_statsSvc->put(path.str(),par).isFailure() ) {
01160       error() << "prepareStats(): Could not register " << name 
01161               << " at " << path << endreq;
01162       delete par;
01163       par = 0;
01164       return StatusCode::FAILURE;
01165     }
01166   }
01167 
01168   // Timestamp: seconds
01169   {
01170     std::string name = "timeSec";
01171     std::ostringstream path;
01172     path << this->getPath(detector) << "/" << name;
01173     TParameter<int>* par = new TParameter<int>(name.c_str(), 
01174                                                context.GetTimeStamp().GetSec());
01175     if( m_statsSvc->put(path.str(),par).isFailure() ) {
01176       error() << "prepareStats(): Could not register " << name 
01177               << " at " << path << endreq;
01178       delete par;
01179       par = 0;
01180       return StatusCode::FAILURE;
01181     }
01182   }
01183 
01184   // Timestamp: nanoseconds
01185   {
01186     std::string name = "timeNanoSec";
01187     std::ostringstream path;
01188     path << this->getPath(detector) << "/" << name;
01189     TParameter<int>* par = 
01190       new TParameter<int>(name.c_str(), context.GetTimeStamp().GetNanoSec());
01191     if( m_statsSvc->put(path.str(),par).isFailure() ) {
01192       error() << "prepareStats(): Could not register " << name 
01193               << " at " << path << endreq;
01194       delete par;
01195       par = 0;
01196       return StatusCode::FAILURE;
01197     }
01198   }
01199 
01200   // Number of Readouts
01201   {
01202     std::string name = "nReadouts";
01203     std::ostringstream path;
01204     path << this->getPath(detector) << "/" << name;
01205     TParameter<int>* par = new TParameter<int>(name.c_str(), 0);
01206     if( m_statsSvc->put(path.str(),par).isFailure() ) {
01207       error() << "prepareStats(): Could not register " << name 
01208               << " at " << path << endreq;
01209       delete par;
01210       par = 0;
01211       return StatusCode::FAILURE;
01212     }
01213   }
01214 
01215   // ADC Sum spectrum
01216   {
01217     std::ostringstream title, path;
01218     std::string name = "adcSum";
01219     title << "ADC Sum for readouts in " << detector.detName(); 
01220     path << this->getPath(detector) << "/" << name;
01221     TH1F* adcSum = new TH1F(name.c_str(),title.str().c_str(),
01222                             2000,0,20000);
01223     adcSum->GetXaxis()->SetTitle("ADC Sum value");
01224     adcSum->GetYaxis()->SetTitle("Entries");
01225     // DEBUG
01226     info() << "name= " << adcSum->GetName() 
01227            << " at path = \"" << path.str() << "\"" << endreq;
01228     if( m_statsSvc->put(path.str(),adcSum).isFailure() ) {
01229       error() << "prepareStats(): Could not register " << path << endreq;
01230       delete adcSum;
01231       return StatusCode::FAILURE;
01232     }
01233   }
01234 
01235   // ADC Sum spectrum with cuts
01236   {
01237     std::ostringstream title, path;
01238     std::string name = "adcSumWithCuts";
01239     title << "ADC Sum for readouts in " << detector.detName(); 
01240     path << this->getPath(detector) << "/" << name;
01241     TH1F* adcSumWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01242                             2000,0,20000);
01243     adcSumWithCuts->GetXaxis()->SetTitle("ADC Sum value With Cuts");
01244     adcSumWithCuts->GetYaxis()->SetTitle("Entries");
01245     // DEBUG
01246     info() << "name= " << adcSumWithCuts->GetName() 
01247            << " at path = \"" << path.str() << "\"" << endreq;
01248     if( m_statsSvc->put(path.str(),adcSumWithCuts).isFailure() ) {
01249       error() << "prepareStats(): Could not register " << path << endreq;
01250       delete adcSumWithCuts;
01251       return StatusCode::FAILURE;
01252     }
01253   }
01254 
01255 
01256   // Gain vs channel spectrum
01257   {
01258     std::ostringstream title, path;
01259     std::string name = "gainChannel";
01260     title << "gainChannel for readouts in " << detector.detName(); 
01261     path << this->getPath(detector) << "/" << name;
01262     TH1F* gainChannel = new TH1F(name.c_str(),title.str().c_str(),
01263                             300,0,300);
01264     gainChannel->GetXaxis()->SetTitle("Channel");
01265     gainChannel->GetYaxis()->SetTitle("gain in ADC bin");
01266     // DEBUG
01267     info() << "name= " << gainChannel->GetName() 
01268            << " at path = \"" << path.str() << "\"" << endreq;
01269     if( m_statsSvc->put(path.str(),gainChannel).isFailure() ) {
01270       error() << "prepareStats(): Could not register " << path << endreq;
01271       delete gainChannel;
01272       return StatusCode::FAILURE;
01273     }
01274   }
01275 
01276   // gain vs channel with cuts
01277   {
01278     std::ostringstream title, path;
01279     std::string name = "gainChannelWithCuts";
01280     title << "gainChannel for readouts in " << detector.detName(); 
01281     path << this->getPath(detector) << "/" << name;
01282     TH1F* gainChannelWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01283                             300,0,300);
01284     gainChannelWithCuts->GetXaxis()->SetTitle("channel");
01285     gainChannelWithCuts->GetYaxis()->SetTitle("gain in ADC bin");
01286     // DEBUG
01287     info() << "name= " << gainChannelWithCuts->GetName() 
01288            << " at path = \"" << path.str() << "\"" << endreq;
01289     if( m_statsSvc->put(path.str(),gainChannelWithCuts).isFailure() ) {
01290       error() << "prepareStats(): Could not register " << path << endreq;
01291       delete gainChannelWithCuts;
01292       return StatusCode::FAILURE;
01293     }
01294   }
01295 
01296   // tdc offset vs channel spectrum
01297   {
01298     std::ostringstream title, path;
01299     std::string name = "tdcOffsetChannel";
01300     title << "tdcOffsetChannel for readouts in " << detector.detName(); 
01301     path << this->getPath(detector) << "/" << name;
01302     TH1F* tdcOffsetChannel = new TH1F(name.c_str(),title.str().c_str(),
01303                             300,0,300);
01304     tdcOffsetChannel->GetXaxis()->SetTitle("Channel");
01305     tdcOffsetChannel->GetYaxis()->SetTitle("tdc offset in TDC bin");
01306     // DEBUG
01307     info() << "name= " << tdcOffsetChannel->GetName() 
01308            << " at path = \"" << path.str() << "\"" << endreq;
01309     if( m_statsSvc->put(path.str(),tdcOffsetChannel).isFailure() ) {
01310       error() << "prepareStats(): Could not register " << path << endreq;
01311       delete tdcOffsetChannel;
01312       return StatusCode::FAILURE;
01313     }
01314   }
01315 
01316   // tdcOffset vs channel with cuts
01317   {
01318     std::ostringstream title, path;
01319     std::string name = "tdcOffsetChannelWithCuts";
01320     title << "tdcOffsetChannel for readouts in " << detector.detName(); 
01321     path << this->getPath(detector) << "/" << name;
01322     TH1F* tdcOffsetChannelWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01323                             300,0,300);
01324     tdcOffsetChannelWithCuts->GetXaxis()->SetTitle("channel");
01325     tdcOffsetChannelWithCuts->GetYaxis()->SetTitle("tdc offset in TDC bin");
01326     // DEBUG
01327     info() << "name= " << tdcOffsetChannelWithCuts->GetName() 
01328            << " at path = \"" << path.str() << "\"" << endreq;
01329     if( m_statsSvc->put(path.str(),tdcOffsetChannelWithCuts).isFailure() ) {
01330       error() << "prepareStats(): Could not register " << path << endreq;
01331       delete tdcOffsetChannelWithCuts;
01332       return StatusCode::FAILURE;
01333     }
01334   }
01335 
01336 
01337   // occupancy vs channel spectrum
01338   {
01339     std::ostringstream title, path;
01340     std::string name = "occupancyChannel";
01341     title << "occupancyChannel for readouts in " << detector.detName(); 
01342     path << this->getPath(detector) << "/" << name;
01343     TH1F* occupancyChannel = new TH1F(name.c_str(),title.str().c_str(),
01344                             300,0,300);
01345     occupancyChannel->GetXaxis()->SetTitle("Channel");
01346     occupancyChannel->GetYaxis()->SetTitle("occupancy");
01347     // DEBUG
01348     info() << "name= " << occupancyChannel->GetName() 
01349            << " at path = \"" << path.str() << "\"" << endreq;
01350     if( m_statsSvc->put(path.str(),occupancyChannel).isFailure() ) {
01351       error() << "prepareStats(): Could not register " << path << endreq;
01352       delete occupancyChannel;
01353       return StatusCode::FAILURE;
01354     }
01355   }
01356 
01357   // occupancy vs channel with cuts
01358   {
01359     std::ostringstream title, path;
01360     std::string name = "occupancyChannelWithCuts";
01361     title << "occupancyChannel for readouts in " << detector.detName(); 
01362     path << this->getPath(detector) << "/" << name;
01363     TH1F* occupancyChannelWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01364                             300,0,300);
01365     occupancyChannelWithCuts->GetXaxis()->SetTitle("channel");
01366     occupancyChannelWithCuts->GetYaxis()->SetTitle("occupancy");
01367     // DEBUG
01368     info() << "name= " << occupancyChannelWithCuts->GetName() 
01369            << " at path = \"" << path.str() << "\"" << endreq;
01370     if( m_statsSvc->put(path.str(),occupancyChannelWithCuts).isFailure() ) {
01371       error() << "prepareStats(): Could not register " << path << endreq;
01372       delete occupancyChannelWithCuts;
01373       return StatusCode::FAILURE;
01374     }
01375   }
01376 
01377 
01378   // TDC Median spectrum
01379   {
01380     std::ostringstream title, path;
01381     std::string name = "tdcMedian";
01382     title << "Median TDC for readouts in " << detector.detName(); 
01383     path << this->getPath(detector) << "/" << name;
01384     TH1F* tdcMedian = new TH1F(name.c_str(),title.str().c_str(),
01385                                800,0,800);
01386     tdcMedian->GetXaxis()->SetTitle("TDC value");
01387     tdcMedian->GetYaxis()->SetTitle("Entries");
01388     if( m_statsSvc->put(path.str(),tdcMedian).isFailure() ) {
01389       error() << "prepareStats(): Could not register " << path << endreq;
01390       delete tdcMedian;
01391       return StatusCode::FAILURE;
01392     }
01393   }
01394 
01395   // TDC Median spectrum with cuts
01396   {
01397     std::ostringstream title, path;
01398     std::string name = "tdcMedianWithCuts";
01399     title << "Median TDC WithCuts for readouts in " << detector.detName(); 
01400     path << this->getPath(detector) << "/" << name;
01401     TH1F* tdcMedianWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01402                                800,0,800);
01403     tdcMedianWithCuts->GetXaxis()->SetTitle("TDC value");
01404     tdcMedianWithCuts->GetYaxis()->SetTitle("Entries");
01405     if( m_statsSvc->put(path.str(),tdcMedianWithCuts).isFailure() ) {
01406       error() << "prepareStats(): Could not register " << path << endreq;
01407       delete tdcMedianWithCuts;
01408       return StatusCode::FAILURE;
01409     }
01410   }
01411 
01412   // Mean TDC Offset by AD ring
01413   {
01414     std::ostringstream title, path;
01415     std::string name = "meanTdcOffset";
01416     title << "Mean TDC offset by ring in " << detector.detName(); 
01417     path << this->getPath(detector) << "/" << name;
01418     TH1F* meanTdcOffset = new TH1F(name.c_str(),title.str().c_str(),
01419                                    8,1,9);
01420     meanTdcOffset->GetXaxis()->SetTitle("AD Ring");
01421     meanTdcOffset->GetYaxis()->SetTitle("Mean TDC Offset");
01422     if( m_statsSvc->put(path.str(),meanTdcOffset).isFailure() ) {
01423       error() << "prepareStats(): Could not register " << path << endreq;
01424       delete meanTdcOffset;
01425       return StatusCode::FAILURE;
01426     }
01427   }
01428 
01429   // Mean TDC Offset WithCuts by AD ring
01430   {
01431     std::ostringstream title, path;
01432     std::string name = "meanTdcOffsetWithCuts";
01433     title << "Mean TDC offset by ring in " << detector.detName(); 
01434     path << this->getPath(detector) << "/" << name;
01435     TH1F* meanTdcOffsetWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01436                                    8,1,9);
01437     meanTdcOffsetWithCuts->GetXaxis()->SetTitle("AD Ring");
01438     meanTdcOffsetWithCuts->GetYaxis()->SetTitle("Mean TDC Offset With Cuts");
01439     if( m_statsSvc->put(path.str(),meanTdcOffsetWithCuts).isFailure() ) {
01440       error() << "prepareStats(): Could not register " << path << endreq;
01441       delete meanTdcOffsetWithCuts;
01442       return StatusCode::FAILURE;
01443     }
01444   }
01445 
01446   // Mean Occupancy by AD ring
01447   {
01448     std::ostringstream title, path;
01449     std::string name = "meanOccupancy";
01450     title << "Mean occupancy by ring in " << detector.detName(); 
01451     path << this->getPath(detector) << "/" << name;
01452     TH1F* meanOccupancy = new TH1F(name.c_str(),title.str().c_str(),
01453                                    8,1,9);
01454     meanOccupancy->GetXaxis()->SetTitle("AD Ring");
01455     meanOccupancy->GetYaxis()->SetTitle("Mean Occupancy");
01456     if( m_statsSvc->put(path.str(),meanOccupancy).isFailure() ) {
01457       error() << "prepareStats(): Could not register " << path << endreq;
01458       delete meanOccupancy;
01459       return StatusCode::FAILURE;
01460     }
01461   }
01462 
01463   // Mean Adc by AD ring
01464   {
01465     std::ostringstream title, path;
01466     std::string name = "meanAdc";
01467     title << "Mean Adc by ring in " << detector.detName(); 
01468     path << this->getPath(detector) << "/" << name;
01469     TH1F* meanAdc = new TH1F(name.c_str(),title.str().c_str(),
01470                                    8,1,9);
01471     meanAdc->GetXaxis()->SetTitle("AD Ring");
01472     meanAdc->GetYaxis()->SetTitle("Mean Adc");
01473     if( m_statsSvc->put(path.str(),meanAdc).isFailure() ) {
01474       error() << "prepareStats(): Could not register " << path << endreq;
01475       delete meanAdc;
01476       return StatusCode::FAILURE;
01477     }
01478   }
01479 
01480   // Mean Occupancy With Cuts by AD ring 
01481   {
01482     std::ostringstream title, path;
01483     std::string name = "meanOccupancyWithCuts";
01484     title << "Mean occupancy With Cuts by ring in " << detector.detName(); 
01485     path << this->getPath(detector) << "/" << name;
01486     TH1F* meanOccupancyWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01487                                    8,1,9);
01488     meanOccupancyWithCuts->GetXaxis()->SetTitle("AD Ring");
01489     meanOccupancyWithCuts->GetYaxis()->SetTitle("Mean Occupancy With Cuts");
01490     if( m_statsSvc->put(path.str(),meanOccupancyWithCuts).isFailure() ) {
01491       error() << "prepareStats(): Could not register " << path << endreq;
01492       delete meanOccupancyWithCuts;
01493       return StatusCode::FAILURE;
01494     }
01495   }
01496 
01497   // Mean Adc With Cuts by AD ring 
01498   {
01499     std::ostringstream title, path;
01500     std::string name = "meanAdcWithCuts";
01501     title << "Mean Adc With Cuts by ring in " << detector.detName(); 
01502     path << this->getPath(detector) << "/" << name;
01503     TH1F* meanAdcWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01504                                    8,1,9);
01505     meanAdcWithCuts->GetXaxis()->SetTitle("AD Ring");
01506     meanAdcWithCuts->GetYaxis()->SetTitle("Mean Adc With Cuts");
01507     if( m_statsSvc->put(path.str(),meanAdcWithCuts).isFailure() ) {
01508       error() << "prepareStats(): Could not register " << path << endreq;
01509       delete meanAdcWithCuts;
01510       return StatusCode::FAILURE;
01511     }
01512   }
01513 
01514 
01515   // Make sure summary histograms for each ring exist in detector statistics
01516   for(int ring = 0; ring <= 8; ring++){
01517     // Occupancy by ring
01518     {
01519       std::ostringstream title, path;
01520       std::string name = "occupancy";
01521       title << "Occupancy for PMTs in " << detector.detName() 
01522             << " ring " << ring; 
01523       path << this->getPath(detector, ring) << "/" << name;
01524       TH1F* occupancy = new TH1F(name.c_str(),title.str().c_str(),
01525                                  24,1,25);
01526       occupancy->GetXaxis()->SetTitle("Column");
01527       occupancy->GetYaxis()->SetTitle("Occupancy");
01528       if( m_statsSvc->put(path.str(),occupancy).isFailure() ) {
01529         error() << "prepareStats(): Could not register " << path << endreq;
01530         delete occupancy;
01531         return StatusCode::FAILURE;
01532       }
01533     }
01534     // TDC offset by ring
01535     {
01536       std::ostringstream title, path;
01537       std::string name = "tdcOffset";
01538       title << "Time Offset for PMTs in " << detector.detName() 
01539             << " ring " << ring; 
01540       path << this->getPath(detector, ring) << "/" << name;
01541       TH1F* tdcOffset = new TH1F(name.c_str(),title.str().c_str(),
01542                            24,1,25);
01543       tdcOffset->GetXaxis()->SetTitle("Column");
01544       tdcOffset->GetYaxis()->SetTitle("Time Offset [tdc]");
01545       if( m_statsSvc->put(path.str(),tdcOffset).isFailure() ) {
01546         error() << "prepareStats(): Could not register " << path << endreq;
01547         delete tdcOffset;
01548         return StatusCode::FAILURE;
01549       }
01550     }
01551     // TDC Raw offset by ring
01552     {
01553       std::ostringstream title, path;
01554       std::string name = "tdcRawOffset";
01555       title << "Time Raw Offset for PMTs in " << detector.detName() 
01556             << " ring " << ring; 
01557       path << this->getPath(detector, ring) << "/" << name;
01558       TH1F* tdcRawOffset = new TH1F(name.c_str(),title.str().c_str(),
01559                            24,1,25);
01560       tdcRawOffset->GetXaxis()->SetTitle("Column");
01561       tdcRawOffset->GetYaxis()->SetTitle("Time Offset [tdc]");
01562       if( m_statsSvc->put(path.str(),tdcRawOffset).isFailure() ) {
01563         error() << "prepareStats(): Could not register " << path << endreq;
01564         delete tdcRawOffset;
01565         return StatusCode::FAILURE;
01566       }
01567     }
01568 
01569     // ADC Median by ring
01570     {
01571       std::ostringstream title, path;
01572       std::string name = "adcMedian";
01573       title << "ADC Median for PMTs in " << detector.detName() 
01574             << " ring " << ring; 
01575       path << this->getPath(detector, ring) << "/" << name;
01576       TH1F* adcMedian = new TH1F(name.c_str(),title.str().c_str(),
01577                                  24,1,25);
01578       adcMedian->GetXaxis()->SetTitle("Column");
01579       adcMedian->GetYaxis()->SetTitle("Median ADC");
01580       if( m_statsSvc->put(path.str(),adcMedian).isFailure() ) {
01581         error() << "prepareStats(): Could not register " << path << endreq;
01582         delete adcMedian;
01583         return StatusCode::FAILURE;
01584       }
01585     }
01586     // ADC Sigma by ring
01587     {
01588       std::ostringstream title, path;
01589       std::string name = "adcSigma";
01590       title << "ADC Sigma for PMTs in " << detector.detName() 
01591             << " ring " << ring;
01592       path << this->getPath(detector, ring) << "/" << name;
01593       TH1F* adcSigma = new TH1F(name.c_str(),title.str().c_str(),
01594                                 24,1,25);
01595       adcSigma->GetXaxis()->SetTitle("Column");
01596       adcSigma->GetYaxis()->SetTitle("ADC Sigma");
01597       if( m_statsSvc->put(path.str(),adcSigma).isFailure() ) {
01598         error() << "prepareStats(): Could not register " << path << endreq;
01599         delete adcSigma;
01600         return StatusCode::FAILURE;
01601       }
01602     }
01603   }
01604 
01605   // Make sure summary histograms for each ring exist in detector statistics with cuts
01606   for(int ring = 0; ring <= 8; ring++){
01607     // Occupancy by ring With Cuts
01608     {
01609       std::ostringstream title, path;
01610       std::string name = "occupancyWithCuts";
01611       title << "Occupancy With Cuts for PMTs in " << detector.detName() 
01612             << " ring " << ring; 
01613       path << this->getPath(detector, ring) << "/" << name;
01614       TH1F* occupancyWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01615                                  24,1,25);
01616       occupancyWithCuts->GetXaxis()->SetTitle("Column");
01617       occupancyWithCuts->GetYaxis()->SetTitle("Occupancy");
01618       if( m_statsSvc->put(path.str(),occupancyWithCuts).isFailure() ) {
01619         error() << "prepareStats(): Could not register " << path << endreq;
01620         delete occupancyWithCuts;
01621         return StatusCode::FAILURE;
01622       }
01623     }
01624     // TDC offset by ring With Cuts
01625     {
01626       std::ostringstream title, path;
01627       std::string name = "tdcOffsetWithCuts";
01628       title << "Time Offset With Cuts for PMTs in " << detector.detName() 
01629             << " ring " << ring; 
01630       path << this->getPath(detector, ring) << "/" << name;
01631       TH1F* tdcOffsetWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01632                            24,1,25);
01633       tdcOffsetWithCuts->GetXaxis()->SetTitle("Column");
01634       tdcOffsetWithCuts->GetYaxis()->SetTitle("Time Offset With Cuts [tdc]");
01635       if( m_statsSvc->put(path.str(),tdcOffsetWithCuts).isFailure() ) {
01636         error() << "prepareStats(): Could not register " << path << endreq;
01637         delete tdcOffsetWithCuts;
01638         return StatusCode::FAILURE;
01639       }
01640     }
01641     // ADC Median by ring With Cuts
01642     {
01643       std::ostringstream title, path;
01644       std::string name = "adcMedianWithCuts";
01645       title << "ADC Median With Cuts for PMTs in " << detector.detName() 
01646             << " ring " << ring; 
01647       path << this->getPath(detector, ring) << "/" << name;
01648       TH1F* adcMedianWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01649                                  24,1,25);
01650       adcMedianWithCuts->GetXaxis()->SetTitle("Column");
01651       adcMedianWithCuts->GetYaxis()->SetTitle("Median ADC With Cuts");
01652       if( m_statsSvc->put(path.str(),adcMedianWithCuts).isFailure() ) {
01653         error() << "prepareStats(): Could not register " << path << endreq;
01654         delete adcMedianWithCuts;
01655         return StatusCode::FAILURE;
01656       }
01657     }
01658     // ADC Sigma by ring With Cuts
01659     {
01660       std::ostringstream title, path;
01661       std::string name = "adcSigmaWithCuts";
01662       title << "ADC Sigma With Cuts for PMTs in " << detector.detName() 
01663             << " ring " << ring;
01664       path << this->getPath(detector, ring) << "/" << name;
01665       TH1F* adcSigmaWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01666                                 24,1,25);
01667       adcSigmaWithCuts->GetXaxis()->SetTitle("Column");
01668       adcSigmaWithCuts->GetYaxis()->SetTitle("ADC Sigma With Cuts");
01669       if( m_statsSvc->put(path.str(),adcSigmaWithCuts).isFailure() ) {
01670         error() << "prepareStats(): Could not register " << path << endreq;
01671         delete adcSigmaWithCuts;
01672         return StatusCode::FAILURE;
01673       }
01674     }
01675   }
01676 
01677 
01678   ServiceMode svcMode(context, 0);
01679 
01680   // Get list of all FEE channels
01681   const std::vector<DayaBay::FeeChannelId>& channelList 
01682     = m_cableSvc->feeChannelIds( svcMode );
01683   std::vector<DayaBay::FeeChannelId>::const_iterator chanIter, 
01684     chanEnd = channelList.end();
01685   // Initialize statistics for each channel
01686   for(chanIter = channelList.begin(); chanIter != chanEnd; chanIter++){
01687     DayaBay::FeeChannelId chanId = *chanIter;
01688 
01689     // Board Number
01690     {
01691       std::string name = "board";
01692       std::ostringstream path;
01693       path << this->getPath(chanId) << "/" << name;
01694       TParameter<int>* par = new TParameter<int>(name.c_str(), chanId.board());
01695       if( m_statsSvc->put(path.str(),par).isFailure() ) {
01696         error() << "prepareStats(): Could not register " << name 
01697                 << " at " << path << endreq;
01698         delete par;
01699         par = 0;
01700         return StatusCode::FAILURE;
01701       }
01702     }
01703 
01704     // Connector Number
01705     {
01706       std::string name = "connector";
01707       std::ostringstream path;
01708       path << this->getPath(chanId) << "/" << name;
01709       TParameter<int>* par = new TParameter<int>(name.c_str(), 
01710                                                  chanId.connector());
01711       if( m_statsSvc->put(path.str(),par).isFailure() ) {
01712         error() << "prepareStats(): Could not register " << name 
01713                 << " at " << path << endreq;
01714         delete par;
01715         par = 0;
01716         return StatusCode::FAILURE;
01717       }
01718     }
01719 
01720     // Number of Hits
01721     {
01722       std::string name = "nHits";
01723       std::ostringstream path;
01724       path << this->getPath(chanId) << "/" << name;
01725       TParameter<int>* par = new TParameter<int>(name.c_str(), 0);
01726       if( m_statsSvc->put(path.str(),par).isFailure() ) {
01727         error() << "prepareStats(): Could not register " << name 
01728                 << " at " << path << endreq;
01729         delete par;
01730         par = 0;
01731         return StatusCode::FAILURE;
01732       }
01733     }
01734 
01735     // Number of Hits with cuts
01736     {
01737       std::string name = "nHitsWithCuts";
01738       std::ostringstream path;
01739       path << this->getPath(chanId) << "/" << name;
01740       TParameter<int>* parWithCuts = new TParameter<int>(name.c_str(), 0);
01741       if( m_statsSvc->put(path.str(),parWithCuts).isFailure() ) {
01742         error() << "prepareStats(): Could not register " << name 
01743                 << " at " << path << endreq;
01744         delete parWithCuts;
01745         parWithCuts = 0;
01746         return StatusCode::FAILURE;
01747       }
01748     }
01749 
01750 
01751     // Raw TDC spectrum
01752     {
01753       std::string name = "tdcRaw";
01754       std::ostringstream title, path;
01755       title << "Raw TDC values for " << detector.detName() 
01756             << " channel " << chanId.board() << "_"
01757             << chanId.connector();
01758       path << this->getPath(chanId) << "/" << name;
01759       TH1F* tdcRaw = new TH1F(name.c_str(),title.str().c_str(),
01760                               300,0,300);
01761       tdcRaw->GetXaxis()->SetTitle("TDC value");
01762       tdcRaw->GetYaxis()->SetTitle("Entries");
01763       if( m_statsSvc->put(path.str(),tdcRaw).isFailure() ) {
01764         error() << "prepareStats(): Could not register " << path << endreq;
01765         delete tdcRaw;
01766         return StatusCode::FAILURE;
01767       }
01768     }
01769 
01770     // TDC spectrum with cuts
01771     {
01772       std::string name = "tdcWithCuts";
01773       std::ostringstream title, path;
01774       title << "TDC values With Cuts for " << detector.detName() 
01775             << " channel " << chanId.board() << "_"
01776             << chanId.connector();
01777       path << this->getPath(chanId) << "/" << name;
01778       TH1F* tdcWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01779                               800,0,800);
01780       tdcWithCuts->GetXaxis()->SetTitle("TDC value");
01781       tdcWithCuts->GetYaxis()->SetTitle("Entries With Cuts");
01782       if( m_statsSvc->put(path.str(),tdcWithCuts).isFailure() ) {
01783         error() << "prepareStats(): Could not register " << path << endreq;
01784         delete tdcWithCuts;
01785         return StatusCode::FAILURE;
01786       }
01787     }
01788 
01789 
01790     // TDC spectrum, median corrected
01791     {
01792       std::string name = "tdcByMedian";
01793       std::ostringstream title, path;
01794       title << "Median-corrected TDC values for " << detector.detName() 
01795             << " channel " << chanId.board() << "_"
01796             << chanId.connector();
01797       path << this->getPath(chanId) << "/" << name;
01798       TH1F* tdcByMedian = new TH1F(name.c_str(),title.str().c_str(),
01799                                    600,-300,300);
01800       tdcByMedian->GetXaxis()->SetTitle("TDC value");
01801       tdcByMedian->GetYaxis()->SetTitle("Entries");
01802       if( m_statsSvc->put(path.str(),tdcByMedian).isFailure() ) {
01803         error() << "prepareStats(): Could not register " << path << endreq;
01804         delete tdcByMedian;
01805         return StatusCode::FAILURE;
01806       }
01807     }
01808 
01809     // TDC spectrum, median corrected with cuts
01810     {
01811       std::string name = "tdcByMedianWithCuts";
01812       std::ostringstream title, path;
01813       title << "Median-corrected TDC values With Cuts for " << detector.detName() 
01814             << " channel " << chanId.board() << "_"
01815             << chanId.connector();
01816       path << this->getPath(chanId) << "/" << name;
01817       TH1F* tdcByMedianWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01818                                    600,-300,300);
01819       tdcByMedianWithCuts->GetXaxis()->SetTitle("TDC value");
01820       tdcByMedianWithCuts->GetYaxis()->SetTitle("Entries With Cuts");
01821       if( m_statsSvc->put(path.str(),tdcByMedianWithCuts).isFailure() ) {
01822         error() << "prepareStats(): Could not register " << path << endreq;
01823         delete tdcByMedianWithCuts;
01824         return StatusCode::FAILURE;
01825       }
01826     }
01827 
01828 
01829     // Raw ADC spectrum
01830     {
01831       std::ostringstream title, path;
01832       std::string name = "adcRaw";
01833       title << "ADC values for " << detector.detName() 
01834             << " channel " << chanId.board() << "_"
01835             << chanId.connector();
01836       path << this->getPath(chanId) << "/" << name;
01837       TH1F* adcRaw = new TH1F(name.c_str(),title.str().c_str(),
01838                               4096,0,4096);
01839       adcRaw->GetXaxis()->SetTitle("ADC value");
01840       adcRaw->GetYaxis()->SetTitle("Entries");
01841       if( m_statsSvc->put(path.str(),adcRaw).isFailure() ) {
01842         error() << "prepareStats(): Could not register " << path << endreq;
01843         delete adcRaw;
01844         return StatusCode::FAILURE;
01845       }
01846     }
01847 
01848     // ADC spectrum with cuts
01849     {
01850       std::ostringstream title, path;
01851       std::string name = "adcRawWithCuts";
01852       title << "ADC values for " << detector.detName() 
01853             << " channel " << chanId.board() << "_"
01854             << chanId.connector();
01855       path << this->getPath(chanId) << "/" << name;
01856       TH1F* adcRawWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01857                               4096,0,4096);
01858       adcRawWithCuts->GetXaxis()->SetTitle("ADC value");
01859       adcRawWithCuts->GetYaxis()->SetTitle("Entries With Cuts");
01860       if( m_statsSvc->put(path.str(),adcRawWithCuts).isFailure() ) {
01861         error() << "prepareStats(): Could not register " << path << endreq;
01862         delete adcRawWithCuts;
01863         return StatusCode::FAILURE;
01864       }
01865     }
01866 
01867     
01868     // ADC spectrum, baseline subtracted
01869     {
01870       std::ostringstream title, path;
01871       std::string name = "adc";
01872       title << "Baseline-subtracted ADC values for " << detector.detName() 
01873             << " channel " << chanId.board() << "_"
01874             << chanId.connector();
01875       path << this->getPath(chanId) << "/" << name;
01876       TH1F* adc = new TH1F(name.c_str(),title.str().c_str(),
01877                            5096,-1000,4096);
01878       adc->GetXaxis()->SetTitle("ADC value");
01879       adc->GetYaxis()->SetTitle("Entries");
01880       if( m_statsSvc->put(path.str(),adc).isFailure() ) {
01881         error() << "prepareStats(): Could not register " << path << endreq;
01882         delete adc;
01883         return StatusCode::FAILURE;
01884       }
01885     }
01886 
01887     // ADC spectrum, baseline subtracted with cuts
01888     {
01889       std::ostringstream title, path;
01890       std::string name = "adcWithCuts";
01891       title << "Baseline-subtracted ADC values With Cuts for " << detector.detName() 
01892             << " channel " << chanId.board() << "_"
01893             << chanId.connector();
01894       path << this->getPath(chanId) << "/" << name;
01895       TH1F* adcWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01896                            5096,-1000,4096);
01897       adcWithCuts->GetXaxis()->SetTitle("ADC value");
01898       adcWithCuts->GetYaxis()->SetTitle("Entries");
01899       if( m_statsSvc->put(path.str(),adcWithCuts).isFailure() ) {
01900         error() << "prepareStats(): Could not register " << path << endreq;
01901         delete adcWithCuts;
01902         return StatusCode::FAILURE;
01903       }
01904     }
01905 
01906 
01907     // ADC spectrum by Clock Cycle
01908     {
01909       std::ostringstream title, path;
01910       std::string name = "adcByClock";
01911       title << "ADC values by clock cycle for " << detector.detName() 
01912             << " channel " << chanId.board() << "_"
01913             << chanId.connector();
01914       path << this->getPath(chanId) << "/" << name;
01915       TH1F* adcByClock = new TH1F(name.c_str(),title.str().c_str(),
01916                                   20,0,20);
01917       adcByClock->GetXaxis()->SetTitle("ADC Clock Cycle");
01918       adcByClock->GetYaxis()->SetTitle("ADC");
01919       if( m_statsSvc->put(path.str(),adcByClock).isFailure() ) {
01920         error() << "prepareStats(): Could not register " << path << endreq;
01921         delete adcByClock;
01922         return StatusCode::FAILURE;
01923       }
01924     }
01925     // ADC spectrum by Clock Cycle with cuts
01926     {
01927       std::ostringstream title, path;
01928       std::string name = "adcByClockWithCuts";
01929       title << "ADC values by clock cycle for " << detector.detName() 
01930             << " channel " << chanId.board() << "_"
01931             << chanId.connector();
01932       path << this->getPath(chanId) << "/" << name;
01933       TH1F* adcByClockWithCuts = new TH1F(name.c_str(),title.str().c_str(),
01934                                   20,0,20);
01935       adcByClockWithCuts->GetXaxis()->SetTitle("ADC Clock Cycle");
01936       adcByClockWithCuts->GetYaxis()->SetTitle("ADC");
01937       if( m_statsSvc->put(path.str(),adcByClockWithCuts).isFailure() ) {
01938         error() << "prepareStats(): Could not register " << path << endreq;
01939         delete adcByClockWithCuts;
01940         return StatusCode::FAILURE;
01941       }
01942     }
01943 
01944 
01945   }  
01946 
01947   m_processedDetectors.push_back(detector);
01948 
01949   return StatusCode::SUCCESS;
01950 }
01951 
01952 std::string PmtCalibLeadingEdgeWithCuts::getPath( const DayaBay::Detector& detector ){
01953   return m_filepath + "/" + detector.detName();
01954 }
01955 
01956 std::string PmtCalibLeadingEdgeWithCuts::getPath(const DayaBay::Detector& detector,
01957                                          int ring){
01958   std::ostringstream path;
01959   path << m_filepath << "/" << detector.detName() << "/ring_" << ring; 
01960   return path.str();
01961 }
01962 
01963 std::string PmtCalibLeadingEdgeWithCuts::getPath(const DayaBay::FeeChannelId& chanId){
01964   std::ostringstream path;
01965   path << m_filepath << "/" << chanId.detName() 
01966        << "/chan_" << chanId.board() << "_" << chanId.connector(); 
01967   return path.str();
01968 }
In This Package:

PmtCalibLeadingEdgeWithCuts.cc
