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