Prototype2/hcalLab/hcalUtil.C

0001 #include <iomanip>
0002 #include <algorithm>
0003 #include <numeric>
0004 #include <iostream>
0005 #include <fstream>
0006 #include <cmath>
0007 #include <ctime>
0008 #include <cstdlib>
0009 #include <stdint.h>
0010 #include <stddef.h>
0011 #include <stdlib.h>
0012 #include <string.h>
0013
0014
0015 // #include <TROOT.h>
0016 // #include <TFile.h>
0017 // #include <TTree.h>
0018 // #include <TGMsgBox.h>
0019 // #include <TCanvas.h>
0020 // #include <TGraph.h>
0021 // #include <TString.h>
0022 // #include <TText.h>
0023 // #include <TApplication.h>
0024 // #include <TStyle.h>
0025 // #include <TSystem.h>
0026 // #include <TGClient.h>
0027 // #include <TGWindow.h>
0028 // #include <TH1.h>
0029 // #include <TH1F.h>
0030 // #include <TH1D.h>
0031 // #include <TProfile.h>
0032 // #include <TF1.h>
0033 // #include <TH2.h>
0034 // #include <TH3.h>
0035 // #include <TMath.h>
0036 // #include <TSpectrum.h>
0037 // #include "TVirtualFitter.h"
0038 // #include <TSystemDirectory.h>
0039 // #include <TSystemFile.h>
0040
0041 //#ifndef __CINT__
0042
0043 #include <fileEventiterator.h>
0044 #include <Event.h>
0045 #include <EventTypes.h>
0046 #include <PHTimeStamp.h>
0047 #include <packetConstants.h>
0048 #include <packet.h>
0049 #include <packet_hbd_fpgashort.h>
0050
0051 //#endif
0052
0053
0054 #include "hcalUtil.h"
0055
0056
0057 hcalUtil    * hcalUtil    :: single = 0;
0058 hLabHelper  * hLabHelper  :: single = 0;
0059 hcalHelper  * hcalHelper  :: single = 0;
0060 hcalTree    * hcalTree    :: single = 0;
0061 tileHelper  * tileHelper  :: single = 0;
0062 tileTree    * tileTree    :: single = 0;
0063
0064 using namespace std;
0065 using namespace TMath;
0066
0067 TCanvas * fitFunc;
0068
0069 TCanvas * scCalibDisplay;
0070
0071
0072 //   ************************************************************************************
0073
0074
0075 //#endif /* _CINT_ */
0076
0077 //  cosmic runs of interest
0078 static const int   cosmicruns = 12;
0079 int   cosmicRN[cosmicruns]    = {212, 215, 216, 218, 223, 337, 338, 339, 405, 406, 229, 240};
0080
0081 //  LED calibration and PEDESTALS
0082 static const int   pedruns    = 1;
0083 int   pedcalibRN[pedruns]     = {646};
0084 static const int   ledruns    = 1;
0085 int   ledcalibRN[ledruns]     = {647};
0086
0087 //  -----------------------------------------------------------------------------------------------
0088
0089 //  Unipolar pulse (0.3164*pow(x,4)*exp(-x*1.5) - Unity integral, peak ~1/3 of integral)
0090 Double_t par0[]       = {     1.,  (Double_t)NSAMPLES/2, 3., 2.,   PEDESTAL,   0. };
0091 Double_t par0Max[]    = {  -100.,  (Double_t)NSAMPLES,   4., 5.,   50000.,     0.3};
0092 //Double_t par0Min[]    = {    1., RISETIME-1,            3., 1.,   1950.,    -0.2};
0093 Double_t par0Min[]    = {-50000.,  0.,                   2., 0.,  -50000.,    -0.3};
0094
0095
0096 TString             runId;
0097
0098 //   The next few functions are to help processing data from HCal Lab recorded via HBD electronics
0099
0100 // **************************************************************************
0101
0102 Double_t ps( Double_t *chadc, Double_t *par)
0103 {
0104   Double_t x = chadc[0];
0105
0106   Double_t mu = par[0];
0107   Double_t g = par[1];
0108   Double_t sigma = par[2];
0109   Double_t norm = par[3];
0110
0111   Double_t sum = 0.0;
0112   for (Double_t n = 1.0; n < 100.0; n++ ) {
0113     sum += TMath::Poisson(n,mu)*TMath::Gaus(x,n*g,sigma,kTRUE);
0114   }
0115
0116   return norm*sum;
0117 }
0118 // **************************************************************************
0119 Int_t  calibPeaks(0), calibMode (0); Double_t spCalib(10.);
0120 // **************************************************************************
0121 //  mode   0  - peaks uncorrelated (unconstrained)
0122 //  mode   1  - peaks equidistant
0123 //  assume background is linear what is most likely not true ....
0124 //  par[0]   background constant
0125 //     [1]   background slope
0126 //     [2]   background quadratic term
0127 //  assume the signal is the sequence of equally spaced gaussian peaks
0128 //  par[3]   peak to peak separation
0129 //     [4]   normalization constant for the first peak (pedestal)
0130 //     [5]   location of the first (pedestal) peak
0131 //     [6]   rms of the first peak - at this stage we will allow it to vary peak-to-peak
0132 //     then 3 parameters per signal peak (normalization, mean and RMS)
0133 Double_t sipmPeaks(Double_t *x, Double_t *par) {
0134    Double_t result  = par[0] + par[1]*x[0] + par[2]*x[0]*x[0];
0135    //  add pedestal peak first
0136    result          += par[4]*TMath::Gaus(x[0],par[5],par[6]);
0137    for (Int_t p=1; p<calibPeaks; p++) {
0138       Double_t norm  = par[4 + p*3];
0139       Double_t mean  = par[5 + p*3];
0140       if(calibMode) {
0141     mean  = (Int_t)((mean+spCalib/2)/spCalib)*par[3];
0142       }
0143       Double_t sigma = par[6 + p*3];
0144       result += norm*TMath::Gaus(x[0],mean,sigma);
0145    }
0146    return result;
0147    ;}
0148
0149
0150 // **************************************************************************
0151 // do single cell calibration using low end data stored in rv_# and fm_# channel histograms
0152 // the assumption - we have beautiful picture with at least 20 peaks !!!!!!
0153 //  mode   0  - peaks uncorrelated (unconstrained)
0154 //  mode   1  - peaks equidistant
0155 Int_t sipmCalib(TString & hName, Int_t mode){
0156   TH1F * data = (TH1F*)(gROOT->FindObject(hName));
0157   if(!data){
0158     cout<<"spmCalib  Data Histogramm "<<hName<<"  not found"<<endl;
0159     return 0;
0160   }
0161   cout<<hName<<"  calibration mode "<<mode<<endl;
0162   calibMode = mode;
0163   Int_t      maxPeaks(20);        //  maximum number of peaks we have chance to resolve
0164   Double_t   range(150.);         //  range of amplitudes for peak search
0165   //  second degree polinomial background
0166   Double_t   bckgr0(100.), bckgr1(-1./range),  bckgr2(0.);
0167   Double_t   pNorm(1000.), pLoc(0.),           pRMS(1.);
0168   Double_t   par[4+maxPeaks*3];   //  backgrownd, pedestal and nPeaks of gaussian peaks with constant peak-to-peak separation
0169   //  initialization
0170   par[0] = bckgr0;
0171   par[1] = bckgr1;
0172   par[2] = bckgr2;
0173   par[3] = spCalib;            //  peak-to-peak separation
0174   par[4] = pNorm;              //  normalization for pedestal peak
0175   par[5] = pLoc;               //  pedestal mean
0176   par[6] = pRMS;               //  pedestal rms
0177   //  all other peaks are initially positioned equidistantly
0178   for (int p=1; p<= maxPeaks; p++) {
0179       par[4 + p*3] = pNorm;                       //  norm of the peak p
0180       par[5 + p*3] = par[3]*p;                    //  mean of the peak p
0181       par[6 + p*3] = pRMS;                        //  rms  of the peak p
0182   }
0183   TH1F * dataCl  = (TH1F*)(data->Clone("dataCl"));
0184   dataCl        -> SetAxisRange(0., range);
0185   TH1F * dataCl2 = (TH1F*)data->Clone("dataCl2");
0186   dataCl2       -> SetAxisRange(0., range);
0187   TString cTitle  = "Peaks from "; cTitle += hName;
0188   TCanvas * peaks = new TCanvas("peaks",cTitle,10,10,1000,900);
0189   peaks         -> Divide(1,2);
0190   peaks         -> cd(1);
0191   dataCl        -> Draw();
0192
0193   //  Use TSpectrum to find the peak candidates
0194   TSpectrum *s = new TSpectrum(maxPeaks);
0195   s->SetResolution(2.);
0196   Int_t nfound = s->Search(dataCl, 1.5, "nobackground    ", 0.10);
0197   printf("TSpectrum:  Found %d candidate peaks to fit\n",nfound);
0198
0199   //Estimate background using TSpectrum::Background
0200   TH1 *hb = s->Background(dataCl, 20, "same");
0201   if (hb)   peaks->Update();
0202   //  if only one or no peaks  found (only pedestals) - we are done
0203   if(nfound<=1) {
0204     cout<<"sipmCalib  One or No peaks were found - EXITING"<<endl;
0205     calibPeaks = nfound;
0206
0207   }  else {
0208     peaks           -> cd(2);
0209
0210     //  function to estimate linear background
0211     //  TF1 *fline       = new TF1("fline","pol1",0.,range);
0212     //  function to estimate quadratic background
0213     TF1 *fline       = new TF1("fline","pol2",0.,range);
0214
0215     dataCl2          -> Fit("fline","qn");
0216     //Loop on all found peaks. Eliminate peaks at the background level
0217     par[0]           = fline->GetParameter(0);
0218     par[1]           = fline->GetParameter(1);
0219     par[2]           = fline->GetParameter(2);
0220     //  cout<<"Line par "<<par[0]<<"   "<<par[1]<<"  "<<par[2]<<endl;
0221     Float_t * xpeaks  = s->GetPositionX();
0222     // peak ordering
0223     Int_t pindex[nfound];
0224     TMath::Sort(nfound, xpeaks, pindex, kFALSE);
0225     //  peaks now  ordered, let's now find the peak to peak separation
0226     Double_t sep[nfound-1];
0227     //    Int_t    sepind[nfound-1];
0228     //  pedestal peak in my definition == what is left of a sygnal after pedestal subtraction - it is definitely shifted
0229     for (int ip = 1; ip<nfound; ip++){
0230       sep[ip-1] = xpeaks[pindex[ip]];
0231       if(ip>1)  sep[ip-1] -= xpeaks[pindex[ip-1]];
0232       cout<<"Separation  "<<ip<<" val "<<sep[ip-1]<<endl;
0233     }
0234     //    TMath::Sort(nfound-1, sep, sepind, kFALSE);
0235     //  zero approximation for peak0-to-peak separation
0236     Int_t cont(0); Double_t avsep(0.);
0237     for (int ip = 0; ip < nfound-1; ip++){
0238       if(ip==0) {
0239     //  avsep = sep[sepind[ip]];
0240     avsep = sep[ip];
0241     cont ++;
0242       } else {
0243     //  if(abs((avsep- sep[sepind[ip]])/(avsep+ sep[sepind[ip]]))<0.2) {
0244     if(abs((avsep- sep[ip])/(avsep+ sep[ip]))<0.2) {
0245       avsep = (avsep*cont+sep[ip])/(cont+1);
0246       cont++;
0247     }
0248       }
0249       cout<<"ip  "<<ip<<" cont "<<cont<<"  avsep "<<avsep<<" var "<<(ip!=0? (avsep- sep[ip])/(avsep+ sep[ip])  : 0.)<<endl;
0250     }
0251     cout<<"Zero approximation to sSingle Pixel Calib = "<<avsep<<"   Contributors "<<cont<<endl;
0252     spCalib    = avsep;
0253     par[3]     = avsep;
0254     calibPeaks = 0;
0255     for (int p=0; p<nfound; p++) {
0256       Double_t xp    = xpeaks[pindex[p]];
0257       Int_t    bin   = dataCl2->GetXaxis()->FindBin(xp);
0258       Double_t yp    = dataCl2->GetBinContent(bin);
0259       cout<<"Testing Peak (LFit) "<<p<<" index "<<pindex[p]<<"  xp "<<xp<<"  bin "<<bin<<"  yP "<<yp<<"  BCKG "<<fline->Eval(xp)<<endl;
0260       //  extract data for this bin from background shape
0261       Double_t bckgr = hb->GetBinContent(bin);
0262       cout<<"Testing Peak (TSp ) "<<p<<"  xp "<<xp<<"  bin "<<bin<<"  Bcgrd "<<bckgr<<" Signal "<<yp-bckgr<<endl;
0263       if (yp-TMath::Sqrt(yp) < bckgr) continue;
0264       par[4 + calibPeaks*3] = yp;
0265       par[5 + calibPeaks*3] = xp;
0266       par[6 + calibPeaks*3] = (p==0?  pRMS  : 2.);
0267       calibPeaks++;
0268     }
0269   }
0270   printf("Found %d useful peaks to fit\n",calibPeaks);
0271   printf("Now fitting: Be patient\n");
0272   TF1 *fit = new TF1("fit",&sipmPeaks,0,range,4+calibPeaks*3);  // first peak is the pedestal
0273   //we may have more than the default 25 parameters
0274   TVirtualFitter::Fitter(dataCl2,10+calibPeaks*3);
0275   fit->SetParameters(par);
0276   //  if mode!=0 we fix all location parameters to what came from TSpectrum
0277   if(mode) for(int ip=0; ip<calibPeaks; ip++) {
0278       fit->FixParameter(5+ip*3, par[5+ip*3]);
0279       fit->SetParLimits(3, 0.5*par[3], 1.5*par[3]);
0280     }  else fit->FixParameter(3, par[3]);
0281   //  in all cases we set max RMS equal to pps/2.
0282   for(int ip=0; ip<calibPeaks; ip++) fit->SetParLimits(6+ip*3, 0.5, par[3]/2.);
0283
0284   fit->SetNpx(1000);
0285   dataCl2->Fit("fit");
0286   fit->GetParameters(par);
0287
0288   char p5[10], p6[10], p3[10];
0289   sprintf(p5, "%.2f ", par[5]);
0290   sprintf(p6, "%.2f ", par[6]);
0291   sprintf(p3, "%.2f ", par[3]);
0292
0293   TString ped = "Pedestal at     "; ped += p5; ped += "  RMS = "; ped += p6;
0294   TString cal = "Pixel Calibtion "; cal += p3; cal += " ADC counts per pixel";
0295   TText label;
0296   Double_t  ymax = data-> GetMaximum();
0297   label.DrawText(30, ymax*0.8,ped);
0298   label.DrawText(30, ymax*0.6,cal);
0299   return nfound;
0300 }
0301
0302 // **************************************************************************
0303 Double_t AsymToX(Double_t ax){
0304   return 12.5*(1 + ax/0.19);
0305 }
0306
0307
0308 // **************************************************************************
0309
0310 Double_t signalShape(Double_t *x, Double_t * par){
0311   Double_t signal(0.), pedestal(0.);
0312   //  cout<<"signalShape x="<<x[0]<<" (0)="<<par[0]<<" (1)="<<par[1]<<" (2)="<<par[2]<<" (3)="<<par[3]<<" (4)="<<par[4]<<" (5)="<<par[5]<<endl;
0313   //  if(x[0]<par[1])  return PEDESTAL;
0314   pedestal = par[4]+x[0]*par[5];
0315   if(x[0]<par[1])   return pedestal;
0316   //  signal contribution
0317   //  cout<<"(0) "<<pow((x[0]-par[1]),par[2])<<endl;
0318   //  cout<<"(1) "<<exp(-(x[0]-par[1])*par[3])<<endl;
0319   // Double_t a0 = par[0];
0320   // Double_t a1 = pow((x[0]-par[1]),par[2]);
0321   // Double_t a2 = exp(-(x[0]-par[1])*par[3]);
0322   // Double_t a  = a0*a1*a2;
0323   // //  cout<<a0<<"   "<<a1<<"  "<<a2<<"  "<<a<<endl;
0324   // signal   = a;
0325   // return a;
0326   signal   = par[0]*pow((x[0]-par[1]),par[2])*exp(-(x[0]-par[1])*par[3]);
0327   //  cout<<"  X  "<<x[0]<<"  signal "<<signal<<"  pedestal  "<<pedestal<<endl;
0328   return   signal+pedestal;
0329 }
0330
0331 // **************************************************************************
0332
0333 // Double_t shape0(Double_t *x, Double_t * par){
0334 //   return par[0]*pow(x[0],4.)*(exp(-x[0]*16./fN));    //  fN is set to (float)NSAMPLES
0335 // }
0336
0337 // **************************************************************************
0338
0339 Double_t erfunc(Double_t *x, Double_t * par){
0340   //  Double_t erfcont =  par[3]*(TMath::Erf((x[0]-par[0])/par[1])+1.);
0341   //  Double_t expcont =  TMath::Exp(par[2]*x[0]);
0342   return par[3]*(TMath::Erf((x[0]-par[0])/par[1])+1.)*TMath::Exp(par[2]*x[0]);
0343 }
0344 // **************************************************************************
0345
0346 Double_t erf_g(Double_t *x, Double_t * par){
0347   float  g = par[0]/(sqrt(6.2830)*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
0348   return g+par[6]*(TMath::Erf((x[0]-par[3])/par[4])+1.)*TMath::Exp(par[5]*x[0]);
0349 }
0350 // **************************************************************************
0351
0352 Double_t powerFun(Double_t *x, Double_t * par){
0353   //  return par[0]*pow((x[0]),par[1]*pow((par[2]-x[0]),par[3]);
0354  return par[0]*pow(x[0],par[1])*pow((1-x[0]),par[2]);
0355  //*pow((0.5-x[0]),par[2]);
0356 }
0357
0358 // **************************************************************************
0359
0360 Double_t power_g(Double_t *x, Double_t * par){
0361   float  g = par[0]/(sqrt(6.2830)*par[2])*TMath::Exp(-(x[0]-par[1])*(x[0]-par[1])/(2.*par[2]*par[2]));
0362   float  p = par[3]*pow(x[0],par[4])*pow((1-x[0]),par[5]);
0363   //*pow((0.5-x[0]),par[5]);
0364   return g+p;
0365 }
0366
0367
0368
0369 //   ====================================================================================
0370 //   ************************************************************************************
0371 //   ====================================================================================
0372
0373
0374 hcalUtil::hcalUtil(){
0375   char * s = getenv("BASE_DIR");
0376   TString basedir = (s==0||strlen(s)==0)? "/home/kistenev/workarea/" : s;
0377 }
0378
0379 //----------------------------------------------------------------------------
0380 // void hcalUtil::setBaseDir(char * bDir){
0381 //   basedir.Resize(0);
0382 //   basedir += bDir;
0383 // }
0384
0385 //----------------------------------------------------------------------------
0386
0387 float hcalUtil::xPeak(TF1* f, float x1, float x2){
0388   //  searches for peak position in the TF1 function - looks into derivatives
0389   float y1 = f->Derivative(x1, 0, 0.002);
0390   float y2 = f->Derivative(x2, 0, 0.002);
0391   if(y1*y2>0) {
0392     return (f->Eval(x1)>f->Eval(x2))? x1 : x2;
0393   }
0394   float xx1 = x1;
0395   float xx2 = x2;
0396   float dx = fabs(xx2-xx1);
0397   float x = 0.;
0398   while(y1*y2<0&&dx>0.001){
0399     x  = (xx1+xx2)/2.;
0400     float y = f->Derivative(x, 0, 0.01);
0401     if(y1*y>0)  xx1 = x; else xx2 = x;
0402     dx = fabs(xx2-xx1);
0403   }
0404   return x;
0405 }
0406
0407 //----------------------------------------------------------------------------
0408
0409 float hcalUtil::width(TF1* f, float x1, float x2, float peak){
0410   //  find width of the function (erfunc only) !!!
0411   float yhm = f->Eval(peak)/2;
0412
0413   float xl = x1;
0414   float xr = peak;
0415
0416   float xm = (xr+xl)/2.;
0417   float ym = f->Eval(xm);
0418   //  cout<<yhm<<" xl "<< xl<<" xr "<<xr<<" xm "<<xm<<" ym "<<ym<<endl;
0419   while(abs(ym-yhm)>0.01*yhm && xr-xl>0.001){
0420     if(ym>yhm) xr = xm; else xl = xm;
0421     xm = (xr+xl)/2.;
0422     ym = f->Eval(xm);
0423     //cout<<"left xl "<<xl<<" xr "<<xr<<" xm "<<xm<<" ym "<<ym<<endl;
0424   }
0425   bool lok = (abs(ym-yhm)>0.1*yhm)? false : true;
0426   float xml = xm;
0427   xl = peak;
0428   xr = x2;
0429   xm = (xr+xl)/2.;
0430   ym = f->Eval(xm);
0431   while(abs(ym-yhm)>0.01*yhm&&xr-xl>0.001){
0432     if(ym>yhm) xl = xm; else xr = xm;
0433     xm = (xr+xl)/2.;
0434     ym = f->Eval(xm);
0435     //cout<<"right xl "<<xl<<" xr "<<xr<<" xm "<<xm<<" ym "<<ym<<endl;
0436   }
0437   float xmr = xm;
0438   bool rok = (abs(ym-yhm)>0.1*yhm)? false : true;
0439   float w   = 0;
0440   if(lok){
0441     if(rok) w = xmr-xml; else  w = 2*(peak-xml);
0442   }  else if (rok){
0443     w = 2*(xmr-peak);
0444   } else w = x2-x1;
0445   return w/2.35;
0446 }
0447
0448 //----------------------------------------------------------------------------
0449
0450 bool hcalUtil::shapeFit(TH1D * proj,  Double_t * pE, float & maxE, Double_t * pG, bool fine, bool MIP){
0451   TF1 * fE;
0452   TF1 * fG;
0453   int    bMax = proj->GetMaximumBin();
0454   float  pMax = proj->GetBinCenter(bMax);
0455   bool   erfOK = true;
0456   TAxis *x;
0457   x = proj->GetXaxis();
0458   double xmin = x->GetXmin();
0459   double xmax = x->GetXmax();
0460   pE[0] = pMax;
0461   pE[1] = proj->GetRMS();
0462   pE[2] = -1./pE[1];
0463   pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
0464   fE    = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
0465   if(fE) delete fE;
0466   fE = new TF1("fE",&erfunc,xmin,xmax,4);
0467   fE->SetParameters(pE);
0468   proj->Fit("fE","NQL0");
0469   fE->GetParameters(pE);
0470   maxE = xPeak(fE,xmin,xmax);
0471   //  cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
0472   if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
0473   float scale = (fine? 1. : 1.5);
0474   //    float scale = (fine? 1. : 1);
0475   pG[0] = proj->GetBinContent(bMax);
0476   if(erfOK){
0477     pG[1] =  maxE;
0478     pG[2] =  width(fE,xmin,xmax,maxE);
0479   } else {
0480     pG[1] = proj->GetMean();
0481     pG[2] = proj->GetRMS();
0482   }
0483   fG = (TF1*)(gROOT->FindObject("fG"));
0484   if(fG) delete fG;
0485   if(MIP) {
0486     fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-0.3*scale),TMath::Min(xmax,pG[1]+0.075*scale));
0487     fG->SetParameters(pG);
0488     fG->SetParLimits(1,TMath::Max(xmin,pG[1]-0.1*scale),TMath::Min(xmax,pG[1]+0.05*scale));
0489   } else {
0490     fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-1.5*scale),TMath::Min(xmax,pG[1]+1.5*scale));
0491     fG->SetParameters(pG);
0492     fG->SetParLimits(1,TMath::Max(xmin,pG[1]-2.*scale),TMath::Min(xmax,pG[1]+2.*scale));
0493   }
0494   proj->Fit("fG","NLQR0");
0495   fG->GetParameters(pG);
0496   return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
0497 }
0498
0499 //----------------------------------------------------------------------------
0500
0501 bool hcalUtil::emcFit(TH1D * proj,  bool fine, bool sing, Double_t * pE, float & maxE, Double_t * pG, bool MIP){
0502   TF1 * fE;
0503   TF1 * fG;
0504   int    bMax = proj->GetMaximumBin();
0505   float  pMax = proj->GetBinCenter(bMax);
0506   bool   erfOK = true;
0507   TAxis *x;
0508   x = proj->GetXaxis();
0509   double xmin = x->GetXmin();
0510   double xmax = x->GetXmax();
0511   pE[0] = pMax;
0512   pE[1] = proj->GetRMS();
0513   pE[2] = -1./pE[1];
0514   pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
0515   fE    = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
0516   if(fE) delete fE;
0517   fE = new TF1("fE",&erfunc,xmin,xmax,4);
0518   fE->SetParameters(pE);
0519   proj->Fit("fE","NQL0");
0520   fE->GetParameters(pE);
0521   maxE = xPeak(fE,xmin,xmax);
0522   //  cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
0523   if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
0524   if(!sing){
0525     float scale = (fine? 1. : 1.5);
0526     //    float scale = (fine? 1. : 1);
0527     pG[0] = proj->GetBinContent(bMax);
0528     if(erfOK){
0529       pG[1] =  maxE;
0530       pG[2] =  width(fE,xmin,xmax,maxE);
0531     } else {
0532       pG[1] = proj->GetMean();
0533       pG[2] = proj->GetRMS();
0534     }
0535     fG = (TF1*)(gROOT->FindObject("fG"));
0536     if(fG) delete fG;
0537     if(MIP) {
0538       fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-0.3*scale),TMath::Min(xmax,pG[1]+0.075*scale));
0539       fG->SetParameters(pG);
0540       fG->SetParLimits(1,TMath::Max(xmin,pG[1]-0.1*scale),TMath::Min(xmax,pG[1]+0.05*scale));
0541     } else {
0542       fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-1.5*scale),TMath::Min(xmax,pG[1]+1.5*scale));
0543       fG->SetParameters(pG);
0544       fG->SetParLimits(1,TMath::Max(xmin,pG[1]-2.*scale),TMath::Min(xmax,pG[1]+2.*scale));
0545     }
0546     proj->Fit("fG","NLQR0");
0547     fG->GetParameters(pG);
0548   }
0549   return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
0550 }
0551
0552 //----------------------------------------------------------------------------
0553
0554 bool hcalUtil::emcFit(TH1D * proj,  bool sing, Double_t * pE, float & maxE, Double_t * pG){
0555   TF1 * fE;
0556   TF1 * fG;
0557   int    bMax = proj->GetMaximumBin();
0558   float  pMax = proj->GetBinCenter(bMax);
0559   bool   erfOK = true;
0560   TAxis *x;
0561   x = proj->GetXaxis();
0562   double xmin = x->GetXmin();
0563   double xmax = x->GetXmax();
0564   pE[0] = pMax;
0565   pE[1] = proj->GetRMS();
0566   pE[2] = -1./pE[1];
0567   pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
0568   fE    = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
0569   if(fE) delete fE;
0570   fE = new TF1("fE",&erfunc,xmin,xmax,4);
0571   fE->SetParameters(pE);
0572   proj->Fit("fE","QNL0");
0573   fE->GetParameters(pE);
0574   maxE = xPeak(fE,xmin,xmax);
0575   //  cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
0576   if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
0577   if(!sing){
0578     pG[0] = proj->GetBinContent(bMax);
0579     if(erfOK){
0580       pG[1] =  maxE;
0581       pG[2] = width(fE,xmin,xmax,maxE) ;
0582     } else {
0583       pG[1] = proj->GetMean();
0584       pG[2] = proj->GetRMS();
0585     }
0586     fG = (TF1*)(gROOT->FindObject("fG"));
0587     if(fG) delete fG;
0588     fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
0589     fG->SetParameters(pG);
0590     fG->SetParLimits(1,TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
0591     proj->Fit("fG","QNLR0");
0592     fG->GetParameters(pG);
0593   }
0594   return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
0595 }
0596 //----------------------------------------------------------------------------
0597
0598 bool hcalUtil::emcFit(TH1D * proj,  bool sing, Double_t * pE, float & maxE, Double_t * pG, float minx, float maxx){
0599   TF1 * fE;
0600   TF1 * fG;
0601   int    bMax = proj->GetMaximumBin();
0602   float  pMax = proj->GetBinCenter(bMax);
0603   bool   erfOK = true;
0604   TAxis *x(NULL);
0605   double xmin(0.);
0606   double xmax = x->GetXmax();
0607   xmin = TMath::Max((double)minx,xmin);
0608   xmax = TMath::Min((double)maxx,xmax);
0609   pE[0] = pMax;
0610   pE[1] = proj->GetRMS();
0611   pE[2] = -1./pE[1];
0612   pE[3] = proj->GetBinContent(bMax)/TMath::Exp(pE[2]*pE[0]);
0613   fE    = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
0614   if(fE) delete fE;
0615   fE = new TF1("fE",&erfunc,xmin,xmax,4);
0616   fE->SetParameters(pE);
0617   proj->Fit("fE","QNL0");
0618   fE->GetParameters(pE);
0619   maxE = xPeak(fE,xmin,xmax);
0620   //  cout<<"Position of Maximum: "<<maxE<<" "<<xmin<<" "<<xmax<<endl;
0621   if(maxE<=xmin||maxE>=xmax||pE[0]<=xmin/*||pE[0]>=xmax*/) erfOK = false;
0622   if(!sing){
0623     pG[0] = proj->GetBinContent(bMax);
0624     if(erfOK){
0625       pG[1] =  maxE;
0626       pG[2] = width(fE,xmin,xmax,maxE) ;
0627       //      cout<<pG[0]<<" "<<pG[1]<<" "<<pG[2]<<endl;
0628     } else {
0629       pG[1] = proj->GetMean();
0630       pG[2] = proj->GetRMS();
0631     }
0632     fG = (TF1*)(gROOT->FindObject("fG"));
0633     if(fG) delete fG;
0634     fG = new TF1("fG","gaus",TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
0635     fG->SetParameters(pG);
0636     fG->SetParLimits(1,TMath::Max(xmin,pG[1]-pG[2]),TMath::Min(xmax,pG[1]+pG[2]));
0637     proj->Fit("fG","QNLR0");
0638     fG->GetParameters(pG);
0639   }
0640   return ((!erfOK||pG[1]<xmin+0.1||pG[1]>xmax-0.05)? false : true);
0641 }
0642
0643 // **************************************************************************
0644
0645 bool hcalUtil::twrMipFit(TH1D * pr,  Double_t * pE, float & maxE, Double_t * pG){
0646   TF1 * fE;
0647   TF1 * fG;
0648   int    bMax = pr->GetMaximumBin();
0649   TAxis *x;
0650   x = pr->GetXaxis();
0651   TH1D * p=(TH1D*)(pr->Clone());
0652   float xmean = p->GetMean();
0653   xmean = 0.3;
0654   int b1 = x->FindBin(xmean-0.07);
0655   int b2 = x->FindBin(xmean+0.07);
0656   for (int bin = b1; bin<=b2; bin++){
0657     p->SetBinError(bin,p->GetBinContent(bin));
0658   }
0659   double xmin = x->GetXmin();
0660   if(xmin<=0.) xmin = 0.05;
0661   //  double xmin = 0.12;
0662   double xmax = x->GetXmax();
0663   pE[0] = pr->Integral()*x->GetBinWidth(0)/(log(xmax)-log(xmin));;
0664   pE[1] = 1.;
0665   pE[2] = 1.;
0666   cout<<"Initial par. values"<<pE[0]<<" "<<pE[1]<<" "<< pE[2]<<endl;
0667   if(pE[0]<=0.) return false;
0668
0669   fE    = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
0670   if(fE) delete fE;
0671   //  cout<<(int)erfunc<<endl;
0672   fE = new TF1("fE",&powerFun,xmin,xmax,3);
0673   fE->SetParameters(pE);
0674   //  p->Fit("fE","NLR0");
0675   p->Fit("fE","NR");
0676   //  p->Draw();
0677   //  fE->Draw("same");
0678   fE->GetParameters(pE);
0679   pG[0] = pr->GetBinContent(bMax)*sqrt(6.28)*0.03;
0680   pG[1] = 0.3;
0681   pG[2] = 0.03;
0682   pG[3] = pE[0];
0683   pG[4] = pE[1];
0684   pG[5] = pE[2];
0685   fG = (TF1*)(gROOT->FindObject("fG"));
0686   if(fG) delete fG;
0687   fG = new TF1("fG",&power_g,xmin,xmax,6);
0688   fG->SetParameters(pG);
0689   fG->SetParLimits(3,pE[0],pE[0]);
0690   fG->SetParLimits(4,pE[1],pE[1]);
0691   fG->SetParLimits(5,pE[2],pE[2]);
0692 //    fG->FixParameter(3,pE[0]);
0693 //    fG->FixParameter(4,pE[1]);
0694 //    fG->FixParameter(5,pE[2]);
0695   //pr->Fit("fG","NLR0");
0696   pr->Fit("fG","NRB");
0697   //  pr->Draw("same");
0698   //  fG->Draw("same");
0699   fG->GetParameters(pG);
0700   maxE = xPeak(fG,pG[1]-0.1,pG[1]+0.1);
0701   cout<<maxE<<endl;
0702   //  delete p;
0703   return ((pG[1]<xmin+0.05||pG[1]>xmax-0.05)? false : true);
0704 }
0705
0706
0707 // **************************************************************************
0708
0709 bool hcalUtil::twrTimeFit(TH1D * pr,  Double_t * pE, float & maxE, Double_t * pG){
0710   TF1 * fE;
0711   TF1 * fG;
0712   int    bMax = pr->GetMaximumBin();
0713   TAxis *x;
0714   x = pr->GetXaxis();
0715   TH1D * p=(TH1D*)(pr->Clone());
0716   float xmean = p->GetMean();
0717   xmean = 0.3;
0718   int b1 = x->FindBin(xmean-0.07);
0719   int b2 = x->FindBin(xmean+0.07);
0720   for (int bin = b1; bin<=b2; bin++){
0721     p->SetBinError(bin,p->GetBinContent(bin));
0722   }
0723   double xmin = x->GetXmin();
0724   ///////  if(xmin<=0.) xmin = 0.05;
0725   //  double xmin = 0.12;
0726   double xmax = x->GetXmax();
0727   pE[0] = pr->Integral()*x->GetBinWidth(0)/(log(xmax)-log(xmin));;
0728   pE[1] = 1.;
0729   pE[2] = 1.;
0730   fE    = (TF1*)(gROOT->GetListOfFunctions()->FindObject("fE"));
0731   if(fE) delete fE;
0732   //  cout<<(int)erfunc<<endl;
0733   fE = new TF1("fE",&powerFun,xmin,xmax,3);
0734   fE->SetParameters(pE);
0735   //  p->Fit("fE","NLR0");
0736   p->Fit("fE","NR");
0737   //  p->Draw();
0738   //  fE->Draw("same");
0739   fE->GetParameters(pE);
0740   pG[0] = pr->GetBinContent(bMax)*sqrt(6.28)*0.03;
0741   pG[1] = 0.3;
0742   pG[2] = 0.03;
0743   pG[3] = pE[0];
0744   pG[4] = pE[1];
0745   pG[5] = pE[2];
0746   fG = (TF1*)(gROOT->FindObject("fG"));
0747   if(fG) delete fG;
0748   fG = new TF1("fG",&power_g,xmin,xmax,6);
0749   fG->SetParameters(pG);
0750   fG->SetParLimits(3,pE[0],pE[0]);
0751   fG->SetParLimits(4,pE[1],pE[1]);
0752   fG->SetParLimits(5,pE[2],pE[2]);
0753 //    fG->FixParameter(3,pE[0]);
0754 //    fG->FixParameter(4,pE[1]);
0755 //    fG->FixParameter(5,pE[2]);
0756   //pr->Fit("fG","NLR0");
0757   pr->Fit("fG","NRB");
0758   //  pr->Draw("same");
0759   //  fG->Draw("same");
0760   fG->GetParameters(pG);
0761   maxE = xPeak(fG,pG[1]-0.1,pG[1]+0.1);
0762   cout<<maxE<<endl;
0763   //  delete p;
0764   return ((pG[1]<xmin+0.05||pG[1]>xmax-0.05)? false : true);
0765 }
0766
0767 //   ====================================================================================
0768 //   ************************************************************************************
0769 //   ====================================================================================
0770
0771 #include "hLabHelper.C"
0772
0773 #include "hcalHelper.C"
0774
0775 #include "hcalTree.C"
0776
0777 #include "tileHelper.C"
0778
0779 #include "tileTree.C"
0780