Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
hcalUtil.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file hcalUtil.C
1 #include <iomanip>
2 #include <algorithm>
3 #include <numeric>
4 #include <iostream>
5 #include <fstream>
6 #include <cmath>
7 #include <ctime>
8 #include <cstdlib>
9 #include <stdint.h>
10 #include <stddef.h>
11 #include <stdlib.h>
12 #include <string.h>
13 
14 
15 // #include <TROOT.h>
16 // #include <TFile.h>
17 // #include <TTree.h>
18 // #include <TGMsgBox.h>
19 // #include <TCanvas.h>
20 // #include <TGraph.h>
21 // #include <TString.h>
22 // #include <TText.h>
23 // #include <TApplication.h>
24 // #include <TStyle.h>
25 // #include <TSystem.h>
26 // #include <TGClient.h>
27 // #include <TGWindow.h>
28 // #include <TH1.h>
29 // #include <TH1F.h>
30 // #include <TH1D.h>
31 // #include <TProfile.h>
32 // #include <TF1.h>
33 // #include <TH2.h>
34 // #include <TH3.h>
35 // #include <TMath.h>
36 // #include <TSpectrum.h>
37 // #include "TVirtualFitter.h"
38 // #include <TSystemDirectory.h>
39 // #include <TSystemFile.h>
40 
41 //#ifndef __CINT__
42 
43 #include <fileEventiterator.h>
44 #include <Event.h>
45 #include <EventTypes.h>
46 #include <PHTimeStamp.h>
47 #include <packetConstants.h>
48 #include <packet.h>
49 #include <packet_hbd_fpgashort.h>
50 
51 //#endif
52 
53 
54 #include "hcalUtil.h"
55 
56 
57 hcalUtil * hcalUtil :: single = 0;
58 hLabHelper * hLabHelper :: single = 0;
59 hcalHelper * hcalHelper :: single = 0;
60 hcalTree * hcalTree :: single = 0;
61 tileHelper * tileHelper :: single = 0;
62 tileTree * tileTree :: single = 0;
63 
64 using namespace std;
65 using namespace TMath;
66 
67 TCanvas * fitFunc;
68 
69 TCanvas * scCalibDisplay;
70 
71 
72 // ************************************************************************************
73 
74 
75 //#endif /* _CINT_ */
76 
77 // cosmic runs of interest
78 static const int cosmicruns = 12;
79 int cosmicRN[cosmicruns] = {212, 215, 216, 218, 223, 337, 338, 339, 405, 406, 229, 240};
80 
81 // LED calibration and PEDESTALS
82 static const int pedruns = 1;
83 int pedcalibRN[pedruns] = {646};
84 static const int ledruns = 1;
85 int ledcalibRN[ledruns] = {647};
86 
87 // -----------------------------------------------------------------------------------------------
88 
89 // Unipolar pulse (0.3164*pow(x,4)*exp(-x*1.5) - Unity integral, peak ~1/3 of integral)
90 Double_t par0[] = { 1., (Double_t)NSAMPLES/2, 3., 2., PEDESTAL, 0. };
91 Double_t par0Max[] = { -100., (Double_t)NSAMPLES, 4., 5., 50000., 0.3};
92 //Double_t par0Min[] = { 1., RISETIME-1, 3., 1., 1950., -0.2};
93 Double_t par0Min[] = {-50000., 0., 2., 0., -50000., -0.3};
94 
95 
96 TString runId;
97 
98 // The next few functions are to help processing data from HCal Lab recorded via HBD electronics
99 
100 // **************************************************************************
101 
102 Double_t ps( Double_t *chadc, Double_t *par)
103 {
104  Double_t x = chadc[0];
105 
106  Double_t mu = par[0];
107  Double_t g = par[1];
108  Double_t sigma = par[2];
109  Double_t norm = par[3];
110 
111  Double_t sum = 0.0;
112  for (Double_t n = 1.0; n < 100.0; n++ ) {
113  sum += TMath::Poisson(n,mu)*TMath::Gaus(x,n*g,sigma,kTRUE);
114  }
115 
116  return norm*sum;
117 }
118 // **************************************************************************
119 Int_t calibPeaks(0), calibMode (0); Double_t spCalib(10.);
120 // **************************************************************************
121 // mode 0 - peaks uncorrelated (unconstrained)
122 // mode 1 - peaks equidistant
123 // assume background is linear what is most likely not true ....
124 // par[0] background constant
125 // [1] background slope
126 // [2] background quadratic term
127 // assume the signal is the sequence of equally spaced gaussian peaks
128 // par[3] peak to peak separation
129 // [4] normalization constant for the first peak (pedestal)
130 // [5] location of the first (pedestal) peak
131 // [6] rms of the first peak - at this stage we will allow it to vary peak-to-peak
132 // then 3 parameters per signal peak (normalization, mean and RMS)
133 Double_t sipmPeaks(Double_t *x, Double_t *par) {
134  Double_t result = par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
135  // add pedestal peak first
136  result += par[4]*TMath::Gaus(x[0],par[5],par[6]);
137  for (Int_t p=1; p<calibPeaks; p++) {
138  Double_t norm = par[4 + p*3];
139  Double_t mean = par[5 + p*3];
140  if(calibMode) {
141  mean = (Int_t)((mean+spCalib/2)/spCalib)*par[3];
142  }
143  Double_t sigma = par[6 + p*3];
144  result += norm*TMath::Gaus(x[0],mean,sigma);
145  }
146  return result;
147  ;}
148 
149 
150 // **************************************************************************
151 // do single cell calibration using low end data stored in rv_# and fm_# channel histograms
152 // the assumption - we have beautiful picture with at least 20 peaks !!!!!!
153 // mode 0 - peaks uncorrelated (unconstrained)
154 // mode 1 - peaks equidistant
155 Int_t sipmCalib(TString & hName, Int_t mode){
156  TH1F * data = (TH1F*)(gROOT->FindObject(hName));
157  if(!data){
158  cout<<"spmCalib Data Histogramm "<<hName<<" not found"<<endl;
159  return 0;
160  }
161  cout<<hName<<" calibration mode "<<mode<<endl;
162  calibMode = mode;
163  Int_t maxPeaks(20); // maximum number of peaks we have chance to resolve
164  Double_t range(150.); // range of amplitudes for peak search
165  // second degree polinomial background
166  Double_t bckgr0(100.), bckgr1(-1./range), bckgr2(0.);
167  Double_t pNorm(1000.), pLoc(0.), pRMS(1.);
168  Double_t par[4+maxPeaks*3]; // backgrownd, pedestal and nPeaks of gaussian peaks with constant peak-to-peak separation
169  // initialization
170  par[0] = bckgr0;
171  par[1] = bckgr1;
172  par[2] = bckgr2;
173  par[3] = spCalib; // peak-to-peak separation
174  par[4] = pNorm; // normalization for pedestal peak
175  par[5] = pLoc; // pedestal mean
176  par[6] = pRMS; // pedestal rms
177  // all other peaks are initially positioned equidistantly
178  for (int p=1; p<= maxPeaks; p++) {
179  par[4 + p*3] = pNorm; // norm of the peak p
180  par[5 + p*3] = par[3]*p; // mean of the peak p
181  par[6 + p*3] = pRMS; // rms of the peak p
182  }
183  TH1F * dataCl = (TH1F*)(data->Clone("dataCl"));
184  dataCl -> SetAxisRange(0., range);
185  TH1F * dataCl2 = (TH1F*)data->Clone("dataCl2");
186  dataCl2 -> SetAxisRange(0., range);
187  TString cTitle = "Peaks from "; cTitle += hName;
188  TCanvas * peaks = new TCanvas("peaks",cTitle,10,10,1000,900);
189  peaks -> Divide(1,2);
190  peaks -> cd(1);
191  dataCl -> Draw();
192 
193  // Use TSpectrum to find the peak candidates
194  TSpectrum *s = new TSpectrum(maxPeaks);
195  s->SetResolution(2.);
196  Int_t nfound = s->Search(dataCl, 1.5, "nobackground ", 0.10);
197  printf("TSpectrum: Found %d candidate peaks to fit\n",nfound);
198 
199  //Estimate background using TSpectrum::Background
200  TH1 *hb = s->Background(dataCl, 20, "same");
201  if (hb) peaks->Update();
202  // if only one or no peaks found (only pedestals) - we are done
203  if(nfound<=1) {
204  cout<<"sipmCalib One or No peaks were found - EXITING"<<endl;
205  calibPeaks = nfound;
206 
207  } else {
208  peaks -> cd(2);
209 
210  // function to estimate linear background
211  // TF1 *fline = new TF1("fline","pol1",0.,range);
212  // function to estimate quadratic background
213  TF1 *fline = new TF1("fline","pol2",0.,range);
214 
215  dataCl2 -> Fit("fline","qn");
216  //Loop on all found peaks. Eliminate peaks at the background level
217  par[0] = fline->GetParameter(0);
218  par[1] = fline->GetParameter(1);
219  par[2] = fline->GetParameter(2);
220  // cout<<"Line par "<<par[0]<<" "<<par[1]<<" "<<par[2]<<endl;
221  Float_t * xpeaks = s->GetPositionX();
222  // peak ordering
223  Int_t pindex[nfound];
224  TMath::Sort(nfound, xpeaks, pindex, kFALSE);
225  // peaks now ordered, let's now find the peak to peak separation
226  Double_t sep[nfound-1];
227  // Int_t sepind[nfound-1];
228  // pedestal peak in my definition == what is left of a sygnal after pedestal subtraction - it is definitely shifted
229  for (int ip = 1; ip<nfound; ip++){
230  sep[ip-1] = xpeaks[pindex[ip]];
231  if(ip>1) sep[ip-1] -= xpeaks[pindex[ip-1]];
232  cout<<"Separation "<<ip<<" val "<<sep[ip-1]<<endl;
233  }
234  // TMath::Sort(nfound-1, sep, sepind, kFALSE);
235  // zero approximation for peak0-to-peak separation
236  Int_t cont(0); Double_t avsep(0.);
237  for (int ip = 0; ip < nfound-1; ip++){
238  if(ip==0) {
239  // avsep = sep[sepind[ip]];
240  avsep = sep[ip];
241  cont ++;
242  } else {
243  // if(abs((avsep- sep[sepind[ip]])/(avsep+ sep[sepind[ip]]))<0.2) {
244  if(abs((avsep- sep[ip])/(avsep+ sep[ip]))<0.2) {
245  avsep = (avsep*cont+sep[ip])/(cont+1);
246  cont++;
247  }
248  }
249  cout<<"ip "<<ip<<" cont "<<cont<<" avsep "<<avsep<<" var "<<(ip!=0? (avsep- sep[ip])/(avsep+ sep[ip]) : 0.)<<endl;
250  }
251  cout<<"Zero approximation to sSingle Pixel Calib = "<<avsep<<" Contributors "<<cont<<endl;
252  spCalib = avsep;
253  par[3] = avsep;
254  calibPeaks = 0;
255  for (int p=0; p<nfound; p++) {
256  Double_t xp = xpeaks[pindex[p]];
257  Int_t bin = dataCl2->GetXaxis()->FindBin(xp);
258  Double_t yp = dataCl2->GetBinContent(bin);
259  cout<<"Testing Peak (LFit) "<<p<<" index "<<pindex[p]<<" xp "<<xp<<" bin "<<bin<<" yP "<<yp<<" BCKG "<<fline->Eval(xp)<<endl;
260  // extract data for this bin from background shape
261  Double_t bckgr = hb->GetBinContent(bin);
262  cout<<"Testing Peak (TSp ) "<<p<<" xp "<<xp<<" bin "<<bin<<" Bcgrd "<<bckgr<<" Signal "<<yp-bckgr<<endl;
263  if (yp-TMath::Sqrt(yp) < bckgr) continue;
264  par[4 + calibPeaks*3] = yp;
265  par[5 + calibPeaks*3] = xp;
266  par[6 + calibPeaks*3] = (p==0? pRMS : 2.);
267  calibPeaks++;
268  }
269  }
270  printf("Found %d useful peaks to fit\n",calibPeaks);
271  printf("Now fitting: Be patient\n");
272  TF1 *fit = new TF1("fit",&sipmPeaks,0,range,4+calibPeaks*3); // first peak is the pedestal
273  //we may have more than the default 25 parameters
274  TVirtualFitter::Fitter(dataCl2,10+calibPeaks*3);
275  fit->SetParameters(par);
276  // if mode!=0 we fix all location parameters to what came from TSpectrum
277  if(mode) for(int ip=0; ip<calibPeaks; ip++) {
278  fit->FixParameter(5+ip*3, par[5+ip*3]);
279  fit->SetParLimits(3, 0.5*par[3], 1.5*par[3]);
280  } else fit->FixParameter(3, par[3]);
281  // in all cases we set max RMS equal to pps/2.
282  for(int ip=0; ip<calibPeaks; ip++) fit->SetParLimits(6+ip*3, 0.5, par[3]/2.);
283 
284  fit->SetNpx(1000);
285  dataCl2->Fit("fit");
286  fit->GetParameters(par);
287 
288  char p5[10], p6[10], p3[10];
289  sprintf(p5, "%.2f ", par[5]);
290  sprintf(p6, "%.2f ", par[6]);
291  sprintf(p3, "%.2f ", par[3]);
292 
293  TString ped = "Pedestal at "; ped += p5; ped += " RMS = "; ped += p6;
294  TString cal = "Pixel Calibtion "; cal += p3; cal += " ADC counts per pixel";
295  TText label;
296  Double_t ymax = data-> GetMaximum();
297  label.DrawText(30, ymax*0.8,ped);
298  label.DrawText(30, ymax*0.6,cal);
299  return nfound;
300 }
301 
302 // **************************************************************************
303 Double_t AsymToX(Double_t ax){
304  return 12.5*(1 + ax/0.19);
305 }
306 
307 
308 // **************************************************************************
309 
310 Double_t signalShape(Double_t *x, Double_t * par){
311  Double_t signal(0.), pedestal(0.);
312  // cout<<"signalShape x="<<x[0]<<" (0)="<<par[0]<<" (1)="<<par[1]<<" (2)="<<par[2]<<" (3)="<<par[3]<<" (4)="<<par[4]<<" (5)="<<par[5]<<endl;
313  // if(x[0]<par[1]) return PEDESTAL;
314  pedestal = par[4]+x[0]*par[5];
315  if(x[0]<par[1]) return pedestal;
316  // signal contribution
317  // cout<<"(0) "<<pow((x[0]-par[1]),par[2])<<endl;
318  // cout<<"(1) "<<exp(-(x[0]-par[1])*par[3])<<endl;
319  // Double_t a0 = par[0];
320  // Double_t a1 = pow((x[0]-par[1]),par[2]);
321  // Double_t a2 = exp(-(x[0]-par[1])*par[3]);
322  // Double_t a = a0*a1*a2;
323  // // cout<<a0<<" "<<a1<<" "<<a2<<" "<<a<<endl;
324  // signal = a;
325  // return a;
326  signal = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
327  // cout<<" X "<<x[0]<<" signal "<<signal<<" pedestal "<<pedestal<<endl;
328  return signal+pedestal;
329 }
330 
331 // **************************************************************************
332 
333 // Double_t shape0(Double_t *x, Double_t * par){
334 // return par[0]*pow(x[0],4.)*(exp(-x[0]*16./fN)); // fN is set to (float)NSAMPLES
335 // }
336 
337 // **************************************************************************
338 
339 Double_t erfunc(Double_t *x, Double_t * par){
340  // Double_t erfcont = par[3]*(TMath::Erf((x[0]-par[0])/par[1])+1.);
341  // Double_t expcont = TMath::Exp(par[2]*x[0]);
342  return par[3]*(TMath::Erf((x[0]-par[0])/par[1])+1.)*TMath::Exp(par[2]*x[0]);
343 }
344 // **************************************************************************
345 
346 Double_t erf_g(Double_t *x, Double_t * par){
347  float g = par[0]/(sqrt(6.2830)*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
348  return g+par[6]*(TMath::Erf((x[0]-par[3])/par[4])+1.)*TMath::Exp(par[5]*x[0]);
349 }
350 // **************************************************************************
351 
352 Double_t powerFun(Double_t *x, Double_t * par){
353  // return par[0]*pow((x[0]),par[1]*pow((par[2]-x[0]),par[3]);
354  return par[0]*pow(x[0],par[1])*pow((1-x[0]),par[2]);
355  //*pow((0.5-x[0]),par[2]);
356 }
357 
358 // **************************************************************************
359 
360 Double_t power_g(Double_t *x, Double_t * par){
361  float g = par[0]/(sqrt(6.2830)*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
362  float p = par[3]*pow(x[0],par[4])*pow((1-x[0]),par[5]);
363  //*pow((0.5-x[0]),par[5]);
364  return g+p;
365 }
366 
367 
368 
369 // ====================================================================================
370 // ************************************************************************************
371 // ====================================================================================
372 
373 
374 hcalUtil::hcalUtil(){
375  char * s = getenv("BASE_DIR");
376  TString basedir = (s==0||strlen(s)==0)? "/home/kistenev/workarea/" : s;
377 }
378 
379 //----------------------------------------------------------------------------
380 // void hcalUtil::setBaseDir(char * bDir){
381 // basedir.Resize(0);
382 // basedir += bDir;
383 // }
384 
385 //----------------------------------------------------------------------------
386 
387 float hcalUtil::xPeak(TF1* f, float x1, float x2){
388  // searches for peak position in the TF1 function - looks into derivatives
389  float y1 = f->Derivative(x1, 0, 0.002);
390  float y2 = f->Derivative(x2, 0, 0.002);
391  if(y1*y2>0) {
392  return (f->Eval(x1)>f->Eval(x2))? x1 : x2;
393  }
394  float xx1 = x1;
395  float xx2 = x2;
396  float dx = fabs(xx2-xx1);
397  float x = 0.;
398  while(y1*y2<0&&dx>0.001){
399  x = (xx1+xx2)/2.;
400  float y = f->Derivative(x, 0, 0.01);
401  if(y1*y>0) xx1 = x; else xx2 = x;
402  dx = fabs(xx2-xx1);
403  }
404  return x;
405 }
406 
407 //----------------------------------------------------------------------------
408 
409 float hcalUtil::width(TF1* f, float x1, float x2, float peak){
410  // find width of the function (erfunc only) !!!
411  float yhm = f->Eval(peak)/2;
412 
413  float xl = x1;
414  float xr = peak;
415 
416  float xm = (xr+xl)/2.;
417  float ym = f->Eval(xm);
418  // cout<<yhm<<" xl "<< xl<<" xr "<<xr<<" xm "<<xm<<" ym "<<ym<<endl;
419  while(abs(ym-yhm)>0.01*yhm && xr-xl>0.001){
420  if(ym>yhm) xr = xm; else xl = xm;
421  xm = (xr+xl)/2.;
422  ym = f->Eval(xm);
423  //cout<<"left xl "<<xl<<" xr "<<xr<<" xm "<<xm<<" ym "<<ym<<endl;
424  }
425  bool lok = (abs(ym-yhm)>0.1*yhm)? false : true;
426  float xml = xm;
427  xl = peak;
428  xr = x2;
429  xm = (xr+xl)/2.;
430  ym = f->Eval(xm);
431  while(abs(ym-yhm)>0.01*yhm&&xr-xl>0.001){
432  if(ym>yhm) xl = xm; else xr = xm;
433  xm = (xr+xl)/2.;
434  ym = f->Eval(xm);
435  //cout<<"right xl "<<xl<<" xr "<<xr<<" xm "<<xm<<" ym "<<ym<<endl;
436  }
437  float xmr = xm;
438  bool rok = (abs(ym-yhm)>0.1*yhm)? false : true;
439  float w = 0;
440  if(lok){
441  if(rok) w = xmr-xml; else w = 2*(peak-xml);
442  } else if (rok){
443  w = 2*(xmr-peak);
444  } else w = x2-x1;
445  return w/2.35;
446 }
447 
448 //----------------------------------------------------------------------------
449 
450 bool hcalUtil::shapeFit(TH1D * proj, Double_t * pE, float & maxE, Double_t * pG, bool fine, bool MIP){
451  TF1 * fE;
452  TF1 * fG;
453  int bMax = proj->GetMaximumBin();
454  float pMax = proj->GetBinCenter(bMax);
455  bool erfOK = true;
456  TAxis *x;
457  x = proj->GetXaxis();
458  double xmin = x->GetXmin();
459  double xmax = x->GetXmax();
460  pE[0] = pMax;
461  pE[1] = proj->GetRMS();
462  pE[2] = -1./pE[1];
463  pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
464  fE = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
465  if(fE) delete fE;
466  fE = new TF1("fE",&erfunc,xmin,xmax,4);
467  fE->SetParameters(pE);
468  proj->Fit("fE","NQL0");
469  fE->GetParameters(pE);
470  maxE = xPeak(fE,xmin,xmax);
471  // cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
472  if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
473  float scale = (fine? 1. : 1.5);
474  // float scale = (fine? 1. : 1);
475  pG[0] = proj->GetBinContent(bMax);
476  if(erfOK){
477  pG[1] = maxE;
478  pG[2] = width(fE,xmin,xmax,maxE);
479  } else {
480  pG[1] = proj->GetMean();
481  pG[2] = proj->GetRMS();
482  }
483  fG = (TF1*)(gROOT->FindObject("fG"));
484  if(fG) delete fG;
485  if(MIP) {
486  fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-0.3*scale),TMath::Min(xmax,pG[1]+0.075*scale));
487  fG->SetParameters(pG);
488  fG->SetParLimits(1,TMath::Max(xmin,pG[1]-0.1*scale),TMath::Min(xmax,pG[1]+0.05*scale));
489  } else {
490  fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-1.5*scale),TMath::Min(xmax,pG[1]+1.5*scale));
491  fG->SetParameters(pG);
492  fG->SetParLimits(1,TMath::Max(xmin,pG[1]-2.*scale),TMath::Min(xmax,pG[1]+2.*scale));
493  }
494  proj->Fit("fG","NLQR0");
495  fG->GetParameters(pG);
496  return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
497 }
498 
499 //----------------------------------------------------------------------------
500 
501 bool hcalUtil::emcFit(TH1D * proj, bool fine, bool sing, Double_t * pE, float & maxE, Double_t * pG, bool MIP){
502  TF1 * fE;
503  TF1 * fG;
504  int bMax = proj->GetMaximumBin();
505  float pMax = proj->GetBinCenter(bMax);
506  bool erfOK = true;
507  TAxis *x;
508  x = proj->GetXaxis();
509  double xmin = x->GetXmin();
510  double xmax = x->GetXmax();
511  pE[0] = pMax;
512  pE[1] = proj->GetRMS();
513  pE[2] = -1./pE[1];
514  pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
515  fE = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
516  if(fE) delete fE;
517  fE = new TF1("fE",&erfunc,xmin,xmax,4);
518  fE->SetParameters(pE);
519  proj->Fit("fE","NQL0");
520  fE->GetParameters(pE);
521  maxE = xPeak(fE,xmin,xmax);
522  // cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
523  if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
524  if(!sing){
525  float scale = (fine? 1. : 1.5);
526  // float scale = (fine? 1. : 1);
527  pG[0] = proj->GetBinContent(bMax);
528  if(erfOK){
529  pG[1] = maxE;
530  pG[2] = width(fE,xmin,xmax,maxE);
531  } else {
532  pG[1] = proj->GetMean();
533  pG[2] = proj->GetRMS();
534  }
535  fG = (TF1*)(gROOT->FindObject("fG"));
536  if(fG) delete fG;
537  if(MIP) {
538  fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-0.3*scale),TMath::Min(xmax,pG[1]+0.075*scale));
539  fG->SetParameters(pG);
540  fG->SetParLimits(1,TMath::Max(xmin,pG[1]-0.1*scale),TMath::Min(xmax,pG[1]+0.05*scale));
541  } else {
542  fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-1.5*scale),TMath::Min(xmax,pG[1]+1.5*scale));
543  fG->SetParameters(pG);
544  fG->SetParLimits(1,TMath::Max(xmin,pG[1]-2.*scale),TMath::Min(xmax,pG[1]+2.*scale));
545  }
546  proj->Fit("fG","NLQR0");
547  fG->GetParameters(pG);
548  }
549  return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
550 }
551 
552 //----------------------------------------------------------------------------
553 
554 bool hcalUtil::emcFit(TH1D * proj, bool sing, Double_t * pE, float & maxE, Double_t * pG){
555  TF1 * fE;
556  TF1 * fG;
557  int bMax = proj->GetMaximumBin();
558  float pMax = proj->GetBinCenter(bMax);
559  bool erfOK = true;
560  TAxis *x;
561  x = proj->GetXaxis();
562  double xmin = x->GetXmin();
563  double xmax = x->GetXmax();
564  pE[0] = pMax;
565  pE[1] = proj->GetRMS();
566  pE[2] = -1./pE[1];
567  pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
568  fE = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
569  if(fE) delete fE;
570  fE = new TF1("fE",&erfunc,xmin,xmax,4);
571  fE->SetParameters(pE);
572  proj->Fit("fE","QNL0");
573  fE->GetParameters(pE);
574  maxE = xPeak(fE,xmin,xmax);
575  // cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
576  if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
577  if(!sing){
578  pG[0] = proj->GetBinContent(bMax);
579  if(erfOK){
580  pG[1] = maxE;
581  pG[2] = width(fE,xmin,xmax,maxE) ;
582  } else {
583  pG[1] = proj->GetMean();
584  pG[2] = proj->GetRMS();
585  }
586  fG = (TF1*)(gROOT->FindObject("fG"));
587  if(fG) delete fG;
588  fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
589  fG->SetParameters(pG);
590  fG->SetParLimits(1,TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
591  proj->Fit("fG","QNLR0");
592  fG->GetParameters(pG);
593  }
594  return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
595 }
596 //----------------------------------------------------------------------------
597 
598 bool hcalUtil::emcFit(TH1D * proj, bool sing, Double_t * pE, float & maxE, Double_t * pG, float minx, float maxx){
599  TF1 * fE;
600  TF1 * fG;
601  int bMax = proj->GetMaximumBin();
602  float pMax = proj->GetBinCenter(bMax);
603  bool erfOK = true;
604  TAxis *x(NULL);
605  double xmin(0.);
606  double xmax = x->GetXmax();
607  xmin = TMath::Max((double)minx,xmin);
608  xmax = TMath::Min((double)maxx,xmax);
609  pE[0] = pMax;
610  pE[1] = proj->GetRMS();
611  pE[2] = -1./pE[1];
612  pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
613  fE = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
614  if(fE) delete fE;
615  fE = new TF1("fE",&erfunc,xmin,xmax,4);
616  fE->SetParameters(pE);
617  proj->Fit("fE","QNL0");
618  fE->GetParameters(pE);
619  maxE = xPeak(fE,xmin,xmax);
620  // cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
621  if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
622  if(!sing){
623  pG[0] = proj->GetBinContent(bMax);
624  if(erfOK){
625  pG[1] = maxE;
626  pG[2] = width(fE,xmin,xmax,maxE) ;
627  // cout<<pG[0]<<" "<<pG[1]<<" "<<pG[2]<<endl;
628  } else {
629  pG[1] = proj->GetMean();
630  pG[2] = proj->GetRMS();
631  }
632  fG = (TF1*)(gROOT->FindObject("fG"));
633  if(fG) delete fG;
634  fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
635  fG->SetParameters(pG);
636  fG->SetParLimits(1,TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
637  proj->Fit("fG","QNLR0");
638  fG->GetParameters(pG);
639  }
640  return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
641 }
642 
643 // **************************************************************************
644 
645 bool hcalUtil::twrMipFit(TH1D * pr, Double_t * pE, float & maxE, Double_t * pG){
646  TF1 * fE;
647  TF1 * fG;
648  int bMax = pr->GetMaximumBin();
649  TAxis *x;
650  x = pr->GetXaxis();
651  TH1D * p=(TH1D*)(pr->Clone());
652  float xmean = p->GetMean();
653  xmean = 0.3;
654  int b1 = x->FindBin(xmean-0.07);
655  int b2 = x->FindBin(xmean+0.07);
656  for (int bin = b1; bin<=b2; bin++){
657  p->SetBinError(bin,p->GetBinContent(bin));
658  }
659  double xmin = x->GetXmin();
660  if(xmin<=0.) xmin = 0.05;
661  // double xmin = 0.12;
662  double xmax = x->GetXmax();
663  pE[0] = pr->Integral()*x->GetBinWidth(0)/(log(xmax)-log(xmin));;
664  pE[1] = 1.;
665  pE[2] = 1.;
666  cout<<"Initial par. values"<<pE[0]<<" "<<pE[1]<<" "<< pE[2]<<endl;
667  if(pE[0]<=0.) return false;
668 
669  fE = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
670  if(fE) delete fE;
671  // cout<<(int)erfunc<<endl;
672  fE = new TF1("fE",&powerFun,xmin,xmax,3);
673  fE->SetParameters(pE);
674  // p->Fit("fE","NLR0");
675  p->Fit("fE","NR");
676  // p->Draw();
677  // fE->Draw("same");
678  fE->GetParameters(pE);
679  pG[0] = pr->GetBinContent(bMax)*sqrt(6.28)*0.03;
680  pG[1] = 0.3;
681  pG[2] = 0.03;
682  pG[3] = pE[0];
683  pG[4] = pE[1];
684  pG[5] = pE[2];
685  fG = (TF1*)(gROOT->FindObject("fG"));
686  if(fG) delete fG;
687  fG = new TF1("fG",&power_g,xmin,xmax,6);
688  fG->SetParameters(pG);
689  fG->SetParLimits(3,pE[0],pE[0]);
690  fG->SetParLimits(4,pE[1],pE[1]);
691  fG->SetParLimits(5,pE[2],pE[2]);
692 // fG->FixParameter(3,pE[0]);
693 // fG->FixParameter(4,pE[1]);
694 // fG->FixParameter(5,pE[2]);
695  //pr->Fit("fG","NLR0");
696  pr->Fit("fG","NRB");
697  // pr->Draw("same");
698  // fG->Draw("same");
699  fG->GetParameters(pG);
700  maxE = xPeak(fG,pG[1]-0.1,pG[1]+0.1);
701  cout<<maxE<<endl;
702  // delete p;
703  return ((pG[1]<xmin+0.05||pG[1]>xmax-0.05)? false : true);
704 }
705 
706 
707 // **************************************************************************
708 
709 bool hcalUtil::twrTimeFit(TH1D * pr, Double_t * pE, float & maxE, Double_t * pG){
710  TF1 * fE;
711  TF1 * fG;
712  int bMax = pr->GetMaximumBin();
713  TAxis *x;
714  x = pr->GetXaxis();
715  TH1D * p=(TH1D*)(pr->Clone());
716  float xmean = p->GetMean();
717  xmean = 0.3;
718  int b1 = x->FindBin(xmean-0.07);
719  int b2 = x->FindBin(xmean+0.07);
720  for (int bin = b1; bin<=b2; bin++){
721  p->SetBinError(bin,p->GetBinContent(bin));
722  }
723  double xmin = x->GetXmin();
725  // double xmin = 0.12;
726  double xmax = x->GetXmax();
727  pE[0] = pr->Integral()*x->GetBinWidth(0)/(log(xmax)-log(xmin));;
728  pE[1] = 1.;
729  pE[2] = 1.;
730  fE = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
731  if(fE) delete fE;
732  // cout<<(int)erfunc<<endl;
733  fE = new TF1("fE",&powerFun,xmin,xmax,3);
734  fE->SetParameters(pE);
735  // p->Fit("fE","NLR0");
736  p->Fit("fE","NR");
737  // p->Draw();
738  // fE->Draw("same");
739  fE->GetParameters(pE);
740  pG[0] = pr->GetBinContent(bMax)*sqrt(6.28)*0.03;
741  pG[1] = 0.3;
742  pG[2] = 0.03;
743  pG[3] = pE[0];
744  pG[4] = pE[1];
745  pG[5] = pE[2];
746  fG = (TF1*)(gROOT->FindObject("fG"));
747  if(fG) delete fG;
748  fG = new TF1("fG",&power_g,xmin,xmax,6);
749  fG->SetParameters(pG);
750  fG->SetParLimits(3,pE[0],pE[0]);
751  fG->SetParLimits(4,pE[1],pE[1]);
752  fG->SetParLimits(5,pE[2],pE[2]);
753 // fG->FixParameter(3,pE[0]);
754 // fG->FixParameter(4,pE[1]);
755 // fG->FixParameter(5,pE[2]);
756  //pr->Fit("fG","NLR0");
757  pr->Fit("fG","NRB");
758  // pr->Draw("same");
759  // fG->Draw("same");
760  fG->GetParameters(pG);
761  maxE = xPeak(fG,pG[1]-0.1,pG[1]+0.1);
762  cout<<maxE<<endl;
763  // delete p;
764  return ((pG[1]<xmin+0.05||pG[1]>xmax-0.05)? false : true);
765 }
766 
767 // ====================================================================================
768 // ************************************************************************************
769 // ====================================================================================
770 
771 #include "hLabHelper.C"
772 
773 #include "hcalHelper.C"
774 
775 #include "hcalTree.C"
776 
777 #include "tileHelper.C"
778 
779 #include "tileTree.C"
780