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
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
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
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
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
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
0146
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
0157
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
0178 }
0179
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
0206 }
0207
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
0234 }
0235
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
0251
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
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
0287
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
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
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
0370 if (flag > 0) {
0371 map<int, float> qiecalib = QieCalibMap[detid];
0372
0373
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
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
0412 char PedSampleNum[20];
0413
0414
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
0435
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
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
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
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
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
0512
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
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
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
0560
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
0598
0599 fRefPedestals = fInputPedestals;
0600 fRefPedestalWidths = fInputPedestalWidths;
0601
0602
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
0616 if (m_nevtsample < 1)
0617 SampleAnalysis();
0618 if (m_nevtsample > 0) {
0619 if (evt % m_nevtsample != 0)
0620 SampleAnalysis();
0621 }
0622
0623
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
0638
0639
0640
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
0655
0656
0657
0658 }
0659
0660
0661
0662 m_file->cd();
0663 m_file->cd("HB");
0664 hbHists.ALLPEDS->Write();
0665 hbHists.PEDRMS->Write();
0666 hbHists.PEDMEAN->Write();
0667
0668 m_file->cd();
0669 m_file->cd("HO");
0670 hoHists.ALLPEDS->Write();
0671 hoHists.PEDRMS->Write();
0672 hoHists.PEDMEAN->Write();
0673
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
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
0718
0719
0720
0721
0722
0723
0724
0725
0726 std::vector<double>::iterator sample_it;
0727
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
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
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
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
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800
0801
0802
0803
0804
0805
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
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
0853
0854
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
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
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
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
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
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
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
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
0972 }
0973
0974 if (erflag == 0)
0975 PedValLog << "HcalPedVal: all pedestals checked OK" << std::endl;
0976
0977
0978
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
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 }