CMS 3D CMS Logo
CastorLedAnalysis.cc
Go to the documentation of this file.
1 
2 #include "CalibFormats/CastorObjects/interface/CastorCoderDb.h"
3 #include "CalibFormats/CastorObjects/interface/CastorCalibrations.h"
4 #include "CalibFormats/CastorObjects/interface/CastorDbService.h"
5 #include "CalibFormats/CastorObjects/interface/CastorDbRecord.h"
6 #include "CondFormats/CastorObjects/interface/CastorQIECoder.h"
7 #include "CondFormats/CastorObjects/interface/CastorPedestals.h"
8 #include "CondFormats/CastorObjects/interface/CastorPedestalWidths.h"
9 
10 #include "CalibCalorimetry/CastorCalib/interface/CastorLedAnalysis.h"
11 #include <TFile.h>
12 #include <cmath>
13 
14 using namespace std;
15 
16 CastorLedAnalysis::CastorLedAnalysis(const edm::ParameterSet& ps) {
17  // init
18  evt = 0;
19  sample = 0;
20  m_file = nullptr;
21  char output[100]{0};
22  // output files
23  for (int k = 0; k < 4; k++)
24  state.push_back(true); // 4 cap-ids (do we care?)
25  m_outputFileText = ps.getUntrackedParameter<string>("outputFileText", "");
26  m_outputFileX = ps.getUntrackedParameter<string>("outputFileXML", "");
27  if (!m_outputFileText.empty()) {
28  cout << "Castor LED results will be saved to " << m_outputFileText.c_str() << endl;
29  m_outFile.open(m_outputFileText.c_str());
30  }
31  m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
32  if (!m_outputFileROOT.empty()) {
33  cout << "Castor LED histograms will be saved to " << m_outputFileROOT.c_str() << endl;
34  }
35 
36  m_nevtsample = ps.getUntrackedParameter<int>("nevtsample", 9999999);
37  if (m_nevtsample < 1)
38  m_nevtsample = 9999999;
39  m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag", 0);
40  if (m_hiSaveflag < 0)
41  m_hiSaveflag = 0;
42  if (m_hiSaveflag > 0)
43  m_hiSaveflag = 1;
44  m_fitflag = ps.getUntrackedParameter<int>("analysisflag", 2);
45  if (m_fitflag < 0)
46  m_fitflag = 0;
47  if (m_fitflag > 4)
48  m_fitflag = 4;
49  m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
50  if (m_startTS < 0)
51  m_startTS = 0;
52  m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
53  m_usecalib = ps.getUntrackedParameter<bool>("usecalib", false);
54  m_logFile.open("CastorLedAnalysis.log");
55 
56  int runNum = ps.getUntrackedParameter<int>("runNumber", 999999);
57 
58  // histogram booking
59  castorHists.ALLLEDS = new TH1F("Castor All LEDs", "HF All Leds", 10, 0, 9);
60  castorHists.LEDRMS = new TH1F("Castor All LED RMS", "HF All LED RMS", 100, 0, 3);
61  castorHists.LEDMEAN = new TH1F("Castor All LED Means", "HF All LED Means", 100, 0, 9);
62  castorHists.CHI2 = new TH1F("Castor Chi2 by ndf for Landau fit", "HF Chi2/ndf Landau", 200, 0., 50.);
63 
64  //XML file header
65  m_outputFileXML.open(m_outputFileX.c_str());
66  snprintf(output, sizeof output, "<?xml version='1.0' encoding='UTF-8'?>");
67  m_outputFileXML << output << endl;
68  snprintf(output, sizeof output, "<ROOT>");
69  m_outputFileXML << output << endl << endl;
70  snprintf(output, sizeof output, " <HEADER>");
71  m_outputFileXML << output << endl;
72  snprintf(output, sizeof output, " <TYPE>");
73  m_outputFileXML << output << endl;
74  snprintf(output, sizeof output, " <EXTENSION_TABLE_NAME>HCAL_LED_TIMING</EXTENSION_TABLE_NAME>");
75  m_outputFileXML << output << endl;
76  snprintf(output, sizeof output, " <NAME>HCAL LED Timing</NAME>");
77  m_outputFileXML << output << endl;
78  snprintf(output, sizeof output, " </TYPE>");
79  m_outputFileXML << output << endl;
80  snprintf(output, sizeof output, " <RUN>");
81  m_outputFileXML << output << endl;
82  snprintf(output, sizeof output, " <RUN_TYPE>hcal-led-timing-test</RUN_TYPE>");
83  m_outputFileXML << output << endl;
84  snprintf(output, sizeof output, " <RUN_NUMBER>%06i</RUN_NUMBER>", runNum);
85  m_outputFileXML << output << endl;
86  snprintf(output, sizeof output, " <RUN_BEGIN_TIMESTAMP>2007-07-09 00:00:00.0</RUN_BEGIN_TIMESTAMP>");
87  m_outputFileXML << output << endl;
88  snprintf(output, sizeof output, " <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>");
89  m_outputFileXML << output << endl;
90  snprintf(output, sizeof output, " </RUN>");
91  m_outputFileXML << output << endl;
92  snprintf(output, sizeof output, " </HEADER>");
93  m_outputFileXML << output << endl;
94  snprintf(output, sizeof output, "<!-- Tags secton -->");
95  m_outputFileXML << output << endl;
96  snprintf(output, sizeof output, " <ELEMENTS>");
97  m_outputFileXML << output << endl;
98  snprintf(output, sizeof output, " <DATA_SET id='-1'/>");
99  m_outputFileXML << output << endl;
100  snprintf(output, sizeof output, " <IOV id='1'>");
101  m_outputFileXML << output << endl;
102  snprintf(output, sizeof output, " <INTERVAL_OF_VALIDITY_BEGIN>2147483647</INTERVAL_OF_VALIDITY_BEGIN>");
103  m_outputFileXML << output << endl;
104  snprintf(output, sizeof output, " <INTERVAL_OF_VALIDITY_END>0</INTERVAL_OF_VALIDITY_END>");
105  m_outputFileXML << output << endl;
106  snprintf(output, sizeof output, " </IOV>");
107  m_outputFileXML << output << endl;
108  snprintf(output, sizeof output, " <TAG id='2' mode='auto'>");
109  m_outputFileXML << output << endl;
110  snprintf(output, sizeof output, " <TAG_NAME>laser_led_%06i<TAG_NAME>", runNum);
111  m_outputFileXML << output << endl;
112  snprintf(output, sizeof output, " <DETECTOR_NAME>HCAL</DETECTOR_NAME>");
113  m_outputFileXML << output << endl;
114  snprintf(output, sizeof output, " <COMMENT_DESCRIPTION></COMMENT_DESCRIPTION>");
115  m_outputFileXML << output << endl;
116  snprintf(output, sizeof output, " </TAG>");
117  m_outputFileXML << output << endl;
118  snprintf(output, sizeof output, " </ELEMENTS>");
119  m_outputFileXML << output << endl;
120  snprintf(output, sizeof output, " <MAPS>");
121  m_outputFileXML << output << endl;
122  snprintf(output, sizeof output, " <TAG idref ='2'>");
123  m_outputFileXML << output << endl;
124  snprintf(output, sizeof output, " <IOV idref='1'>");
125  m_outputFileXML << output << endl;
126  snprintf(output, sizeof output, " <DATA_SET idref='-1' />");
127  m_outputFileXML << output << endl;
128  snprintf(output, sizeof output, " </IOV>");
129  m_outputFileXML << output << endl;
130  snprintf(output, sizeof output, " </TAG>");
131  m_outputFileXML << output << endl;
132  snprintf(output, sizeof output, " </MAPS>");
133  m_outputFileXML << output << endl;
134 }
135 
136 //-----------------------------------------------------------------------------
137 CastorLedAnalysis::~CastorLedAnalysis() {
139 
140  for (_meol = castorHists.LEDTRENDS.begin(); _meol != castorHists.LEDTRENDS.end(); _meol++) {
141  for (int i = 0; i < 15; i++)
142  _meol->second[i].first->Delete();
143  }
144 
145  castorHists.ALLLEDS->Delete();
146  castorHists.LEDRMS->Delete();
147  castorHists.LEDMEAN->Delete();
148  castorHists.CHI2->Delete();
149 }
150 
151 //-----------------------------------------------------------------------------
152 void CastorLedAnalysis::LedSetup(const std::string& m_outputFileROOT) {
153  // open the histogram file, create directories within
154  m_file = new TFile(m_outputFileROOT.c_str(), "RECREATE");
155  m_file->mkdir("Castor");
156  m_file->cd();
157 }
158 
159 //-----------------------------------------------------------------------------
160 void CastorLedAnalysis::GetLedConst(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
161  double time2 = 0;
162  double time1 = 0;
163  double time3 = 0;
164  double time4 = 0;
165  double dtime2 = 0;
166  double dtime1 = 0;
167  double dtime3 = 0;
168  double dtime4 = 0;
169  char output[100]{0};
170 
171  if (!m_outputFileText.empty()) {
172  if (m_fitflag == 0 || m_fitflag == 2)
173  m_outFile << "Det Eta,Phi,D Mean Error" << std::endl;
174  else if (m_fitflag == 1 || m_fitflag == 3)
175  m_outFile << "Det Eta,Phi,D Peak Error" << std::endl;
176  else if (m_fitflag == 4)
177  m_outFile << "Det Eta,Phi,D Mean Error Peak Error MeanEv Error PeakEv Error"
178  << std::endl;
179  }
180  for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
181  // scale the LED pulse to 1 event
182  _meol->second[10].first->Scale(1. / evt_curr);
183  if (m_fitflag == 0 || m_fitflag == 4) {
184  time1 = _meol->second[10].first->GetMean();
185  dtime1 = _meol->second[10].first->GetRMS() / sqrt((float)evt_curr * (m_endTS - m_startTS + 1));
186  }
187  if (m_fitflag == 1 || m_fitflag == 4) {
188  // put proper errors
189  for (int j = 0; j < 10; j++)
190  _meol->second[10].first->SetBinError(j + 1, _meol->second[j].first->GetRMS() / sqrt((float)evt_curr));
191  }
192  if (m_fitflag == 1 || m_fitflag == 3 || m_fitflag == 4) {
193  _meol->second[10].first->Fit("landau", "Q");
194  // _meol->second[10].first->Fit("gaus","Q");
195  TF1* fit = _meol->second[10].first->GetFunction("landau");
196  // TF1 *fit = _meol->second[10].first->GetFunction("gaus");
197  time2 = fit->GetParameter(1);
198  dtime2 = fit->GetParError(1);
199  }
200  if (m_fitflag == 2 || m_fitflag == 4) {
201  time3 = _meol->second[12].first->GetMean();
202  dtime3 = _meol->second[12].first->GetRMS() / sqrt((float)_meol->second[12].first->GetEntries());
203  }
204  if (m_fitflag == 3 || m_fitflag == 4) {
205  time4 = _meol->second[13].first->GetMean();
206  dtime4 = _meol->second[13].first->GetRMS() / sqrt((float)_meol->second[13].first->GetEntries());
207  }
208  for (int i = 0; i < 10; i++) {
209  _meol->second[i].first->GetXaxis()->SetTitle("Pulse height (fC)");
210  _meol->second[i].first->GetYaxis()->SetTitle("Counts");
211  // if(m_hiSaveflag>0)_meol->second[i].first->Write();
212  }
213  _meol->second[10].first->GetXaxis()->SetTitle("Time slice");
214  _meol->second[10].first->GetYaxis()->SetTitle("Averaged pulse (fC)");
215  if (m_hiSaveflag > 0)
216  _meol->second[10].first->Write();
217  _meol->second[10].second.first[0].push_back(time1);
218  _meol->second[10].second.first[1].push_back(dtime1);
219  _meol->second[11].second.first[0].push_back(time2);
220  _meol->second[11].second.first[1].push_back(dtime2);
221  _meol->second[12].first->GetXaxis()->SetTitle("Mean TS");
222  _meol->second[12].first->GetYaxis()->SetTitle("Counts");
223  if (m_fitflag == 2 && m_hiSaveflag > 0)
224  _meol->second[12].first->Write();
225  _meol->second[12].second.first[0].push_back(time3);
226  _meol->second[12].second.first[1].push_back(dtime3);
227  _meol->second[13].first->GetXaxis()->SetTitle("Peak TS");
228  _meol->second[13].first->GetYaxis()->SetTitle("Counts");
229  if (m_fitflag > 2 && m_hiSaveflag > 0)
230  _meol->second[13].first->Write();
231  _meol->second[13].second.first[0].push_back(time4);
232  _meol->second[13].second.first[1].push_back(dtime4);
233  _meol->second[14].first->GetXaxis()->SetTitle("Peak TS error");
234  _meol->second[14].first->GetYaxis()->SetTitle("Counts");
235  if (m_fitflag > 2 && m_hiSaveflag > 0)
236  _meol->second[14].first->Write();
237  _meol->second[15].first->GetXaxis()->SetTitle("Chi2/NDF");
238  _meol->second[15].first->GetYaxis()->SetTitle("Counts");
239  if (m_fitflag > 2 && m_hiSaveflag > 0)
240  _meol->second[15].first->Write();
241  _meol->second[16].first->GetXaxis()->SetTitle("Integrated Signal");
242  _meol->second[16].first->Write();
243 
244  // Ascii printout (need to modify to include new info)
245  HcalDetId detid = _meol->first;
246 
247  if (!m_outputFileText.empty()) {
248  if (m_fitflag == 0) {
249  m_outFile << detid << " " << time1 << " " << dtime1 << std::endl;
250  snprintf(output, sizeof output, " <DATA_SET>");
251  m_outputFileXML << output << endl;
252  snprintf(output, sizeof output, " <VERSION>version:1</VERSION>");
253  m_outputFileXML << output << endl;
254  snprintf(output, sizeof output, " <CHANNEL>");
255  m_outputFileXML << output << endl;
256  snprintf(output, sizeof output, " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
257  m_outputFileXML << output << endl;
258  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
259  m_outputFileXML << output << endl;
260  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
261  m_outputFileXML << output << endl;
262  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
263  m_outputFileXML << output << endl;
264  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
265  m_outputFileXML << output << endl;
266  if (detid.subdet() == 1)
267  snprintf(output, sizeof output, " <DETECTOR_NAME>HB</DETECTOR_NAME>");
268  if (detid.subdet() == 2)
269  snprintf(output, sizeof output, " <DETECTOR_NAME>HE</DETECTOR_NAME>");
270  if (detid.subdet() == 3)
271  snprintf(output, sizeof output, " <DETECTOR_NAME>HO</DETECTOR_NAME>");
272  if (detid.subdet() == 4)
273  snprintf(output, sizeof output, " <DETECTOR_NAME>HF</DETECTOR_NAME>");
274  m_outputFileXML << output << endl;
275  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
276  m_outputFileXML << output << endl;
277  snprintf(output, sizeof output, " </CHANNEL>");
278  m_outputFileXML << output << endl;
279  snprintf(output, sizeof output, " <DATA>");
280  m_outputFileXML << output << endl;
281  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time1);
282  m_outputFileXML << output << endl;
283  snprintf(output, sizeof output, " <OFFSET_TIME> 0</OFFSET_TIME>");
284  m_outputFileXML << output << endl;
285  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime1);
286  m_outputFileXML << output << endl;
287  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
288  m_outputFileXML << output << endl;
289  snprintf(output, sizeof output, " <STATUS_WORD> 0</STATUS_WORD>");
290  m_outputFileXML << output << endl;
291  snprintf(output, sizeof output, " </DATA>");
292  m_outputFileXML << output << endl;
293  snprintf(output, sizeof output, " </DATA_SET>");
294  m_outputFileXML << output << endl;
295 
296  } else if (m_fitflag == 1) {
297  m_outFile << detid << " " << time2 << " " << dtime2 << std::endl;
298  snprintf(output, sizeof output, " <DATA_SET>");
299  m_outputFileXML << output << endl;
300  snprintf(output, sizeof output, " <VERSION>version:1</VERSION>");
301  m_outputFileXML << output << endl;
302  snprintf(output, sizeof output, " <CHANNEL>");
303  m_outputFileXML << output << endl;
304  snprintf(output, sizeof output, " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
305  m_outputFileXML << output << endl;
306  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
307  m_outputFileXML << output << endl;
308  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
309  m_outputFileXML << output << endl;
310  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
311  m_outputFileXML << output << endl;
312  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
313  m_outputFileXML << output << endl;
314  if (detid.subdet() == 1)
315  snprintf(output, sizeof output, " <DETECTOR_NAME>HB</DETECTOR_NAME>");
316  if (detid.subdet() == 2)
317  snprintf(output, sizeof output, " <DETECTOR_NAME>HE</DETECTOR_NAME>");
318  if (detid.subdet() == 3)
319  snprintf(output, sizeof output, " <DETECTOR_NAME>HO</DETECTOR_NAME>");
320  if (detid.subdet() == 4)
321  snprintf(output, sizeof output, " <DETECTOR_NAME>HF</DETECTOR_NAME>");
322  m_outputFileXML << output << endl;
323  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
324  m_outputFileXML << output << endl;
325  snprintf(output, sizeof output, " </CHANNEL>");
326  m_outputFileXML << output << endl;
327  snprintf(output, sizeof output, " <DATA>");
328  m_outputFileXML << output << endl;
329  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time2);
330  m_outputFileXML << output << endl;
331  snprintf(output, sizeof output, " <OFFSET_TIME> 0</OFFSET_TIME>");
332  m_outputFileXML << output << endl;
333  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime2);
334  m_outputFileXML << output << endl;
335  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
336  m_outputFileXML << output << endl;
337  snprintf(output, sizeof output, " <STATUS_WORD> 0</STATUS_WORD>");
338  m_outputFileXML << output << endl;
339  snprintf(output, sizeof output, " </DATA>");
340  m_outputFileXML << output << endl;
341  snprintf(output, sizeof output, " </DATA_SET>");
342  m_outputFileXML << output << endl;
343  }
344 
345  else if (m_fitflag == 2) {
346  m_outFile << detid << " " << time3 << " " << dtime3 << std::endl;
347  snprintf(output, sizeof output, " <DATA_SET>");
348  m_outputFileXML << output << endl;
349  snprintf(output, sizeof output, " <VERSION>version:1</VERSION>");
350  m_outputFileXML << output << endl;
351  snprintf(output, sizeof output, " <CHANNEL>");
352  m_outputFileXML << output << endl;
353  snprintf(output, sizeof output, " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
354  m_outputFileXML << output << endl;
355  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
356  m_outputFileXML << output << endl;
357  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
358  m_outputFileXML << output << endl;
359  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
360  m_outputFileXML << output << endl;
361  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
362  m_outputFileXML << output << endl;
363  if (detid.subdet() == 1)
364  snprintf(output, sizeof output, " <DETECTOR_NAME>HB</DETECTOR_NAME>");
365  if (detid.subdet() == 2)
366  snprintf(output, sizeof output, " <DETECTOR_NAME>HE</DETECTOR_NAME>");
367  if (detid.subdet() == 3)
368  snprintf(output, sizeof output, " <DETECTOR_NAME>HO</DETECTOR_NAME>");
369  if (detid.subdet() == 4)
370  snprintf(output, sizeof output, " <DETECTOR_NAME>HF</DETECTOR_NAME>");
371  m_outputFileXML << output << endl;
372  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
373  m_outputFileXML << output << endl;
374  snprintf(output, sizeof output, " </CHANNEL>");
375  m_outputFileXML << output << endl;
376  snprintf(output, sizeof output, " <DATA>");
377  m_outputFileXML << output << endl;
378  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time3);
379  m_outputFileXML << output << endl;
380  snprintf(output, sizeof output, " <OFFSET_TIME> 0</OFFSET_TIME>");
381  m_outputFileXML << output << endl;
382  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime3);
383  m_outputFileXML << output << endl;
384  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
385  m_outputFileXML << output << endl;
386  snprintf(output, sizeof output, " <STATUS_WORD> 0</STATUS_WORD>");
387  m_outputFileXML << output << endl;
388  snprintf(output, sizeof output, " </DATA>");
389  m_outputFileXML << output << endl;
390  snprintf(output, sizeof output, " </DATA_SET>");
391  m_outputFileXML << output << endl;
392  } else if (m_fitflag == 3) {
393  m_outFile << detid << " " << time4 << " " << dtime4 << std::endl;
394  snprintf(output, sizeof output, " <DATA_SET>");
395  m_outputFileXML << output << endl;
396  snprintf(output, sizeof output, " <VERSION>version:1</VERSION>");
397  m_outputFileXML << output << endl;
398  snprintf(output, sizeof output, " <CHANNEL>");
399  m_outputFileXML << output << endl;
400  snprintf(output, sizeof output, " <EXTENSION_TABLE_NAME>HCAL_CHANNELS</EXTENSION_TABLE_NAME>");
401  m_outputFileXML << output << endl;
402  snprintf(output, sizeof output, " <ETA>%2i</ETA>", detid.ietaAbs());
403  m_outputFileXML << output << endl;
404  snprintf(output, sizeof output, " <PHI>%2i</PHI>", detid.iphi());
405  m_outputFileXML << output << endl;
406  snprintf(output, sizeof output, " <DEPTH>%2i</DEPTH>", detid.depth());
407  m_outputFileXML << output << endl;
408  snprintf(output, sizeof output, " <Z>%2i</Z>", detid.zside());
409  m_outputFileXML << output << endl;
410  if (detid.subdet() == 1)
411  snprintf(output, sizeof output, " <DETECTOR_NAME>HB</DETECTOR_NAME>");
412  if (detid.subdet() == 2)
413  snprintf(output, sizeof output, " <DETECTOR_NAME>HE</DETECTOR_NAME>");
414  if (detid.subdet() == 3)
415  snprintf(output, sizeof output, " <DETECTOR_NAME>HO</DETECTOR_NAME>");
416  if (detid.subdet() == 4)
417  snprintf(output, sizeof output, " <DETECTOR_NAME>HF</DETECTOR_NAME>");
418  m_outputFileXML << output << endl;
419  snprintf(output, sizeof output, " <HCAL_CHANNEL_ID>%10i</HCAL_CHANNEL_ID>", detid.rawId());
420  m_outputFileXML << output << endl;
421  snprintf(output, sizeof output, " </CHANNEL>");
422  m_outputFileXML << output << endl;
423  snprintf(output, sizeof output, " <DATA>");
424  m_outputFileXML << output << endl;
425  snprintf(output, sizeof output, " <MEAN_TIME>%7f</MEAN_TIME>", time4);
426  m_outputFileXML << output << endl;
427  snprintf(output, sizeof output, " <OFFSET_TIME> 0</OFFSET_TIME>");
428  m_outputFileXML << output << endl;
429  snprintf(output, sizeof output, " <ERROR_STAT>%7f</ERROR_STAT>", dtime4);
430  m_outputFileXML << output << endl;
431  snprintf(output, sizeof output, " <ANALYSIS_FLAG>%2i</ANALYSIS_FLAG>", m_fitflag + 1);
432  m_outputFileXML << output << endl;
433  snprintf(output, sizeof output, " <STATUS_WORD> 0</STATUS_WORD>");
434  m_outputFileXML << output << endl;
435  snprintf(output, sizeof output, " </DATA>");
436  m_outputFileXML << output << endl;
437  snprintf(output, sizeof output, " </DATA_SET>");
438  m_outputFileXML << output << endl;
439  }
440 
441  else if (m_fitflag == 4) {
442  m_outFile << detid << " " << time1 << " " << dtime1 << " " << time2 << " " << dtime2 << " " << time3
443  << " " << dtime3 << " " << time4 << " " << dtime4 << std::endl;
444  }
445  }
446  }
447 }
448 
449 //-----------------------------------------------------------------------------
450 void CastorLedAnalysis::LedSampleAnalysis() {
451  // it is called every m_nevtsample events (a sample) and the end of run
452  char LedSampleNum[20];
453 
454  sprintf(LedSampleNum, "LedSample_%d", sample);
455  m_file->cd();
456  m_file->mkdir(LedSampleNum);
457  m_file->cd(LedSampleNum);
458 
459  // Compute LED constants
460  GetLedConst(castorHists.LEDTRENDS);
461 }
462 
463 //-----------------------------------------------------------------------------
464 void CastorLedAnalysis::LedTrendings(map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
465  for (_meol = toolT.begin(); _meol != toolT.end(); _meol++) {
466  char name[1024];
467  HcalDetId detid = _meol->first;
468  sprintf(name, "LED timing trend, eta=%d phi=%d depth=%d", detid.ieta(), detid.iphi(), detid.depth());
469  int bins = _meol->second[10 + m_fitflag].second.first[0].size();
470  float lo = 0.5;
471  float hi = (float)bins + 0.5;
472  _meol->second[10 + m_fitflag].second.second.push_back(new TH1F(name, name, bins, lo, hi));
473 
474  std::vector<double>::iterator sample_it;
475  // LED timing - put content and errors
476  int j = 0;
477  for (sample_it = _meol->second[10 + m_fitflag].second.first[0].begin();
478  sample_it != _meol->second[10 + m_fitflag].second.first[0].end();
479  ++sample_it) {
480  _meol->second[10 + m_fitflag].second.second[0]->SetBinContent(++j, *sample_it);
481  }
482  j = 0;
483  for (sample_it = _meol->second[10 + m_fitflag].second.first[1].begin();
484  sample_it != _meol->second[10 + m_fitflag].second.first[1].end();
485  ++sample_it) {
486  _meol->second[10 + m_fitflag].second.second[0]->SetBinError(++j, *sample_it);
487  }
488  sprintf(name, "Sample (%d events)", m_nevtsample);
489  _meol->second[10 + m_fitflag].second.second[0]->GetXaxis()->SetTitle(name);
490  _meol->second[10 + m_fitflag].second.second[0]->GetYaxis()->SetTitle("Peak position");
491  _meol->second[10 + m_fitflag].second.second[0]->Write();
492  }
493 }
494 
495 //-----------------------------------------------------------------------------
496 void CastorLedAnalysis::LedDone() {
497  // First process the last sample (remaining events).
498  if (evt % m_nevtsample != 0)
499  LedSampleAnalysis();
500 
501  // Now do the end of run analysis: trending histos
502  if (sample > 1 && m_fitflag != 4) {
503  m_file->cd();
504  m_file->cd("Castor");
505  LedTrendings(castorHists.LEDTRENDS);
506  }
507 
508  // Write other histograms.
509  m_file->cd();
510  m_file->cd("Castor");
511  castorHists.ALLLEDS->Write();
512  castorHists.LEDRMS->Write();
513  castorHists.LEDMEAN->Write();
514 
515  // Write the histo file and close it
516  // m_file->Write();
517  m_file->Close();
518  cout << "Castor histograms written to " << m_outputFileROOT.c_str() << endl;
519 }
520 
521 //-----------------------------------------------------------------------------
522 void CastorLedAnalysis::processLedEvent(const CastorDigiCollection& castor, const CastorDbService& cond) {
523  evt++;
524  sample = (evt - 1) / m_nevtsample + 1;
525  evt_curr = evt % m_nevtsample;
526  if (evt_curr == 0)
527  evt_curr = m_nevtsample;
528 
529  // HF/Castor
530  try {
531  if (castor.empty())
532  throw(int) castor.size();
533  for (CastorDigiCollection::const_iterator j = castor.begin(); j != castor.end(); ++j) {
534  const CastorDataFrame digi = (const CastorDataFrame)(*j);
535  _meol = castorHists.LEDTRENDS.find(digi.id());
536  if (_meol == castorHists.LEDTRENDS.end()) {
537  SetupLEDHists(2, digi.id(), castorHists.LEDTRENDS);
538  }
539  LedCastorHists(digi.id(), digi, castorHists.LEDTRENDS, cond);
540  }
541  } catch (int i) {
542  // m_logFile << "Event with " << i<<" Castor Digis passed." << std::endl;
543  }
544 
545  // Call the function every m_nevtsample events
546  if (evt % m_nevtsample == 0)
547  LedSampleAnalysis();
548 }
549 //----------------------------------------------------------------------------
550 void CastorLedAnalysis::SetupLEDHists(int id, const HcalDetId detid, map<HcalDetId, map<int, LEDBUNCH> >& toolT) {
551  string type = "HBHE";
552  if (id == 1)
553  type = "HO";
554  if (id == 2)
555  type = "HF";
556 
557  _meol = toolT.find(detid);
558  if (_meol == toolT.end()) {
559  // if histos for this channel do not exist, create them
560  map<int, LEDBUNCH> insert;
561  char name[1024];
562  for (int i = 0; i < 10; i++) {
563  sprintf(name,
564  "%s Pulse height, eta=%d phi=%d depth=%d TS=%d",
565  type.c_str(),
566  detid.ieta(),
567  detid.iphi(),
568  detid.depth(),
569  i);
570  insert[i].first = new TH1F(name, name, 200, 0., 2000.);
571  }
572  sprintf(name, "%s LED Mean pulse, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
573  insert[10].first = new TH1F(name, name, 10, -0.5, 9.5);
574  sprintf(name, "%s LED Pulse, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
575  insert[11].first = new TH1F(name, name, 10, -0.5, 9.5);
576  sprintf(name, "%s Mean TS, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
577  insert[12].first = new TH1F(name, name, 200, 0., 10.);
578  sprintf(name, "%s Peak TS, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
579  insert[13].first = new TH1F(name, name, 200, 0., 10.);
580  sprintf(name, "%s Peak TS error, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
581  insert[14].first = new TH1F(name, name, 200, 0., 0.05);
582  sprintf(name, "%s Fit chi2, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
583  insert[15].first = new TH1F(name, name, 100, 0., 50.);
584  sprintf(
585  name, "%s Integrated Signal, eta=%d phi=%d depth=%d", type.c_str(), detid.ieta(), detid.iphi(), detid.depth());
586  insert[16].first = new TH1F(name, name, 500, 0., 5000.);
587 
588  toolT[detid] = insert;
589  }
590 }
591 //-----------------------------------------------------------------------------
592 void CastorLedAnalysis::LedCastorHists(const HcalDetId& detid,
593  const CastorDataFrame& ledDigi,
594  map<HcalDetId, map<int, LEDBUNCH> >& toolT,
595  const CastorDbService& cond) {
596  map<int, LEDBUNCH> _mei;
597  _meol = toolT.find(detid);
598  _mei = _meol->second;
599  // Rest the histos if we're at the end of a 'bunch'
600  if ((evt - 1) % m_nevtsample == 0 && state[0]) {
601  for (int k = 0; k < (int)state.size(); k++)
602  state[k] = false;
603  for (int i = 0; i < 16; i++)
604  _mei[i].first->Reset();
605  }
606 
607  // now we have the signal in fC, let's get rid of that darn pedestal
608  // Most of this is borrowed from CastorSimpleReconstructor, so thanks Jeremy/Phil
609 
610  float max_fC = 0;
611  float ta = 0;
612  m_coder = cond.getCastorCoder(detid);
613  m_ped = cond.getPedestal(detid);
614  m_shape = cond.getCastorShape();
615  //cout << "New Digi!!!!!!!!!!!!!!!!!!!!!!" << endl;
616  for (int TS = m_startTS; TS < m_endTS && TS < ledDigi.size(); TS++) {
617  int capid = ledDigi[TS].capid();
618  // BE CAREFUL: this is assuming peds are stored in ADCs
619  int adc = (int)(ledDigi[TS].adc() - m_ped->getValue(capid));
620  if (adc < 0) {
621  adc = 0;
622  } // to prevent negative adcs after ped subtraction, which should really only happen
623  // if you're using the wrong peds.
624  double fC = m_coder->charge(*m_shape, adc, capid);
625  //ta = (fC - m_ped->getValue(capid));
626  ta = fC;
627  //cout << "DetID: " << detid << " CapID: " << capid << " ADC: " << adc << " Ped: " << m_ped->getValue(capid) << " fC: " << fC << endl;
628  _mei[TS].first->Fill(ta);
629  _mei[10].first->AddBinContent(TS + 1, ta); // This is average pulse, could probably do better (Profile?)
630  if (m_fitflag > 1) {
631  if (TS == m_startTS)
632  _mei[11].first->Reset();
633  _mei[11].first->SetBinContent(TS + 1, ta);
634  }
635 
636  // keep track of max TS and max amplitude (in fC)
637  if (ta > max_fC) {
638  max_fC = ta;
639  }
640  }
641 
642  // Now we have a sample with pedestals subtracted and in units of fC
643  // If we are using a weighted mean (m_fitflag = 2) to extraxt timing
644  // we now want to use Phil's timing correction. This is not necessary
645  // if we are performing a Landau fit (m_fitflag = 3)
646 
647  float sum = 0.;
648  for (int i = 0; i < 10; i++)
649  sum = sum + _mei[11].first->GetBinContent(i + 1);
650  if (sum > 100) {
651  if (m_fitflag == 2 || m_fitflag == 4) {
652  float timmean = _mei[11].first->GetMean(); // let's use Phil's way instead
653  float timmeancorr = BinsizeCorr(timmean);
654  _mei[12].first->Fill(timmeancorr);
655  }
656  _mei[16].first->Fill(
657  _mei[11].first->Integral()); // Integrated charge (may be more usfull to convert to Energy first?)
658  if (m_fitflag == 3 || m_fitflag == 4) {
659  _mei[11].first->Fit("landau", "Q");
660  TF1* fit = _mei[11].first->GetFunction("landau");
661  _mei[13].first->Fill(fit->GetParameter(1));
662  _mei[14].first->Fill(fit->GetParError(1));
663  _mei[15].first->Fill(fit->GetChisquare() / fit->GetNDF());
664  }
665  }
666 }
667 
668 //-----------------------------------------------------------------------------
669 float CastorLedAnalysis::BinsizeCorr(float time) {
670  // this is the bin size correction to be applied for laser data (from Andy),
671  // it comes from a pulse shape measured from TB04 data (from Jordan)
672  // This should eventually be replaced with the more thorough treatment from Phil
673 
674  float corrtime = 0.;
675  static const float tstrue[32] = {0.003, 0.03425, 0.06548, 0.09675, 0.128, 0.15925, 0.1905, 0.22175,
676  0.253, 0.28425, 0.3155, 0.34675, 0.378, 0.40925, 0.4405, 0.47175,
677  0.503, 0.53425, 0.5655, 0.59675, 0.628, 0.65925, 0.6905, 0.72175,
678  0.753, 0.78425, 0.8155, 0.84675, 0.878, 0.90925, 0.9405, 0.97175};
679  static const float tsreco[32] = {-0.00422, 0.01815, 0.04409, 0.07346, 0.09799, 0.12192, 0.15072, 0.18158,
680  0.21397, 0.24865, 0.28448, 0.31973, 0.35449, 0.39208, 0.43282, 0.47244,
681  0.5105, 0.55008, 0.58827, 0.62828, 0.6717, 0.70966, 0.74086, 0.77496,
682  0.80843, 0.83472, 0.86044, 0.8843, 0.90674, 0.92982, 0.95072, 0.9726};
683 
684  int inttime = (int)time;
685  float restime = time - inttime;
686  for (int i = 0; i <= 32; i++) {
687  float lolim = 0.;
688  float uplim = 1.;
689  float tsdown;
690  float tsup;
691  if (i > 0) {
692  lolim = tsreco[i - 1];
693  tsdown = tstrue[i - 1];
694  } else
695  tsdown = tstrue[31] - 1.;
696  if (i < 32) {
697  uplim = tsreco[i];
698  tsup = tstrue[i];
699  } else
700  tsup = tstrue[0] + 1.;
701  if (restime >= lolim && restime < uplim) {
702  corrtime = (tsdown * (uplim - restime) + tsup * (restime - lolim)) / (uplim - lolim);
703  }
704  }
705  corrtime += inttime;
706 
707  return corrtime;
708 }
709 //-----------------------------------------------------------------------------
type
type
Definition: HCALResponse.h:21
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
CastorLedAnalysis::~CastorLedAnalysis
~CastorLedAnalysis()
Destructor.
Definition: CastorLedAnalysis.cc:137
mps_fire.i
i
Definition: mps_fire.py:338
HcalDetId::subdet
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:146
CastorDbService::getPedestal
const CastorPedestal * getPedestal(const HcalGenericDetId &fId) const
Definition: CastorDbService.cc:141
simplePhotonAnalyzer_cfi.sample
sample
Definition: simplePhotonAnalyzer_cfi.py:12
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
HcalDetId::zside
int zside() const
get the z-side of the cell (1/-1)
Definition: HcalDetId.h:149
CastorLedAnalysis::LedSetup
void LedSetup(const std::string &m_outputFileROOT)
Definition: CastorLedAnalysis.cc:152
DetId::rawId
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
edm::SortedCollection::const_iterator
std::vector< T >::const_iterator const_iterator
Definition: SortedCollection.h:82
convertSQLitetoXML_cfg.output
output
Definition: convertSQLitetoXML_cfg.py:32
CastorLedAnalysis::LedTrendings
void LedTrendings(std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT)
Definition: CastorLedAnalysis.cc:464
CastorLedAnalysis::GetLedConst
void GetLedConst(std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT)
Definition: CastorLedAnalysis.cc:160
CastorLedAnalysis::CastorLedAnalysis
CastorLedAnalysis(const edm::ParameterSet &ps)
Constructor.
Definition: CastorLedAnalysis.cc:16
HcalDetId::depth
int depth() const
get the tower depth
Definition: HcalDetId.h:166
CastorLedAnalysis::LedCastorHists
void LedCastorHists(const HcalDetId &detid, const CastorDataFrame &ledDigi, std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT, const CastorDbService &cond)
Definition: CastorLedAnalysis.cc:592
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:18
CastorDbService::getCastorShape
const CastorQIEShape * getCastorShape() const
Definition: CastorDbService.cc:176
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:159
edm::eventsetup::heterocontainer::insert
bool insert(Storage &iStorage, ItemType *iItem, const IdTag &iIdTag)
Definition: HCMethods.h:50
ecalMGPA::adc
constexpr int adc(sample_type sample)
get the ADC sample (12 bits)
Definition: EcalMGPASample.h:11
CastorLedAnalysis::BinsizeCorr
float BinsizeCorr(float time)
Definition: CastorLedAnalysis.cc:669
create_public_pileup_plots.bins
bins
Definition: create_public_pileup_plots.py:231
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:79
HcalDetId::ietaAbs
int ietaAbs() const
get the absolute value of the cell ieta
Definition: HcalDetId.h:154
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:310
HcalDetId::iphi
int iphi() const
get the cell iphi
Definition: HcalDetId.h:161
genParticles_cff.map
map
Definition: genParticles_cff.py:11
protons_cff.time
time
Definition: protons_cff.py:35
castorDigitizer_cfi.castor
castor
Definition: castorDigitizer_cfi.py:9
CastorLedAnalysis::SetupLEDHists
void SetupLEDHists(int id, const HcalDetId detid, std::map< HcalDetId, std::map< int, LEDBUNCH > > &toolT)
Definition: CastorLedAnalysis.cc:550
CastorLedAnalysis::processLedEvent
void processLedEvent(const CastorDigiCollection &castor, const CastorDbService &cond)
Definition: CastorLedAnalysis.cc:522
cond
Definition: plugin.cc:24
edm::SortedCollection::size
size_type size() const
Definition: SortedCollection.h:236
gather_cfg.cout
cout
Definition: gather_cfg.py:144
CastorDataFrame::id
const HcalCastorDetId & id() const
Definition: CastorDataFrame.h:22
CastorDataFrame::size
int size() const
total number of samples in the digi
Definition: CastorDataFrame.h:26
edm::SortedCollection::begin
const_iterator begin() const
Definition: SortedCollection.h:303
CastorDbService::getCastorCoder
const CastorQIECoder * getCastorCoder(const HcalGenericDetId &fId) const
Definition: CastorDbService.cc:169
Generated for CMSSW Reference Manual by doxygen 1.8.11