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