CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
CastorPedestalAnalysis.cc
Go to the documentation of this file.
8 
10 #include <TFile.h>
11 #include <math.h>
12 
13 using namespace std;
14 
16  : fRefPedestals (0),
17  fRefPedestalWidths (0),
18  fRawPedestals (0),
19  fRawPedestalWidths (0),
20  fValPedestals (0),
21  fValPedestalWidths (0)
22 {
23  evt=0;
24  sample=0;
25  m_file=0;
26  m_AllPedsOK=0;
27  for(int i=0; i<4; i++) m_stat[i]=0;
28  for(int k=0;k<4;k++) state.push_back(true);
29 
30 // user cfg parameters
31  m_outputFileMean = ps.getUntrackedParameter<string>("outputFileMeans", "");
32  if ( m_outputFileMean.size() != 0 ) {
33  cout << "Castor pedestal means will be saved to " << m_outputFileMean.c_str() << endl;
34  }
35  m_outputFileWidth = ps.getUntrackedParameter<string>("outputFileWidths", "");
36  if ( m_outputFileWidth.size() != 0 ) {
37  cout << "Castor pedestal widths will be saved to " << m_outputFileWidth.c_str() << endl;
38  }
39  m_outputFileROOT = ps.getUntrackedParameter<string>("outputFileHist", "");
40  if ( m_outputFileROOT.size() != 0 ) {
41  cout << "Castor pedestal histograms will be saved to " << m_outputFileROOT.c_str() << endl;
42  }
43  m_nevtsample = ps.getUntrackedParameter<int>("nevtsample",0);
44 // for compatibility with previous versions
45  if(m_nevtsample==9999999) m_nevtsample=0;
46  m_pedsinADC = ps.getUntrackedParameter<int>("pedsinADC",0);
47  m_hiSaveflag = ps.getUntrackedParameter<int>("hiSaveflag",0);
48  m_pedValflag = ps.getUntrackedParameter<int>("pedValflag",0);
49  if(m_pedValflag<0) m_pedValflag=0;
50  if (m_nevtsample>0 && m_pedValflag>0) {
51  cout<<"WARNING - incompatible cfg options: nevtsample = "<<m_nevtsample<<", pedValflag = "<<m_pedValflag<<endl;
52  cout<<"Setting pedValflag = 0"<<endl;
53  m_pedValflag=0;
54  }
55  if(m_pedValflag>1) m_pedValflag=1;
56  m_startTS = ps.getUntrackedParameter<int>("firstTS", 0);
57  if(m_startTS<0) m_startTS=0;
58  m_endTS = ps.getUntrackedParameter<int>("lastTS", 9);
59 
60 // m_logFile.open("CastorPedestalAnalysis.log");
61 
62  castorHists.ALLPEDS = new TH1F("Castor All Pedestals","HF All Peds",10,0,9);
63  castorHists.PEDRMS= new TH1F("Castor All Pedestal Widths","HF All Pedestal RMS",100,0,3);
64  castorHists.PEDMEAN= new TH1F("Castor All Pedestal Means","HF All Pedestal Means",100,0,9);
65  castorHists.CHI2= new TH1F("Castor Chi2/ndf for whole range Gauss fit","HF Chi2/ndf Gauss",200,0.,50.);
66 }
67 
68 //-----------------------------------------------------------------------------
70 
71  for(_meot=castorHists.PEDTRENDS.begin(); _meot!=castorHists.PEDTRENDS.end(); _meot++){
72  for(int i=0; i<16; i++) _meot->second[i].first->Delete();
73  }
74 
75  castorHists.ALLPEDS->Delete();
76  castorHists.PEDRMS->Delete();
77  castorHists.PEDMEAN->Delete();
78  castorHists.CHI2->Delete();
79 }
80 
81 //-----------------------------------------------------------------------------
82 void CastorPedestalAnalysis::setup(const std::string& m_outputFileROOT) {
83  // open the histogram file, create directories within
84  m_file=new TFile(m_outputFileROOT.c_str(),"RECREATE");
85  m_file->mkdir("Castor");
86  m_file->cd();
87 }
88 
89 //-----------------------------------------------------------------------------
91  const CastorDbService& cond)
92 {
93  evt++;
94  sample=1;
95  evt_curr=evt;
96  if(m_nevtsample>0) {
97  sample = (evt-1)/m_nevtsample +1;
100  }
101 
102  m_shape = cond.getCastorShape();
103  // HF
104  try{
105  if(!castor.size()) throw (int)castor.size();
106  for (CastorDigiCollection::const_iterator j=castor.begin(); j!=castor.end(); j++){
107  const CastorDataFrame digi = (const CastorDataFrame)(*j);
108  m_coder = cond.getCastorCoder(digi.id());
109  for (int i=m_startTS; i<digi.size() && i<=m_endTS; i++) {
110  for(int flag=0; flag<4; flag++){
111  if(i+flag<digi.size() && i+flag<=m_endTS){
112  per2CapsHists(flag,2,digi.id(),digi.sample(i),digi.sample(i+flag),castorHists.PEDTRENDS,cond);
113  }
114  }
115  }
116  if(m_startTS==0 && m_endTS>4){
117  AllChanHists(digi.id(),digi.sample(0),digi.sample(1),digi.sample(2),digi.sample(3),digi.sample(4),digi.sample(5),castorHists.PEDTRENDS);
118  }
119  }
120  }
121  catch (int i ) {
122 // m_logFile << "Event with " << i<<" Castor Digis passed." << std::endl;
123  }
124  // Call the function every m_nevtsample events
125  if(m_nevtsample>0) {
127  }
128 }
129 
130 //-----------------------------------------------------------------------------
131 void CastorPedestalAnalysis::per2CapsHists(int flag, int id, const HcalDetId detid, const HcalQIESample& qie1, const HcalQIESample& qie2, map<HcalDetId, map<int,PEDBUNCH> > &toolT, const CastorDbService& cond) {
132 
133 // this function is due to be called for every time slice, it fills either a charge
134 // histo for a single capID (flag=0) or a product histo for two capIDs (flag>0)
135 
136  static const int bins=10;
137  static const int bins2=100;
138  float lo=-0.5; float hi=9.5;
139  map<int,PEDBUNCH> _mei;
140  static map<HcalDetId, map<int,float> > QieCalibMap;
141  string type = "Castor";
142 
143  /*
144  if(id==0){
145  if(detid.ieta()<16) type = "HB";
146  if(detid.ieta()>16) type = "HE";
147  if(detid.ieta()==16){
148  if(detid.depth()<3) type = "HB";
149  if(detid.depth()==3) type = "HE";
150  }
151  }
152  else if(id==1) type = "HO";
153  else if(id==2) type = "HF";
154  */
155 
156  _meot = toolT.find(detid);
157 
158 // if histos for the current channel do not exist, first create them,
159  if (_meot==toolT.end()){
160  map<int,PEDBUNCH> insert;
161  map<int,float> qiecalib;
162  char name[1024];
163  for(int i=0; i<4; i++){
164  lo=-0.5;
165  // fix from Andy: if you convert to fC and then bin in units of 1, you may 'skip' a bin while
166  // filling, since the ADCs are quantized
167  if (m_pedsinADC) hi=9.5;
168  else hi = 11.5;
169  sprintf(name,"%s Pedestal, eta=%d phi=%d d=%d cap=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth(),i);
170  insert[i].first = new TH1F(name,name,bins,lo,hi);
171  sprintf(name,"%s Product, eta=%d phi=%d d=%d caps=%d*%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth(),i,(i+1)%4);
172  insert[4+i].first = new TH1F(name,name,bins2,0.,100.);
173  sprintf(name,"%s Product, eta=%d phi=%d d=%d caps=%d*%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth(),i,(i+2)%4);
174  insert[8+i].first = new TH1F(name,name,bins2,0.,100.);
175  sprintf(name,"%s Product, eta=%d phi=%d d=%d caps=%d*%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth(),i,(i+3)%4);
176  insert[12+i].first = new TH1F(name,name,bins2,0.,100.);
177  }
178  sprintf(name,"%s Signal in TS 4+5, eta=%d phi=%d d=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
179  insert[16].first = new TH1F(name,name,21,-0.5,20.5);
180  sprintf(name,"%s Signal in TS 4+5-2-3, eta=%d phi=%d d=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
181  insert[17].first = new TH1F(name,name,21,-10.5,10.5);
182  sprintf(name,"%s Signal in TS 4+5-(0+1+2+3)/2., eta=%d phi=%d d=%d",type.c_str(),detid.ieta(),detid.iphi(),detid.depth());
183  insert[18].first = new TH1F(name,name,21,-10.5,10.5);
184  toolT[detid] = insert;
185  _meot = toolT.find(detid);
186 // store QIE calibrations in a map for later reuse
187  QieCalibMap[detid]=qiecalib;
188  }
189 
190  _mei = _meot->second;
191 
192  const CastorQIECoder* coder = cond.getCastorCoder(detid);
193  const CastorQIEShape* shape = cond.getCastorShape();
194  float charge1 = coder->charge(*shape,qie1.adc(),qie1.capid());
195  float charge2 = coder->charge(*shape,qie2.adc(),qie2.capid());
196 
197 // fill single capID histo
198  if(flag==0){
199  if(m_nevtsample>0) {
200  if((evt-1)%m_nevtsample==0 && state[qie1.capid()]){
201  state[qie1.capid()]=false;
202  _mei[qie1.capid()].first->Reset();
203  _mei[qie1.capid()+4].first->Reset();
204  _mei[qie1.capid()+8].first->Reset();
205  _mei[qie1.capid()+12].first->Reset();
206  }
207  }
208  if (qie1.adc()<bins){
209  if (m_pedsinADC) _mei[qie1.capid()].first->Fill(qie1.adc());
210  else _mei[qie1.capid()].first->Fill(charge1);
211  }
212  else if(qie1.adc()>=bins){
213  _mei[qie1.capid()].first->AddBinContent(bins+1,1);
214  }
215  }
216 
217 // fill 2 capID histo
218  if(flag>0){
219  map<int,float> qiecalib = QieCalibMap[detid];
220  //float charge1=(qie1.adc()-qiecalib[qie1.capid()+4])/qiecalib[qie1.capid()];
221  //float charge2=(qie2.adc()-qiecalib[qie2.capid()+4])/qiecalib[qie2.capid()];
222  if (charge1*charge2<bins2){
223  _mei[qie1.capid()+4*flag].first->Fill(charge1*charge2);
224  }
225  else{
226  _mei[qie1.capid()+4*flag].first->Fill(bins2);
227  }
228  }
229 
230  if(flag==0){
231  // if(id==0) hbHists.ALLPEDS->Fill(qie1.adc());
232  // else if(id==1) hoHists.ALLPEDS->Fill(qie1.adc());
233  // else if(id==2) castorHists.ALLPEDS->Fill(qie1.adc());
234  castorHists.ALLPEDS->Fill(qie1.adc());
235  }
236 }
237 
238 //-----------------------------------------------------------------------------
239 void CastorPedestalAnalysis::AllChanHists(const HcalDetId detid, const HcalQIESample& qie0, const HcalQIESample& qie1, const HcalQIESample& qie2, const HcalQIESample& qie3, const HcalQIESample& qie4, const HcalQIESample& qie5, map<HcalDetId, map<int,PEDBUNCH> > &toolT) {
240 
241 // this function is due to be called for every channel
242 
243  _meot = toolT.find(detid);
244  map<int,PEDBUNCH> _mei = _meot->second;
245  _mei[16].first->Fill(qie4.adc()+qie5.adc()-1.);
246  _mei[17].first->Fill(qie4.adc()+qie5.adc()-qie2.adc()-qie3.adc());
247  _mei[18].first->Fill(qie4.adc()+qie5.adc()-(qie0.adc()+qie1.adc()+qie2.adc()+qie3.adc())/2.);
248 }
249 
250 //-----------------------------------------------------------------------------
252  // it is called every m_nevtsample events (a sample) and the end of run
253  char PedSampleNum[20];
254 
255 // Compute pedestal constants for each HBHE, HO, HF
256  sprintf(PedSampleNum,"Castor_Sample%d",sample);
257  m_file->cd();
258  m_file->mkdir(PedSampleNum);
259  m_file->cd(PedSampleNum);
260  GetPedConst(castorHists.PEDTRENDS,castorHists.PEDMEAN,castorHists.PEDRMS);
261 }
262 
263 //-----------------------------------------------------------------------------
264 void CastorPedestalAnalysis::GetPedConst(map<HcalDetId, map<int,PEDBUNCH> > &toolT, TH1F* PedMeans, TH1F* PedWidths)
265 {
266 // Completely rewritten version oct 2006
267 // Compute pedestal constants and fill into CastorPedestals and CastorPedestalWidths objects
268  float cap[4]; float sig[4][4]; float dcap[4]; float dsig[4][4]; float chi2[4];
269 
270  for(_meot=toolT.begin(); _meot!=toolT.end(); _meot++){
271  HcalDetId detid = _meot->first;
272 
273 // take mean and width from a Gaussian fit or directly from the histo
274  if(fitflag>0){
275  for (int i=0; i<4; i++) {
276  TF1 *fit = _meot->second[i].first->GetFunction("gaus");
277  chi2[i]=0;
278  if(fit->GetNDF()!=0) chi2[i]=fit->GetChisquare()/fit->GetNDF();
279  cap[i]=fit->GetParameter(1);
280  sig[i][i]=fit->GetParameter(2);
281  dcap[i]=fit->GetParError(1);
282  dsig[i][i]=fit->GetParError(2);
283  }
284  }
285  else{
286  for (int i=0; i<4; i++) {
287  cap[i]=_meot->second[i].first->GetMean();
288  sig[i][i]=_meot->second[i].first->GetRMS();
289  m_stat[i]=0;
290 
291  for(int j=m_startTS; j<m_endTS+1; j++){
292  m_stat[i]+=_meot->second[i].first->GetBinContent(j+1);
293  }
294  dcap[i] = sig[i][i]/sqrt(m_stat[i]);
295 // dsig[i][i] = dcap[i]*sig[i][i]/cap[i];
296  dsig[i][i] = sig[i][i]/sqrt(2.*m_stat[i]);
297  chi2[i]=0.;
298  }
299  }
300 
301  for (int i=0; i<4; i++) {
302  if(m_hiSaveflag>0) {
303  if (m_pedsinADC)
304  _meot->second[i].first->GetXaxis()->SetTitle("ADC");
305  else _meot->second[i].first->GetXaxis()->SetTitle("Charge, fC");
306  _meot->second[i].first->GetYaxis()->SetTitle("CapID samplings");
307  _meot->second[i].first->Write();
308  }
309  if(m_nevtsample>0) {
310  _meot->second[i].second.first[0].push_back(cap[i]);
311  _meot->second[i].second.first[1].push_back(dcap[i]);
312  _meot->second[i].second.first[2].push_back(sig[i][i]);
313  _meot->second[i].second.first[3].push_back(dsig[i][i]);
314  _meot->second[i].second.first[4].push_back(chi2[i]);
315  }
316  PedMeans->Fill(cap[i]);
317  PedWidths->Fill(sig[i][i]);
318  }
319 
320 // special histos for Shuichi
321  if(m_hiSaveflag==-100){
322  for(int i=16; i<19; i++){
323  if (m_pedsinADC)
324  _meot->second[i].first->GetXaxis()->SetTitle("ADC");
325  else _meot->second[i].first->GetXaxis()->SetTitle("Charge, fC");
326  _meot->second[i].first->GetYaxis()->SetTitle("Events");
327  _meot->second[i].first->Write();
328  }
329  }
330 
331 // diagonal sigma is width squared
332  sig[0][0]=sig[0][0]*sig[0][0];
333  sig[1][1]=sig[1][1]*sig[1][1];
334  sig[2][2]=sig[2][2]*sig[2][2];
335  sig[3][3]=sig[3][3]*sig[3][3];
336 
337 // off diagonal sigmas (correlations) are computed from 3 histograms
338 // here we still have all 4*3=12 combinations
339  sig[0][1]= _meot->second[4].first->GetMean()-cap[0]*cap[1];
340  sig[0][2]= _meot->second[8].first->GetMean()-cap[0]*cap[2];
341  sig[1][2]= _meot->second[5].first->GetMean()-cap[1]*cap[2];
342  sig[1][3]= _meot->second[9].first->GetMean()-cap[1]*cap[3];
343  sig[2][3]= _meot->second[6].first->GetMean()-cap[2]*cap[3];
344  sig[0][3]= _meot->second[12].first->GetMean()-cap[0]*cap[3];
345  sig[1][0]= _meot->second[13].first->GetMean()-cap[1]*cap[0];
346  sig[2][0]= _meot->second[10].first->GetMean()-cap[2]*cap[0];
347  sig[2][1]= _meot->second[14].first->GetMean()-cap[2]*cap[1];
348  sig[3][1]= _meot->second[11].first->GetMean()-cap[3]*cap[1];
349  sig[3][2]= _meot->second[15].first->GetMean()-cap[3]*cap[2];
350  sig[3][0]= _meot->second[7].first->GetMean()-cap[3]*cap[0];
351 
352 // there is no proper error calculation for the correlation coefficients
353  for(int i=0; i<4; i++){
354  if(m_nevtsample>0) {
355  _meot->second[i].second.first[5].push_back(sig[i][(i+1)%4]);
356  _meot->second[i].second.first[6].push_back(2*sig[i][i]*dsig[i][i]);
357  _meot->second[i].second.first[7].push_back(sig[i][(i+2)%4]);
358  _meot->second[i].second.first[8].push_back(2*sig[i][i]*dsig[i][i]);
359  _meot->second[i].second.first[9].push_back(sig[i][(i+3)%4]);
360  _meot->second[i].second.first[10].push_back(2*sig[i][i]*dsig[i][i]);
361  }
362 // save product histos if desired
363  if(m_hiSaveflag>10) {
364  if (m_pedsinADC)
365  _meot->second[i+4].first->GetXaxis()->SetTitle("ADC^2");
366  else _meot->second[i+4].first->GetXaxis()->SetTitle("Charge^2, fC^2");
367  _meot->second[i+4].first->GetYaxis()->SetTitle("2-CapID samplings");
368  _meot->second[i+4].first->Write();
369  if (m_pedsinADC)
370  _meot->second[i+8].first->GetXaxis()->SetTitle("ADC^2");
371  else _meot->second[i+8].first->GetXaxis()->SetTitle("Charge^2, fC^2");
372  _meot->second[i+8].first->GetYaxis()->SetTitle("2-CapID samplings");
373  _meot->second[i+8].first->Write();
374  if (m_pedsinADC)
375  _meot->second[i+12].first->GetXaxis()->SetTitle("ADC^2");
376  else _meot->second[i+12].first->GetXaxis()->SetTitle("Charge^2, fC^2");
377  _meot->second[i+12].first->GetYaxis()->SetTitle("2-CapID samplings");
378  _meot->second[i+12].first->Write();
379  }
380  }
381 
382 // fill the objects - at this point only close and medium correlations are stored
383 // and the matrix is assumed symmetric
384  if (m_nevtsample<1) {
385  sig[1][0]=sig[0][1];
386  sig[2][0]=sig[0][2];
387  sig[2][1]=sig[1][2];
388  sig[3][1]=sig[1][3];
389  sig[3][2]=sig[2][3];
390  sig[0][3]=sig[3][0];
391  if (fRawPedestals) {
392  CastorPedestal item(detid,cap[0],cap[1],cap[2],cap[3]);
393  fRawPedestals->addValues(item);
394  }
395  if (fRawPedestalWidths) {
396  CastorPedestalWidth widthsp(detid);
397  widthsp.setSigma(0,0,sig[0][0]);
398  widthsp.setSigma(0,1,sig[0][1]);
399  widthsp.setSigma(0,2,sig[0][2]);
400  widthsp.setSigma(1,1,sig[1][1]);
401  widthsp.setSigma(1,2,sig[1][2]);
402  widthsp.setSigma(1,3,sig[1][3]);
403  widthsp.setSigma(2,2,sig[2][2]);
404  widthsp.setSigma(2,3,sig[2][3]);
405  widthsp.setSigma(3,3,sig[3][3]);
406  widthsp.setSigma(3,0,sig[0][3]);
407  fRawPedestalWidths->addValues(widthsp);
408  }
409  }
410  }
411 }
412 
413 //-----------------------------------------------------------------------------
414 int CastorPedestalAnalysis::done(const CastorPedestals* fInputPedestals,
415  const CastorPedestalWidths* fInputPedestalWidths,
416  CastorPedestals* fOutputPedestals,
417  CastorPedestalWidths* fOutputPedestalWidths)
418 {
419  int nstat[4];
420 
421 // Pedestal objects
422  // inputs...
423  fRefPedestals = fInputPedestals;
424  fRefPedestalWidths = fInputPedestalWidths;
425 
426  // outputs...
427  if(m_pedValflag>0) {
428  fValPedestals = fOutputPedestals;
429  fValPedestalWidths = fOutputPedestalWidths;
432  }
433  else {
434  fRawPedestals = fOutputPedestals;
435  fRawPedestalWidths = fOutputPedestalWidths;
438  }
439 
440 // compute pedestal constants
442  if(m_nevtsample>0) {
444  }
445 
446 // trending histos
447  if(m_nevtsample>0){
448  m_file->cd();
449  m_file->cd("Castor");
450  Trendings(castorHists.PEDTRENDS,castorHists.CHI2,castorHists.CAPID_AVERAGE,castorHists.CAPID_CHI2);
451  }
452 
453  if (m_nevtsample<1) {
454 
455 // pedestal validation: m_AllPedsOK=-1 means not validated,
456 // 0 everything OK,
457 // N>0 : mod(N,100000) drifts + width changes
458 // int(N/100000) missing channels
459  m_AllPedsOK=-1;
460  if(m_pedValflag>0) {
461  for (int i=0; i<4; i++) nstat[i]=(int)m_stat[i];
462  int NPedErrors=CastorPedVal(nstat,fRefPedestals,fRefPedestalWidths,
465  m_AllPedsOK=NPedErrors;
466  }
467 // setting m_AllPedsOK=-2 will inhibit writing pedestals out
468 // if(m_pedValflag==1){
469 // if(evt<100)m_AllPedsOK=-2;
470 // }
471 
472  }
473 
474  // Write other histograms.
475 
476  // Castor
477  m_file->cd();
478  m_file->cd("Castor");
479  castorHists.ALLPEDS->Write();
480  castorHists.PEDRMS->Write();
481  castorHists.PEDMEAN->Write();
482 
483  m_file->Close();
484  cout << "Hcal/Castor histograms written to " << m_outputFileROOT.c_str() << endl;
485  return (int)m_AllPedsOK;
486 }
487 
488 //-----------------------------------------------------------------------------
489 void CastorPedestalAnalysis::Trendings(map<HcalDetId, map<int,PEDBUNCH> > &toolT, TH1F* Chi2, TH1F* CapidAverage, TH1F* CapidChi2){
490 
491 // check stability of pedestal constants in a single long run
492 
493  map<int, std::vector<double> > AverageValues;
494 
495  for(_meot=toolT.begin(); _meot!=toolT.end(); _meot++){
496  for(int i=0; i<4; i++){
497  char name[1024];
498  HcalDetId detid = _meot->first;
499  sprintf(name,"Pedestal trend, eta=%d phi=%d d=%d cap=%d",detid.ieta(),detid.iphi(),detid.depth(),i);
500  int bins = _meot->second[i].second.first[0].size();
501  float lo =0.5;
502  float hi = (float)bins+0.5;
503  _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi));
504  sprintf(name,"Width trend, eta=%d phi=%d d=%d cap=%d",detid.ieta(),detid.iphi(),detid.depth(),i);
505  bins = _meot->second[i].second.first[2].size();
506  hi = (float)bins+0.5;
507  _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi));
508  sprintf(name,"Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",detid.ieta(),detid.iphi(),detid.depth(),i,(i+1)%4);
509  bins = _meot->second[i].second.first[5].size();
510  hi = (float)bins+0.5;
511  _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi));
512 /* sprintf(name,"Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",detid.ieta(),detid.iphi(),detid.depth(),i,(i+2)%4);
513  bins = _meot->second[i].second.first[7].size();
514  hi = (float)bins+0.5;
515  _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi));
516  sprintf(name,"Correlation trend, eta=%d phi=%d d=%d caps=%d*%d",detid.ieta(),detid.iphi(),detid.depth(),i,(i+3)%4);
517  bins = _meot->second[i].second.first[9].size();
518  hi = (float)bins+0.5;
519  _meot->second[i].second.second.push_back(new TH1F(name,name,bins,lo,hi)); */
520 
521  std::vector<double>::iterator sample_it;
522  // Pedestal mean - put content and errors
523  int j=0;
524  for(sample_it=_meot->second[i].second.first[0].begin();
525  sample_it!=_meot->second[i].second.first[0].end();sample_it++){
526  _meot->second[i].second.second[0]->SetBinContent(++j,*sample_it);
527  }
528  j=0;
529  for(sample_it=_meot->second[i].second.first[1].begin();
530  sample_it!=_meot->second[i].second.first[1].end();sample_it++){
531  _meot->second[i].second.second[0]->SetBinError(++j,*sample_it);
532  }
533  // fit with a constant - extract parameters
534  _meot->second[i].second.second[0]->Fit("pol0","Q");
535  TF1 *fit = _meot->second[i].second.second[0]->GetFunction("pol0");
536  AverageValues[0].push_back(fit->GetParameter(0));
537  AverageValues[1].push_back(fit->GetParError(0));
538  if(sample>1)
539  AverageValues[2].push_back(fit->GetChisquare()/fit->GetNDF());
540  else
541  AverageValues[2].push_back(fit->GetChisquare());
542  sprintf(name,"Sample (%d events)",m_nevtsample);
543  _meot->second[i].second.second[0]->GetXaxis()->SetTitle(name);
544  _meot->second[i].second.second[0]->GetYaxis()->SetTitle("Pedestal value");
545  _meot->second[i].second.second[0]->Write();
546  // Pedestal width - put content and errors
547  j=0;
548  for(sample_it=_meot->second[i].second.first[2].begin();
549  sample_it!=_meot->second[i].second.first[2].end();sample_it++){
550  _meot->second[i].second.second[1]->SetBinContent(++j,*sample_it);
551  }
552  j=0;
553  for(sample_it=_meot->second[i].second.first[3].begin();
554  sample_it!=_meot->second[i].second.first[3].end();sample_it++){
555  _meot->second[i].second.second[1]->SetBinError(++j,*sample_it);
556  }
557  _meot->second[i].second.second[1]->GetXaxis()->SetTitle(name);
558  _meot->second[i].second.second[1]->GetYaxis()->SetTitle("Pedestal width");
559  _meot->second[i].second.second[1]->Write();
560  // Correlation coeffs - put contents and errors
561  j=0;
562  for(sample_it=_meot->second[i].second.first[5].begin();
563  sample_it!=_meot->second[i].second.first[5].end();sample_it++){
564  _meot->second[i].second.second[2]->SetBinContent(++j,*sample_it);
565  }
566  j=0;
567  for(sample_it=_meot->second[i].second.first[6].begin();
568  sample_it!=_meot->second[i].second.first[6].end();sample_it++){
569  _meot->second[i].second.second[2]->SetBinError(++j,*sample_it);
570  }
571  _meot->second[i].second.second[2]->GetXaxis()->SetTitle(name);
572  _meot->second[i].second.second[2]->GetYaxis()->SetTitle("Close correlation");
573  _meot->second[i].second.second[2]->Write();
574  /* j=0;
575  for(sample_it=_meot->second[i].second.first[7].begin();
576  sample_it!=_meot->second[i].second.first[7].end();sample_it++){
577  _meot->second[i].second.second[3]->SetBinContent(++j,*sample_it);
578  }
579  j=0;
580  for(sample_it=_meot->second[i].second.first[8].begin();
581  sample_it!=_meot->second[i].second.first[8].end();sample_it++){
582  _meot->second[i].second.second[3]->SetBinError(++j,*sample_it);
583  }
584  _meot->second[i].second.second[3]->GetXaxis()->SetTitle(name);
585  _meot->second[i].second.second[3]->GetYaxis()->SetTitle("Intermediate correlation");
586  _meot->second[i].second.second[3]->Write();
587  j=0;
588  for(sample_it=_meot->second[i].second.first[9].begin();
589  sample_it!=_meot->second[i].second.first[9].end();sample_it++){
590  _meot->second[i].second.second[4]->SetBinContent(++j,*sample_it);
591  }
592  j=0;
593  for(sample_it=_meot->second[i].second.first[10].begin();
594  sample_it!=_meot->second[i].second.first[10].end();sample_it++){
595  _meot->second[i].second.second[4]->SetBinError(++j,*sample_it);
596  }
597  _meot->second[i].second.second[4]->GetXaxis()->SetTitle(name);
598  _meot->second[i].second.second[4]->GetYaxis()->SetTitle("Distant correlation");
599  _meot->second[i].second.second[4]->Write(); */
600  // chi2
601  j=0;
602  for(sample_it=_meot->second[i].second.first[4].begin();
603  sample_it!=_meot->second[i].second.first[4].end();sample_it++){
604  Chi2->Fill(*sample_it);
605  }
606  }
607  }
608  CapidAverage= new TH1F("Constant fit: Pedestal Values",
609  "Constant fit: Pedestal Values",
610  AverageValues[0].size(),0.,AverageValues[0].size());
611  std::vector<double>::iterator sample_it;
612  int j=0;
613  for(sample_it=AverageValues[0].begin();
614  sample_it!=AverageValues[0].end();sample_it++){
615  CapidAverage->SetBinContent(++j,*sample_it);
616  }
617  j=0;
618  for(sample_it=AverageValues[1].begin();
619  sample_it!=AverageValues[1].end();sample_it++){
620  CapidAverage->SetBinError(++j,*sample_it);
621  }
622  CapidChi2= new TH1F("Constant fit: Chi2/ndf",
623  "Constant fit: Chi2/ndf",
624  AverageValues[2].size(),0.,AverageValues[2].size());
625  j=0;
626  for(sample_it=AverageValues[2].begin();
627  sample_it!=AverageValues[2].end();sample_it++){
628  CapidChi2->SetBinContent(++j,*sample_it);
629  //CapidChi2->SetBinError(++j,0);
630  }
631  Chi2->GetXaxis()->SetTitle("Chi2/ndf");
632  Chi2->GetYaxis()->SetTitle("50 x [(16+2) x 4 x 4] `events`");
633  Chi2->Write();
634  CapidAverage->GetYaxis()->SetTitle("Pedestal value");
635  CapidAverage->GetXaxis()->SetTitle("(16+2) x 4 x 4 `events`");
636  CapidAverage->Write();
637  CapidChi2->GetYaxis()->SetTitle("Chi2/ndf");
638  CapidChi2->GetXaxis()->SetTitle("(16+2) x 4 x 4 `events`");
639  CapidChi2->Write();
640 
641 }
642 
643 //-----------------------------------------------------------------------------
644 int CastorPedestalAnalysis::CastorPedVal(int nstat[4], const CastorPedestals* fRefPedestals,
645  const CastorPedestalWidths* fRefPedestalWidths,
646  CastorPedestals* fRawPedestals,
647  CastorPedestalWidths* fRawPedestalWidths,
648  CastorPedestals* fValPedestals,
649  CastorPedestalWidths* fValPedestalWidths)
650 {
651 // new version of pedestal validation - it is designed to be as independent of
652 // all the rest as possible - you only need to provide valid pedestal objects
653 // and a vector of statistics per capID to use this as standalone code
655  float RefPedVals[4]; float RefPedSigs[4][4];
656  float RawPedVals[4]; float RawPedSigs[4][4];
657  map<HcalDetId,bool> isinRaw;
658  map<HcalDetId,bool> isinRef;
659  std::vector<DetId> RefChanns=fRefPedestals->getAllChannels();
660  std::vector<DetId> RawChanns=fRawPedestals->getAllChannels();
661  std::ofstream PedValLog;
662  PedValLog.open("CastorPedVal.log");
663 
664  if(nstat[0]+nstat[1]+nstat[2]+nstat[3]<2500) PedValLog<<"CastorPedVal: warning - low statistics"<<std::endl;
665 // find complete list of channels in current data and reference
666  for (int i=0; i<(int)RawChanns.size(); i++){
667  isinRef[HcalDetId(RawChanns[i])]=false;
668  }
669  for (int i=0; i<(int)RefChanns.size(); i++){
670  detid=HcalDetId(RefChanns[i]);
671  isinRaw[detid]=false;
672  isinRef[detid]=true;
673  }
674  for (int i=0; i<(int)RawChanns.size(); i++){
675  detid=HcalDetId(RawChanns[i]);
676  isinRaw[detid]=true;
677  if (isinRef[detid]==false) {
678  PedValLog<<"CastorPedVal: channel "<<detid<<" not found in reference set"<<std::endl;
679  std::cerr<<"CastorPedVal: channel "<<detid<<" not found in reference set"<<std::endl;
680  }
681  }
682 
683 // main loop over channels
684  int erflag=0;
685  for (int i=0; i<(int)RefChanns.size(); i++){
686  detid=HcalDetId(RefChanns[i]);
687  for (int icap=0; icap<4; icap++) {
688  RefPedVals[icap]=fRefPedestals->getValues(detid)->getValue(icap);
689  for (int icap2=icap; icap2<4; icap2++) {
690  RefPedSigs[icap][icap2]=fRefPedestalWidths->getValues(detid)->getSigma(icap,icap2);
691  if(icap2!=icap)RefPedSigs[icap2][icap]=RefPedSigs[icap][icap2];
692  }
693  }
694 
695 // read new raw values
696  if(isinRaw[detid]) {
697  for (int icap=0; icap<4; icap++) {
698  RawPedVals[icap]=fRawPedestals->getValues(detid)->getValue(icap);
699  for (int icap2=icap; icap2<4; icap2++) {
700  RawPedSigs[icap][icap2]=fRawPedestalWidths->getValues(detid)->getSigma(icap,icap2);
701  if(icap2!=icap)RawPedSigs[icap2][icap]=RawPedSigs[icap][icap2];
702  }
703  }
704 
705 // first quick check if raw values make sense: if not, the channel is treated like absent
706  for (int icap=0; icap<4; icap++) {
707  if(RawPedVals[icap]<1. || RawPedSigs[icap][icap]<0.01) isinRaw[detid]=false;
708  for (int icap2=icap; icap2<4; icap2++){
709  if(fabs(RawPedSigs[icap][icap2]/sqrt(RawPedSigs[icap][icap]*RawPedSigs[icap2][icap2]))>1.) isinRaw[detid]=false;
710  }
711  }
712  }
713 
714 // check raw values against reference
715  if(isinRaw[detid]) {
716  for (int icap=0; icap<4; icap++) {
717  int icap2=(icap+1)%4;
718  float width=sqrt(RawPedSigs[icap][icap]);
719  float erof1=width/sqrt((float)nstat[icap]);
720  float erof2=sqrt(erof1*erof1+RawPedSigs[icap][icap]/(float)nstat[icap]);
721  float erofwidth=width/sqrt(2.*nstat[icap]);
722  float diffof1=RawPedVals[icap]-RefPedVals[icap];
723  float diffof2=RawPedVals[icap]+RawPedVals[icap2]-RefPedVals[icap]-RefPedVals[icap2];
724  float diffofw=width-sqrt(RefPedSigs[icap][icap]);
725 
726 // validation in 2 TS for HB, HE, HO, in 1 TS for HF
727  int nTS=2;
728  if(detid.subdet()==HcalForward) nTS=1;
729  if(nTS==1 && fabs(diffof1)>0.5+erof1) {
730  erflag+=1;
731  PedValLog<<"HcalPedVal: drift in channel "<<detid<<" cap "<<icap<<": "<<RawPedVals[icap]<<" - "<<RefPedVals[icap]<<" = "<<diffof1<<std::endl;
732  }
733  if(nTS==2 && fabs(diffof2)>0.5+erof2) {
734  erflag+=1;
735  PedValLog<<"HcalPedVal: drift in channel "<<detid<<" caps "<<icap<<"+"<<icap2<<": "<<RawPedVals[icap]<<"+"<<RawPedVals[icap2]<<" - "<<RefPedVals[icap]<<"+"<<RefPedVals[icap2]<<" = "<<diffof2<<std::endl;
736  }
737  if(fabs(diffofw)>0.15*width+erofwidth) {
738  erflag+=1;
739  PedValLog<<"HcalPedVal: width changed in channel "<<detid<<" cap "<<icap<<": "<<width<<" - "<<sqrt(RefPedSigs[icap][icap])<<" = "<<diffofw<<std::endl;
740  }
741  }
742  }
743 
744 // for disconnected/bad channels restore reference values
745  else {
746  PedValLog<<"HcalPedVal: no valid data from channel "<<detid<<std::endl;
747  erflag+=100000;
748  CastorPedestal item(detid,RefPedVals[0],RefPedVals[1],RefPedVals[2],RefPedVals[3]);
749  fValPedestals->addValues(item);
750  CastorPedestalWidth widthsp(detid);
751  for (int icap=0; icap<4; icap++) {
752  for (int icap2=icap; icap2<4; icap2++) widthsp.setSigma(icap2,icap,RefPedSigs[icap2][icap]);
753  }
754  fValPedestalWidths->addValues(widthsp);
755  }
756 
757 // end of channel loop
758  }
759 
760  if(erflag==0) PedValLog<<"HcalPedVal: all pedestals checked OK"<<std::endl;
761 
762 // now construct the remaining part of the validated objects
763 // if nothing changed outside tolerance, validated set = reference set
764  if(erflag%100000 == 0) {
765  for (int i=0; i<(int)RefChanns.size(); i++){
766  detid=HcalDetId(RefChanns[i]);
767  if (isinRaw[detid]) {
768  CastorPedestalWidth widthsp(detid);
769  for (int icap=0; icap<4; icap++) {
770  RefPedVals[icap]=fRefPedestals->getValues(detid)->getValue(icap);
771  for (int icap2=icap; icap2<4; icap2++) {
772  RefPedSigs[icap][icap2]=fRefPedestalWidths->getValues(detid)->getSigma(icap,icap2);
773  if(icap2!=icap)RefPedSigs[icap2][icap]=RefPedSigs[icap][icap2];
774  widthsp.setSigma(icap2,icap,RefPedSigs[icap2][icap]);
775  }
776  }
777  fValPedestalWidths->addValues(widthsp);
778  CastorPedestal item(detid,RefPedVals[0],RefPedVals[1],RefPedVals[2],RefPedVals[3]);
779  fValPedestals->addValues(item);
780  }
781  }
782  }
783 
784 // if anything changed, validated set = raw set + reference for missing/bad channels
785  else {
786  for (int i=0; i<(int)RawChanns.size(); i++){
787  detid=HcalDetId(RawChanns[i]);
788  if (isinRaw[detid]) {
789  CastorPedestalWidth widthsp(detid);
790  for (int icap=0; icap<4; icap++) {
791  RawPedVals[icap]=fRawPedestals->getValues(detid)->getValue(icap);
792  for (int icap2=icap; icap2<4; icap2++) {
793  RawPedSigs[icap][icap2]=fRawPedestalWidths->getValues(detid)->getSigma(icap,icap2);
794  if(icap2!=icap)RawPedSigs[icap2][icap]=RawPedSigs[icap][icap2];
795  widthsp.setSigma(icap2,icap,RawPedSigs[icap2][icap]);
796  }
797  }
798  fValPedestalWidths->addValues(widthsp);
799  CastorPedestal item(detid,RawPedVals[0],RawPedVals[1],RawPedVals[2],RawPedVals[3]);
800  fValPedestals->addValues(item);
801  }
802  }
803  }
804  return erflag;
805 }
806 
type
Definition: HCALResponse.h:21
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
float getValue(int fCapId) const
get value for capId = 0..3
static const double cap
Definition: HcalTimeSlew.cc:7
void setSigma(int fCapId1, int fCapId2, float fSigma)
CastorPedestalWidths * fRawPedestalWidths
std::map< HcalDetId, std::map< int, PEDBUNCH > >::iterator _meot
struct CastorPedestalAnalysis::@41 castorHists
const CastorQIECoder * m_coder
int adc() const
get the ADC sample
Definition: HcalQIESample.h:22
const HcalQIESample & sample(int i) const
access a sample
std::vector< CastorDataFrame >::const_iterator const_iterator
std::vector< DetId > getAllChannels() const
static int CastorPedVal(int nstat[4], const CastorPedestals *fRefPedestals, const CastorPedestalWidths *fRefPedestalWidths, CastorPedestals *fRawPedestals, CastorPedestalWidths *fRawPedestalWidths, CastorPedestals *fValPedestals, CastorPedestalWidths *fValPedestalWidths)
const Item * getValues(DetId fId, bool throwOnFail=true) const
void setup(const std::string &m_outputFileROOT)
void processEvent(const CastorDigiCollection &castor, const CastorDbService &cond)
int depth() const
get the tower depth
Definition: HcalDetId.h:55
void Trendings(std::map< HcalDetId, std::map< int, PEDBUNCH > > &toolT, TH1F *Chi2, TH1F *CapidAverage, TH1F *CapidChi2)
CastorPedestalWidths * fValPedestalWidths
void AllChanHists(const HcalDetId detid, const HcalQIESample &qie0, const HcalQIESample &qie1, const HcalQIESample &qie2, const HcalQIESample &qie3, const HcalQIESample &qie4, const HcalQIESample &qie5, std::map< HcalDetId, std::map< int, PEDBUNCH > > &toolT)
T sqrt(T t)
Definition: SSEVec.h:48
float getSigma(int fCapId1, int fCapId2) const
get correlation element for capId1/2 = 0..3
const CastorQIEShape * getCastorShape() const
int ieta() const
get the cell ieta
Definition: HcalDetId.h:51
float charge(const CastorQIEShape &fShape, unsigned fAdc, unsigned fCapId) const
ADC [0..127] + capid [0..3] -&gt; fC conversion.
int j
Definition: DBlmapReader.cc:9
bool insert(Storage &iStorage, ItemType *iItem, const IdTag &iIdTag)
Definition: HCMethods.h:49
const_iterator end() const
int iphi() const
get the cell iphi
Definition: HcalDetId.h:53
int capid() const
get the Capacitor id
Definition: HcalQIESample.h:26
const CastorQIEShape * m_shape
int done(const CastorPedestals *fInputPedestals, const CastorPedestalWidths *fInputWidths, CastorPedestals *fOutputPedestals, CastorPedestalWidths *fOutputWidths)
bool addValues(const Item &myItem)
#define begin
Definition: vmac.h:30
const CastorPedestalWidths * fRefPedestalWidths
CastorPedestalAnalysis(const edm::ParameterSet &ps)
Constructor.
size_type size() const
void GetPedConst(std::map< HcalDetId, std::map< int, PEDBUNCH > > &toolT, TH1F *PedMeans, TH1F *PedWidths)
tuple cout
Definition: gather_cfg.py:121
Definition: Chi2.h:17
const HcalCastorDetId & id() const
int size() const
total number of samples in the digi
tuple size
Write out results.
void per2CapsHists(int flag, int id, const HcalDetId detid, const HcalQIESample &qie1, const HcalQIESample &qie2, std::map< HcalDetId, std::map< int, PEDBUNCH > > &toolT, const CastorDbService &cond)
const_iterator begin() const
const CastorQIECoder * getCastorCoder(const HcalGenericDetId &fId) const
const CastorPedestals * fRefPedestals