Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-02-14 12:47:15 

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