Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 12:46:28 

0001 #include "CalibFormats/CastorObjects/interface/CastorCoderDb.h"
0002 #include "CalibFormats/CastorObjects/interface/CastorCalibrations.h"
0003 #include "CalibFormats/CastorObjects/interface/CastorDbService.h"
0004 #include "CalibFormats/CastorObjects/interface/CastorDbRecord.h"
0005 #include "CondFormats/CastorObjects/interface/CastorQIECoder.h"
0006 #include "CondFormats/CastorObjects/interface/CastorPedestals.h"
0007 #include "CondFormats/CastorObjects/interface/CastorPedestalWidths.h"
0008 
0009 #include "CalibCalorimetry/CastorCalib/interface/CastorPedestalAnalysis.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include <TFile.h>
0012 #include <cmath>
0013 
0014 CastorPedestalAnalysis::CastorPedestalAnalysis(const edm::ParameterSet& ps)
0015     : fRefPedestals(nullptr),
0016       fRefPedestalWidths(nullptr),
0017       fRawPedestals(nullptr),
0018       fRawPedestalWidths(nullptr),
0019       fValPedestals(nullptr),
0020       fValPedestalWidths(nullptr) {
0021   evt = 0;
0022   sample = 0;
0023   m_file = nullptr;
0024   m_AllPedsOK = 0;
0025   for (int i = 0; i < 4; i++)
0026     m_stat[i] = 0;
0027   for (int k = 0; k < 4; k++)
0028     state.push_back(true);
0029 
0030   // user cfg parameters
0031   m_outputFileMean = ps.getUntrackedParameter<std::string>("outputFileMeans", "");
0032   if (!m_outputFileMean.empty()) {
0033     edm::LogWarning("Castor") << "Castor pedestal means will be saved to " << m_outputFileMean.c_str();
0034   }
0035   m_outputFileWidth = ps.getUntrackedParameter<std::string>("outputFileWidths", "");
0036   if (!m_outputFileWidth.empty()) {
0037     edm::LogWarning("Castor") << "Castor pedestal widths will be saved to " << m_outputFileWidth.c_str();
0038   }
0039   m_outputFileROOT = ps.getUntrackedParameter<std::string>("outputFileHist", "");
0040   if (!m_outputFileROOT.empty()) {
0041     edm::LogWarning("Castor") << "Castor pedestal histograms will be saved to " << m_outputFileROOT.c_str();
0042   }
0043   m_nevtsample = ps.getUntrackedParameter<int>("nevtsample", 0);
0044   // for compatibility with previous versions
0045   if (m_nevtsample == 9999999)
0046     m_nevtsample = 0;
0047   m_pedsinADC = ps.getUntrackedParameter<int>("pedsinADC", 0);
0048   m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag", 0);
0049   m_pedValflag = ps.getUntrackedParameter<int>("pedValflag", 0);
0050   if (m_pedValflag < 0)
0051     m_pedValflag = 0;
0052   if (m_nevtsample > 0 && m_pedValflag > 0) {
0053     edm::LogWarning("Castor") << "WARNING - incompatible cfg options: nevtsample = " << m_nevtsample
0054                               << ", pedValflag = " << m_pedValflag;
0055     edm::LogWarning("Castor") << "Setting pedValflag = 0";
0056     m_pedValflag = 0;
0057   }
0058   if (m_pedValflag > 1)
0059     m_pedValflag = 1;
0060   m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
0061   if (m_startTS < 0)
0062     m_startTS = 0;
0063   m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
0064 
0065   //  m_logFile.open("CastorPedestalAnalysis.log");
0066 
0067   castorHists.ALLPEDS = new TH1F("Castor All Pedestals", "HF All Peds", 10, 0, 9);
0068   castorHists.PEDRMS = new TH1F("Castor All Pedestal Widths", "HF All Pedestal RMS", 100, 0, 3);
0069   castorHists.PEDMEAN = new TH1F("Castor All Pedestal Means", "HF All Pedestal Means", 100, 0, 9);
0070   castorHists.CHI2 = new TH1F("Castor Chi2/ndf for whole range Gauss fit", "HF Chi2/ndf Gauss", 200, 0., 50.);
0071 }
0072 
0073 //-----------------------------------------------------------------------------
0074 CastorPedestalAnalysis::~CastorPedestalAnalysis() {
0075   for (_meot = castorHists.PEDTRENDS.begin(); _meot != castorHists.PEDTRENDS.end(); _meot++) {
0076     for (int i = 0; i < 16; i++)
0077       _meot->second[i].first->Delete();
0078   }
0079 
0080   castorHists.ALLPEDS->Delete();
0081   castorHists.PEDRMS->Delete();
0082   castorHists.PEDMEAN->Delete();
0083   castorHists.CHI2->Delete();
0084 }
0085 
0086 //-----------------------------------------------------------------------------
0087 void CastorPedestalAnalysis::setup(const std::string& m_outputFileROOT) {
0088   // open the histogram file, create directories within
0089   m_file = new TFile(m_outputFileROOT.c_str(), "RECREATE");
0090   m_file->mkdir("Castor");
0091   m_file->cd();
0092 }
0093 
0094 //-----------------------------------------------------------------------------
0095 void CastorPedestalAnalysis::processEvent(const CastorDigiCollection& castor, const CastorDbService& cond) {
0096   evt++;
0097   sample = 1;
0098   evt_curr = evt;
0099   if (m_nevtsample > 0) {
0100     sample = (evt - 1) / m_nevtsample + 1;
0101     evt_curr = evt % m_nevtsample;
0102     if (evt_curr == 0)
0103       evt_curr = m_nevtsample;
0104   }
0105 
0106   m_shape = cond.getCastorShape();
0107   // HF
0108   try {
0109     if (castor.empty())
0110       throw(int) castor.size();
0111     for (CastorDigiCollection::const_iterator j = castor.begin(); j != castor.end(); ++j) {
0112       const CastorDataFrame digi = (const CastorDataFrame)(*j);
0113       m_coder = cond.getCastorCoder(digi.id());
0114       for (int i = m_startTS; i < digi.size() && i <= m_endTS; i++) {
0115         for (int flag = 0; flag < 4; flag++) {
0116           if (i + flag < digi.size() && i + flag <= m_endTS) {
0117             per2CapsHists(flag, 2, digi.id(), digi.sample(i), digi.sample(i + flag), castorHists.PEDTRENDS, cond);
0118           }
0119         }
0120       }
0121       if (m_startTS == 0 && m_endTS > 4) {
0122         AllChanHists(digi.id(),
0123                      digi.sample(0),
0124                      digi.sample(1),
0125                      digi.sample(2),
0126                      digi.sample(3),
0127                      digi.sample(4),
0128                      digi.sample(5),
0129                      castorHists.PEDTRENDS);
0130       }
0131     }
0132   } catch (int i) {
0133     //    m_logFile << "Event with " << i<<" Castor Digis passed." << std::endl;
0134   }
0135   // Call the function every m_nevtsample events
0136   if (m_nevtsample > 0) {
0137     if (evt % m_nevtsample == 0)
0138       SampleAnalysis();
0139   }
0140 }
0141 
0142 //-----------------------------------------------------------------------------
0143 void CastorPedestalAnalysis::per2CapsHists(int flag,
0144                                            int id,
0145                                            const HcalDetId detid,
0146                                            const HcalQIESample& qie1,
0147                                            const HcalQIESample& qie2,
0148                                            std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT,
0149                                            const CastorDbService& cond) {
0150   // this function is due to be called for every time slice, it fills either a charge
0151   // histo for a single capID (flag=0) or a product histo for two capIDs (flag>0)
0152 
0153   static const int bins = 10;
0154   static const int bins2 = 100;
0155   float lo = -0.5;
0156   float hi = 9.5;
0157   std::map<int, PEDBUNCH> _mei;
0158   static std::map<HcalDetId, std::map<int, float> > QieCalibMap;
0159   std::string type = "Castor";
0160 
0161   /*
0162   if(id==0){
0163     if(detid.ieta()<16) type = "HB";
0164     if(detid.ieta()>16) type = "HE";
0165     if(detid.ieta()==16){
0166       if(detid.depth()<3) type = "HB";
0167       if(detid.depth()==3) type = "HE";
0168     }
0169   } 
0170   else if(id==1) type = "HO";
0171   else if(id==2) type = "HF"; 
0172   */
0173 
0174   _meot = toolT.find(detid);
0175 
0176   // if histos for the current channel do not exist, first create them,
0177   if (_meot == toolT.end()) {
0178     std::map<int, PEDBUNCH> insert;
0179     std::map<int, float> qiecalib;
0180     char name[1024];
0181     for (int i = 0; i < 4; i++) {
0182       lo = -0.5;
0183       // fix from Andy: if you convert to fC and then bin in units of 1, you may 'skip' a bin while
0184       // filling, since the ADCs are quantized
0185       if (m_pedsinADC)
0186         hi = 9.5;
0187       else
0188         hi = 11.5;
0189       sprintf(
0190           name, "%s Pedestal, eta=%d phi=%d d=%d cap=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth(), i);
0191       insert[i].first = new TH1F(name, name, bins, lo, hi);
0192       sprintf(name,
0193               "%s Product, eta=%d phi=%d d=%d caps=%d*%d",
0194               type.c_str(),
0195               detid.ieta(),
0196               detid.iphi(),
0197               detid.depth(),
0198               i,
0199               (i + 1) % 4);
0200       insert[4 + i].first = new TH1F(name, name, bins2, 0., 100.);
0201       sprintf(name,
0202               "%s Product, eta=%d phi=%d d=%d caps=%d*%d",
0203               type.c_str(),
0204               detid.ieta(),
0205               detid.iphi(),
0206               detid.depth(),
0207               i,
0208               (i + 2) % 4);
0209       insert[8 + i].first = new TH1F(name, name, bins2, 0., 100.);
0210       sprintf(name,
0211               "%s Product, eta=%d phi=%d d=%d caps=%d*%d",
0212               type.c_str(),
0213               detid.ieta(),
0214               detid.iphi(),
0215               detid.depth(),
0216               i,
0217               (i + 3) % 4);
0218       insert[12 + i].first = new TH1F(name, name, bins2, 0., 100.);
0219     }
0220     sprintf(name, "%s Signal in TS 4+5, eta=%d phi=%d d=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0221     insert[16].first = new TH1F(name, name, 21, -0.5, 20.5);
0222     sprintf(
0223         name, "%s Signal in TS 4+5-2-3, eta=%d phi=%d d=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0224     insert[17].first = new TH1F(name, name, 21, -10.5, 10.5);
0225     sprintf(name,
0226             "%s Signal in TS 4+5-(0+1+2+3)/2., eta=%d phi=%d d=%d",
0227             type.c_str(),
0228             detid.ieta(),
0229             detid.iphi(),
0230             detid.depth());
0231     insert[18].first = new TH1F(name, name, 21, -10.5, 10.5);
0232     toolT[detid] = insert;
0233     _meot = toolT.find(detid);
0234     // store QIE calibrations in a map for later reuse
0235     QieCalibMap[detid] = qiecalib;
0236   }
0237 
0238   _mei = _meot->second;
0239 
0240   const CastorQIECoder* coder = cond.getCastorCoder(detid);
0241   const CastorQIEShape* shape = cond.getCastorShape();
0242   float charge1 = coder->charge(*shape, qie1.adc(), qie1.capid());
0243   float charge2 = coder->charge(*shape, qie2.adc(), qie2.capid());
0244 
0245   // fill single capID histo
0246   if (flag == 0) {
0247     if (m_nevtsample > 0) {
0248       if ((evt - 1) % m_nevtsample == 0 && state[qie1.capid()]) {
0249         state[qie1.capid()] = false;
0250         _mei[qie1.capid()].first->Reset();
0251         _mei[qie1.capid() + 4].first->Reset();
0252         _mei[qie1.capid() + 8].first->Reset();
0253         _mei[qie1.capid() + 12].first->Reset();
0254       }
0255     }
0256     if (qie1.adc() < bins) {
0257       if (m_pedsinADC)
0258         _mei[qie1.capid()].first->Fill(qie1.adc());
0259       else
0260         _mei[qie1.capid()].first->Fill(charge1);
0261     } else if (qie1.adc() >= bins) {
0262       _mei[qie1.capid()].first->AddBinContent(bins + 1, 1);
0263     }
0264   }
0265 
0266   // fill 2 capID histo
0267   if (flag > 0) {
0268     std::map<int, float> qiecalib = QieCalibMap[detid];
0269     //float charge1=(qie1.adc()-qiecalib[qie1.capid()+4])/qiecalib[qie1.capid()];
0270     //float charge2=(qie2.adc()-qiecalib[qie2.capid()+4])/qiecalib[qie2.capid()];
0271     if (charge1 * charge2 < bins2) {
0272       _mei[qie1.capid() + 4 * flag].first->Fill(charge1 * charge2);
0273     } else {
0274       _mei[qie1.capid() + 4 * flag].first->Fill(bins2);
0275     }
0276   }
0277 
0278   if (flag == 0) {
0279     //    if(id==0) hbHists.ALLPEDS->Fill(qie1.adc());
0280     //   else if(id==1) hoHists.ALLPEDS->Fill(qie1.adc());
0281     //   else if(id==2) castorHists.ALLPEDS->Fill(qie1.adc());
0282     castorHists.ALLPEDS->Fill(qie1.adc());
0283   }
0284 }
0285 
0286 //-----------------------------------------------------------------------------
0287 void CastorPedestalAnalysis::AllChanHists(const HcalDetId detid,
0288                                           const HcalQIESample& qie0,
0289                                           const HcalQIESample& qie1,
0290                                           const HcalQIESample& qie2,
0291                                           const HcalQIESample& qie3,
0292                                           const HcalQIESample& qie4,
0293                                           const HcalQIESample& qie5,
0294                                           std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT) {
0295   // this function is due to be called for every channel
0296 
0297   _meot = toolT.find(detid);
0298   std::map<int, PEDBUNCH> _mei = _meot->second;
0299   _mei[16].first->Fill(qie4.adc() + qie5.adc() - 1.);
0300   _mei[17].first->Fill(qie4.adc() + qie5.adc() - qie2.adc() - qie3.adc());
0301   _mei[18].first->Fill(qie4.adc() + qie5.adc() - (qie0.adc() + qie1.adc() + qie2.adc() + qie3.adc()) / 2.);
0302 }
0303 
0304 //-----------------------------------------------------------------------------
0305 void CastorPedestalAnalysis::SampleAnalysis() {
0306   // it is called every m_nevtsample events (a sample) and the end of run
0307   char PedSampleNum[20];
0308 
0309   // Compute pedestal constants for each HBHE, HO, HF
0310   sprintf(PedSampleNum, "Castor_Sample%d", sample);
0311   m_file->cd();
0312   m_file->mkdir(PedSampleNum);
0313   m_file->cd(PedSampleNum);
0314   GetPedConst(castorHists.PEDTRENDS, castorHists.PEDMEAN, castorHists.PEDRMS);
0315 }
0316 
0317 //-----------------------------------------------------------------------------
0318 void CastorPedestalAnalysis::GetPedConst(std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT,
0319                                          TH1F* PedMeans,
0320                                          TH1F* PedWidths) {
0321   // Completely rewritten version oct 2006
0322   // Compute pedestal constants and fill into CastorPedestals and CastorPedestalWidths objects
0323   float cap[4];
0324   float sig[4][4];
0325   float dcap[4];
0326   float dsig[4][4];
0327   float chi2[4];
0328 
0329   for (_meot = toolT.begin(); _meot != toolT.end(); _meot++) {
0330     HcalDetId detid = _meot->first;
0331 
0332     // take mean and width from a Gaussian fit or directly from the histo
0333     if (fitflag > 0) {
0334       for (int i = 0; i < 4; i++) {
0335         TF1* fit = _meot->second[i].first->GetFunction("gaus");
0336         chi2[i] = 0;
0337         if (fit->GetNDF() != 0)
0338           chi2[i] = fit->GetChisquare() / fit->GetNDF();
0339         cap[i] = fit->GetParameter(1);
0340         sig[i][i] = fit->GetParameter(2);
0341         dcap[i] = fit->GetParError(1);
0342         dsig[i][i] = fit->GetParError(2);
0343       }
0344     } else {
0345       for (int i = 0; i < 4; i++) {
0346         cap[i] = _meot->second[i].first->GetMean();
0347         sig[i][i] = _meot->second[i].first->GetRMS();
0348         m_stat[i] = 0;
0349 
0350         for (int j = m_startTS; j < m_endTS + 1; j++) {
0351           m_stat[i] += _meot->second[i].first->GetBinContent(j + 1);
0352         }
0353         dcap[i] = sig[i][i] / sqrt(m_stat[i]);
0354         //        dsig[i][i] = dcap[i]*sig[i][i]/cap[i];
0355         dsig[i][i] = sig[i][i] / sqrt(2. * m_stat[i]);
0356         chi2[i] = 0.;
0357       }
0358     }
0359 
0360     for (int i = 0; i < 4; i++) {
0361       if (m_hiSaveflag > 0) {
0362         if (m_pedsinADC)
0363           _meot->second[i].first->GetXaxis()->SetTitle("ADC");
0364         else
0365           _meot->second[i].first->GetXaxis()->SetTitle("Charge, fC");
0366         _meot->second[i].first->GetYaxis()->SetTitle("CapID samplings");
0367         _meot->second[i].first->Write();
0368       }
0369       if (m_nevtsample > 0) {
0370         _meot->second[i].second.first[0].push_back(cap[i]);
0371         _meot->second[i].second.first[1].push_back(dcap[i]);
0372         _meot->second[i].second.first[2].push_back(sig[i][i]);
0373         _meot->second[i].second.first[3].push_back(dsig[i][i]);
0374         _meot->second[i].second.first[4].push_back(chi2[i]);
0375       }
0376       PedMeans->Fill(cap[i]);
0377       PedWidths->Fill(sig[i][i]);
0378     }
0379 
0380     // special histos for Shuichi
0381     if (m_hiSaveflag == -100) {
0382       for (int i = 16; i < 19; i++) {
0383         if (m_pedsinADC)
0384           _meot->second[i].first->GetXaxis()->SetTitle("ADC");
0385         else
0386           _meot->second[i].first->GetXaxis()->SetTitle("Charge, fC");
0387         _meot->second[i].first->GetYaxis()->SetTitle("Events");
0388         _meot->second[i].first->Write();
0389       }
0390     }
0391 
0392     // diagonal sigma is width squared
0393     sig[0][0] = sig[0][0] * sig[0][0];
0394     sig[1][1] = sig[1][1] * sig[1][1];
0395     sig[2][2] = sig[2][2] * sig[2][2];
0396     sig[3][3] = sig[3][3] * sig[3][3];
0397 
0398     // off diagonal sigmas (correlations) are computed from 3 histograms
0399     // here we still have all 4*3=12 combinations
0400     sig[0][1] = _meot->second[4].first->GetMean() - cap[0] * cap[1];
0401     sig[0][2] = _meot->second[8].first->GetMean() - cap[0] * cap[2];
0402     sig[1][2] = _meot->second[5].first->GetMean() - cap[1] * cap[2];
0403     sig[1][3] = _meot->second[9].first->GetMean() - cap[1] * cap[3];
0404     sig[2][3] = _meot->second[6].first->GetMean() - cap[2] * cap[3];
0405     sig[0][3] = _meot->second[12].first->GetMean() - cap[0] * cap[3];
0406     sig[1][0] = _meot->second[13].first->GetMean() - cap[1] * cap[0];
0407     sig[2][0] = _meot->second[10].first->GetMean() - cap[2] * cap[0];
0408     sig[2][1] = _meot->second[14].first->GetMean() - cap[2] * cap[1];
0409     sig[3][1] = _meot->second[11].first->GetMean() - cap[3] * cap[1];
0410     sig[3][2] = _meot->second[15].first->GetMean() - cap[3] * cap[2];
0411     sig[3][0] = _meot->second[7].first->GetMean() - cap[3] * cap[0];
0412 
0413     // there is no proper error calculation for the correlation coefficients
0414     for (int i = 0; i < 4; i++) {
0415       if (m_nevtsample > 0) {
0416         _meot->second[i].second.first[5].push_back(sig[i][(i + 1) % 4]);
0417         _meot->second[i].second.first[6].push_back(2 * sig[i][i] * dsig[i][i]);
0418         _meot->second[i].second.first[7].push_back(sig[i][(i + 2) % 4]);
0419         _meot->second[i].second.first[8].push_back(2 * sig[i][i] * dsig[i][i]);
0420         _meot->second[i].second.first[9].push_back(sig[i][(i + 3) % 4]);
0421         _meot->second[i].second.first[10].push_back(2 * sig[i][i] * dsig[i][i]);
0422       }
0423       // save product histos if desired
0424       if (m_hiSaveflag > 10) {
0425         if (m_pedsinADC)
0426           _meot->second[i + 4].first->GetXaxis()->SetTitle("ADC^2");
0427         else
0428           _meot->second[i + 4].first->GetXaxis()->SetTitle("Charge^2, fC^2");
0429         _meot->second[i + 4].first->GetYaxis()->SetTitle("2-CapID samplings");
0430         _meot->second[i + 4].first->Write();
0431         if (m_pedsinADC)
0432           _meot->second[i + 8].first->GetXaxis()->SetTitle("ADC^2");
0433         else
0434           _meot->second[i + 8].first->GetXaxis()->SetTitle("Charge^2, fC^2");
0435         _meot->second[i + 8].first->GetYaxis()->SetTitle("2-CapID samplings");
0436         _meot->second[i + 8].first->Write();
0437         if (m_pedsinADC)
0438           _meot->second[i + 12].first->GetXaxis()->SetTitle("ADC^2");
0439         else
0440           _meot->second[i + 12].first->GetXaxis()->SetTitle("Charge^2, fC^2");
0441         _meot->second[i + 12].first->GetYaxis()->SetTitle("2-CapID samplings");
0442         _meot->second[i + 12].first->Write();
0443       }
0444     }
0445 
0446     // fill the objects - at this point only close and medium correlations are stored
0447     // and the matrix is assumed symmetric
0448     if (m_nevtsample < 1) {
0449       sig[1][0] = sig[0][1];
0450       sig[2][0] = sig[0][2];
0451       sig[2][1] = sig[1][2];
0452       sig[3][1] = sig[1][3];
0453       sig[3][2] = sig[2][3];
0454       sig[0][3] = sig[3][0];
0455       if (fRawPedestals) {
0456         CastorPedestal item(detid, cap[0], cap[1], cap[2], cap[3]);
0457         fRawPedestals->addValues(item);
0458       }
0459       if (fRawPedestalWidths) {
0460         CastorPedestalWidth widthsp(detid);
0461         widthsp.setSigma(0, 0, sig[0][0]);
0462         widthsp.setSigma(0, 1, sig[0][1]);
0463         widthsp.setSigma(0, 2, sig[0][2]);
0464         widthsp.setSigma(1, 1, sig[1][1]);
0465         widthsp.setSigma(1, 2, sig[1][2]);
0466         widthsp.setSigma(1, 3, sig[1][3]);
0467         widthsp.setSigma(2, 2, sig[2][2]);
0468         widthsp.setSigma(2, 3, sig[2][3]);
0469         widthsp.setSigma(3, 3, sig[3][3]);
0470         widthsp.setSigma(3, 0, sig[0][3]);
0471         fRawPedestalWidths->addValues(widthsp);
0472       }
0473     }
0474   }
0475 }
0476 
0477 //-----------------------------------------------------------------------------
0478 int CastorPedestalAnalysis::done(const CastorPedestals* fInputPedestals,
0479                                  const CastorPedestalWidths* fInputPedestalWidths,
0480                                  CastorPedestals* fOutputPedestals,
0481                                  CastorPedestalWidths* fOutputPedestalWidths) {
0482   int nstat[4];
0483 
0484   // Pedestal objects
0485   // inputs...
0486   fRefPedestals = fInputPedestals;
0487   fRefPedestalWidths = fInputPedestalWidths;
0488 
0489   // outputs...
0490   if (m_pedValflag > 0) {
0491     fValPedestals = fOutputPedestals;
0492     fValPedestalWidths = fOutputPedestalWidths;
0493     fRawPedestals = new CastorPedestals();
0494     fRawPedestalWidths = new CastorPedestalWidths();
0495   } else {
0496     fRawPedestals = fOutputPedestals;
0497     fRawPedestalWidths = fOutputPedestalWidths;
0498     fValPedestals = new CastorPedestals();
0499     fValPedestalWidths = new CastorPedestalWidths();
0500   }
0501 
0502   // compute pedestal constants
0503   if (m_nevtsample < 1)
0504     SampleAnalysis();
0505   if (m_nevtsample > 0) {
0506     if (evt % m_nevtsample != 0)
0507       SampleAnalysis();
0508   }
0509 
0510   // trending histos
0511   if (m_nevtsample > 0) {
0512     m_file->cd();
0513     m_file->cd("Castor");
0514     Trendings(castorHists.PEDTRENDS, castorHists.CHI2, castorHists.CAPID_AVERAGE, castorHists.CAPID_CHI2);
0515   }
0516 
0517   if (m_nevtsample < 1) {
0518     // pedestal validation: m_AllPedsOK=-1 means not validated,
0519     //                                   0 everything OK,
0520     //                                   N>0 : mod(N,100000) drifts + width changes
0521     //                                         int(N/100000) missing channels
0522     m_AllPedsOK = -1;
0523     if (m_pedValflag > 0) {
0524       for (int i = 0; i < 4; i++)
0525         nstat[i] = (int)m_stat[i];
0526       int NPedErrors = CastorPedVal(nstat,
0527                                     fRefPedestals,
0528                                     fRefPedestalWidths,
0529                                     fRawPedestals,
0530                                     fRawPedestalWidths,
0531                                     fValPedestals,
0532                                     fValPedestalWidths);
0533       m_AllPedsOK = NPedErrors;
0534     }
0535     // setting m_AllPedsOK=-2 will inhibit writing pedestals out
0536     //    if(m_pedValflag==1){
0537     //      if(evt<100)m_AllPedsOK=-2;
0538     //    }
0539   }
0540 
0541   // Write other histograms.
0542 
0543   // Castor
0544   m_file->cd();
0545   m_file->cd("Castor");
0546   castorHists.ALLPEDS->Write();
0547   castorHists.PEDRMS->Write();
0548   castorHists.PEDMEAN->Write();
0549 
0550   m_file->Close();
0551   edm::LogWarning("Castor") << "Hcal/Castor histograms written to " << m_outputFileROOT.c_str();
0552   return (int)m_AllPedsOK;
0553 }
0554 
0555 //-----------------------------------------------------------------------------
0556 void CastorPedestalAnalysis::Trendings(std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT,
0557                                        TH1F* Chi2,
0558                                        TH1F* CapidAverage,
0559                                        TH1F* CapidChi2) {
0560   // check stability of pedestal constants in a single long run
0561 
0562   std::map<int, std::vector<double> > AverageValues;
0563 
0564   for (_meot = toolT.begin(); _meot != toolT.end(); _meot++) {
0565     for (int i = 0; i < 4; i++) {
0566       char name[1024];
0567       HcalDetId detid = _meot->first;
0568       sprintf(name, "Pedestal trend, eta=%d phi=%d d=%d cap=%d", detid.ieta(), detid.iphi(), detid.depth(), i);
0569       int bins = _meot->second[i].second.first[0].size();
0570       float lo = 0.5;
0571       float hi = (float)bins + 0.5;
0572       _meot->second[i].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0573       sprintf(name, "Width trend, eta=%d phi=%d d=%d cap=%d", detid.ieta(), detid.iphi(), detid.depth(), i);
0574       bins = _meot->second[i].second.first[2].size();
0575       hi = (float)bins + 0.5;
0576       _meot->second[i].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0577       sprintf(name,
0578               "Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",
0579               detid.ieta(),
0580               detid.iphi(),
0581               detid.depth(),
0582               i,
0583               (i + 1) % 4);
0584       bins = _meot->second[i].second.first[5].size();
0585       hi = (float)bins + 0.5;
0586       _meot->second[i].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0587       /*      sprintf(name,"Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",detid.ieta(),detid.iphi(),detid.depth(),i,(i+2)%4);
0588       bins = _meot->second[i].second.first[7].size();
0589       hi = (float)bins+0.5;
0590       _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi));
0591       sprintf(name,"Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",detid.ieta(),detid.iphi(),detid.depth(),i,(i+3)%4);
0592       bins = _meot->second[i].second.first[9].size();
0593       hi = (float)bins+0.5;
0594       _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi)); */
0595 
0596       std::vector<double>::iterator sample_it;
0597       // Pedestal mean - put content and errors
0598       int j = 0;
0599       for (sample_it = _meot->second[i].second.first[0].begin(); sample_it != _meot->second[i].second.first[0].end();
0600            ++sample_it) {
0601         _meot->second[i].second.second[0]->SetBinContent(++j, *sample_it);
0602       }
0603       j = 0;
0604       for (sample_it = _meot->second[i].second.first[1].begin(); sample_it != _meot->second[i].second.first[1].end();
0605            ++sample_it) {
0606         _meot->second[i].second.second[0]->SetBinError(++j, *sample_it);
0607       }
0608       // fit with a constant - extract parameters
0609       _meot->second[i].second.second[0]->Fit("pol0", "Q");
0610       TF1* fit = _meot->second[i].second.second[0]->GetFunction("pol0");
0611       AverageValues[0].push_back(fit->GetParameter(0));
0612       AverageValues[1].push_back(fit->GetParError(0));
0613       if (sample > 1)
0614         AverageValues[2].push_back(fit->GetChisquare() / fit->GetNDF());
0615       else
0616         AverageValues[2].push_back(fit->GetChisquare());
0617       sprintf(name, "Sample (%d events)", m_nevtsample);
0618       _meot->second[i].second.second[0]->GetXaxis()->SetTitle(name);
0619       _meot->second[i].second.second[0]->GetYaxis()->SetTitle("Pedestal value");
0620       _meot->second[i].second.second[0]->Write();
0621       // Pedestal width - put content and errors
0622       j = 0;
0623       for (sample_it = _meot->second[i].second.first[2].begin(); sample_it != _meot->second[i].second.first[2].end();
0624            ++sample_it) {
0625         _meot->second[i].second.second[1]->SetBinContent(++j, *sample_it);
0626       }
0627       j = 0;
0628       for (sample_it = _meot->second[i].second.first[3].begin(); sample_it != _meot->second[i].second.first[3].end();
0629            ++sample_it) {
0630         _meot->second[i].second.second[1]->SetBinError(++j, *sample_it);
0631       }
0632       _meot->second[i].second.second[1]->GetXaxis()->SetTitle(name);
0633       _meot->second[i].second.second[1]->GetYaxis()->SetTitle("Pedestal width");
0634       _meot->second[i].second.second[1]->Write();
0635       // Correlation coeffs - put contents and errors
0636       j = 0;
0637       for (sample_it = _meot->second[i].second.first[5].begin(); sample_it != _meot->second[i].second.first[5].end();
0638            ++sample_it) {
0639         _meot->second[i].second.second[2]->SetBinContent(++j, *sample_it);
0640       }
0641       j = 0;
0642       for (sample_it = _meot->second[i].second.first[6].begin(); sample_it != _meot->second[i].second.first[6].end();
0643            ++sample_it) {
0644         _meot->second[i].second.second[2]->SetBinError(++j, *sample_it);
0645       }
0646       _meot->second[i].second.second[2]->GetXaxis()->SetTitle(name);
0647       _meot->second[i].second.second[2]->GetYaxis()->SetTitle("Close correlation");
0648       _meot->second[i].second.second[2]->Write();
0649       /*     j=0;
0650       for(sample_it=_meot->second[i].second.first[7].begin();
0651           sample_it!=_meot->second[i].second.first[7].end();sample_it++){
0652         _meot->second[i].second.second[3]->SetBinContent(++j,*sample_it);
0653       }
0654       j=0;
0655       for(sample_it=_meot->second[i].second.first[8].begin();
0656           sample_it!=_meot->second[i].second.first[8].end();sample_it++){
0657         _meot->second[i].second.second[3]->SetBinError(++j,*sample_it);
0658       }
0659       _meot->second[i].second.second[3]->GetXaxis()->SetTitle(name);
0660       _meot->second[i].second.second[3]->GetYaxis()->SetTitle("Intermediate correlation");
0661       _meot->second[i].second.second[3]->Write();
0662       j=0;
0663       for(sample_it=_meot->second[i].second.first[9].begin();
0664           sample_it!=_meot->second[i].second.first[9].end();sample_it++){
0665         _meot->second[i].second.second[4]->SetBinContent(++j,*sample_it);
0666       }
0667       j=0;
0668       for(sample_it=_meot->second[i].second.first[10].begin();
0669           sample_it!=_meot->second[i].second.first[10].end();sample_it++){
0670         _meot->second[i].second.second[4]->SetBinError(++j,*sample_it);
0671       }
0672       _meot->second[i].second.second[4]->GetXaxis()->SetTitle(name);
0673       _meot->second[i].second.second[4]->GetYaxis()->SetTitle("Distant correlation");
0674       _meot->second[i].second.second[4]->Write(); */
0675       // chi2
0676       j = 0;
0677       for (sample_it = _meot->second[i].second.first[4].begin(); sample_it != _meot->second[i].second.first[4].end();
0678            ++sample_it) {
0679         Chi2->Fill(*sample_it);
0680       }
0681     }
0682   }
0683   CapidAverage = new TH1F("Constant fit: Pedestal Values",
0684                           "Constant fit: Pedestal Values",
0685                           AverageValues[0].size(),
0686                           0.,
0687                           AverageValues[0].size());
0688   std::vector<double>::iterator sample_it;
0689   int j = 0;
0690   for (sample_it = AverageValues[0].begin(); sample_it != AverageValues[0].end(); ++sample_it) {
0691     CapidAverage->SetBinContent(++j, *sample_it);
0692   }
0693   j = 0;
0694   for (sample_it = AverageValues[1].begin(); sample_it != AverageValues[1].end(); ++sample_it) {
0695     CapidAverage->SetBinError(++j, *sample_it);
0696   }
0697   CapidChi2 = new TH1F(
0698       "Constant fit: Chi2/ndf", "Constant fit: Chi2/ndf", AverageValues[2].size(), 0., AverageValues[2].size());
0699   j = 0;
0700   for (sample_it = AverageValues[2].begin(); sample_it != AverageValues[2].end(); ++sample_it) {
0701     CapidChi2->SetBinContent(++j, *sample_it);
0702     //CapidChi2->SetBinError(++j,0);
0703   }
0704   Chi2->GetXaxis()->SetTitle("Chi2/ndf");
0705   Chi2->GetYaxis()->SetTitle("50 x [(16+2) x 4 x 4] `events`");
0706   Chi2->Write();
0707   CapidAverage->GetYaxis()->SetTitle("Pedestal value");
0708   CapidAverage->GetXaxis()->SetTitle("(16+2) x 4 x 4 `events`");
0709   CapidAverage->Write();
0710   CapidChi2->GetYaxis()->SetTitle("Chi2/ndf");
0711   CapidChi2->GetXaxis()->SetTitle("(16+2) x 4 x 4 `events`");
0712   CapidChi2->Write();
0713 }
0714 
0715 //-----------------------------------------------------------------------------
0716 int CastorPedestalAnalysis::CastorPedVal(int nstat[4],
0717                                          const CastorPedestals* fRefPedestals,
0718                                          const CastorPedestalWidths* fRefPedestalWidths,
0719                                          CastorPedestals* fRawPedestals,
0720                                          CastorPedestalWidths* fRawPedestalWidths,
0721                                          CastorPedestals* fValPedestals,
0722                                          CastorPedestalWidths* fValPedestalWidths) {
0723   // new version of pedestal validation - it is designed to be as independent of
0724   // all the rest as possible - you only need to provide valid pedestal objects
0725   // and a vector of statistics per capID to use this as standalone code
0726   HcalDetId detid;
0727   float RefPedVals[4];
0728   float RefPedSigs[4][4];
0729   float RawPedVals[4];
0730   float RawPedSigs[4][4];
0731   std::map<HcalDetId, bool> isinRaw;
0732   std::map<HcalDetId, bool> isinRef;
0733   std::vector<DetId> RefChanns = fRefPedestals->getAllChannels();
0734   std::vector<DetId> RawChanns = fRawPedestals->getAllChannels();
0735   std::ofstream PedValLog;
0736   PedValLog.open("CastorPedVal.log");
0737 
0738   if (nstat[0] + nstat[1] + nstat[2] + nstat[3] < 2500)
0739     PedValLog << "CastorPedVal: warning - low statistics" << std::endl;
0740   // find complete list of channels in current data and reference
0741   for (int i = 0; i < (int)RawChanns.size(); i++) {
0742     isinRef[HcalDetId(RawChanns[i])] = false;
0743   }
0744   for (int i = 0; i < (int)RefChanns.size(); i++) {
0745     detid = HcalDetId(RefChanns[i]);
0746     isinRaw[detid] = false;
0747     isinRef[detid] = true;
0748   }
0749   for (int i = 0; i < (int)RawChanns.size(); i++) {
0750     detid = HcalDetId(RawChanns[i]);
0751     isinRaw[detid] = true;
0752     if (isinRef[detid] == false) {
0753       PedValLog << "CastorPedVal: channel " << detid << " not found in reference set" << std::endl;
0754       std::cerr << "CastorPedVal: channel " << detid << " not found in reference set" << std::endl;
0755     }
0756   }
0757 
0758   // main loop over channels
0759   int erflag = 0;
0760   for (int i = 0; i < (int)RefChanns.size(); i++) {
0761     detid = HcalDetId(RefChanns[i]);
0762     for (int icap = 0; icap < 4; icap++) {
0763       RefPedVals[icap] = fRefPedestals->getValues(detid)->getValue(icap);
0764       for (int icap2 = icap; icap2 < 4; icap2++) {
0765         RefPedSigs[icap][icap2] = fRefPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0766         if (icap2 != icap)
0767           RefPedSigs[icap2][icap] = RefPedSigs[icap][icap2];
0768       }
0769     }
0770 
0771     // read new raw values
0772     if (isinRaw[detid]) {
0773       for (int icap = 0; icap < 4; icap++) {
0774         RawPedVals[icap] = fRawPedestals->getValues(detid)->getValue(icap);
0775         for (int icap2 = icap; icap2 < 4; icap2++) {
0776           RawPedSigs[icap][icap2] = fRawPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0777           if (icap2 != icap)
0778             RawPedSigs[icap2][icap] = RawPedSigs[icap][icap2];
0779         }
0780       }
0781 
0782       // first quick check if raw values make sense: if not, the channel is treated like absent
0783       for (int icap = 0; icap < 4; icap++) {
0784         if (RawPedVals[icap] < 1. || RawPedSigs[icap][icap] < 0.01)
0785           isinRaw[detid] = false;
0786         for (int icap2 = icap; icap2 < 4; icap2++) {
0787           if (fabs(RawPedSigs[icap][icap2] / sqrt(RawPedSigs[icap][icap] * RawPedSigs[icap2][icap2])) > 1.)
0788             isinRaw[detid] = false;
0789         }
0790       }
0791     }
0792 
0793     // check raw values against reference
0794     if (isinRaw[detid]) {
0795       for (int icap = 0; icap < 4; icap++) {
0796         int icap2 = (icap + 1) % 4;
0797         float width = sqrt(RawPedSigs[icap][icap]);
0798         float erof1 = width / sqrt((float)nstat[icap]);
0799         float erof2 = sqrt(erof1 * erof1 + RawPedSigs[icap][icap] / (float)nstat[icap]);
0800         float erofwidth = width / sqrt(2. * nstat[icap]);
0801         float diffof1 = RawPedVals[icap] - RefPedVals[icap];
0802         float diffof2 = RawPedVals[icap] + RawPedVals[icap2] - RefPedVals[icap] - RefPedVals[icap2];
0803         float diffofw = width - sqrt(RefPedSigs[icap][icap]);
0804 
0805         // validation in 2 TS for HB, HE, HO, in 1 TS for HF
0806         int nTS = 2;
0807         if (detid.subdet() == HcalForward)
0808           nTS = 1;
0809         if (nTS == 1 && fabs(diffof1) > 0.5 + erof1) {
0810           erflag += 1;
0811           PedValLog << "HcalPedVal: drift in channel " << detid << " cap " << icap << ": " << RawPedVals[icap] << " - "
0812                     << RefPedVals[icap] << " = " << diffof1 << std::endl;
0813         }
0814         if (nTS == 2 && fabs(diffof2) > 0.5 + erof2) {
0815           erflag += 1;
0816           PedValLog << "HcalPedVal: drift in channel " << detid << " caps " << icap << "+" << icap2 << ": "
0817                     << RawPedVals[icap] << "+" << RawPedVals[icap2] << " - " << RefPedVals[icap] << "+"
0818                     << RefPedVals[icap2] << " = " << diffof2 << std::endl;
0819         }
0820         if (fabs(diffofw) > 0.15 * width + erofwidth) {
0821           erflag += 1;
0822           PedValLog << "HcalPedVal: width changed in channel " << detid << " cap " << icap << ": " << width << " - "
0823                     << sqrt(RefPedSigs[icap][icap]) << " = " << diffofw << std::endl;
0824         }
0825       }
0826     }
0827 
0828     // for disconnected/bad channels restore reference values
0829     else {
0830       PedValLog << "HcalPedVal: no valid data from channel " << detid << std::endl;
0831       erflag += 100000;
0832       CastorPedestal item(detid, RefPedVals[0], RefPedVals[1], RefPedVals[2], RefPedVals[3]);
0833       fValPedestals->addValues(item);
0834       CastorPedestalWidth widthsp(detid);
0835       for (int icap = 0; icap < 4; icap++) {
0836         for (int icap2 = icap; icap2 < 4; icap2++)
0837           widthsp.setSigma(icap2, icap, RefPedSigs[icap2][icap]);
0838       }
0839       fValPedestalWidths->addValues(widthsp);
0840     }
0841 
0842     // end of channel loop
0843   }
0844 
0845   if (erflag == 0)
0846     PedValLog << "HcalPedVal: all pedestals checked OK" << std::endl;
0847 
0848   // now construct the remaining part of the validated objects
0849   // if nothing changed outside tolerance, validated set = reference set
0850   if (erflag % 100000 == 0) {
0851     for (int i = 0; i < (int)RefChanns.size(); i++) {
0852       detid = HcalDetId(RefChanns[i]);
0853       if (isinRaw[detid]) {
0854         CastorPedestalWidth widthsp(detid);
0855         for (int icap = 0; icap < 4; icap++) {
0856           RefPedVals[icap] = fRefPedestals->getValues(detid)->getValue(icap);
0857           for (int icap2 = icap; icap2 < 4; icap2++) {
0858             RefPedSigs[icap][icap2] = fRefPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0859             if (icap2 != icap)
0860               RefPedSigs[icap2][icap] = RefPedSigs[icap][icap2];
0861             widthsp.setSigma(icap2, icap, RefPedSigs[icap2][icap]);
0862           }
0863         }
0864         fValPedestalWidths->addValues(widthsp);
0865         CastorPedestal item(detid, RefPedVals[0], RefPedVals[1], RefPedVals[2], RefPedVals[3]);
0866         fValPedestals->addValues(item);
0867       }
0868     }
0869   }
0870 
0871   // if anything changed, validated set = raw set + reference for missing/bad channels
0872   else {
0873     for (int i = 0; i < (int)RawChanns.size(); i++) {
0874       detid = HcalDetId(RawChanns[i]);
0875       if (isinRaw[detid]) {
0876         CastorPedestalWidth widthsp(detid);
0877         for (int icap = 0; icap < 4; icap++) {
0878           RawPedVals[icap] = fRawPedestals->getValues(detid)->getValue(icap);
0879           for (int icap2 = icap; icap2 < 4; icap2++) {
0880             RawPedSigs[icap][icap2] = fRawPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0881             if (icap2 != icap)
0882               RawPedSigs[icap2][icap] = RawPedSigs[icap][icap2];
0883             widthsp.setSigma(icap2, icap, RawPedSigs[icap2][icap]);
0884           }
0885         }
0886         fValPedestalWidths->addValues(widthsp);
0887         CastorPedestal item(detid, RawPedVals[0], RawPedVals[1], RawPedVals[2], RawPedVals[3]);
0888         fValPedestals->addValues(item);
0889       }
0890     }
0891   }
0892   return erflag;
0893 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!