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