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