CMS 3D CMS Logo
HcalLedAnalysis.cc
Go to the documentation of this file.
1 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
2 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
3 #include "CondFormats/HcalObjects/interface/HcalPedestals.h"
4 #include "FWCore/MessageLogger/interface/MessageLogger.h"
5 #include "CalibCalorimetry/HcalAlgos/interface/HcalLedAnalysis.h"
6 #include "TFile.h"
7 #include <cmath>
8 using namespace std;
9 
10 HcalLedAnalysis::HcalLedAnalysis(const edm::ParameterSet& ps) {
11  // init
12 
13  m_coder = nullptr;
14  m_ped = nullptr;
15  m_shape = nullptr;
16  evt = 0;
17  sample = 0;
18  m_file = nullptr;
19  // output files
20  for (int k = 0; k < 4; k++)
21  state.push_back(true); // 4 cap-ids (do we care?)
22  m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
23  m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML", "");
24  if (!m_outputFileText.empty()) {
25  edm::LogInfo("HcalLedAnalysis") << "Hcal LED results will be saved to " << m_outputFileText.c_str() << endl;
26  m_outFile.open(m_outputFileText.c_str());
27  }
28  m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
29  if (!m_outputFileROOT.empty()) {
30  edm::LogInfo("HcalLedAnalysis") << "Hcal LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
31  }
32 
33  m_nevtsample = ps.getUntrackedParameter<int>("nevtsample", 9999999);
34  if (m_nevtsample < 1)
35  m_nevtsample = 9999999;
36  m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag", 0);
37  if (m_hiSaveflag < 0)
38  m_hiSaveflag = 0;
39  if (m_hiSaveflag > 0)
40  m_hiSaveflag = 1;
41  m_fitflag = ps.getUntrackedParameter<int>("analysisflag", 2);
42  if (m_fitflag < 0)
43  m_fitflag = 0;
44  if (m_fitflag > 4)
45  m_fitflag = 4;
46  m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
47  if (m_startTS < 0)
48  m_startTS = 0;
49  m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
50  m_usecalib = ps.getUntrackedParameter<bool>("usecalib", false);
51 
52  int runNum = ps.getUntrackedParameter<int>("runNumber", 999999);
53 
54  // histogram booking
55  hbHists.ALLLEDS = new TH1F("HBHE All LEDs", "HB/HE All Leds", 10, 0, 9);
56  hbHists.LEDRMS = new TH1F("HBHE All LED RMS", "HB/HE All LED RMS", 100, 0, 3);
57  hbHists.LEDMEAN = new TH1F("HBHE All LED Means", "HB/HE All LED Means", 100, 0, 9);
58  hbHists.CHI2 = new TH1F("HBHE Chi2 by ndf for Landau fit", "HB/HE Chi2/ndf Landau", 200, 0., 50.);
59 
60  hoHists.ALLLEDS = new TH1F("HO All LEDs", "HO All Leds", 10, 0, 9);
61  hoHists.LEDRMS = new TH1F("HO All LED RMS", "HO All LED RMS", 100, 0, 3);
62  hoHists.LEDMEAN = new TH1F("HO All LED Means", "HO All LED Means", 100, 0, 9);
63  hoHists.CHI2 = new TH1F("HO Chi2 by ndf for Landau fit", "HO Chi2/ndf Landau", 200, 0., 50.);
64 
65  hfHists.ALLLEDS = new TH1F("HF All LEDs", "HF All Leds", 10, 0, 9);
66  hfHists.LEDRMS = new TH1F("HF All LED RMS", "HF All LED RMS", 100, 0, 3);
67  hfHists.LEDMEAN = new TH1F("HF All LED Means", "HF All LED Means", 100, 0, 9);
68  hfHists.CHI2 = new TH1F("HF Chi2 by ndf for Landau fit", "HF Chi2/ndf Landau", 200, 0., 50.);
69 
70  char output[100]{0};
71 
72  //XML file header
73  m_outputFileXML.open(m_outputFileX.c_str());
74 
75  m_outputFileXML << "<?xml version='1.0' encoding='UTF-8'?>" << endl;
76 
77  m_outputFileXML << "<ROOT>" << endl;
78 
79  m_outputFileXML << " <HEADER>" << endl;
80 
81  m_outputFileXML << " <TYPE>" << endl;
82 
83  m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>" << endl;
84 
85  m_outputFileXML << " <NAME>HCAL LED Timing</NAME>" << endl;
86 
87  m_outputFileXML << " </TYPE>" << endl;
88 
89  m_outputFileXML << " <RUN>" << endl;
90 
91  m_outputFileXML << " <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>" << endl;
92 
93  snprintf(output, sizeof output, " <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
94  m_outputFileXML << output << endl;
95 
96  m_outputFileXML << " <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>" << endl;
97 
98  m_outputFileXML << " <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
99 
100  m_outputFileXML << " </RUN>" << endl;
101 
102  m_outputFileXML << " </HEADER>" << endl;
103 
104  m_outputFileXML << "<!-- Tags secton -->" << endl;
105 
106  m_outputFileXML << " <ELEMENTS>" << endl;
107 
108  m_outputFileXML << " <DATA_SET id='-1'/>" << endl;
109 
110  m_outputFileXML << " <IOV id='1'>" << endl;
111 
112  m_outputFileXML << " <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>" << endl;
113 
114  m_outputFileXML << " <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>" << endl;
115 
116  m_outputFileXML << " </IOV>" << endl;
117 
118  m_outputFileXML << " <TAG id='2' mode='auto'>" << endl;
119 
120  snprintf(output, sizeof output, " <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
121  m_outputFileXML << output << endl;
122 
123  m_outputFileXML << " <DETECTOR_NAME>HCAL</DETECTOR_NAME>" << endl;
124 
125  m_outputFileXML << " <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>" << endl;
126 
127  m_outputFileXML << " </TAG>" << endl;
128 
129  m_outputFileXML << " </ELEMENTS>" << endl;
130 
131  m_outputFileXML << " <MAPS>" << endl;
132 
133  m_outputFileXML << " <TAG idref ='2'>" << endl;
134 
135  m_outputFileXML << " <IOV idref='1'>" << endl;
136 
137  m_outputFileXML << " <DATA_SET idref='-1' />" << endl;
138 
139  m_outputFileXML << " </IOV>" << endl;
140 
141  m_outputFileXML << " </TAG>" << endl;
142 
143  m_outputFileXML << " </MAPS>" << endl;
144 }
145 
146 //-----------------------------------------------------------------------------
147 HcalLedAnalysis::~HcalLedAnalysis() {
149  for (_meol = hbHists.LEDTRENDS.begin(); _meol != hbHists.LEDTRENDS.end(); _meol++) {
150  for (int i = 0; i < 15; i++)
151  _meol->second[i].first->Delete();
152  }
153  for (_meol = hoHists.LEDTRENDS.begin(); _meol != hoHists.LEDTRENDS.end(); _meol++) {
154  for (int i = 0; i < 15; i++)
155  _meol->second[i].first->Delete();
156  }
157  for (_meol = hfHists.LEDTRENDS.begin(); _meol != hfHists.LEDTRENDS.end(); _meol++) {
158  for (int i = 0; i < 15; i++)
159  _meol->second[i].first->Delete();
160  }
161  hbHists.ALLLEDS->Delete();
162  hbHists.LEDRMS->Delete();
163  hbHists.LEDMEAN->Delete();
164  hbHists.CHI2->Delete();
165 
166  hoHists.ALLLEDS->Delete();
167  hoHists.LEDRMS->Delete();
168  hoHists.LEDMEAN->Delete();
169  hoHists.CHI2->Delete();
170 
171  hfHists.ALLLEDS->Delete();
172  hfHists.LEDRMS->Delete();
173  hfHists.LEDMEAN->Delete();
174  hfHists.CHI2->Delete();
175 }
176 
177 //-----------------------------------------------------------------------------
178 void HcalLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
179  // open the histogram file, create directories within
180  m_file = new TFile(m_outputFileROOT.c_str(), "RECREATE");
181  m_file->mkdir("HBHE");
182  m_file->cd();
183  m_file->mkdir("HO");
184  m_file->cd();
185  m_file->mkdir("HF");
186  m_file->cd();
187  m_file->mkdir("Calib");
188  m_file->cd();
189 }
190 
191 //-----------------------------------------------------------------------------
192 /*
193 void HcalLedAnalysis::doPeds(const HcalPedestal* fInputPedestals){
194 // put all pedestals in a map m_AllPedVals, to be used in processLedEvent -
195 // sorry, this is the only way I was able to implement pedestal subtraction
196 
197 // DEPRECATED
198 // This is no longer useful, better ways of doing it -A
199  HcalDetId detid;
200  map<int,float> PedVals;
201  pedCan = fInputPedestals;
202  if(pedCan){
203  std::vector<DetId> Channs=pedCan->getAllChannels();
204  for (int i=0; i<(int)Channs.size(); i++){
205  detid=HcalDetId (Channs[i]);
206  for (int icap=0; icap<4; icap++) PedVals[icap]=pedCan->getValue(detid,icap);
207  m_AllPedVals[detid]=PedVals;
208  }
209  }
210 }
211 */
212 //-----------------------------------------------------------------------------
213 void HcalLedAnalysis::GetLedConst(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
214  double time2 = 0;
215  double time1 = 0;
216  double time3 = 0;
217  double time4 = 0;
218  double dtime2 = 0;
219  double dtime1 = 0;
220  double dtime3 = 0;
221  double dtime4 = 0;
222  char output[256]{0};
223 
224  if (!m_outputFileText.empty()) {
225  if (m_fitflag == 0 || m_fitflag == 2)
226  m_outFile << "Det Eta,Phi,D Mean Error" << std::endl;
227  else if (m_fitflag == 1 || m_fitflag == 3)
228  m_outFile << "Det Eta,Phi,D Peak Error" << std::endl;
229  else if (m_fitflag == 4)
230  m_outFile << "Det Eta,Phi,D Mean Error Peak Error MeanEv Error PeakEv Error"
231  << std::endl;
232  }
233  for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
234  // scale the LED pulse to 1 event
235  _meol->second[10].first->Scale(1. / evt_curr);
236  if (m_fitflag == 0 || m_fitflag == 4) {
237  time1 = _meol->second[10].first->GetMean();
238  dtime1 = _meol->second[10].first->GetRMS() / sqrt((float)evt_curr * (m_endTS - m_startTS + 1));
239  }
240  if (m_fitflag == 1 || m_fitflag == 4) {
241  // put proper errors
242  for (int j = 0; j < 10; j++)
243  _meol->second[10].first->SetBinError(j + 1, _meol->second[j].first->GetRMS() / sqrt((float)evt_curr));
244  }
245  if (m_fitflag == 1 || m_fitflag == 3 || m_fitflag == 4) {
246  _meol->second[10].first->Fit("landau", "Q");
247  // _meol->second[10].first->Fit("gaus","Q");
248  TF1* fit = _meol->second[10].first->GetFunction("landau");
249  // TF1 *fit = _meol->second[10].first->GetFunction("gaus");
250  time2 = fit->GetParameter(1);
251  dtime2 = fit->GetParError(1);
252  }
253  if (m_fitflag == 2 || m_fitflag == 4) {
254  time3 = _meol->second[12].first->GetMean();
255  dtime3 = _meol->second[12].first->GetRMS() / sqrt((float)_meol->second[12].first->GetEntries());
256  }
257  if (m_fitflag == 3 || m_fitflag == 4) {
258  time4 = _meol->second[13].first->GetMean();
259  dtime4 = _meol->second[13].first->GetRMS() / sqrt((float)_meol->second[13].first->GetEntries());
260  }
261  for (int i = 0; i < 10; i++) {
262  _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
263  _meol->second[i].first->GetYaxis()->SetTitle("Counts");
264  // if(m_hiSaveflag>0)_meol->second[i].first->Write();
265  }
266  _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
267  _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
268  if (m_hiSaveflag > 0)
269  _meol->second[10].first->Write();
270  _meol->second[10].second.first[0].push_back(time1);
271  _meol->second[10].second.first[1].push_back(dtime1);
272  _meol->second[11].second.first[0].push_back(time2);
273  _meol->second[11].second.first[1].push_back(dtime2);
274  _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
275  _meol->second[12].first->GetYaxis()->SetTitle("Counts");
276  if (m_fitflag == 2 && m_hiSaveflag > 0)
277  _meol->second[12].first->Write();
278  _meol->second[12].second.first[0].push_back(time3);
279  _meol->second[12].second.first[1].push_back(dtime3);
280  _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
281  _meol->second[13].first->GetYaxis()->SetTitle("Counts");
282  if (m_fitflag > 2 && m_hiSaveflag > 0)
283  _meol->second[13].first->Write();
284  _meol->second[13].second.first[0].push_back(time4);
285  _meol->second[13].second.first[1].push_back(dtime4);
286  _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
287  _meol->second[14].first->GetYaxis()->SetTitle("Counts");
288  if (m_fitflag > 2 && m_hiSaveflag > 0)
289  _meol->second[14].first->Write();
290  _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
291  _meol->second[15].first->GetYaxis()->SetTitle("Counts");
292  if (m_fitflag > 2 && m_hiSaveflag > 0)
293  _meol->second[15].first->Write();
294  _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
295  _meol->second[16].first->Write();
296 
297  // Ascii printout (need to modify to include new info)
298  HcalDetId detid = _meol->first;
299 
300  if (!m_outputFileText.empty()) {
301  if (m_fitflag == 0) {
302  m_outFile << detid << " " << time1 << " " << dtime1 << std::endl;
303  m_outputFileXML << " <DATA_SET>" << endl;
304  m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
305  m_outputFileXML << " <CHANNEL>" << endl;
306  m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
307  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
308  m_outputFileXML << output << endl;
309  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
310  m_outputFileXML << output << endl;
311  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
312  m_outputFileXML << output << endl;
313  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
314  m_outputFileXML << output << endl;
315 
316  if (detid.subdet() == 1)
317  m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
318  if (detid.subdet() == 2)
319  m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
320  if (detid.subdet() == 3)
321  m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
322  if (detid.subdet() == 4)
323  m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
324  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
325  m_outputFileXML << output << endl;
326  m_outputFileXML << " </CHANNEL>" << endl;
327  m_outputFileXML << " <DATA>" << endl;
328  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time1);
329  m_outputFileXML << output << endl;
330  m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
331  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
332  m_outputFileXML << output << endl;
333  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
334  m_outputFileXML << output << endl;
335  m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
336  m_outputFileXML << " </DATA>" << endl;
337  m_outputFileXML << " </DATA_SET>" << endl;
338 
339  } else if (m_fitflag == 1) {
340  m_outFile << detid << " " << time2 << " " << dtime2 << std::endl;
341  m_outputFileXML << " <DATA_SET>" << endl;
342  m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
343  m_outputFileXML << " <CHANNEL>" << endl;
344  m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
345  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
346  m_outputFileXML << output << endl;
347  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
348  m_outputFileXML << output << endl;
349  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
350  m_outputFileXML << output << endl;
351  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
352  m_outputFileXML << output << endl;
353  if (detid.subdet() == 1)
354  m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
355  if (detid.subdet() == 2)
356  m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
357  if (detid.subdet() == 3)
358  m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
359  if (detid.subdet() == 4)
360  m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
361  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
362  m_outputFileXML << output << endl;
363  m_outputFileXML << " </CHANNEL>" << endl;
364  m_outputFileXML << " <DATA>" << endl;
365  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time2);
366  m_outputFileXML << output << endl;
367  m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
368  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime2);
369  m_outputFileXML << output << endl;
370  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
371  m_outputFileXML << output << endl;
372  m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
373  m_outputFileXML << " </DATA>" << endl;
374  m_outputFileXML << " </DATA_SET>" << endl;
375  }
376 
377  else if (m_fitflag == 2) {
378  m_outFile << detid << " " << time3 << " " << dtime3 << std::endl;
379  m_outputFileXML << " <DATA_SET>" << endl;
380  m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
381  m_outputFileXML << " <CHANNEL>" << endl;
382  m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
383  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
384  m_outputFileXML << output << endl;
385  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
386  m_outputFileXML << output << endl;
387  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
388  m_outputFileXML << output << endl;
389  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
390  m_outputFileXML << output << endl;
391  if (detid.subdet() == 1)
392  m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
393  if (detid.subdet() == 2)
394  m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
395  if (detid.subdet() == 3)
396  m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
397  if (detid.subdet() == 4)
398  m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
399  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
400  m_outputFileXML << output << endl;
401  m_outputFileXML << " </CHANNEL>" << endl;
402  m_outputFileXML << " <DATA>" << endl;
403  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time3);
404  m_outputFileXML << output << endl;
405  m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
406  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime3);
407  m_outputFileXML << output << endl;
408  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
409  m_outputFileXML << output << endl;
410  m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
411  m_outputFileXML << " </DATA>" << endl;
412  m_outputFileXML << " </DATA_SET>" << endl;
413  } else if (m_fitflag == 3) {
414  m_outFile << detid << " " << time4 << " " << dtime4 << std::endl;
415  m_outputFileXML << " <DATA_SET>" << endl;
416  m_outputFileXML << " <VERSION>version:1</VERSION>" << endl;
417  m_outputFileXML << " <CHANNEL>" << endl;
418  m_outputFileXML << " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>" << endl;
419  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
420  m_outputFileXML << output << endl;
421  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
422  m_outputFileXML << output << endl;
423  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
424  m_outputFileXML << output << endl;
425  sprintf(output, " <Z>%2i</Z>", detid.zside());
426  m_outputFileXML << output << endl;
427  if (detid.subdet() == 1)
428  m_outputFileXML << " <DETECTOR_NAME>HB</DETECTOR_NAME>" << endl;
429  if (detid.subdet() == 2)
430  m_outputFileXML << " <DETECTOR_NAME>HE</DETECTOR_NAME>" << endl;
431  if (detid.subdet() == 3)
432  m_outputFileXML << " <DETECTOR_NAME>HO</DETECTOR_NAME>" << endl;
433  if (detid.subdet() == 4)
434  m_outputFileXML << " <DETECTOR_NAME>HF</DETECTOR_NAME>" << endl;
435  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
436  m_outputFileXML << output << endl;
437  m_outputFileXML << " </CHANNEL>" << endl;
438  m_outputFileXML << " <DATA>" << endl;
439  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time4);
440  m_outputFileXML << output << endl;
441  m_outputFileXML << " <OFFSET_TIME> 0</OFFSET_TIME>" << endl;
442  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
443  m_outputFileXML << output << endl;
444  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
445  m_outputFileXML << output << endl;
446  m_outputFileXML << " <STATUS_WORD> 0</STATUS_WORD>" << endl;
447  m_outputFileXML << " </DATA>" << endl;
448  m_outputFileXML << " </DATA_SET>" << endl;
449  }
450 
451  else if (m_fitflag == 4) {
452  m_outFile << detid << " " << time1 << " " << dtime1 << " " << time2 << " " << dtime2 << " " << time3
453  << " " << dtime3 << " " << time4 << " " << dtime4 << std::endl;
454  }
455  }
456  }
457 }
458 
459 //-----------------------------------------------------------------------------
460 void HcalLedAnalysis::LedSampleAnalysis() {
461  // it is called every m_nevtsample events (a sample) and the end of run
462  char LedSampleNum[20];
463 
464  snprintf(LedSampleNum, sizeof LedSampleNum, "LedSample_%d", sample);
465  m_file->cd();
466  m_file->mkdir(LedSampleNum);
467  m_file->cd(LedSampleNum);
468 
469  // Compute LED constants for each HB/HE, HO, HF
470  GetLedConst(hbHists.LEDTRENDS);
471  GetLedConst(hoHists.LEDTRENDS);
472  GetLedConst(hfHists.LEDTRENDS);
473 }
474 
475 //-----------------------------------------------------------------------------
476 void HcalLedAnalysis::LedTrendings(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
477  for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
478  char name[1024];
479  HcalDetId detid = _meol->first;
480  snprintf(name, sizeof name, "LED timing trend, eta=%d phi=%d depth=%d", detid.ieta(), detid.iphi(), detid.depth());
481  int bins = _meol->second[10 + m_fitflag].second.first[0].size();
482  float lo = 0.5;
483  float hi = (float)bins + 0.5;
484  _meol->second[10 + m_fitflag].second.second.push_back(new TH1F(name, name, bins, lo, hi));
485 
486  std::vector<double>::iterator sample_it;
487  // LED timing - put content and errors
488  int j = 0;
489  for (sample_it = _meol->second[10 + m_fitflag].second.first[0].begin();
490  sample_it != _meol->second[10 + m_fitflag].second.first[0].end();
491  ++sample_it) {
492  _meol->second[10 + m_fitflag].second.second[0]->SetBinContent(++j, *sample_it);
493  }
494  j = 0;
495  for (sample_it = _meol->second[10 + m_fitflag].second.first[1].begin();
496  sample_it != _meol->second[10 + m_fitflag].second.first[1].end();
497  ++sample_it) {
498  _meol->second[10 + m_fitflag].second.second[0]->SetBinError(++j, *sample_it);
499  }
500  snprintf(name, sizeof name, "Sample (%d events)", m_nevtsample);
501  _meol->second[10 + m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
502  _meol->second[10 + m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
503  _meol->second[10 + m_fitflag].second.second[0]->Write();
504  }
505 }
506 
507 //-----------------------------------------------------------------------------
508 void HcalLedAnalysis::LedDone() {
509  // First process the last sample (remaining events).
510  if (evt % m_nevtsample != 0)
511  LedSampleAnalysis();
512 
513  // Now do the end of run analysis: trending histos
514  if (sample > 1 && m_fitflag != 4) {
515  m_file->cd();
516  m_file->cd("HBHE");
517  LedTrendings(hbHists.LEDTRENDS);
518  m_file->cd();
519  m_file->cd("HO");
520  LedTrendings(hoHists.LEDTRENDS);
521  m_file->cd();
522  m_file->cd("HF");
523  LedTrendings(hfHists.LEDTRENDS);
524  }
525 
526  // Write other histograms.
527  // HB
528  m_file->cd();
529  m_file->cd("HBHE");
530  hbHists.ALLLEDS->Write();
531  hbHists.LEDRMS->Write();
532  hbHists.LEDMEAN->Write();
533  // HO
534  m_file->cd();
535  m_file->cd("HO");
536  hoHists.ALLLEDS->Write();
537  hoHists.LEDRMS->Write();
538  hoHists.LEDMEAN->Write();
539  // HF
540  m_file->cd();
541  m_file->cd("HF");
542  hfHists.ALLLEDS->Write();
543  hfHists.LEDRMS->Write();
544  hfHists.LEDMEAN->Write();
545  // Calib
546  m_file->cd();
547  m_file->cd("Calib");
548  for (_meca = calibHists.begin(); _meca != calibHists.end(); _meca++) {
549  _meca->second.avePulse->Write();
550  _meca->second.integPulse->Write();
551  }
552 
553  // Write the histo file and close it
554  // m_file->Write();
555  m_file->Close();
556  edm::LogInfo("HcalLedAnalysis") << "Hcal histograms written to " << m_outputFileROOT.c_str() << endl;
557 }
558 
559 //-----------------------------------------------------------------------------
560 void HcalLedAnalysis::processLedEvent(const HBHEDigiCollection& hbhe,
561  const HODigiCollection& ho,
562  const HFDigiCollection& hf,
563  const HcalCalibDigiCollection& calib,
564  const HcalDbService& cond) {
565  evt++;
566  sample = (evt - 1) / m_nevtsample + 1;
567  evt_curr = evt % m_nevtsample;
568  if (evt_curr == 0)
569  evt_curr = m_nevtsample;
570 
571  // Calib
572 
573  if (m_usecalib) {
574  if (calib.empty()) {
575  edm::LogError("HcalLedAnalysis") << "Event with " << (int)calib.size() << " Calib Digis passed.";
576  return;
577  }
578  // this is effectively a loop over electronic channels
579  for (HcalCalibDigiCollection::const_iterator j = calib.begin(); j != calib.end(); ++j) {
580  const HcalCalibDataFrame digi = (const HcalCalibDataFrame)(*j);
581  HcalElectronicsId elecId = digi.elecId();
582  HcalCalibDetId calibId = digi.id();
583  ProcessCalibEvent(elecId.fiberChanId(),
584  calibId,
585  digi); //Shouldn't depend on anything in elecId but not sure how else to do it
586  }
587  }
588 
589  // HB + HE
590  if (hbhe.empty()) {
591  edm::LogError("HcalLedAnalysis") << "Event with " << (int)hbhe.size() << " HBHE Digis passed.";
592  return;
593  }
594  // this is effectively a loop over electronic channels
595  for (HBHEDigiCollection::const_iterator j = hbhe.begin(); j != hbhe.end(); ++j) {
596  const HBHEDataFrame digi = (const HBHEDataFrame)(*j);
597  for (int k = 0; k < (int)state.size(); k++)
598  state[k] = true;
599  // See if histos exist for this channel, and if not, create them
600  _meol = hbHists.LEDTRENDS.find(digi.id());
601  if (_meol == hbHists.LEDTRENDS.end()) {
602  SetupLEDHists(0, digi.id(), hbHists.LEDTRENDS);
603  }
604  LedHBHEHists(digi.id(), digi, hbHists.LEDTRENDS, cond);
605  }
606 
607  // HO
608  if (ho.empty()) {
609  edm::LogError("HcalLedAnalysis") << "Event with " << (int)ho.size() << " HO Digis passed.";
610  return;
611  }
612  for (HODigiCollection::const_iterator j = ho.begin(); j != ho.end(); ++j) {
613  const HODataFrame digi = (const HODataFrame)(*j);
614  _meol = hoHists.LEDTRENDS.find(digi.id());
615  if (_meol == hoHists.LEDTRENDS.end()) {
616  SetupLEDHists(1, digi.id(), hoHists.LEDTRENDS);
617  }
618  LedHOHists(digi.id(), digi, hoHists.LEDTRENDS, cond);
619  }
620 
621  // HF
622  if (hf.empty()) {
623  edm::LogError("HcalLedAnalysis") << "Event with " << (int)hf.size() << " HF Digis passed.";
624  return;
625  }
626  for (HFDigiCollection::const_iterator j = hf.begin(); j != hf.end(); ++j) {
627  const HFDataFrame digi = (const HFDataFrame)(*j);
628  _meol = hfHists.LEDTRENDS.find(digi.id());
629  if (_meol == hfHists.LEDTRENDS.end()) {
630  SetupLEDHists(2, digi.id(), hfHists.LEDTRENDS);
631  }
632  LedHFHists(digi.id(), digi, hfHists.LEDTRENDS, cond);
633  }
634 
635  // Call the function every m_nevtsample events
636  if (evt % m_nevtsample == 0)
637  LedSampleAnalysis();
638 }
639 //----------------------------------------------------------------------------
640 void HcalLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
641  string type = "HBHE";
642  if (id == 1)
643  type = "HO";
644  if (id == 2)
645  type = "HF";
646 
647  _meol = toolT.find(detid);
648  if (_meol == toolT.end()) {
649  // if histos for this channel do not exist, create them
650  map<int, LEDBUNCH> insert;
651  char name[1024];
652  for (int i = 0; i < 10; i++) {
653  snprintf(name,
654  sizeof name,
655  "%s Pulse height, eta=%d phi=%d depth=%d TS=%d",
656  type.c_str(),
657  detid.ieta(),
658  detid.iphi(),
659  detid.depth(),
660  i);
661  insert[i].first = new TH1F(name, name, 200, 0., 2000.);
662  }
663  snprintf(name,
664  sizeof name,
665  "%s LED Mean pulse, eta=%d phi=%d depth=%d",
666  type.c_str(),
667  detid.ieta(),
668  detid.iphi(),
669  detid.depth());
670  insert[10].first = new TH1F(name, name, 10, -0.5, 9.5);
671  snprintf(name,
672  sizeof name,
673  "%s LED Pulse, eta=%d phi=%d depth=%d",
674  type.c_str(),
675  detid.ieta(),
676  detid.iphi(),
677  detid.depth());
678  insert[11].first = new TH1F(name, name, 10, -0.5, 9.5);
679  snprintf(name,
680  sizeof name,
681  "%s Mean TS, eta=%d phi=%d depth=%d",
682  type.c_str(),
683  detid.ieta(),
684  detid.iphi(),
685  detid.depth());
686  insert[12].first = new TH1F(name, name, 200, 0., 10.);
687  snprintf(name,
688  sizeof name,
689  "%s Peak TS, eta=%d phi=%d depth=%d",
690  type.c_str(),
691  detid.ieta(),
692  detid.iphi(),
693  detid.depth());
694  insert[13].first = new TH1F(name, name, 200, 0., 10.);
695  snprintf(name,
696  sizeof name,
697  "%s Peak TS error, eta=%d phi=%d depth=%d",
698  type.c_str(),
699  detid.ieta(),
700  detid.iphi(),
701  detid.depth());
702  insert[14].first = new TH1F(name, name, 200, 0., 0.05);
703  snprintf(name,
704  sizeof name,
705  "%s Fit chi2, eta=%d phi=%d depth=%d",
706  type.c_str(),
707  detid.ieta(),
708  detid.iphi(),
709  detid.depth());
710  insert[15].first = new TH1F(name, name, 100, 0., 50.);
711  snprintf(name,
712  sizeof name,
713  "%s Integrated Signal, eta=%d phi=%d depth=%d",
714  type.c_str(),
715  detid.ieta(),
716  detid.iphi(),
717  detid.depth());
718  insert[16].first = new TH1F(name, name, 500, 0., 5000.);
719 
720  toolT[detid] = insert;
721  }
722 }
723 //-----------------------------------------------------------------------------
724 void HcalLedAnalysis::LedHBHEHists(const HcalDetId& detid,
725  const HBHEDataFrame& ledDigi,
726  map<HcalDetId, map<int, LEDBUNCH> >& toolT,
727  const HcalDbService& cond) {
728  map<int, LEDBUNCH> _mei;
729  _meol = toolT.find(detid);
730  _mei = _meol->second;
731 
732  // Reset the histos if we're at the end of a 'bunch'
733  if ((evt - 1) % m_nevtsample == 0 && state[0]) {
734  for (int k = 0; k < (int)state.size(); k++)
735  state[k] = false;
736  for (int i = 0; i < 16; i++)
737  _mei[i].first->Reset();
738  }
739 
740  // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
741 
742  // int maxTS = -1;
743  float max_fC = 0;
744  float ta = 0;
745  m_coder = cond.getHcalCoder(detid);
746  m_ped = cond.getPedestal(detid);
747  m_shape = cond.getHcalShape(m_coder);
748  for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
749  int capid = ledDigi[TS].capid();
750  int adc = ledDigi[TS].adc();
751  double fC = m_coder->charge(*m_shape, adc, capid);
752  ta = (fC - m_ped->getValue(capid));
753  //cout << "DetID: " << detid << " CapID: " << capid << " ADC: " << adc << " fC: " << fC << endl;
754  _mei[TS].first->Fill(ta);
755  _mei[10].first->AddBinContent(TS + 1, ta); // This is average pulse, could probably do better (Profile?)
756  if (m_fitflag > 1) {
757  if (TS == m_startTS)
758  _mei[11].first->Reset();
759  _mei[11].first->SetBinContent(TS + 1, ta);
760  }
761  // keep track of max TS and max amplitude (in fC)
762  if (ta > max_fC) {
763  max_fC = ta;
764  // maxTS = TS;
765  }
766  }
767 
768  // Now we have a sample with pedestals subtracted and in units of fC
769  // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
770  // we now want to use Phil's timing correction. This is not necessary
771  // if we are performing a Landau fit (m_fitflag = 3)
772 
773  float sum = 0.;
774  for (int i = 0; i < 10; i++)
775  sum = sum + _mei[11].first->GetBinContent(i + 1);
776  if (sum > 100) {
777  if (m_fitflag == 2 || m_fitflag == 4) {
778  float timmean = _mei[11].first->GetMean(); // let's use Phil's way instead
779  float timmeancorr = BinsizeCorr(timmean);
780  _mei[12].first->Fill(timmeancorr);
781  }
782  _mei[16].first->Fill(
783  _mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
784  if (m_fitflag == 3 || m_fitflag == 4) {
785  _mei[11].first->Fit("landau", "Q");
786  TF1* fit = _mei[11].first->GetFunction("landau");
787  _mei[13].first->Fill(fit->GetParameter(1));
788  _mei[14].first->Fill(fit->GetParError(1));
789  _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
790  }
791  }
792 }
793 
794 //-----------------------------------------------------------------------------
795 void HcalLedAnalysis::LedHOHists(const HcalDetId& detid,
796  const HODataFrame& ledDigi,
797  map<HcalDetId, map<int, LEDBUNCH> >& toolT,
798  const HcalDbService& cond) {
799  map<int, LEDBUNCH> _mei;
800  _meol = toolT.find(detid);
801  _mei = _meol->second;
802  // Rest the histos if we're at the end of a 'bunch'
803  if ((evt - 1) % m_nevtsample == 0 && state[0]) {
804  for (int k = 0; k < (int)state.size(); k++)
805  state[k] = false;
806  for (int i = 0; i < 16; i++)
807  _mei[i].first->Reset();
808  }
809 
810  // now we have the signal in fC, let's get rid of that darn pedestal
811  // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
812 
813  // int maxTS = -1;
814  float max_fC = 0;
815  float ta = 0;
816  m_coder = cond.getHcalCoder(detid);
817  m_ped = cond.getPedestal(detid);
818  m_shape = cond.getHcalShape(m_coder);
819  for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
820  int capid = ledDigi[TS].capid();
821  int adc = ledDigi[TS].adc();
822  double fC = m_coder->charge(*m_shape, adc, capid);
823  ta = (fC - m_ped->getValue(capid));
824  _mei[TS].first->Fill(ta);
825  _mei[10].first->AddBinContent(TS + 1, ta); // This is average pulse, could probably do better (Profile?)
826  if (m_fitflag > 1) {
827  if (TS == m_startTS)
828  _mei[11].first->Reset();
829  _mei[11].first->SetBinContent(TS + 1, ta);
830  }
831  // keep track of max TS and max amplitude (in fC)
832  if (ta > max_fC) {
833  max_fC = ta;
834  // maxTS = TS;
835  }
836  }
837 
838  // Now we have a sample with pedestals subtracted and in units of fC
839  // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
840  // we now want to use Phil's timing correction. This is not necessary
841  // if we are performing a Landau fit (m_fitflag = 3)
842 
843  float sum = 0.;
844  for (int i = 0; i < 10; i++)
845  sum = sum + _mei[11].first->GetBinContent(i + 1);
846  if (sum > 100) {
847  if (m_fitflag == 2 || m_fitflag == 4) {
848  float timmean = _mei[11].first->GetMean(); // let's use Phil's way instead
849  float timmeancorr = BinsizeCorr(timmean);
850  _mei[12].first->Fill(timmeancorr);
851  }
852  _mei[16].first->Fill(
853  _mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
854  if (m_fitflag == 3 || m_fitflag == 4) {
855  _mei[11].first->Fit("landau", "Q");
856  TF1* fit = _mei[11].first->GetFunction("landau");
857  _mei[13].first->Fill(fit->GetParameter(1));
858  _mei[14].first->Fill(fit->GetParError(1));
859  _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
860  }
861  }
862 }
863 
864 //-----------------------------------------------------------------------------
865 void HcalLedAnalysis::LedHFHists(const HcalDetId& detid,
866  const HFDataFrame& ledDigi,
867  map<HcalDetId, map<int, LEDBUNCH> >& toolT,
868  const HcalDbService& cond) {
869  map<int, LEDBUNCH> _mei;
870  _meol = toolT.find(detid);
871  _mei = _meol->second;
872  // Rest the histos if we're at the end of a 'bunch'
873  if ((evt - 1) % m_nevtsample == 0 && state[0]) {
874  for (int k = 0; k < (int)state.size(); k++)
875  state[k] = false;
876  for (int i = 0; i < 16; i++)
877  _mei[i].first->Reset();
878  }
879 
880  // now we have the signal in fC, let's get rid of that darn pedestal
881  // Most of this is borrowed from HcalSimpleReconstructor, so thanks Jeremy/Phil
882 
883  // int maxTS = -1;
884  float max_fC = 0;
885  float ta = 0;
886  m_coder = cond.getHcalCoder(detid);
887  m_ped = cond.getPedestal(detid);
888  m_shape = cond.getHcalShape(m_coder);
889  //cout << "New Digi!!!!!!!!!!!!!!!!!!!!!!" << endl;
890  for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
891  int capid = ledDigi[TS].capid();
892  // BE CAREFUL: this is assuming peds are stored in ADCs
893  int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
894  if (adc < 0) {
895  adc = 0;
896  } // to prevent negative adcs after ped subtraction, which should really only happen
897  // if you're using the wrong peds.
898  double fC = m_coder->charge(*m_shape, adc, capid);
899  //ta = (fC - m_ped->getValue(capid));
900  ta = fC;
901  //cout << "DetID: " << detid << " CapID: " << capid << " ADC: " << adc << " Ped: " << m_ped->getValue(capid) << " fC: " << fC << endl;
902  _mei[TS].first->Fill(ta);
903  _mei[10].first->AddBinContent(TS + 1, ta); // This is average pulse, could probably do better (Profile?)
904  if (m_fitflag > 1) {
905  if (TS == m_startTS)
906  _mei[11].first->Reset();
907  _mei[11].first->SetBinContent(TS + 1, ta);
908  }
909 
910  // keep track of max TS and max amplitude (in fC)
911  if (ta > max_fC) {
912  max_fC = ta;
913  // maxTS = TS;
914  }
915  }
916 
917  // Now we have a sample with pedestals subtracted and in units of fC
918  // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
919  // we now want to use Phil's timing correction. This is not necessary
920  // if we are performing a Landau fit (m_fitflag = 3)
921 
922  float sum = 0.;
923  for (int i = 0; i < 10; i++)
924  sum = sum + _mei[11].first->GetBinContent(i + 1);
925  if (sum > 100) {
926  if (m_fitflag == 2 || m_fitflag == 4) {
927  float timmean = _mei[11].first->GetMean(); // let's use Phil's way instead
928  float timmeancorr = BinsizeCorr(timmean);
929  _mei[12].first->Fill(timmeancorr);
930  }
931  _mei[16].first->Fill(
932  _mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
933  if (m_fitflag == 3 || m_fitflag == 4) {
934  _mei[11].first->Fit("landau", "Q");
935  TF1* fit = _mei[11].first->GetFunction("landau");
936  _mei[13].first->Fill(fit->GetParameter(1));
937  _mei[14].first->Fill(fit->GetParError(1));
938  _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
939  }
940  }
941 }
942 
943 //-----------------------------------------------------------------------------
944 float HcalLedAnalysis::BinsizeCorr(float time) {
945  // this is the bin size correction to be applied for laser data (from Andy),
946  // it comes from a pulse shape measured from TB04 data (from Jordan)
947  // This should eventually be replaced with the more thorough treatment from Phil
948 
949  float corrtime = 0.;
950  static const float tstrue[32] = {0.003, 0.03425, 0.06548, 0.09675, 0.128, 0.15925, 0.1905, 0.22175,
951  0.253, 0.28425, 0.3155, 0.34675, 0.378, 0.40925, 0.4405, 0.47175,
952  0.503, 0.53425, 0.5655, 0.59675, 0.628, 0.65925, 0.6905, 0.72175,
953  0.753, 0.78425, 0.8155, 0.84675, 0.878, 0.90925, 0.9405, 0.97175};
954  static const float tsreco[32] = {-0.00422, 0.01815, 0.04409, 0.07346, 0.09799, 0.12192, 0.15072, 0.18158,
955  0.21397, 0.24865, 0.28448, 0.31973, 0.35449, 0.39208, 0.43282, 0.47244,
956  0.5105, 0.55008, 0.58827, 0.62828, 0.6717, 0.70966, 0.74086, 0.77496,
957  0.80843, 0.83472, 0.86044, 0.8843, 0.90674, 0.92982, 0.95072, 0.9726};
958 
959  int inttime = (int)time;
960  float restime = time - inttime;
961  for (int i = 0; i <= 32; i++) {
962  float lolim = 0.;
963  float uplim = 1.;
964  float tsdown;
965  float tsup;
966  if (i > 0) {
967  lolim = tsreco[i - 1];
968  tsdown = tstrue[i - 1];
969  } else
970  tsdown = tstrue[31] - 1.;
971  if (i < 32) {
972  uplim = tsreco[i];
973  tsup = tstrue[i];
974  } else
975  tsup = tstrue[0] + 1.;
976  if (restime >= lolim && restime < uplim) {
977  corrtime = (tsdown * (uplim - restime) + tsup * (restime - lolim)) / (uplim - lolim);
978  }
979  }
980  corrtime += inttime;
981 
982  return corrtime;
983 }
984 //-----------------------------------------------------------------------------
985 
986 // Will try to implement Phil's time slew correction here at some point
987 
988 //-----------------------------------------------------------------------------
989 void HcalLedAnalysis::ProcessCalibEvent(int fiberChan, HcalCalibDetId calibId, const HcalCalibDataFrame& digi) {
990  _meca = calibHists.find(calibId);
991  if (_meca == calibHists.end()) {
992  // if histos for this channel do not exist, first create them
993  char name[1024];
994  std::string prefix;
995  if (calibId.calibFlavor() == HcalCalibDetId::CalibrationBox) {
996  std::string sector = (calibId.hcalSubdet() == HcalBarrel) ? ("HB")
997  : (calibId.hcalSubdet() == HcalEndcap) ? ("HE")
998  : (calibId.hcalSubdet() == HcalOuter) ? ("HO")
999  : (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 }
HFDataFrame::id
constexpr HcalDetId const & id() const
Definition: HFDataFrame.h:23
HcalLedAnalysis::LedHBHEHists
void LedHBHEHists(const HcalDetId &detid, const HBHEDataFrame &ledDigi, std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT, const HcalDbService &cond)
Definition: HcalLedAnalysis.cc:724
mps_fire.i
i
Definition: mps_fire.py:429
hcallasereventfilter2012_cfi.prefix
prefix
Definition: hcallasereventfilter2012_cfi.py:10
HcalDetId::zside
constexpr int zside() const
get the z-side of the cell (1/-1)
Definition: HcalDetId.h:141
HcalCalibDetId::cboxChannelString
std::string cboxChannelString() const
get the calibration box channel as a string (if relevant)
Definition: HcalCalibDetId.cc:118
HcalDetId::ietaAbs
constexpr int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:148
HcalLedAnalysis::LedSetup
void LedSetup(const std::string &m_outputFileROOT)
Definition: HcalLedAnalysis.cc:178
HcalLedAnalysis::processLedEvent
void processLedEvent(const HBHEDigiCollection &hbhe, const HODigiCollection &ho, const HFDigiCollection &hf, const HcalCalibDigiCollection &calib, const HcalDbService &cond)
Definition: HcalLedAnalysis.cc:560
edm::SortedCollection::const_iterator
std::vector< T >::const_iterator const_iterator
Definition: SortedCollection.h:80
HcalLedAnalysis::~HcalLedAnalysis
~HcalLedAnalysis()
Destructor.
Definition: HcalLedAnalysis.cc:147
HcalCalibDataFrame::size
int size() const
total number of samples in the digi
Definition: HcalCalibDataFrame.h:26
photonIsolationHIProducer_cfi.ho
ho
Definition: photonIsolationHIProducer_cfi.py:10
HcalCalibDetId::ieta
int ieta() const
Definition: HcalCalibDetId.cc:94
edm::LogError
Log< level::Error, false > LogError
Definition: MessageLogger.h:129
type
type
Definition: SiPixelVCal_PayloadInspector.cc:39
HcalCalibDetId::hcalSubdet
HcalSubdetector hcalSubdet() const
get the HcalSubdetector (if relevant)
Definition: HcalCalibDetId.cc:90
hcalRecHitTable_cff.time
time
Definition: hcalRecHitTable_cff.py:14
HcalLedAnalysis::HcalLedAnalysis
HcalLedAnalysis(const edm::ParameterSet &ps)
Constructor.
Definition: HcalLedAnalysis.cc:10
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
ecalGpuTask_cfi.sample
sample
Definition: ecalGpuTask_cfi.py:48
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
trigObjTnPSource_cfi.bins
bins
Definition: trigObjTnPSource_cfi.py:20
photonIsolationHIProducer_cfi.hbhe
hbhe
Definition: photonIsolationHIProducer_cfi.py:8
HcalDetId::subdet
constexpr HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:138
HcalCalibDataFrame::sample
const HcalQIESample & sample(int i) const
access a sample
Definition: HcalCalibDataFrame.h:39
HcalCalibDataFrame::elecId
const HcalElectronicsId & elecId() const
Definition: HcalCalibDataFrame.h:23
HcalLedAnalysis::SetupLEDHists
void SetupLEDHists(int id, const HcalDetId detid, std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT)
Definition: HcalLedAnalysis.cc:640
nano_mu_digi_cff.float
float
Definition: nano_mu_digi_cff.py:14
HcalLedAnalysis::GetLedConst
void GetLedConst(std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT)
Definition: HcalLedAnalysis.cc:213
HcalCalibDataFrame::id
const HcalCalibDetId & id() const
Definition: HcalCalibDataFrame.h:22
HcalElectronicsId::fiberChanId
constexpr int32_t fiberChanId() const
get the fiber channel id (which of channels on a fiber)
Definition: HcalElectronicsId.h:78
hi
Definition: EPCuts.h:4
HcalLedAnalysis::ProcessCalibEvent
void ProcessCalibEvent(int fiberChan, HcalCalibDetId calibId, const HcalCalibDataFrame &digi)
Definition: HcalLedAnalysis.cc:989
HcalDetId::ieta
constexpr int ieta() const
get the cell ieta
Definition: HcalDetId.h:155
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
HcalCalibDetId::iphi
int iphi() const
get the low-edge iphi (if relevant)
Definition: HcalCalibDetId.cc:104
edm::eventsetup::heterocontainer::insert
bool insert(Storage &iStorage, ItemType *iItem, const IdTag &iIdTag)
Definition: HCMethods.h:50
HFDataFrame::size
constexpr int size() const
total number of samples in the digi
Definition: HFDataFrame.h:27
HBHEDataFrame::id
constexpr const HcalDetId & id() const
Definition: HBHEDataFrame.h:23
photonIsolationHIProducer_cfi.hf
hf
Definition: photonIsolationHIProducer_cfi.py:9
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:131
genParticles_cff.map
map
Definition: genParticles_cff.py:11
HcalLedAnalysis::BinsizeCorr
float BinsizeCorr(float time)
Definition: HcalLedAnalysis.cc:944
HcalLedAnalysis::LedHFHists
void LedHFHists(const HcalDetId &detid, const HFDataFrame &ledDigi, std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT, const HcalDbService &cond)
Definition: HcalLedAnalysis.cc:865
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:57
HcalCalibDetId::calibFlavor
CalibDetType calibFlavor() const
get the flavor of this calibration detid
Definition: HcalCalibDetId.h:77
HODataFrame::size
constexpr int size() const
total number of samples in the digi
Definition: HODataFrame.h:27
HcalQIESample::adc
constexpr int adc() const
get the ADC sample
Definition: HcalQIESample.h:43
HcalLedAnalysis::LedTrendings
void LedTrendings(std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT)
Definition: HcalLedAnalysis.cc:476
output
Definition: output.py:1
nano_mu_digi_cff.sector
sector
Definition: nano_mu_digi_cff.py:21
HcalLedAnalysis::LedHOHists
void LedHOHists(const HcalDetId &detid, const HODataFrame &ledDigi, std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT, const HcalDbService &cond)
Definition: HcalLedAnalysis.cc:795
HBHEDataFrame::size
constexpr int size() const
total number of samples in the digi
Definition: HBHEDataFrame.h:27
HcalElectronicsId
Readout chain identification for Hcal.
Definition: HcalElectronicsId.h:32
HcalDetId::iphi
constexpr int iphi() const
get the cell iphi
Definition: HcalDetId.h:157
HODataFrame::id
constexpr HcalDetId const & id() const
Definition: HODataFrame.h:23
gpuClustering::adc
uint16_t *__restrict__ uint16_t const *__restrict__ adc
Definition: gpuClusterChargeCut.h:19
DeadROC_duringRun.runNum
runNum
Definition: DeadROC_duringRun.py:21
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
HcalDetId::depth
constexpr int depth() const
get the tower depth
Definition: HcalDetId.h:164
Generated for CMSSW Reference Manual by doxygen 1.8.14