CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
HFLightCalRand.cc
Go to the documentation of this file.
1 // Analysis of HF LED/Laser run:
2 // Case when signal has random time position in TS
3 // SPE calibration for low light intensity or raw SPE calibration for high light intensity
4 // and HF performance based on this analysis
5 //
6 // Igor Vodopiyanov. Oct-2007 .... update Sept-2008
7 // Thanks G.Safronov, M.Mohammadi, F.Ratnikov
8 //
9 #include <memory>
10 #include <string>
11 #include <iostream>
12 
13 #include "TH1F.h"
14 #include "TH2F.h"
15 #include "TFile.h"
16 #include "math.h"
17 #include "TMath.h"
18 #include "TF1.h"
19 
20 #include "CalibCalorimetry/HcalStandardModules/interface/HFLightCalRand.h"
21 
22 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
23 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
24 
25 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
26 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
27 #include "CondFormats/HcalObjects/interface/HcalQIEShape.h"
28 #include "CondFormats/HcalObjects/interface/HcalQIECoder.h"
29 
30 #include "FWCore/Framework/interface/Event.h"
31 #include "FWCore/ParameterSet/interface/ParameterSet.h"
32 #include "FWCore/Framework/interface/EventSetup.h"
33 #include "FWCore/Framework/interface/ESHandle.h"
34 
35 using namespace std;
36 Int_t NEvents, runNumb=0, EventN=0;
37 
38 namespace {
39  //bool verbose = true;
40  bool verbose = false;
41 }
42 
43 HFLightCalRand::HFLightCalRand (const edm::ParameterSet& fConfiguration) {
44  //std::string histfile = fConfiguration.getUntrackedParameter<string>("rootFile");
45  histfile = fConfiguration.getUntrackedParameter<string>("rootFile");
46  textfile = fConfiguration.getUntrackedParameter<string>("textFile");
47 }
48 
49 HFLightCalRand::~HFLightCalRand () {
50  //delete mFile;
51 }
52 
53 void HFLightCalRand::beginJob() {
54 
55  char htit[64];
56  mFile = new TFile (histfile.c_str(),"RECREATE");
57  if ((tFile = fopen(textfile.c_str(),"w"))==NULL)printf("\nNo textfile open\n\n");
58  // General Histos
59  htmax = new TH1F("htmax","Max TS",10,-0.5,9.5);
60  htmean = new TH1F("htmean","Mean signal TS",100,0,10);
61  hsignalmean = new TH1F("hsignalmean","Mean ADC 4maxTS",1201,-25,30000);
62  hsignalrms = new TH1F("hsignalrms","RMS ADC 4maxTS",500,0,500);
63  hpedmean = new TH1F("hpedmean","Mean ADC 4lowTS",200,-10,90);
64  hpedrms = new TH1F("hpedrms","RMS ADC 4lowTS",200,0,100);
65  hspes = new TH1F("hspes","SPE if measured",200,0,40);
66  hnpevar = new TH1F("hnpevar","~N PE input",500,0,500);
67  hsignalmapP = new TH2F("hsignalmapP","Mean(Response) - Mean(Pedestal) HFP;#eta;#phi",26,28.5,41.5,36,0,72);
68  hsignalRMSmapP = new TH2F("hsignalRMSmapP","RMS Response HFP;#eta;#phi",26,28.5,41.5,36,0,72);
69  hnpemapP = new TH2F("hnpemapP","~N PE input HFP;#eta;#phi",26,28.5,41.5,36,0,72);
70  hnpemapP->SetOption("COLZ");hsignalmapP->SetOption("COLZ");hsignalRMSmapP->SetOption("COLZ");
71  hsignalmapM = new TH2F("hsignalmapM","Mean(Response) - Mean(Pedestal) HFM;#eta;#phi",26,-41.5,-28.5,36,0,72);
72  hsignalRMSmapM = new TH2F("hsignalRMSmapM","RMS Response HFM;#eta;#phi",26,-41.5,-28.5,36,0,72);
73  hnpemapM = new TH2F("hnpemapM","~N PE input HFM;#eta;#phi",26,-41.5,-28.5,36,0,72);
74  hnpemapM->SetOption("COLZ");hsignalmapM->SetOption("COLZ");hsignalRMSmapM->SetOption("COLZ");
75  // Channel-by-channel histos
76  for (int i=0;i<13;i++) for (int j=0;j<36;j++) for (int k=0;k<2;k++) {
77  if (i>10 && j%2==0) continue;
78  sprintf(htit,"ts_+%d_%d_%d",i+29,j*2+1,k+1);
79  hts[i][j][k] = new TH1F(htit,htit,10,-0.5,9.5); // TimeSlices (pulse shape)
80  sprintf(htit,"tsmean_+%d_%d_%d",i+29,j*2+1,k+1);
81  htsm[i][j][k] = new TH1F(htit,htit,100,0,10); // Mean signal time estimated from TS
82  sprintf(htit,"sp_+%d_%d_%d",i+29,j*2+1,k+1);
83  hsp[i][j][k] = new TH1F(htit,htit,1201,-25,30000); // Big-scale spectrum (linear ADC)
84  sprintf(htit,"spe_+%d_%d_%d",i+29,j*2+1,k+1);
85  hspe[i][j][k] = new TH1F(htit,htit,200,-9.5,190.5); // Small-scale spectrum (linear ADC)
86  sprintf(htit,"ped_+%d_%d_%d",i+29,j*2+1,k+1);
87  hped[i][j][k] = new TH1F(htit,htit,200,-9.5,190.5); // Pedestal spectrum
88  sprintf(htit,"ts_-%d_%d_%d",i+29,j*2+1,k+1);
89  hts[i+13][j][k] = new TH1F(htit,htit,10,-0.5,9.5);
90  sprintf(htit,"tsmean_-%d_%d_%d",i+29,j*2+1,k+1);
91  htsm[i+13][j][k] = new TH1F(htit,htit,100,0,10);
92  sprintf(htit,"sp_-%d_%d_%d",i+29,j*2+1,k+1);
93  hsp[i+13][j][k] = new TH1F(htit,htit,1201,-25,30000);
94  sprintf(htit,"spe_-%d_%d_%d",i+29,j*2+1,k+1);
95  hspe[i+13][j][k] = new TH1F(htit,htit,200,-9.5,190.5);
96  sprintf(htit,"ped_-%d_%d_%d",i+29,j*2+1,k+1);
97  hped[i+13][j][k] = new TH1F(htit,htit,200,-9.5,190.5);
98  }
99  // PIN-diodes histos
100  for (int i=0;i<4;i++) for (int j=0;j<3;j++) {
101  sprintf(htit,"ts_PIN%d_+Q%d",j+1,i+1);
102  htspin[i][j] = new TH1F(htit,htit,10,-0.5,9.5);
103  sprintf(htit,"sp_PIN%d_+Q%d",j+1,i+1);
104  hsppin[i][j] = new TH1F(htit,htit,1601,-25,40000);
105  sprintf(htit,"spe_PIN%d_+Q%d",j+1,i+1);
106  hspepin[i][j] = new TH1F(htit,htit,200,-9.5,190.5);
107  sprintf(htit,"ped_PIN%d_+Q%d",j+1,i+1);
108  hpedpin[i][j] = new TH1F(htit,htit,200,-9.5,190.5);
109  sprintf(htit,"tsmean_PIN%d_+Q%d",j+1,i+1);
110  htsmpin[i][j] = new TH1F(htit,htit,100,0,10);
111  sprintf(htit,"ts_PIN%d_-Q%d",j+1,i+1);
112  htspin[i+4][j] = new TH1F(htit,htit,10,-0.5,9.5);
113  sprintf(htit,"sp_PIN%d_-Q%d",j+1,i+1);
114  hsppin[i+4][j] = new TH1F(htit,htit,1601,-25,40000);
115  sprintf(htit,"spe_PIN%d_-Q%d",j+1,i+1);
116  hspepin[i+4][j] = new TH1F(htit,htit,200,-9.5,190.5);
117  sprintf(htit,"ped_PIN%d_-Q%d",j+1,i+1);
118  hpedpin[i+4][j] = new TH1F(htit,htit,200,-9.5,190.5);
119  sprintf(htit,"tsmean_PIN%d_-Q%d",j+1,i+1);
120  htsmpin[i+4][j] = new TH1F(htit,htit,100,0,10);
121  }
122  std::cout<<std::endl<<"beginJob: histfile="<<histfile.c_str()<<" textfile="<<textfile.c_str()<<std::endl;
123  return;
124 }
125 
126 void HistSpec(TH1F* hist, Double_t &mean, Double_t &rms, Double_t range=4) {
127  Double_t xmin,xmax;
128  mean=hist->GetMean();
129  rms=hist->GetRMS();
130  xmin=hist->GetXaxis()->GetXmin();
131  xmax=hist->GetXaxis()->GetXmax();
132  hist->SetAxisRange(mean-range*rms-100,mean+range*rms+100);
133  mean=hist->GetMean();
134  rms=hist->GetRMS();
135  hist->SetAxisRange(mean-range*rms-100,mean+range*rms+100);
136  mean=hist->GetMean();
137  rms=hist->GetRMS();
138  hist->SetAxisRange(xmin,xmax);
139  return;
140 }
141 
142 Double_t Fit3Peak(Double_t *x, Double_t *par) {
143 // Spectra fit function: Pedestal Gaussian + asymmetric 1PE + 2PE +3PE peaks
144 
145  Double_t sum,xx,A0,C0,r0,sigma0,mean1,sigma1,A1,C1,r1,mean2,sigma2,A2,C2,r2,mean3,sigma3,A3,C3,r3;
146 
147  const Double_t k0=2.0,k1=1.6, k2=2.0;
148 
149  xx=x[0];
150  sigma0 = par[2];
151  A0 = 2*NEvents/(2+2*par[0]+par[0]*par[0]+par[0]*par[0]*par[0]/3);
152  r0 = ((xx-par[1])/sigma0);
153  C0 = 1/(sigma0* TMath::Exp(-k0*k0/2)/k0 +
154  sigma0*sqrt(2*3.14159)*0.5*(1+TMath::Erf(k0/1.41421)));
155  //sum = 1/(sqrt(2*3.14159)*par[2])*A0*TMath::Exp(-0.5*r0*r0);
156  if(r0 < k0) sum = C0*A0*TMath::Exp(-0.5*r0*r0);
157  else sum = C0*A0*TMath::Exp(0.5*k0*k0-k0*r0);
158 
159  mean1 = par[1]+par[3];
160  //sigma1 = par[4];
161  sigma1 = 1.547+0.125*par[3]+0.004042*par[3]*par[3];
162  sigma1 = (sigma1+(9.1347e-3+3.845e-2*par[3])*par[4]*2.0)*par[2];
163  A1 = A0*par[0];
164  C1 = 1/(sigma1* TMath::Exp(-k1*k1/2)/k1 +
165  sigma1*sqrt(2*3.14159)*0.5*(1+TMath::Erf(k1/1.41421)));
166  r1 = ((xx-mean1)/sigma1);
167  if(r1 < k1) sum += C1*A1*TMath::Exp(-0.5*r1*r1);
168  else sum += C1*A1*TMath::Exp(0.5*k1*k1-k1*r1);
169 
170  mean2 = 2*par[3]+par[1];
171  sigma2 = sqrt(2*sigma1*sigma1 - par[2]*par[2]);
172  //A2 = A0*par[5]*par[0]*par[0]/2;
173  A2 = A0*par[0]*par[0]/2;
174  C2 = 1/(sigma2* TMath::Exp(-k2*k2/2)/k2 +
175  sigma2*sqrt(2*3.14159)*0.5*(1+TMath::Erf(k2/1.41421)));
176  r2 = ((xx-mean2)/sigma2);
177  if(r2 < k2) sum += C2*A2*TMath::Exp(-0.5*r2*r2);
178  else sum += C2*A2*TMath::Exp(0.5*k2*k2-k2*r2);
179 
180  mean3 = 3*par[3]+par[1];
181  sigma3 = sqrt(3*sigma1*sigma1 - 2*par[2]*par[2]);
182  A3 = A0*par[0]*par[0]*par[0]/6;
183  C3 = 1/(sigma3*sqrt(2*3.14159));
184  r3 = ((xx-mean3)/sigma3);
185  sum += C3*A3*TMath::Exp(-0.5*r3*r3);
186 
187  return sum;
188 }
189 
190 void HFLightCalRand::endJob(void)
191 {
192  Double_t mean,rms,meanped,rmsped,npevar;
193  Double_t par[5],parm[5],dspe=0,dnpe;
194  Int_t tsmax,intspe;
195  fprintf(tFile,"#RunN %d Events processed %d",runNumb,EventN);
196  std::cout<<"endJob: histos processing..."<<std::endl;
197  std::cout<<"RunN= "<<runNumb<<" Events processed= "<<EventN<<std::endl;
198 
199  for (int i=0;i<26;i++) for (int j=0;j<36;j++) for (int k=0;k<2;k++) {
200  if (i>10 && i<13 && j%2==0) continue;
201  if (i>23 && j%2==0) continue;
202  meanped=rmsped=mean=rms=0;
203  if (hsp[i][j][k]->Integral()>0) {
204  HistSpec(hped[i][j][k],meanped,rmsped);
205  HistSpec(hsp[i][j][k],mean,rms);
206  if (hspe[i][j][k]->Integral()>hsp[i][j][k]->Integral()*0.9 || mean<100) {
207  HistSpec(hspe[i][j][k],mean,rms);
208  }
209  hsignalmean->Fill(mean);
210  hsignalrms->Fill(rms);
211  hpedmean->Fill(meanped);
212  hpedrms->Fill(rmsped);
213  }
214  }
215 
216  meanped=hpedmean->GetMean();
217  rmsped=hpedrms->GetMean();
218  mean=hsignalmean->GetMean();
219  rms=hsignalrms->GetMean();
220  fprintf(tFile," MeanInput=<%.2f> [linADCcount] RMS=<%.2f>\n",mean,rms);
221  fprintf(tFile,"#eta/phi/depth sum4maxTS RMS ~N_PE sum4lowTS RMS maxTS SPE +/- Err Comment\n");
222  TF1* fPed = new TF1("fPed","gaus",0,120);
223  fPed->SetNpx(200);
224  TF1 *fTot = new TF1("fTot",Fit3Peak ,0,200,5);
225  fTot->SetNpx(800);
226  for (int i=0;i<26;i++) for (int j=0;j<36;j++) for (int k=0;k<2;k++) {
227  if (i>10 && i<13 && j%2==0) continue;
228  if (i>23 && j%2==0) continue;
229  HistSpec(hped[i][j][k],meanped,rmsped);
230  HistSpec(hsp[i][j][k],mean,rms);
231  par[3]=0;
232  if (hspe[i][j][k]->Integral()>hsp[i][j][k]->Integral()*0.9 || mean<100) {
233  HistSpec(hspe[i][j][k],mean,rms);
234  if (hspe[i][j][k]->Integral(1,(int) (meanped+3*rmsped+12))/NEvents>0.1) {
235  //if (hspe[i][j][k]->Integral()>100 && mean-meanped<100) {
236  if (mean+rms*3-meanped-rmsped*3>1 && rmsped>0) { // SPE fit if low intensity>0
237  par[1] = meanped;
238  par[2] = rmsped;
239  par[0] = hped[i][j][k]->GetMaximum();
240  fPed->SetParameters(par);
241  hped[i][j][k]->Fit(fPed,"BQ0");
242  fPed->GetParameters(&par[0]);
243  hped[i][j][k]->Fit(fPed,"B0Q","",par[1]-par[2]*3,par[1]+par[2]*3);
244  fPed->GetParameters(par);
245  hped[i][j][k]->Fit(fPed,"BLIQ","",par[1]-par[2]*3,par[1]+par[2]*3);
246  fPed->GetParameters(&par[0]);
247  fPed->GetParameters(&parm[0]);
248  NEvents = (int) hspe[i][j][k]->Integral();
249  par[0]=0.1;
250  par[3]=10;
251  par[4]=6;
252  par[2]=par[2]/0.93;
253  fTot->SetParameters(par);
254  fTot->SetParLimits(0,0,2);
255  //fTot->FixParameter(1,par[1]);
256  //fTot->FixParameter(2,par[2]);
257  fTot->SetParLimits(1,par[1]-4,par[1]+4);
258  fTot->SetParLimits(2,par[2]*0.9,par[2]);
259  fTot->SetParLimits(3,1.2,50);
260  //fTot->SetParLimits(4,par[2]*1.1,100);
261  par[4]=0;
262  fTot->SetParLimits(4,-1.64,1.64);
263  hspe[i][j][k]->Fit(fTot,"BLEQ","");
264  fTot->GetParameters(par);
265  dspe=fTot->GetParError(3);
266  dnpe=fTot->GetParError(0);
267  if (par[3]<1.21 || dnpe>par[0]) par[3]=-1;
268  else if (par[0]>1.96 || par[3]>49) par[3]=0;
269  else {
270  hspes->Fill(par[3]);
271  }
272  }
273  }
274  }
275 
276  // NPE
277  npevar=0;
278  if (par[3]>0) npevar=par[0]; // NPE from SPE fit
279  else { // NPE from high intensity signal
280  if (hspe[i][j][k]->Integral()>hsp[i][j][k]->Integral()*0.98) {
281  HistSpec(hspe[i][j][k],mean,rms,3);
282  }
283  else {
284  HistSpec(hsp[i][j][k],mean,rms,3);
285  }
286  if (rmsped>0) {
287  if ((rms*rms-rmsped*rmsped/0.8836)>1 && mean>meanped+2) {
288  npevar=(mean-meanped-2)*(mean-meanped-2)/(rms*rms-rmsped*rmsped/0.8836);
289  }
290  else if (mean<100) {
291  intspe=int(hspe[i][j][k]->Integral());
292  hspe[i][j][k]->SetAxisRange(meanped+2+rmsped*4,300);
293  npevar=hspe[i][j][k]->Integral()/intspe;
294  if (npevar>0.01) npevar=-1;
295  else npevar=0;
296  hspe[i][j][k]->SetAxisRange(-20,300);
297  }
298  }
299  }
300  if (npevar>5.0e-5) hnpevar->Fill(npevar);
301 
302  if (i<13) {
303  hsignalmapP->Fill(i+28.6+k/2.0,j*2+1,mean-meanped-1.8);
304  hsignalRMSmapP->Fill(i+28.6+k/2.0,j*2+1,rms);
305  if (npevar>0) hnpemapP->Fill(i+28.6+k/2.0,j*2+1,npevar);
306  fprintf(tFile,"%3d%4d%5d %11.2f%8.2f",i+29,j*2+1,k+1,mean,rms);
307  }
308  else {
309  fprintf(tFile,"%3d%4d%5d %11.2f%8.2f",13-i-29,j*2+1,k+1,mean,rms);
310  hsignalmapM->Fill(13-i-28.6-k/2.0,j*2+1,mean-meanped-1.8);
311  hsignalRMSmapM->Fill(13-i-28.6-k/2.0,j*2+1,rms);
312  if (npevar>0) hnpemapM->Fill(13-i-28.6-k/2.0,j*2+1,npevar);
313  }
314  if (npevar>0) fprintf(tFile," %9.4f",npevar);
315  else fprintf(tFile," 0 ");
316  fprintf(tFile," %8.2f%8.2f",meanped,rmsped);
317  tsmax=hts[i][j][k]->GetMaximumBin()-1;
318  fprintf(tFile," %4d",tsmax);
319  if (par[3]>0 && par[3]<99) fprintf(tFile,"%8.2f%7.2f",par[3]+1,dspe);
320  else if (npevar>0) fprintf(tFile,"%8.2f 0 ",(mean-meanped-1)/npevar);
321  else fprintf(tFile," 0 0 ");
322 
323  // Diagnostics
324  fprintf(tFile," ");
325  if (hsp[i][j][k]->GetEntries()<=0) fprintf(tFile,"NoSignal\n");
326  else if (hsp[i][j][k]->GetEntries()<=10) fprintf(tFile,"Nev<10\n");
327  else {
328  if (hsp[i][j][k]->Integral()<=10 || mean>12000) fprintf(tFile,"SignalOffRange\n");
329  else {
330  if (hsp[i][j][k]->Integral()<100) fprintf(tFile,"Nev<100/");
331  if (npevar>0 && par[3]>0 && (npevar*NEvents<10 || npevar<0.001))
332  fprintf(tFile,"LowSignal/");
333  else if (npevar==0 && fabs(mean-meanped)<3) fprintf(tFile,"LowSignal/");
334  if (par[3]<0) fprintf(tFile,"BadFit/");
335  else if (par[3]==0) fprintf(tFile,"NoSPEFit/");
336  else if (par[3]>0 && npevar>1) fprintf(tFile,"NPE>1/");
337  if (npevar<0) fprintf(tFile,"Problem/");
338  if (mean<2) fprintf(tFile,"LowMean/");
339  if (rms<0.5) fprintf(tFile,"LowRMS/");
340  if (meanped<-1) fprintf(tFile,"LowPed/");
341  else if (meanped>25) fprintf(tFile,"HighPed/");
342  if (rmsped<0.5 && rmsped>0) fprintf(tFile,"NarrowPed/");
343  else if (rmsped>10) fprintf(tFile,"WidePed/");
344  if (hped[i][j][k]->GetBinContent(201)>10) fprintf(tFile,"PedOffRange");
345  fprintf(tFile,"-\n");
346  }
347  }
348  }
349 
350  for (int i=0;i<8;i++) for (int j=0;j<3;j++) {
351  HistSpec(hpedpin[i][j],meanped,rmsped);
352  HistSpec(hsppin[i][j],mean,rms);
353  if (hspepin[i][j]->Integral()>hsppin[i][j]->Integral()*0.9 || mean<100) {
354  HistSpec(hspepin[i][j],mean,rms);
355  }
356  if (i<4) fprintf(tFile," PIN%d +Q%d %12.2f %6.2f",j+1,i+1,mean,rms);
357  else fprintf(tFile," PIN%d -Q%d %12.2f %6.2f",j+1,i-3,mean,rms);
358  fprintf(tFile," %15.2f %6.2f",meanped,rmsped);
359  tsmax=htspin[i][j]->GetMaximumBin()-1;
360  fprintf(tFile," %4d\n",tsmax);
361  }
362 
363  mFile->Write();
364  mFile->Close();
365  fclose(tFile);
366  std::cout<<std::endl<<" --endJob-- done"<<std::endl;
367  return;
368 }
369 
370 void HFLightCalRand::analyze(const edm::Event& fEvent, const edm::EventSetup& fSetup) {
371 
372  // event ID
373  edm::EventID eventId = fEvent.id();
374  int runNumber = eventId.run ();
375  int eventNumber = eventId.event ();
376  if (runNumb==0) runNumb=runNumber;
377  EventN++;
378  if (verbose) std::cout << "========================================="<<std::endl
379  << "run/event: "<<runNumber<<'/'<<eventNumber<<std::endl;
380 
381  Double_t buf[20];
382  Double_t maxADC,signal,ped=0,meant;
383  Int_t maxisample=0,i1=3,i2=6;
384 
385  // HF PIN-diodes
386  edm::Handle<HcalCalibDigiCollection> calib;
387  fEvent.getByType(calib);
388  if (verbose) std::cout<<"Analysis-> total CAL digis= "<<calib->size()<<std::endl;
389  /* COMMENTED OUT by J. Mans (7-28-2008) as major changes needed with new Calib DetId
390  re-commented out by R.Ofierzynski (11.Nov.2008) - to be able to provide a consistent code for CMSSW_3_0_0_pre3:
391  major changes are needed for the new Calib DetId which does not have the old methods any more
392 
393  for (unsigned j = 0; j < calib->size (); ++j) {
394  const HcalCalibDataFrame digi = (*calib)[j];
395  HcalElectronicsId elecId = digi.elecId();
396  HcalCalibDetId calibId = digi.id();
397  if (verbose) std::cout<<calibId.sectorString().c_str()<<" "<<calibId.rbx()<<" "<<elecId.fiberChanId()<<std::endl;
398  int isector = calibId.rbx()-1;
399  int ipin = elecId.fiberChanId();
400  int iside = -1;
401  if (calibId.sectorString() == "HFP") iside = 0;
402  else if (calibId.sectorString() == "HFM") iside = 4;
403 
404  if (iside != -1) {
405  maxADC=-99;
406  for (int isample = 0; isample < digi.size(); ++isample) {
407  int adc = digi[isample].adc();
408  int capid = digi[isample].capid ();
409  double linear_ADC = digi[isample].nominal_fC();
410  if (verbose) std::cout<<"PIN linear_ADC = "<<linear_ADC<<std::endl;
411  htspin[isector+iside][ipin]->Fill(isample,linear_ADC);
412  buf[isample]=linear_ADC;
413  if (maxADC<linear_ADC) {
414  maxADC=linear_ADC;
415  maxisample=isample;
416  }
417  }
418  i1=maxisample-1;
419  i2=maxisample+2;
420  if (i1<0) {i1=0;i2=3;}
421  else if (i2>9) {i1=6;i2=9;}
422  if (i1==0) ped=buf[8]+buf[9]+buf[6]+buf[7];
423  else if (i1==1) ped=buf[8]+buf[9]+buf[6]+buf[7];
424  else if (i1==2) ped=buf[0]+buf[1]+buf[6]+buf[7];
425  else if (i1==3) ped=buf[0]+buf[1]+buf[2]+buf[7];
426  else if (i1>=4) ped=buf[0]+buf[1]+buf[2]+buf[3];
427  signal=0;
428  for (ii=0;ii<4;ii++) signal+=TMath::Max(buf[ii+i1],ped/4.0);
429  hsppin[isector+iside][ipin]->Fill(signal);
430  hspepin[isector+iside][ipin]->Fill(signal);
431  hpedpin[isector+iside][ipin]->Fill(ped);
432 
433  // Mean PIN signal time estimation
434  ped=ped/4;
435  meant=0;
436  for (ii=0;ii<4;ii++) meant+=(TMath::Max(buf[ii+i1],ped)-ped)*(ii+i1);
437  if (signal-ped*4>0) meant/=(signal-ped*4);
438  else meant=i1+1;
439  htsmpin[isector+iside][ipin]->Fill(meant);
440  }
441  }
442  */
443  // HF
444  edm::Handle<HFDigiCollection> hf_digi;
445  fEvent.getByType(hf_digi);
446  if (verbose) std::cout<<"Analysis-> total HF digis= "<<hf_digi->size()<<std::endl;
447 
448  for (unsigned ihit = 0; ihit < hf_digi->size (); ++ihit) {
449  const HFDataFrame& frame = (*hf_digi)[ihit];
450  HcalDetId detId = frame.id();
451  int ieta = detId.ieta();
452  int iphi = detId.iphi();
453  int depth = detId.depth();
454  if (verbose) std::cout <<"HF digi # " <<ihit<<": eta/phi/depth: "
455  <<ieta<<'/'<<iphi<<'/'<< depth << std::endl;
456 
457  if (ieta>0) ieta = ieta-29;
458  else ieta = 13-ieta-29;
459 
460  maxADC=-99;
461  for (int isample = 0; isample < frame.size(); ++isample) {
462  int adc = frame[isample].adc();
463  int capid = frame[isample].capid ();
464  double linear_ADC = frame[isample].nominal_fC();
465  double nominal_fC = detId.subdet () == HcalForward ? 2.6 * linear_ADC : linear_ADC;
466 
467  if (verbose) std::cout << "Analysis-> HF sample # " << isample
468  << ", capid=" << capid
469  << ": ADC=" << adc
470  << ", linearized ADC=" << linear_ADC
471  << ", nominal fC=" << nominal_fC << std::endl;
472 
473  hts[ieta][(iphi-1)/2][depth-1]->Fill(isample,linear_ADC);
474  buf[isample]=linear_ADC;
475  if (maxADC<linear_ADC) {
476  maxADC=linear_ADC;
477  maxisample=isample;
478  }
479  }
480 
481  // Signal is four capIDs around maxTS, Pedestal is four capID off the signal
482  maxADC=-99;
483  for (int ibf=0; ibf<9; ibf++) {
484  Double_t sumbuf=0;
485  for (int jbf=0; jbf<2; jbf++) {
486  sumbuf += buf[ibf+jbf];
487  if (ibf+jbf<2) sumbuf -= (buf[ibf+jbf+4]+buf[ibf+jbf+8])/2.0;
488  else if (ibf+jbf<4) sumbuf -= buf[ibf+jbf+4];
489  else if (ibf+jbf<6) sumbuf -= (buf[ibf+jbf-4]+buf[ibf+jbf+4])/2.0;
490  else if (ibf+jbf<8) sumbuf -= buf[ibf+jbf-4];
491  else if (ibf+jbf<10) sumbuf -= (buf[ibf+jbf-4]+buf[ibf+jbf-8])/2.0;
492  }
493  if (sumbuf>maxADC) {
494  maxADC=sumbuf;
495  if (buf[ibf]>buf[ibf+1]) maxisample=ibf;
496  else maxisample=ibf+1;
497  //maxisample=ibf;
498  }
499  }
500  htmax->Fill(maxisample);
501  i1=maxisample-1;
502  i2=maxisample+2;
503  if (i1<0) {i1=0;i2=3;}
504  else if (i2>9) {i1=6;i2=9;}
505  else if (i2<9) {
506  if (buf[i1]<buf[i2+1]) {i1=i1+1;i2=i2+1;}
507  }
508  signal=buf[i1]+buf[i1+1]+buf[i1+2]+buf[i1+3];
509  hsp[ieta][(iphi-1)/2][depth-1]->Fill(signal);
510  hspe[ieta][(iphi-1)/2][depth-1]->Fill(signal);
511 
512  if (i1==0) ped=(buf[4]+buf[8])/2.0+(buf[5]+buf[9])/2.0+buf[6]+buf[7];
513  else if (i1==1) ped=(buf[0]+buf[8])/2.0+(buf[5]+buf[9])/2.0+buf[6]+buf[7];
514  else if (i1==2) ped=(buf[0]+buf[8])/2.0+(buf[1]+buf[9])/2.0+buf[6]+buf[7];
515  else if (i1==3) ped=(buf[0]+buf[8])/2.0+(buf[1]+buf[9])/2.0+buf[2]+buf[7];
516  else if (i1==4) ped=(buf[0]+buf[8])/2.0+(buf[1]+buf[9])/2.0+buf[2]+buf[3];
517  else if (i1==5) ped=(buf[0]+buf[4])/2.0+(buf[1]+buf[9])/2.0+buf[2]+buf[3];
518  else if (i1==6) ped=(buf[0]+buf[4])/2.0+(buf[1]+buf[5])/2.0+buf[2]+buf[3];
519  /*
520  if (i1==0) ped=buf[8]+buf[9]+buf[6]+buf[7];
521  else if (i1==1) ped=buf[0]+buf[9]+buf[6]+buf[7];
522  else if (i1==2) ped=buf[0]+buf[1]+buf[6]+buf[7];
523  else if (i1==3) ped=buf[0]+buf[1]+buf[2]+buf[7];
524  else if (i1>=4) ped=buf[0]+buf[1]+buf[2]+buf[3];
525  */
526  hped[ieta][(iphi-1)/2][depth-1]->Fill(ped);
527 
528  // Mean signal time estimation
529  ped=ped/4;
530  meant=(buf[i1]-ped)*i1+(buf[i1+1]-ped)*(i1+1)+(buf[i1+2]-ped)*(i1+2)+(buf[i1+3]-ped)*(i1+3);
531  meant /= (buf[i1]-ped)+(buf[i1+1]-ped)+(buf[i1+2]-ped)+(buf[i1+3]-ped);
532  htmean->Fill(meant);
533  htsm[ieta][(iphi-1)/2][depth-1]->Fill(meant);
534  }
535 }
536 
ecalMGPA::adc
int adc(sample_type sample)
get the ADC sample (12 bits)
Definition: EcalMGPASample.h:11
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:42
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:44
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
Fit3Peak
Double_t Fit3Peak(Double_t *x, Double_t *par)
Definition: HFLightCalRand.cc:142
HcalDetId::subdet
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:32
NULL
#define NULL
Definition: scimark2.h:8
edm::Event::getByType
bool getByType(Handle< PROD > &result) const
Definition: Event.h:398
HFLightCalRand::beginJob
virtual void beginJob()
Definition: HFLightCalRand.cc:53
runNumb
Int_t runNumb
Definition: HFLightCalRand.cc:36
HcalDetId::depth
int depth() const
get the tower depth
Definition: HcalDetId.h:42
calib
MVATrainerComputer * calib
Definition: MVATrainer.cc:64
EventN
Int_t EventN
Definition: HFLightCalRand.cc:36
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
HcalDetId::ieta
int ieta() const
get the cell ieta
Definition: HcalDetId.h:38
j
int j
Definition: DBlmapReader.cc:9
HFLightCalRand::HFLightCalRand
HFLightCalRand(const edm::ParameterSet &fConfiguration)
Definition: HFLightCalRand.cc:43
HFLightCalRand::endJob
virtual void endJob(void)
Definition: HFLightCalRand.cc:190
gen::k
int k[5][pyjets_maxn]
Definition: Cascade2Hadronizer.cc:81
HcalDetId::iphi
int iphi() const
get the cell iphi
Definition: HcalDetId.h:40
HFLightCalRand::analyze
virtual void analyze(const edm::Event &fEvent, const edm::EventSetup &fSetup)
Definition: HFLightCalRand.cc:370
HFDataFrame::size
int size() const
total number of samples in the digi
Definition: HFDataFrame.h:26
HFLightCalRand::~HFLightCalRand
virtual ~HFLightCalRand()
Definition: HFLightCalRand.cc:49
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
NEvents
Int_t NEvents
Definition: HFLightCalRand.cc:36
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
HFDataFrame::id
const HcalDetId & id() const
Definition: HFDataFrame.h:22
x
Definition: DDAxes.h:10
HistSpec
void HistSpec(TH1F *hist, Double_t &mean, Double_t &rms, Double_t range=4)
Definition: HFLightCalRand.cc:126
Generated for CMSSW Reference Manual by doxygen 1.8.5