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