Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-03-23 15:57:14 

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