File indexing completed on 2024-09-07 04:34:53
0001 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0002 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
0003 #include "CondFormats/HcalObjects/interface/HcalPedestals.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "CalibCalorimetry/HcalAlgos/interface/HcalLedAnalysis.h"
0006 #include "TFile.h"
0007 #include <cmath>
0008 using namespace std;
0009
0010 HcalLedAnalysis::HcalLedAnalysis(const edm::ParameterSet& ps) {
0011
0012
0013 m_coder = nullptr;
0014 m_ped = nullptr;
0015 m_shape = nullptr;
0016 evt = 0;
0017 sample = 0;
0018 m_file = nullptr;
0019
0020 for (int k = 0; k < 4; k++)
0021 state.push_back(true);
0022 m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
0023 m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML", "");
0024 if (!m_outputFileText.empty()) {
0025 edm::LogInfo("HcalLedAnalysis") << "Hcal LED results will be saved to " << m_outputFileText.c_str() << endl;
0026 m_outFile.open(m_outputFileText.c_str());
0027 }
0028 m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
0029 if (!m_outputFileROOT.empty()) {
0030 edm::LogInfo("HcalLedAnalysis") << "Hcal LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
0031 }
0032
0033 m_nevtsample = ps.getUntrackedParameter<int>("nevtsample", 9999999);
0034 if (m_nevtsample < 1)
0035 m_nevtsample = 9999999;
0036 m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag", 0);
0037 if (m_hiSaveflag < 0)
0038 m_hiSaveflag = 0;
0039 if (m_hiSaveflag > 0)
0040 m_hiSaveflag = 1;
0041 m_fitflag = ps.getUntrackedParameter<int>("analysisflag", 2);
0042 if (m_fitflag < 0)
0043 m_fitflag = 0;
0044 if (m_fitflag > 4)
0045 m_fitflag = 4;
0046 m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
0047 if (m_startTS < 0)
0048 m_startTS = 0;
0049 m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
0050 m_usecalib = ps.getUntrackedParameter<bool>("usecalib", false);
0051
0052 int runNum = ps.getUntrackedParameter<int>("runNumber", 999999);
0053
0054
0055 hbHists.ALLLEDS = new TH1F("HBHE All LEDs", "HB/HE All Leds", 10, 0, 9);
0056 hbHists.LEDRMS = new TH1F("HBHE All LED RMS", "HB/HE All LED RMS", 100, 0, 3);
0057 hbHists.LEDMEAN = new TH1F("HBHE All LED Means", "HB/HE All LED Means", 100, 0, 9);
0058 hbHists.CHI2 = new TH1F("HBHE Chi2 by ndf for Landau fit", "HB/HE Chi2/ndf Landau", 200, 0., 50.);
0059
0060 hoHists.ALLLEDS = new TH1F("HO All LEDs", "HO All Leds", 10, 0, 9);
0061 hoHists.LEDRMS = new TH1F("HO All LED RMS", "HO All LED RMS", 100, 0, 3);
0062 hoHists.LEDMEAN = new TH1F("HO All LED Means", "HO All LED Means", 100, 0, 9);
0063 hoHists.CHI2 = new TH1F("HO Chi2 by ndf for Landau fit", "HO Chi2/ndf Landau", 200, 0., 50.);
0064
0065 hfHists.ALLLEDS = new TH1F("HF All LEDs", "HF All Leds", 10, 0, 9);
0066 hfHists.LEDRMS = new TH1F("HF All LED RMS", "HF All LED RMS", 100, 0, 3);
0067 hfHists.LEDMEAN = new TH1F("HF All LED Means", "HF All LED Means", 100, 0, 9);
0068 hfHists.CHI2 = new TH1F("HF Chi2 by ndf for Landau fit", "HF Chi2/ndf Landau", 200, 0., 50.);
0069
0070 char output[100]{0};
0071
0072
0073 m_outputFileXML.open(m_outputFileX.c_str());
0074
0075 m_outputFileXML << "<?xml version='1.0' encoding='UTF-8'?>" << endl;
0076
0077 m_outputFileXML << "<ROOT>" << endl;
0078
0079 m_outputFileXML << " <HEADER>" << endl;
0080
0081 m_outputFileXML << " <TYPE>" << endl;
0082
0083 m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>" << endl;
0084
0085 m_outputFileXML << " <NAME>HCAL LED Timing</NAME>" << endl;
0086
0087 m_outputFileXML << " </TYPE>" << endl;
0088
0089 m_outputFileXML << " <RUN>" << endl;
0090
0091 m_outputFileXML << " <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>" << endl;
0092
0093 snprintf(output, sizeof output, " <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
0094 m_outputFileXML << output << endl;
0095
0096 m_outputFileXML << " <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>" << endl;
0097
0098 m_outputFileXML << " <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
0099
0100 m_outputFileXML << " </RUN>" << endl;
0101
0102 m_outputFileXML << " </HEADER>" << endl;
0103
0104 m_outputFileXML << "<!-- Tags secton -->" << endl;
0105
0106 m_outputFileXML << " <ELEMENTS>" << endl;
0107
0108 m_outputFileXML << " <DATA_SET id='-1'/>" << endl;
0109
0110 m_outputFileXML << " <IOV id='1'>" << endl;
0111
0112 m_outputFileXML << " <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>" << endl;
0113
0114 m_outputFileXML << " <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>" << endl;
0115
0116 m_outputFileXML << " </IOV>" << endl;
0117
0118 m_outputFileXML << " <TAG id='2' mode='auto'>" << endl;
0119
0120 snprintf(output, sizeof output, " <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
0121 m_outputFileXML << output << endl;
0122
0123 m_outputFileXML << " <DETECTOR_NAME>HCAL</DETECTOR_NAME>" << endl;
0124
0125 m_outputFileXML << " <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
0126
0127 m_outputFileXML << " </TAG>" << endl;
0128
0129 m_outputFileXML << " </ELEMENTS>" << endl;
0130
0131 m_outputFileXML << " <MAPS>" << endl;
0132
0133 m_outputFileXML << " <TAG idref ='2'>" << endl;
0134
0135 m_outputFileXML << " <IOV idref='1'>" << endl;
0136
0137 m_outputFileXML << " <DATA_SET idref='-1' />" << endl;
0138
0139 m_outputFileXML << " </IOV>" << endl;
0140
0141 m_outputFileXML << " </TAG>" << endl;
0142
0143 m_outputFileXML << " </MAPS>" << endl;
0144 }
0145
0146
0147 HcalLedAnalysis::~HcalLedAnalysis() {
0148
0149 for (_meol = hbHists.LEDTRENDS.begin(); _meol != hbHists.LEDTRENDS.end(); _meol++) {
0150 for (int i = 0; i < 15; i++)
0151 _meol->second[i].first->Delete();
0152 }
0153 for (_meol = hoHists.LEDTRENDS.begin(); _meol != hoHists.LEDTRENDS.end(); _meol++) {
0154 for (int i = 0; i < 15; i++)
0155 _meol->second[i].first->Delete();
0156 }
0157 for (_meol = hfHists.LEDTRENDS.begin(); _meol != hfHists.LEDTRENDS.end(); _meol++) {
0158 for (int i = 0; i < 15; i++)
0159 _meol->second[i].first->Delete();
0160 }
0161 hbHists.ALLLEDS->Delete();
0162 hbHists.LEDRMS->Delete();
0163 hbHists.LEDMEAN->Delete();
0164 hbHists.CHI2->Delete();
0165
0166 hoHists.ALLLEDS->Delete();
0167 hoHists.LEDRMS->Delete();
0168 hoHists.LEDMEAN->Delete();
0169 hoHists.CHI2->Delete();
0170
0171 hfHists.ALLLEDS->Delete();
0172 hfHists.LEDRMS->Delete();
0173 hfHists.LEDMEAN->Delete();
0174 hfHists.CHI2->Delete();
0175 }
0176
0177
0178 void HcalLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
0179
0180 m_file = new TFile(m_outputFileROOT.c_str(), "RECREATE");
0181 m_file->mkdir("HBHE");
0182 m_file->cd();
0183 m_file->mkdir("HO");
0184 m_file->cd();
0185 m_file->mkdir("HF");
0186 m_file->cd();
0187 m_file->mkdir("Calib");
0188 m_file->cd();
0189 }
0190
0191
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201
0202
0203
0204
0205
0206
0207
0208
0209
0210
0211
0212
0213 void HcalLedAnalysis::GetLedConst(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0214 double time2 = 0;
0215 double time1 = 0;
0216 double time3 = 0;
0217 double time4 = 0;
0218 double dtime2 = 0;
0219 double dtime1 = 0;
0220 double dtime3 = 0;
0221 double dtime4 = 0;
0222 char output[256]{0};
0223
0224 if (!m_outputFileText.empty()) {
0225 if (m_fitflag == 0 || m_fitflag == 2)
0226 m_outFile << "Det Eta,Phi,D Mean Error" << std::endl;
0227 else if (m_fitflag == 1 || m_fitflag == 3)
0228 m_outFile << "Det Eta,Phi,D Peak Error" << std::endl;
0229 else if (m_fitflag == 4)
0230 m_outFile << "Det Eta,Phi,D Mean Error Peak Error MeanEv Error PeakEv Error"
0231 << std::endl;
0232 }
0233 for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
0234
0235 _meol->second[10].first->Scale(1. / evt_curr);
0236 if (m_fitflag == 0 || m_fitflag == 4) {
0237 time1 = _meol->second[10].first->GetMean();
0238 dtime1 = _meol->second[10].first->GetRMS() / sqrt((float)evt_curr * (m_endTS - m_startTS + 1));
0239 }
0240 if (m_fitflag == 1 || m_fitflag == 4) {
0241
0242 for (int j = 0; j < 10; j++)
0243 _meol->second[10].first->SetBinError(j + 1, _meol->second[j].first->GetRMS() / sqrt((float)evt_curr));
0244 }
0245 if (m_fitflag == 1 || m_fitflag == 3 || m_fitflag == 4) {
0246 _meol->second[10].first->Fit("landau", "Q");
0247
0248 TF1* fit = _meol->second[10].first->GetFunction("landau");
0249
0250 time2 = fit->GetParameter(1);
0251 dtime2 = fit->GetParError(1);
0252 }
0253 if (m_fitflag == 2 || m_fitflag == 4) {
0254 time3 = _meol->second[12].first->GetMean();
0255 dtime3 = _meol->second[12].first->GetRMS() / sqrt((float)_meol->second[12].first->GetEntries());
0256 }
0257 if (m_fitflag == 3 || m_fitflag == 4) {
0258 time4 = _meol->second[13].first->GetMean();
0259 dtime4 = _meol->second[13].first->GetRMS() / sqrt((float)_meol->second[13].first->GetEntries());
0260 }
0261 for (int i = 0; i < 10; i++) {
0262 _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
0263 _meol->second[i].first->GetYaxis()->SetTitle("Counts");
0264
0265 }
0266 _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
0267 _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
0268 if (m_hiSaveflag > 0)
0269 _meol->second[10].first->Write();
0270 _meol->second[10].second.first[0].push_back(time1);
0271 _meol->second[10].second.first[1].push_back(dtime1);
0272 _meol->second[11].second.first[0].push_back(time2);
0273 _meol->second[11].second.first[1].push_back(dtime2);
0274 _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
0275 _meol->second[12].first->GetYaxis()->SetTitle("Counts");
0276 if (m_fitflag == 2 && m_hiSaveflag > 0)
0277 _meol->second[12].first->Write();
0278 _meol->second[12].second.first[0].push_back(time3);
0279 _meol->second[12].second.first[1].push_back(dtime3);
0280 _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
0281 _meol->second[13].first->GetYaxis()->SetTitle("Counts");
0282 if (m_fitflag > 2 && m_hiSaveflag > 0)
0283 _meol->second[13].first->Write();
0284 _meol->second[13].second.first[0].push_back(time4);
0285 _meol->second[13].second.first[1].push_back(dtime4);
0286 _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
0287 _meol->second[14].first->GetYaxis()->SetTitle("Counts");
0288 if (m_fitflag > 2 && m_hiSaveflag > 0)
0289 _meol->second[14].first->Write();
0290 _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
0291 _meol->second[15].first->GetYaxis()->SetTitle("Counts");
0292 if (m_fitflag > 2 && m_hiSaveflag > 0)
0293 _meol->second[15].first->Write();
0294 _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
0295 _meol->second[16].first->Write();
0296
0297
0298 HcalDetId detid = _meol->first;
0299
0300 if (!m_outputFileText.empty()) {
0301 if (m_fitflag == 0) {
0302 m_outFile << detid << " " << time1 << " " << dtime1 << std::endl;
0303 m_outputFileXML << " <DATA_SET>" << endl;
0304 m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
0305 m_outputFileXML << " <CHANNEL>" << endl;
0306 m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0307 snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
0308 m_outputFileXML << output << endl;
0309 snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
0310 m_outputFileXML << output << endl;
0311 snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
0312 m_outputFileXML << output << endl;
0313 snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
0314 m_outputFileXML << output << endl;
0315
0316 if (detid.subdet() == 1)
0317 m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0318 if (detid.subdet() == 2)
0319 m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0320 if (detid.subdet() == 3)
0321 m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0322 if (detid.subdet() == 4)
0323 m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0324 snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0325 m_outputFileXML << output << endl;
0326 m_outputFileXML << " </CHANNEL>" << endl;
0327 m_outputFileXML << " <DATA>" << endl;
0328 snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time1);
0329 m_outputFileXML << output << endl;
0330 m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0331 snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
0332 m_outputFileXML << output << endl;
0333 snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0334 m_outputFileXML << output << endl;
0335 m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
0336 m_outputFileXML << " </DATA>" << endl;
0337 m_outputFileXML << " </DATA_SET>" << endl;
0338
0339 } else if (m_fitflag == 1) {
0340 m_outFile << detid << " " << time2 << " " << dtime2 << std::endl;
0341 m_outputFileXML << " <DATA_SET>" << endl;
0342 m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
0343 m_outputFileXML << " <CHANNEL>" << endl;
0344 m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0345 snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
0346 m_outputFileXML << output << endl;
0347 snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
0348 m_outputFileXML << output << endl;
0349 snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
0350 m_outputFileXML << output << endl;
0351 snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
0352 m_outputFileXML << output << endl;
0353 if (detid.subdet() == 1)
0354 m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0355 if (detid.subdet() == 2)
0356 m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0357 if (detid.subdet() == 3)
0358 m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0359 if (detid.subdet() == 4)
0360 m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0361 snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0362 m_outputFileXML << output << endl;
0363 m_outputFileXML << " </CHANNEL>" << endl;
0364 m_outputFileXML << " <DATA>" << endl;
0365 snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time2);
0366 m_outputFileXML << output << endl;
0367 m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0368 snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime2);
0369 m_outputFileXML << output << endl;
0370 snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0371 m_outputFileXML << output << endl;
0372 m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
0373 m_outputFileXML << " </DATA>" << endl;
0374 m_outputFileXML << " </DATA_SET>" << endl;
0375 }
0376
0377 else if (m_fitflag == 2) {
0378 m_outFile << detid << " " << time3 << " " << dtime3 << std::endl;
0379 m_outputFileXML << " <DATA_SET>" << endl;
0380 m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
0381 m_outputFileXML << " <CHANNEL>" << endl;
0382 m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0383 snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
0384 m_outputFileXML << output << endl;
0385 snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
0386 m_outputFileXML << output << endl;
0387 snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
0388 m_outputFileXML << output << endl;
0389 snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
0390 m_outputFileXML << output << endl;
0391 if (detid.subdet() == 1)
0392 m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0393 if (detid.subdet() == 2)
0394 m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0395 if (detid.subdet() == 3)
0396 m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0397 if (detid.subdet() == 4)
0398 m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0399 snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0400 m_outputFileXML << output << endl;
0401 m_outputFileXML << " </CHANNEL>" << endl;
0402 m_outputFileXML << " <DATA>" << endl;
0403 snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time3);
0404 m_outputFileXML << output << endl;
0405 m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0406 snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime3);
0407 m_outputFileXML << output << endl;
0408 snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0409 m_outputFileXML << output << endl;
0410 m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
0411 m_outputFileXML << " </DATA>" << endl;
0412 m_outputFileXML << " </DATA_SET>" << endl;
0413 } else if (m_fitflag == 3) {
0414 m_outFile << detid << " " << time4 << " " << dtime4 << std::endl;
0415 m_outputFileXML << " <DATA_SET>" << endl;
0416 m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
0417 m_outputFileXML << " <CHANNEL>" << endl;
0418 m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
0419 snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
0420 m_outputFileXML << output << endl;
0421 snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
0422 m_outputFileXML << output << endl;
0423 snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
0424 m_outputFileXML << output << endl;
0425 sprintf(output, " <Z>%2i</Z>", detid.zside());
0426 m_outputFileXML << output << endl;
0427 if (detid.subdet() == 1)
0428 m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
0429 if (detid.subdet() == 2)
0430 m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
0431 if (detid.subdet() == 3)
0432 m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
0433 if (detid.subdet() == 4)
0434 m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
0435 snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
0436 m_outputFileXML << output << endl;
0437 m_outputFileXML << " </CHANNEL>" << endl;
0438 m_outputFileXML << " <DATA>" << endl;
0439 snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time4);
0440 m_outputFileXML << output << endl;
0441 m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
0442 snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
0443 m_outputFileXML << output << endl;
0444 snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
0445 m_outputFileXML << output << endl;
0446 m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
0447 m_outputFileXML << " </DATA>" << endl;
0448 m_outputFileXML << " </DATA_SET>" << endl;
0449 }
0450
0451 else if (m_fitflag == 4) {
0452 m_outFile << detid << " " << time1 << " " << dtime1 << " " << time2 << " " << dtime2 << " " << time3
0453 << " " << dtime3 << " " << time4 << " " << dtime4 << std::endl;
0454 }
0455 }
0456 }
0457 }
0458
0459
0460 void HcalLedAnalysis::LedSampleAnalysis() {
0461
0462 char LedSampleNum[20];
0463
0464 snprintf(LedSampleNum, sizeof LedSampleNum, "LedSample_%d", sample);
0465 m_file->cd();
0466 m_file->mkdir(LedSampleNum);
0467 m_file->cd(LedSampleNum);
0468
0469
0470 GetLedConst(hbHists.LEDTRENDS);
0471 GetLedConst(hoHists.LEDTRENDS);
0472 GetLedConst(hfHists.LEDTRENDS);
0473 }
0474
0475
0476 void HcalLedAnalysis::LedTrendings(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0477 for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
0478 char name[1024];
0479 HcalDetId detid = _meol->first;
0480 snprintf(name, sizeof name, "LED timing trend, eta=%d phi=%d depth=%d", detid.ieta(), detid.iphi(), detid.depth());
0481 int bins = _meol->second[10 + m_fitflag].second.first[0].size();
0482 float lo = 0.5;
0483 float hi = (float)bins + 0.5;
0484 _meol->second[10 + m_fitflag].second.second.push_back(new TH1F(name, name, bins, lo, hi));
0485
0486 std::vector<double>::iterator sample_it;
0487
0488 int j = 0;
0489 for (sample_it = _meol->second[10 + m_fitflag].second.first[0].begin();
0490 sample_it != _meol->second[10 + m_fitflag].second.first[0].end();
0491 ++sample_it) {
0492 _meol->second[10 + m_fitflag].second.second[0]->SetBinContent(++j, *sample_it);
0493 }
0494 j = 0;
0495 for (sample_it = _meol->second[10 + m_fitflag].second.first[1].begin();
0496 sample_it != _meol->second[10 + m_fitflag].second.first[1].end();
0497 ++sample_it) {
0498 _meol->second[10 + m_fitflag].second.second[0]->SetBinError(++j, *sample_it);
0499 }
0500 snprintf(name, sizeof name, "Sample (%d events)", m_nevtsample);
0501 _meol->second[10 + m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
0502 _meol->second[10 + m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
0503 _meol->second[10 + m_fitflag].second.second[0]->Write();
0504 }
0505 }
0506
0507
0508 void HcalLedAnalysis::LedDone() {
0509
0510 if (evt % m_nevtsample != 0)
0511 LedSampleAnalysis();
0512
0513
0514 if (sample > 1 && m_fitflag != 4) {
0515 m_file->cd();
0516 m_file->cd("HBHE");
0517 LedTrendings(hbHists.LEDTRENDS);
0518 m_file->cd();
0519 m_file->cd("HO");
0520 LedTrendings(hoHists.LEDTRENDS);
0521 m_file->cd();
0522 m_file->cd("HF");
0523 LedTrendings(hfHists.LEDTRENDS);
0524 }
0525
0526
0527
0528 m_file->cd();
0529 m_file->cd("HBHE");
0530 hbHists.ALLLEDS->Write();
0531 hbHists.LEDRMS->Write();
0532 hbHists.LEDMEAN->Write();
0533
0534 m_file->cd();
0535 m_file->cd("HO");
0536 hoHists.ALLLEDS->Write();
0537 hoHists.LEDRMS->Write();
0538 hoHists.LEDMEAN->Write();
0539
0540 m_file->cd();
0541 m_file->cd("HF");
0542 hfHists.ALLLEDS->Write();
0543 hfHists.LEDRMS->Write();
0544 hfHists.LEDMEAN->Write();
0545
0546 m_file->cd();
0547 m_file->cd("Calib");
0548 for (_meca = calibHists.begin(); _meca != calibHists.end(); _meca++) {
0549 _meca->second.avePulse->Write();
0550 _meca->second.integPulse->Write();
0551 }
0552
0553
0554
0555 m_file->Close();
0556 edm::LogInfo("HcalLedAnalysis") << "Hcal histograms written to " << m_outputFileROOT.c_str() << endl;
0557 }
0558
0559
0560 void HcalLedAnalysis::processLedEvent(const HBHEDigiCollection& hbhe,
0561 const HODigiCollection& ho,
0562 const HFDigiCollection& hf,
0563 const HcalCalibDigiCollection& calib,
0564 const HcalDbService& cond) {
0565 evt++;
0566 sample = (evt - 1) / m_nevtsample + 1;
0567 evt_curr = evt % m_nevtsample;
0568 if (evt_curr == 0)
0569 evt_curr = m_nevtsample;
0570
0571
0572
0573 if (m_usecalib) {
0574 if (calib.empty()) {
0575 edm::LogError("HcalLedAnalysis") << "Event with " << (int)calib.size() << " Calib Digis passed.";
0576 return;
0577 }
0578
0579 for (HcalCalibDigiCollection::const_iterator j = calib.begin(); j != calib.end(); ++j) {
0580 const HcalCalibDataFrame digi = (const HcalCalibDataFrame)(*j);
0581 HcalElectronicsId elecId = digi.elecId();
0582 HcalCalibDetId calibId = digi.id();
0583 ProcessCalibEvent(elecId.fiberChanId(),
0584 calibId,
0585 digi);
0586 }
0587 }
0588
0589
0590 if (hbhe.empty()) {
0591 edm::LogError("HcalLedAnalysis") << "Event with " << (int)hbhe.size() << " HBHE Digis passed.";
0592 return;
0593 }
0594
0595 for (HBHEDigiCollection::const_iterator j = hbhe.begin(); j != hbhe.end(); ++j) {
0596 const HBHEDataFrame digi = (const HBHEDataFrame)(*j);
0597 for (int k = 0; k < (int)state.size(); k++)
0598 state[k] = true;
0599
0600 _meol = hbHists.LEDTRENDS.find(digi.id());
0601 if (_meol == hbHists.LEDTRENDS.end()) {
0602 SetupLEDHists(0, digi.id(), hbHists.LEDTRENDS);
0603 }
0604 LedHBHEHists(digi.id(), digi, hbHists.LEDTRENDS, cond);
0605 }
0606
0607
0608 if (ho.empty()) {
0609 edm::LogError("HcalLedAnalysis") << "Event with " << (int)ho.size() << " HO Digis passed.";
0610 return;
0611 }
0612 for (HODigiCollection::const_iterator j = ho.begin(); j != ho.end(); ++j) {
0613 const HODataFrame digi = (const HODataFrame)(*j);
0614 _meol = hoHists.LEDTRENDS.find(digi.id());
0615 if (_meol == hoHists.LEDTRENDS.end()) {
0616 SetupLEDHists(1, digi.id(), hoHists.LEDTRENDS);
0617 }
0618 LedHOHists(digi.id(), digi, hoHists.LEDTRENDS, cond);
0619 }
0620
0621
0622 if (hf.empty()) {
0623 edm::LogError("HcalLedAnalysis") << "Event with " << (int)hf.size() << " HF Digis passed.";
0624 return;
0625 }
0626 for (HFDigiCollection::const_iterator j = hf.begin(); j != hf.end(); ++j) {
0627 const HFDataFrame digi = (const HFDataFrame)(*j);
0628 _meol = hfHists.LEDTRENDS.find(digi.id());
0629 if (_meol == hfHists.LEDTRENDS.end()) {
0630 SetupLEDHists(2, digi.id(), hfHists.LEDTRENDS);
0631 }
0632 LedHFHists(digi.id(), digi, hfHists.LEDTRENDS, cond);
0633 }
0634
0635
0636 if (evt % m_nevtsample == 0)
0637 LedSampleAnalysis();
0638 }
0639
0640 void HcalLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
0641 string type = "HBHE";
0642 if (id == 1)
0643 type = "HO";
0644 if (id == 2)
0645 type = "HF";
0646
0647 _meol = toolT.find(detid);
0648 if (_meol == toolT.end()) {
0649
0650 map<int, LEDBUNCH> insert;
0651 char name[1024];
0652 for (int i = 0; i < 10; i++) {
0653 snprintf(name,
0654 sizeof name,
0655 "%s Pulse height, eta=%d phi=%d depth=%d TS=%d",
0656 type.c_str(),
0657 detid.ieta(),
0658 detid.iphi(),
0659 detid.depth(),
0660 i);
0661 insert[i].first = new TH1F(name, name, 200, 0., 2000.);
0662 }
0663 snprintf(name,
0664 sizeof name,
0665 "%s LED Mean pulse, eta=%d phi=%d depth=%d",
0666 type.c_str(),
0667 detid.ieta(),
0668 detid.iphi(),
0669 detid.depth());
0670 insert[10].first = new TH1F(name, name, 10, -0.5, 9.5);
0671 snprintf(name,
0672 sizeof name,
0673 "%s LED Pulse, eta=%d phi=%d depth=%d",
0674 type.c_str(),
0675 detid.ieta(),
0676 detid.iphi(),
0677 detid.depth());
0678 insert[11].first = new TH1F(name, name, 10, -0.5, 9.5);
0679 snprintf(name,
0680 sizeof name,
0681 "%s Mean TS, eta=%d phi=%d depth=%d",
0682 type.c_str(),
0683 detid.ieta(),
0684 detid.iphi(),
0685 detid.depth());
0686 insert[12].first = new TH1F(name, name, 200, 0., 10.);
0687 snprintf(name,
0688 sizeof name,
0689 "%s Peak TS, eta=%d phi=%d depth=%d",
0690 type.c_str(),
0691 detid.ieta(),
0692 detid.iphi(),
0693 detid.depth());
0694 insert[13].first = new TH1F(name, name, 200, 0., 10.);
0695 snprintf(name,
0696 sizeof name,
0697 "%s Peak TS error, eta=%d phi=%d depth=%d",
0698 type.c_str(),
0699 detid.ieta(),
0700 detid.iphi(),
0701 detid.depth());
0702 insert[14].first = new TH1F(name, name, 200, 0., 0.05);
0703 snprintf(name,
0704 sizeof name,
0705 "%s Fit chi2, eta=%d phi=%d depth=%d",
0706 type.c_str(),
0707 detid.ieta(),
0708 detid.iphi(),
0709 detid.depth());
0710 insert[15].first = new TH1F(name, name, 100, 0., 50.);
0711 snprintf(name,
0712 sizeof name,
0713 "%s Integrated Signal, eta=%d phi=%d depth=%d",
0714 type.c_str(),
0715 detid.ieta(),
0716 detid.iphi(),
0717 detid.depth());
0718 insert[16].first = new TH1F(name, name, 500, 0., 5000.);
0719
0720 toolT[detid] = insert;
0721 }
0722 }
0723
0724 void HcalLedAnalysis::LedHBHEHists(const HcalDetId& detid,
0725 const HBHEDataFrame& ledDigi,
0726 map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0727 const HcalDbService& cond) {
0728 map<int, LEDBUNCH> _mei;
0729 _meol = toolT.find(detid);
0730 _mei = _meol->second;
0731
0732
0733 if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0734 for (int k = 0; k < (int)state.size(); k++)
0735 state[k] = false;
0736 for (int i = 0; i < 16; i++)
0737 _mei[i].first->Reset();
0738 }
0739
0740
0741
0742
0743 float max_fC = 0;
0744 float ta = 0;
0745 m_coder = cond.getHcalCoder(detid);
0746 m_ped = cond.getPedestal(detid);
0747 m_shape = cond.getHcalShape(m_coder);
0748 for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0749 int capid = ledDigi[TS].capid();
0750 int adc = ledDigi[TS].adc();
0751 double fC = m_coder->charge(*m_shape, adc, capid);
0752 ta = (fC - m_ped->getValue(capid));
0753
0754 _mei[TS].first->Fill(ta);
0755 _mei[10].first->AddBinContent(TS + 1, ta);
0756 if (m_fitflag > 1) {
0757 if (TS == m_startTS)
0758 _mei[11].first->Reset();
0759 _mei[11].first->SetBinContent(TS + 1, ta);
0760 }
0761
0762 if (ta > max_fC) {
0763 max_fC = ta;
0764
0765 }
0766 }
0767
0768
0769
0770
0771
0772
0773 float sum = 0.;
0774 for (int i = 0; i < 10; i++)
0775 sum = sum + _mei[11].first->GetBinContent(i + 1);
0776 if (sum > 100) {
0777 if (m_fitflag == 2 || m_fitflag == 4) {
0778 float timmean = _mei[11].first->GetMean();
0779 float timmeancorr = BinsizeCorr(timmean);
0780 _mei[12].first->Fill(timmeancorr);
0781 }
0782 _mei[16].first->Fill(
0783 _mei[11].first->Integral());
0784 if (m_fitflag == 3 || m_fitflag == 4) {
0785 _mei[11].first->Fit("landau", "Q");
0786 TF1* fit = _mei[11].first->GetFunction("landau");
0787 _mei[13].first->Fill(fit->GetParameter(1));
0788 _mei[14].first->Fill(fit->GetParError(1));
0789 _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0790 }
0791 }
0792 }
0793
0794
0795 void HcalLedAnalysis::LedHOHists(const HcalDetId& detid,
0796 const HODataFrame& ledDigi,
0797 map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0798 const HcalDbService& cond) {
0799 map<int, LEDBUNCH> _mei;
0800 _meol = toolT.find(detid);
0801 _mei = _meol->second;
0802
0803 if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0804 for (int k = 0; k < (int)state.size(); k++)
0805 state[k] = false;
0806 for (int i = 0; i < 16; i++)
0807 _mei[i].first->Reset();
0808 }
0809
0810
0811
0812
0813
0814 float max_fC = 0;
0815 float ta = 0;
0816 m_coder = cond.getHcalCoder(detid);
0817 m_ped = cond.getPedestal(detid);
0818 m_shape = cond.getHcalShape(m_coder);
0819 for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0820 int capid = ledDigi[TS].capid();
0821 int adc = ledDigi[TS].adc();
0822 double fC = m_coder->charge(*m_shape, adc, capid);
0823 ta = (fC - m_ped->getValue(capid));
0824 _mei[TS].first->Fill(ta);
0825 _mei[10].first->AddBinContent(TS + 1, ta);
0826 if (m_fitflag > 1) {
0827 if (TS == m_startTS)
0828 _mei[11].first->Reset();
0829 _mei[11].first->SetBinContent(TS + 1, ta);
0830 }
0831
0832 if (ta > max_fC) {
0833 max_fC = ta;
0834
0835 }
0836 }
0837
0838
0839
0840
0841
0842
0843 float sum = 0.;
0844 for (int i = 0; i < 10; i++)
0845 sum = sum + _mei[11].first->GetBinContent(i + 1);
0846 if (sum > 100) {
0847 if (m_fitflag == 2 || m_fitflag == 4) {
0848 float timmean = _mei[11].first->GetMean();
0849 float timmeancorr = BinsizeCorr(timmean);
0850 _mei[12].first->Fill(timmeancorr);
0851 }
0852 _mei[16].first->Fill(
0853 _mei[11].first->Integral());
0854 if (m_fitflag == 3 || m_fitflag == 4) {
0855 _mei[11].first->Fit("landau", "Q");
0856 TF1* fit = _mei[11].first->GetFunction("landau");
0857 _mei[13].first->Fill(fit->GetParameter(1));
0858 _mei[14].first->Fill(fit->GetParError(1));
0859 _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0860 }
0861 }
0862 }
0863
0864
0865 void HcalLedAnalysis::LedHFHists(const HcalDetId& detid,
0866 const HFDataFrame& ledDigi,
0867 map<HcalDetId, map<int, LEDBUNCH> >& toolT,
0868 const HcalDbService& cond) {
0869 map<int, LEDBUNCH> _mei;
0870 _meol = toolT.find(detid);
0871 _mei = _meol->second;
0872
0873 if ((evt - 1) % m_nevtsample == 0 && state[0]) {
0874 for (int k = 0; k < (int)state.size(); k++)
0875 state[k] = false;
0876 for (int i = 0; i < 16; i++)
0877 _mei[i].first->Reset();
0878 }
0879
0880
0881
0882
0883
0884 float max_fC = 0;
0885 float ta = 0;
0886 m_coder = cond.getHcalCoder(detid);
0887 m_ped = cond.getPedestal(detid);
0888 m_shape = cond.getHcalShape(m_coder);
0889
0890 for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
0891 int capid = ledDigi[TS].capid();
0892
0893 int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
0894 if (adc < 0) {
0895 adc = 0;
0896 }
0897
0898 double fC = m_coder->charge(*m_shape, adc, capid);
0899
0900 ta = fC;
0901
0902 _mei[TS].first->Fill(ta);
0903 _mei[10].first->AddBinContent(TS + 1, ta);
0904 if (m_fitflag > 1) {
0905 if (TS == m_startTS)
0906 _mei[11].first->Reset();
0907 _mei[11].first->SetBinContent(TS + 1, ta);
0908 }
0909
0910
0911 if (ta > max_fC) {
0912 max_fC = ta;
0913
0914 }
0915 }
0916
0917
0918
0919
0920
0921
0922 float sum = 0.;
0923 for (int i = 0; i < 10; i++)
0924 sum = sum + _mei[11].first->GetBinContent(i + 1);
0925 if (sum > 100) {
0926 if (m_fitflag == 2 || m_fitflag == 4) {
0927 float timmean = _mei[11].first->GetMean();
0928 float timmeancorr = BinsizeCorr(timmean);
0929 _mei[12].first->Fill(timmeancorr);
0930 }
0931 _mei[16].first->Fill(
0932 _mei[11].first->Integral());
0933 if (m_fitflag == 3 || m_fitflag == 4) {
0934 _mei[11].first->Fit("landau", "Q");
0935 TF1* fit = _mei[11].first->GetFunction("landau");
0936 _mei[13].first->Fill(fit->GetParameter(1));
0937 _mei[14].first->Fill(fit->GetParError(1));
0938 _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
0939 }
0940 }
0941 }
0942
0943
0944 float HcalLedAnalysis::BinsizeCorr(float time) {
0945
0946
0947
0948
0949 float corrtime = 0.;
0950 static const float tstrue[32] = {0.003, 0.03425, 0.06548, 0.09675, 0.128, 0.15925, 0.1905, 0.22175,
0951 0.253, 0.28425, 0.3155, 0.34675, 0.378, 0.40925, 0.4405, 0.47175,
0952 0.503, 0.53425, 0.5655, 0.59675, 0.628, 0.65925, 0.6905, 0.72175,
0953 0.753, 0.78425, 0.8155, 0.84675, 0.878, 0.90925, 0.9405, 0.97175};
0954 static const float tsreco[32] = {-0.00422, 0.01815, 0.04409, 0.07346, 0.09799, 0.12192, 0.15072, 0.18158,
0955 0.21397, 0.24865, 0.28448, 0.31973, 0.35449, 0.39208, 0.43282, 0.47244,
0956 0.5105, 0.55008, 0.58827, 0.62828, 0.6717, 0.70966, 0.74086, 0.77496,
0957 0.80843, 0.83472, 0.86044, 0.8843, 0.90674, 0.92982, 0.95072, 0.9726};
0958
0959 int inttime = (int)time;
0960 float restime = time - inttime;
0961 for (int i = 0; i <= 32; i++) {
0962 float lolim = 0.;
0963 float uplim = 1.;
0964 float tsdown;
0965 float tsup;
0966 if (i > 0) {
0967 lolim = tsreco[i - 1];
0968 tsdown = tstrue[i - 1];
0969 } else
0970 tsdown = tstrue[31] - 1.;
0971 if (i < 32) {
0972 uplim = tsreco[i];
0973 tsup = tstrue[i];
0974 } else
0975 tsup = tstrue[0] + 1.;
0976 if (restime >= lolim && restime < uplim) {
0977 corrtime = (tsdown * (uplim - restime) + tsup * (restime - lolim)) / (uplim - lolim);
0978 }
0979 }
0980 corrtime += inttime;
0981
0982 return corrtime;
0983 }
0984
0985
0986
0987
0988
0989 void HcalLedAnalysis::ProcessCalibEvent(int fiberChan, HcalCalibDetId calibId, const HcalCalibDataFrame& digi) {
0990 _meca = calibHists.find(calibId);
0991 if (_meca == calibHists.end()) {
0992
0993 char name[1024];
0994 std::string prefix;
0995 if (calibId.calibFlavor() == HcalCalibDetId::CalibrationBox) {
0996 std::string sector = (calibId.hcalSubdet() == HcalBarrel) ? ("HB")
0997 : (calibId.hcalSubdet() == HcalEndcap) ? ("HE")
0998 : (calibId.hcalSubdet() == HcalOuter) ? ("HO")
0999 : (calibId.hcalSubdet() == HcalForward) ? ("HF")
1000 : "";
1001 snprintf(name,
1002 sizeof name,
1003 "%s %+d iphi=%d %s",
1004 sector.c_str(),
1005 calibId.ieta(),
1006 calibId.iphi(),
1007 calibId.cboxChannelString().c_str());
1008 prefix = name;
1009 }
1010
1011 snprintf(name, sizeof name, "%s Pin Diode Mean", prefix.c_str());
1012 calibHists[calibId].avePulse = new TProfile(name, name, 10, -0.5, 9.5, 0, 1000);
1013 snprintf(name, sizeof name, "%s Pin Diode Current Pulse", prefix.c_str());
1014 calibHists[calibId].thisPulse = new TH1F(name, name, 10, -0.5, 9.5);
1015 snprintf(name, sizeof name, "%s Pin Diode Integrated Pulse", prefix.c_str());
1016 calibHists[calibId].integPulse = new TH1F(name, name, 200, 0, 500);
1017 } else {
1018 for (int i = m_startTS; i < digi.size() && i <= m_endTS; i++) {
1019 calibHists[calibId].avePulse->Fill(i, digi.sample(i).adc());
1020 calibHists[calibId].thisPulse->SetBinContent(i + 1, digi.sample(i).adc());
1021 }
1022 calibHists[calibId].integPulse->Fill(calibHists[calibId].thisPulse->Integral());
1023 }
1024 }