Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:58:07 

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