File indexing completed on 2021-02-14 12:46:28
0001 #include "CalibFormats/CastorObjects/interface/CastorCoderDb.h"
0002 #include "CalibFormats/CastorObjects/interface/CastorCalibrations.h"
0003 #include "CalibFormats/CastorObjects/interface/CastorDbService.h"
0004 #include "CalibFormats/CastorObjects/interface/CastorDbRecord.h"
0005 #include "CondFormats/CastorObjects/interface/CastorQIECoder.h"
0006 #include "CondFormats/CastorObjects/interface/CastorPedestals.h"
0007 #include "CondFormats/CastorObjects/interface/CastorPedestalWidths.h"
0008
0009 #include "CalibCalorimetry/CastorCalib/interface/CastorPedestalAnalysis.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include <TFile.h>
0012 #include <cmath>
0013
0014 CastorPedestalAnalysis::CastorPedestalAnalysis(const edm::ParameterSet& ps)
0015 : fRefPedestals(nullptr),
0016 fRefPedestalWidths(nullptr),
0017 fRawPedestals(nullptr),
0018 fRawPedestalWidths(nullptr),
0019 fValPedestals(nullptr),
0020 fValPedestalWidths(nullptr) {
0021 evt = 0;
0022 sample = 0;
0023 m_file = nullptr;
0024 m_AllPedsOK = 0;
0025 for (int i = 0; i < 4; i++)
0026 m_stat[i] = 0;
0027 for (int k = 0; k < 4; k++)
0028 state.push_back(true);
0029
0030
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
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
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
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
0108 try {
0109 if (castor.empty())
0110 throw(int) castor.size();
0111 for (CastorDigiCollection::const_iterator j = castor.begin(); j != castor.end(); ++j) {
0112 const CastorDataFrame digi = (const CastorDataFrame)(*j);
0113 m_coder = cond.getCastorCoder(digi.id());
0114 for (int i = m_startTS; i < digi.size() && i <= m_endTS; i++) {
0115 for (int flag = 0; flag < 4; flag++) {
0116 if (i + flag < digi.size() && i + flag <= m_endTS) {
0117 per2CapsHists(flag, 2, digi.id(), digi.sample(i), digi.sample(i + flag), castorHists.PEDTRENDS, cond);
0118 }
0119 }
0120 }
0121 if (m_startTS == 0 && m_endTS > 4) {
0122 AllChanHists(digi.id(),
0123 digi.sample(0),
0124 digi.sample(1),
0125 digi.sample(2),
0126 digi.sample(3),
0127 digi.sample(4),
0128 digi.sample(5),
0129 castorHists.PEDTRENDS);
0130 }
0131 }
0132 } catch (int i) {
0133
0134 }
0135
0136 if (m_nevtsample > 0) {
0137 if (evt % m_nevtsample == 0)
0138 SampleAnalysis();
0139 }
0140 }
0141
0142
0143 void CastorPedestalAnalysis::per2CapsHists(int flag,
0144 int id,
0145 const HcalDetId detid,
0146 const HcalQIESample& qie1,
0147 const HcalQIESample& qie2,
0148 std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT,
0149 const CastorDbService& cond) {
0150
0151
0152
0153 static const int bins = 10;
0154 static const int bins2 = 100;
0155 float lo = -0.5;
0156 float hi = 9.5;
0157 std::map<int, PEDBUNCH> _mei;
0158 static std::map<HcalDetId, std::map<int, float> > QieCalibMap;
0159 std::string type = "Castor";
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170
0171
0172
0173
0174 _meot = toolT.find(detid);
0175
0176
0177 if (_meot == toolT.end()) {
0178 std::map<int, PEDBUNCH> insert;
0179 std::map<int, float> qiecalib;
0180 char name[1024];
0181 for (int i = 0; i < 4; i++) {
0182 lo = -0.5;
0183
0184
0185 if (m_pedsinADC)
0186 hi = 9.5;
0187 else
0188 hi = 11.5;
0189 sprintf(
0190 name, "%s Pedestal, eta=%d phi=%d d=%d cap=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth(), i);
0191 insert[i].first = new TH1F(name, name, bins, lo, hi);
0192 sprintf(name,
0193 "%s Product, eta=%d phi=%d d=%d caps=%d*%d",
0194 type.c_str(),
0195 detid.ieta(),
0196 detid.iphi(),
0197 detid.depth(),
0198 i,
0199 (i + 1) % 4);
0200 insert[4 + i].first = new TH1F(name, name, bins2, 0., 100.);
0201 sprintf(name,
0202 "%s Product, eta=%d phi=%d d=%d caps=%d*%d",
0203 type.c_str(),
0204 detid.ieta(),
0205 detid.iphi(),
0206 detid.depth(),
0207 i,
0208 (i + 2) % 4);
0209 insert[8 + i].first = new TH1F(name, name, bins2, 0., 100.);
0210 sprintf(name,
0211 "%s Product, eta=%d phi=%d d=%d caps=%d*%d",
0212 type.c_str(),
0213 detid.ieta(),
0214 detid.iphi(),
0215 detid.depth(),
0216 i,
0217 (i + 3) % 4);
0218 insert[12 + i].first = new TH1F(name, name, bins2, 0., 100.);
0219 }
0220 sprintf(name, "%s Signal in TS 4+5, eta=%d phi=%d d=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0221 insert[16].first = new TH1F(name, name, 21, -0.5, 20.5);
0222 sprintf(
0223 name, "%s Signal in TS 4+5-2-3, eta=%d phi=%d d=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
0224 insert[17].first = new TH1F(name, name, 21, -10.5, 10.5);
0225 sprintf(name,
0226 "%s Signal in TS 4+5-(0+1+2+3)/2., eta=%d phi=%d d=%d",
0227 type.c_str(),
0228 detid.ieta(),
0229 detid.iphi(),
0230 detid.depth());
0231 insert[18].first = new TH1F(name, name, 21, -10.5, 10.5);
0232 toolT[detid] = insert;
0233 _meot = toolT.find(detid);
0234
0235 QieCalibMap[detid] = qiecalib;
0236 }
0237
0238 _mei = _meot->second;
0239
0240 const CastorQIECoder* coder = cond.getCastorCoder(detid);
0241 const CastorQIEShape* shape = cond.getCastorShape();
0242 float charge1 = coder->charge(*shape, qie1.adc(), qie1.capid());
0243 float charge2 = coder->charge(*shape, qie2.adc(), qie2.capid());
0244
0245
0246 if (flag == 0) {
0247 if (m_nevtsample > 0) {
0248 if ((evt - 1) % m_nevtsample == 0 && state[qie1.capid()]) {
0249 state[qie1.capid()] = false;
0250 _mei[qie1.capid()].first->Reset();
0251 _mei[qie1.capid() + 4].first->Reset();
0252 _mei[qie1.capid() + 8].first->Reset();
0253 _mei[qie1.capid() + 12].first->Reset();
0254 }
0255 }
0256 if (qie1.adc() < bins) {
0257 if (m_pedsinADC)
0258 _mei[qie1.capid()].first->Fill(qie1.adc());
0259 else
0260 _mei[qie1.capid()].first->Fill(charge1);
0261 } else if (qie1.adc() >= bins) {
0262 _mei[qie1.capid()].first->AddBinContent(bins + 1, 1);
0263 }
0264 }
0265
0266
0267 if (flag > 0) {
0268 std::map<int, float> qiecalib = QieCalibMap[detid];
0269
0270
0271 if (charge1 * charge2 < bins2) {
0272 _mei[qie1.capid() + 4 * flag].first->Fill(charge1 * charge2);
0273 } else {
0274 _mei[qie1.capid() + 4 * flag].first->Fill(bins2);
0275 }
0276 }
0277
0278 if (flag == 0) {
0279
0280
0281
0282 castorHists.ALLPEDS->Fill(qie1.adc());
0283 }
0284 }
0285
0286
0287 void CastorPedestalAnalysis::AllChanHists(const HcalDetId detid,
0288 const HcalQIESample& qie0,
0289 const HcalQIESample& qie1,
0290 const HcalQIESample& qie2,
0291 const HcalQIESample& qie3,
0292 const HcalQIESample& qie4,
0293 const HcalQIESample& qie5,
0294 std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT) {
0295
0296
0297 _meot = toolT.find(detid);
0298 std::map<int, PEDBUNCH> _mei = _meot->second;
0299 _mei[16].first->Fill(qie4.adc() + qie5.adc() - 1.);
0300 _mei[17].first->Fill(qie4.adc() + qie5.adc() - qie2.adc() - qie3.adc());
0301 _mei[18].first->Fill(qie4.adc() + qie5.adc() - (qie0.adc() + qie1.adc() + qie2.adc() + qie3.adc()) / 2.);
0302 }
0303
0304
0305 void CastorPedestalAnalysis::SampleAnalysis() {
0306
0307 char PedSampleNum[20];
0308
0309
0310 sprintf(PedSampleNum, "Castor_Sample%d", sample);
0311 m_file->cd();
0312 m_file->mkdir(PedSampleNum);
0313 m_file->cd(PedSampleNum);
0314 GetPedConst(castorHists.PEDTRENDS, castorHists.PEDMEAN, castorHists.PEDRMS);
0315 }
0316
0317
0318 void CastorPedestalAnalysis::GetPedConst(std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT,
0319 TH1F* PedMeans,
0320 TH1F* PedWidths) {
0321
0322
0323 float cap[4];
0324 float sig[4][4];
0325 float dcap[4];
0326 float dsig[4][4];
0327 float chi2[4];
0328
0329 for (_meot = toolT.begin(); _meot != toolT.end(); _meot++) {
0330 HcalDetId detid = _meot->first;
0331
0332
0333 if (fitflag > 0) {
0334 for (int i = 0; i < 4; i++) {
0335 TF1* fit = _meot->second[i].first->GetFunction("gaus");
0336 chi2[i] = 0;
0337 if (fit->GetNDF() != 0)
0338 chi2[i] = fit->GetChisquare() / fit->GetNDF();
0339 cap[i] = fit->GetParameter(1);
0340 sig[i][i] = fit->GetParameter(2);
0341 dcap[i] = fit->GetParError(1);
0342 dsig[i][i] = fit->GetParError(2);
0343 }
0344 } else {
0345 for (int i = 0; i < 4; i++) {
0346 cap[i] = _meot->second[i].first->GetMean();
0347 sig[i][i] = _meot->second[i].first->GetRMS();
0348 m_stat[i] = 0;
0349
0350 for (int j = m_startTS; j < m_endTS + 1; j++) {
0351 m_stat[i] += _meot->second[i].first->GetBinContent(j + 1);
0352 }
0353 dcap[i] = sig[i][i] / sqrt(m_stat[i]);
0354
0355 dsig[i][i] = sig[i][i] / sqrt(2. * m_stat[i]);
0356 chi2[i] = 0.;
0357 }
0358 }
0359
0360 for (int i = 0; i < 4; i++) {
0361 if (m_hiSaveflag > 0) {
0362 if (m_pedsinADC)
0363 _meot->second[i].first->GetXaxis()->SetTitle("ADC");
0364 else
0365 _meot->second[i].first->GetXaxis()->SetTitle("Charge, fC");
0366 _meot->second[i].first->GetYaxis()->SetTitle("CapID samplings");
0367 _meot->second[i].first->Write();
0368 }
0369 if (m_nevtsample > 0) {
0370 _meot->second[i].second.first[0].push_back(cap[i]);
0371 _meot->second[i].second.first[1].push_back(dcap[i]);
0372 _meot->second[i].second.first[2].push_back(sig[i][i]);
0373 _meot->second[i].second.first[3].push_back(dsig[i][i]);
0374 _meot->second[i].second.first[4].push_back(chi2[i]);
0375 }
0376 PedMeans->Fill(cap[i]);
0377 PedWidths->Fill(sig[i][i]);
0378 }
0379
0380
0381 if (m_hiSaveflag == -100) {
0382 for (int i = 16; i < 19; i++) {
0383 if (m_pedsinADC)
0384 _meot->second[i].first->GetXaxis()->SetTitle("ADC");
0385 else
0386 _meot->second[i].first->GetXaxis()->SetTitle("Charge, fC");
0387 _meot->second[i].first->GetYaxis()->SetTitle("Events");
0388 _meot->second[i].first->Write();
0389 }
0390 }
0391
0392
0393 sig[0][0] = sig[0][0] * sig[0][0];
0394 sig[1][1] = sig[1][1] * sig[1][1];
0395 sig[2][2] = sig[2][2] * sig[2][2];
0396 sig[3][3] = sig[3][3] * sig[3][3];
0397
0398
0399
0400 sig[0][1] = _meot->second[4].first->GetMean() - cap[0] * cap[1];
0401 sig[0][2] = _meot->second[8].first->GetMean() - cap[0] * cap[2];
0402 sig[1][2] = _meot->second[5].first->GetMean() - cap[1] * cap[2];
0403 sig[1][3] = _meot->second[9].first->GetMean() - cap[1] * cap[3];
0404 sig[2][3] = _meot->second[6].first->GetMean() - cap[2] * cap[3];
0405 sig[0][3] = _meot->second[12].first->GetMean() - cap[0] * cap[3];
0406 sig[1][0] = _meot->second[13].first->GetMean() - cap[1] * cap[0];
0407 sig[2][0] = _meot->second[10].first->GetMean() - cap[2] * cap[0];
0408 sig[2][1] = _meot->second[14].first->GetMean() - cap[2] * cap[1];
0409 sig[3][1] = _meot->second[11].first->GetMean() - cap[3] * cap[1];
0410 sig[3][2] = _meot->second[15].first->GetMean() - cap[3] * cap[2];
0411 sig[3][0] = _meot->second[7].first->GetMean() - cap[3] * cap[0];
0412
0413
0414 for (int i = 0; i < 4; i++) {
0415 if (m_nevtsample > 0) {
0416 _meot->second[i].second.first[5].push_back(sig[i][(i + 1) % 4]);
0417 _meot->second[i].second.first[6].push_back(2 * sig[i][i] * dsig[i][i]);
0418 _meot->second[i].second.first[7].push_back(sig[i][(i + 2) % 4]);
0419 _meot->second[i].second.first[8].push_back(2 * sig[i][i] * dsig[i][i]);
0420 _meot->second[i].second.first[9].push_back(sig[i][(i + 3) % 4]);
0421 _meot->second[i].second.first[10].push_back(2 * sig[i][i] * dsig[i][i]);
0422 }
0423
0424 if (m_hiSaveflag > 10) {
0425 if (m_pedsinADC)
0426 _meot->second[i + 4].first->GetXaxis()->SetTitle("ADC^2");
0427 else
0428 _meot->second[i + 4].first->GetXaxis()->SetTitle("Charge^2, fC^2");
0429 _meot->second[i + 4].first->GetYaxis()->SetTitle("2-CapID samplings");
0430 _meot->second[i + 4].first->Write();
0431 if (m_pedsinADC)
0432 _meot->second[i + 8].first->GetXaxis()->SetTitle("ADC^2");
0433 else
0434 _meot->second[i + 8].first->GetXaxis()->SetTitle("Charge^2, fC^2");
0435 _meot->second[i + 8].first->GetYaxis()->SetTitle("2-CapID samplings");
0436 _meot->second[i + 8].first->Write();
0437 if (m_pedsinADC)
0438 _meot->second[i + 12].first->GetXaxis()->SetTitle("ADC^2");
0439 else
0440 _meot->second[i + 12].first->GetXaxis()->SetTitle("Charge^2, fC^2");
0441 _meot->second[i + 12].first->GetYaxis()->SetTitle("2-CapID samplings");
0442 _meot->second[i + 12].first->Write();
0443 }
0444 }
0445
0446
0447
0448 if (m_nevtsample < 1) {
0449 sig[1][0] = sig[0][1];
0450 sig[2][0] = sig[0][2];
0451 sig[2][1] = sig[1][2];
0452 sig[3][1] = sig[1][3];
0453 sig[3][2] = sig[2][3];
0454 sig[0][3] = sig[3][0];
0455 if (fRawPedestals) {
0456 CastorPedestal item(detid, cap[0], cap[1], cap[2], cap[3]);
0457 fRawPedestals->addValues(item);
0458 }
0459 if (fRawPedestalWidths) {
0460 CastorPedestalWidth widthsp(detid);
0461 widthsp.setSigma(0, 0, sig[0][0]);
0462 widthsp.setSigma(0, 1, sig[0][1]);
0463 widthsp.setSigma(0, 2, sig[0][2]);
0464 widthsp.setSigma(1, 1, sig[1][1]);
0465 widthsp.setSigma(1, 2, sig[1][2]);
0466 widthsp.setSigma(1, 3, sig[1][3]);
0467 widthsp.setSigma(2, 2, sig[2][2]);
0468 widthsp.setSigma(2, 3, sig[2][3]);
0469 widthsp.setSigma(3, 3, sig[3][3]);
0470 widthsp.setSigma(3, 0, sig[0][3]);
0471 fRawPedestalWidths->addValues(widthsp);
0472 }
0473 }
0474 }
0475 }
0476
0477
0478 int CastorPedestalAnalysis::done(const CastorPedestals* fInputPedestals,
0479 const CastorPedestalWidths* fInputPedestalWidths,
0480 CastorPedestals* fOutputPedestals,
0481 CastorPedestalWidths* fOutputPedestalWidths) {
0482 int nstat[4];
0483
0484
0485
0486 fRefPedestals = fInputPedestals;
0487 fRefPedestalWidths = fInputPedestalWidths;
0488
0489
0490 if (m_pedValflag > 0) {
0491 fValPedestals = fOutputPedestals;
0492 fValPedestalWidths = fOutputPedestalWidths;
0493 fRawPedestals = new CastorPedestals();
0494 fRawPedestalWidths = new CastorPedestalWidths();
0495 } else {
0496 fRawPedestals = fOutputPedestals;
0497 fRawPedestalWidths = fOutputPedestalWidths;
0498 fValPedestals = new CastorPedestals();
0499 fValPedestalWidths = new CastorPedestalWidths();
0500 }
0501
0502
0503 if (m_nevtsample < 1)
0504 SampleAnalysis();
0505 if (m_nevtsample > 0) {
0506 if (evt % m_nevtsample != 0)
0507 SampleAnalysis();
0508 }
0509
0510
0511 if (m_nevtsample > 0) {
0512 m_file->cd();
0513 m_file->cd("Castor");
0514 Trendings(castorHists.PEDTRENDS, castorHists.CHI2, castorHists.CAPID_AVERAGE, castorHists.CAPID_CHI2);
0515 }
0516
0517 if (m_nevtsample < 1) {
0518
0519
0520
0521
0522 m_AllPedsOK = -1;
0523 if (m_pedValflag > 0) {
0524 for (int i = 0; i < 4; i++)
0525 nstat[i] = (int)m_stat[i];
0526 int NPedErrors = CastorPedVal(nstat,
0527 fRefPedestals,
0528 fRefPedestalWidths,
0529 fRawPedestals,
0530 fRawPedestalWidths,
0531 fValPedestals,
0532 fValPedestalWidths);
0533 m_AllPedsOK = NPedErrors;
0534 }
0535
0536
0537
0538
0539 }
0540
0541
0542
0543
0544 m_file->cd();
0545 m_file->cd("Castor");
0546 castorHists.ALLPEDS->Write();
0547 castorHists.PEDRMS->Write();
0548 castorHists.PEDMEAN->Write();
0549
0550 m_file->Close();
0551 edm::LogWarning("Castor") << "Hcal/Castor histograms written to " << m_outputFileROOT.c_str();
0552 return (int)m_AllPedsOK;
0553 }
0554
0555
0556 void CastorPedestalAnalysis::Trendings(std::map<HcalDetId, std::map<int, PEDBUNCH> >& toolT,
0557 TH1F* Chi2,
0558 TH1F* CapidAverage,
0559 TH1F* CapidChi2) {
0560
0561
0562 std::map<int, std::vector<double> > AverageValues;
0563
0564 for (_meot = toolT.begin(); _meot != toolT.end(); _meot++) {
0565 for (int i = 0; i < 4; i++) {
0566 char name[1024];
0567 HcalDetId detid = _meot->first;
0568 sprintf(name, "Pedestal trend, eta=%d phi=%d d=%d cap=%d", detid.ieta(), detid.iphi(), detid.depth(), i);
0569 int bins = _meot->second[i].second.first[0].size();
0570 float lo = 0.5;
0571 float hi = (float)bins + 0.5;
0572 _meot->second[i].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0573 sprintf(name, "Width trend, eta=%d phi=%d d=%d cap=%d", detid.ieta(), detid.iphi(), detid.depth(), i);
0574 bins = _meot->second[i].second.first[2].size();
0575 hi = (float)bins + 0.5;
0576 _meot->second[i].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0577 sprintf(name,
0578 "Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",
0579 detid.ieta(),
0580 detid.iphi(),
0581 detid.depth(),
0582 i,
0583 (i + 1) % 4);
0584 bins = _meot->second[i].second.first[5].size();
0585 hi = (float)bins + 0.5;
0586 _meot->second[i].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0587
0588
0589
0590
0591
0592
0593
0594
0595
0596 std::vector<double>::iterator sample_it;
0597
0598 int j = 0;
0599 for (sample_it = _meot->second[i].second.first[0].begin(); sample_it != _meot->second[i].second.first[0].end();
0600 ++sample_it) {
0601 _meot->second[i].second.second[0]->SetBinContent(++j, *sample_it);
0602 }
0603 j = 0;
0604 for (sample_it = _meot->second[i].second.first[1].begin(); sample_it != _meot->second[i].second.first[1].end();
0605 ++sample_it) {
0606 _meot->second[i].second.second[0]->SetBinError(++j, *sample_it);
0607 }
0608
0609 _meot->second[i].second.second[0]->Fit("pol0", "Q");
0610 TF1* fit = _meot->second[i].second.second[0]->GetFunction("pol0");
0611 AverageValues[0].push_back(fit->GetParameter(0));
0612 AverageValues[1].push_back(fit->GetParError(0));
0613 if (sample > 1)
0614 AverageValues[2].push_back(fit->GetChisquare() / fit->GetNDF());
0615 else
0616 AverageValues[2].push_back(fit->GetChisquare());
0617 sprintf(name, "Sample (%d events)", m_nevtsample);
0618 _meot->second[i].second.second[0]->GetXaxis()->SetTitle(name);
0619 _meot->second[i].second.second[0]->GetYaxis()->SetTitle("Pedestal value");
0620 _meot->second[i].second.second[0]->Write();
0621
0622 j = 0;
0623 for (sample_it = _meot->second[i].second.first[2].begin(); sample_it != _meot->second[i].second.first[2].end();
0624 ++sample_it) {
0625 _meot->second[i].second.second[1]->SetBinContent(++j, *sample_it);
0626 }
0627 j = 0;
0628 for (sample_it = _meot->second[i].second.first[3].begin(); sample_it != _meot->second[i].second.first[3].end();
0629 ++sample_it) {
0630 _meot->second[i].second.second[1]->SetBinError(++j, *sample_it);
0631 }
0632 _meot->second[i].second.second[1]->GetXaxis()->SetTitle(name);
0633 _meot->second[i].second.second[1]->GetYaxis()->SetTitle("Pedestal width");
0634 _meot->second[i].second.second[1]->Write();
0635
0636 j = 0;
0637 for (sample_it = _meot->second[i].second.first[5].begin(); sample_it != _meot->second[i].second.first[5].end();
0638 ++sample_it) {
0639 _meot->second[i].second.second[2]->SetBinContent(++j, *sample_it);
0640 }
0641 j = 0;
0642 for (sample_it = _meot->second[i].second.first[6].begin(); sample_it != _meot->second[i].second.first[6].end();
0643 ++sample_it) {
0644 _meot->second[i].second.second[2]->SetBinError(++j, *sample_it);
0645 }
0646 _meot->second[i].second.second[2]->GetXaxis()->SetTitle(name);
0647 _meot->second[i].second.second[2]->GetYaxis()->SetTitle("Close correlation");
0648 _meot->second[i].second.second[2]->Write();
0649
0650
0651
0652
0653
0654
0655
0656
0657
0658
0659
0660
0661
0662
0663
0664
0665
0666
0667
0668
0669
0670
0671
0672
0673
0674
0675
0676 j = 0;
0677 for (sample_it = _meot->second[i].second.first[4].begin(); sample_it != _meot->second[i].second.first[4].end();
0678 ++sample_it) {
0679 Chi2->Fill(*sample_it);
0680 }
0681 }
0682 }
0683 CapidAverage = new TH1F("Constant fit: Pedestal Values",
0684 "Constant fit: Pedestal Values",
0685 AverageValues[0].size(),
0686 0.,
0687 AverageValues[0].size());
0688 std::vector<double>::iterator sample_it;
0689 int j = 0;
0690 for (sample_it = AverageValues[0].begin(); sample_it != AverageValues[0].end(); ++sample_it) {
0691 CapidAverage->SetBinContent(++j, *sample_it);
0692 }
0693 j = 0;
0694 for (sample_it = AverageValues[1].begin(); sample_it != AverageValues[1].end(); ++sample_it) {
0695 CapidAverage->SetBinError(++j, *sample_it);
0696 }
0697 CapidChi2 = new TH1F(
0698 "Constant fit: Chi2/ndf", "Constant fit: Chi2/ndf", AverageValues[2].size(), 0., AverageValues[2].size());
0699 j = 0;
0700 for (sample_it = AverageValues[2].begin(); sample_it != AverageValues[2].end(); ++sample_it) {
0701 CapidChi2->SetBinContent(++j, *sample_it);
0702
0703 }
0704 Chi2->GetXaxis()->SetTitle("Chi2/ndf");
0705 Chi2->GetYaxis()->SetTitle("50 x [(16+2) x 4 x 4] `events`");
0706 Chi2->Write();
0707 CapidAverage->GetYaxis()->SetTitle("Pedestal value");
0708 CapidAverage->GetXaxis()->SetTitle("(16+2) x 4 x 4 `events`");
0709 CapidAverage->Write();
0710 CapidChi2->GetYaxis()->SetTitle("Chi2/ndf");
0711 CapidChi2->GetXaxis()->SetTitle("(16+2) x 4 x 4 `events`");
0712 CapidChi2->Write();
0713 }
0714
0715
0716 int CastorPedestalAnalysis::CastorPedVal(int nstat[4],
0717 const CastorPedestals* fRefPedestals,
0718 const CastorPedestalWidths* fRefPedestalWidths,
0719 CastorPedestals* fRawPedestals,
0720 CastorPedestalWidths* fRawPedestalWidths,
0721 CastorPedestals* fValPedestals,
0722 CastorPedestalWidths* fValPedestalWidths) {
0723
0724
0725
0726 HcalDetId detid;
0727 float RefPedVals[4];
0728 float RefPedSigs[4][4];
0729 float RawPedVals[4];
0730 float RawPedSigs[4][4];
0731 std::map<HcalDetId, bool> isinRaw;
0732 std::map<HcalDetId, bool> isinRef;
0733 std::vector<DetId> RefChanns = fRefPedestals->getAllChannels();
0734 std::vector<DetId> RawChanns = fRawPedestals->getAllChannels();
0735 std::ofstream PedValLog;
0736 PedValLog.open("CastorPedVal.log");
0737
0738 if (nstat[0] + nstat[1] + nstat[2] + nstat[3] < 2500)
0739 PedValLog << "CastorPedVal: warning - low statistics" << std::endl;
0740
0741 for (int i = 0; i < (int)RawChanns.size(); i++) {
0742 isinRef[HcalDetId(RawChanns[i])] = false;
0743 }
0744 for (int i = 0; i < (int)RefChanns.size(); i++) {
0745 detid = HcalDetId(RefChanns[i]);
0746 isinRaw[detid] = false;
0747 isinRef[detid] = true;
0748 }
0749 for (int i = 0; i < (int)RawChanns.size(); i++) {
0750 detid = HcalDetId(RawChanns[i]);
0751 isinRaw[detid] = true;
0752 if (isinRef[detid] == false) {
0753 PedValLog << "CastorPedVal: channel " << detid << " not found in reference set" << std::endl;
0754 std::cerr << "CastorPedVal: channel " << detid << " not found in reference set" << std::endl;
0755 }
0756 }
0757
0758
0759 int erflag = 0;
0760 for (int i = 0; i < (int)RefChanns.size(); i++) {
0761 detid = HcalDetId(RefChanns[i]);
0762 for (int icap = 0; icap < 4; icap++) {
0763 RefPedVals[icap] = fRefPedestals->getValues(detid)->getValue(icap);
0764 for (int icap2 = icap; icap2 < 4; icap2++) {
0765 RefPedSigs[icap][icap2] = fRefPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0766 if (icap2 != icap)
0767 RefPedSigs[icap2][icap] = RefPedSigs[icap][icap2];
0768 }
0769 }
0770
0771
0772 if (isinRaw[detid]) {
0773 for (int icap = 0; icap < 4; icap++) {
0774 RawPedVals[icap] = fRawPedestals->getValues(detid)->getValue(icap);
0775 for (int icap2 = icap; icap2 < 4; icap2++) {
0776 RawPedSigs[icap][icap2] = fRawPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0777 if (icap2 != icap)
0778 RawPedSigs[icap2][icap] = RawPedSigs[icap][icap2];
0779 }
0780 }
0781
0782
0783 for (int icap = 0; icap < 4; icap++) {
0784 if (RawPedVals[icap] < 1. || RawPedSigs[icap][icap] < 0.01)
0785 isinRaw[detid] = false;
0786 for (int icap2 = icap; icap2 < 4; icap2++) {
0787 if (fabs(RawPedSigs[icap][icap2] / sqrt(RawPedSigs[icap][icap] * RawPedSigs[icap2][icap2])) > 1.)
0788 isinRaw[detid] = false;
0789 }
0790 }
0791 }
0792
0793
0794 if (isinRaw[detid]) {
0795 for (int icap = 0; icap < 4; icap++) {
0796 int icap2 = (icap + 1) % 4;
0797 float width = sqrt(RawPedSigs[icap][icap]);
0798 float erof1 = width / sqrt((float)nstat[icap]);
0799 float erof2 = sqrt(erof1 * erof1 + RawPedSigs[icap][icap] / (float)nstat[icap]);
0800 float erofwidth = width / sqrt(2. * nstat[icap]);
0801 float diffof1 = RawPedVals[icap] - RefPedVals[icap];
0802 float diffof2 = RawPedVals[icap] + RawPedVals[icap2] - RefPedVals[icap] - RefPedVals[icap2];
0803 float diffofw = width - sqrt(RefPedSigs[icap][icap]);
0804
0805
0806 int nTS = 2;
0807 if (detid.subdet() == HcalForward)
0808 nTS = 1;
0809 if (nTS == 1 && fabs(diffof1) > 0.5 + erof1) {
0810 erflag += 1;
0811 PedValLog << "HcalPedVal: drift in channel " << detid << " cap " << icap << ": " << RawPedVals[icap] << " - "
0812 << RefPedVals[icap] << " = " << diffof1 << std::endl;
0813 }
0814 if (nTS == 2 && fabs(diffof2) > 0.5 + erof2) {
0815 erflag += 1;
0816 PedValLog << "HcalPedVal: drift in channel " << detid << " caps " << icap << "+" << icap2 << ": "
0817 << RawPedVals[icap] << "+" << RawPedVals[icap2] << " - " << RefPedVals[icap] << "+"
0818 << RefPedVals[icap2] << " = " << diffof2 << std::endl;
0819 }
0820 if (fabs(diffofw) > 0.15 * width + erofwidth) {
0821 erflag += 1;
0822 PedValLog << "HcalPedVal: width changed in channel " << detid << " cap " << icap << ": " << width << " - "
0823 << sqrt(RefPedSigs[icap][icap]) << " = " << diffofw << std::endl;
0824 }
0825 }
0826 }
0827
0828
0829 else {
0830 PedValLog << "HcalPedVal: no valid data from channel " << detid << std::endl;
0831 erflag += 100000;
0832 CastorPedestal item(detid, RefPedVals[0], RefPedVals[1], RefPedVals[2], RefPedVals[3]);
0833 fValPedestals->addValues(item);
0834 CastorPedestalWidth widthsp(detid);
0835 for (int icap = 0; icap < 4; icap++) {
0836 for (int icap2 = icap; icap2 < 4; icap2++)
0837 widthsp.setSigma(icap2, icap, RefPedSigs[icap2][icap]);
0838 }
0839 fValPedestalWidths->addValues(widthsp);
0840 }
0841
0842
0843 }
0844
0845 if (erflag == 0)
0846 PedValLog << "HcalPedVal: all pedestals checked OK" << std::endl;
0847
0848
0849
0850 if (erflag % 100000 == 0) {
0851 for (int i = 0; i < (int)RefChanns.size(); i++) {
0852 detid = HcalDetId(RefChanns[i]);
0853 if (isinRaw[detid]) {
0854 CastorPedestalWidth widthsp(detid);
0855 for (int icap = 0; icap < 4; icap++) {
0856 RefPedVals[icap] = fRefPedestals->getValues(detid)->getValue(icap);
0857 for (int icap2 = icap; icap2 < 4; icap2++) {
0858 RefPedSigs[icap][icap2] = fRefPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0859 if (icap2 != icap)
0860 RefPedSigs[icap2][icap] = RefPedSigs[icap][icap2];
0861 widthsp.setSigma(icap2, icap, RefPedSigs[icap2][icap]);
0862 }
0863 }
0864 fValPedestalWidths->addValues(widthsp);
0865 CastorPedestal item(detid, RefPedVals[0], RefPedVals[1], RefPedVals[2], RefPedVals[3]);
0866 fValPedestals->addValues(item);
0867 }
0868 }
0869 }
0870
0871
0872 else {
0873 for (int i = 0; i < (int)RawChanns.size(); i++) {
0874 detid = HcalDetId(RawChanns[i]);
0875 if (isinRaw[detid]) {
0876 CastorPedestalWidth widthsp(detid);
0877 for (int icap = 0; icap < 4; icap++) {
0878 RawPedVals[icap] = fRawPedestals->getValues(detid)->getValue(icap);
0879 for (int icap2 = icap; icap2 < 4; icap2++) {
0880 RawPedSigs[icap][icap2] = fRawPedestalWidths->getValues(detid)->getSigma(icap, icap2);
0881 if (icap2 != icap)
0882 RawPedSigs[icap2][icap] = RawPedSigs[icap][icap2];
0883 widthsp.setSigma(icap2, icap, RawPedSigs[icap2][icap]);
0884 }
0885 }
0886 fValPedestalWidths->addValues(widthsp);
0887 CastorPedestal item(detid, RawPedVals[0], RawPedVals[1], RawPedVals[2], RawPedVals[3]);
0888 fValPedestals->addValues(item);
0889 }
0890 }
0891 }
0892 return erflag;
0893 }