epic-lfhcal-tbana/OldStructure/FittingHeader.h

0001 /***********************************************************************************************
0002 *** provided by,                                                                          ******
0003 ***     Friederike Bock, fbock@cern.ch                                                    ******
0004 ************************************************************************************************/
0005
0006 #ifndef FITTINGGENERAL
0007 #define FITTINGGENERAL
0008
0009   Int_t nPeaksMultGauss = 0;
0010
0011   //__________________________________________________________________________________________________________
0012   // multiple Gauss fit
0013   //__________________________________________________________________________________________________________
0014   Double_t multGauss(Double_t * x, Double_t * par){
0015     Double_t result = par[0] + par[1] * x[0];
0016     for (Int_t p = 0; p < nPeaksMultGauss; p++){
0017       Double_t norm = par[3 * p +2]; // "height" or "area"
0018       Double_t mean = par[3 * p + 3];
0019       Double_t sigma = par[3 * p + 4];
0020       #if defined(__PEAKS_C_FIT_AREAS__)
0021         norm /= sigma * (TMath::Sqrt(TMath::TwoPi())); // "area"
0022       #endif                                               /* defined(__PEAKS_C_FIT_AREAS__) */
0023         result += norm * TMath::Gaus(x[0], mean, sigma);
0024     }
0025     return result;
0026   }
0027
0028   //__________________________________________________________________________________________________________
0029   // Fit noise nicely
0030   //__________________________________________________________________________________________________________
0031   Bool_t FitNoise (TH1D* histo, TF1* &fit, Double_t &mean, Double_t &meanErr, Double_t &sigma, Double_t &sigmaErr, Int_t cb, Int_t cc, TString baseName, TString nameGain, Int_t verbosity = 0){
0032
0033     if (!histo) return kFALSE;
0034     if (!(histo->GetEntries() > 0)) return kFALSE;
0035
0036     fit = new TF1(Form("%s_B%d_C%02d",baseName.Data(),cb,cc), "gaus", 0, 100);
0037     // preset parames
0038     fit->SetParameter(0, histo->GetMaximum()); // Amplitude
0039     fit->SetParameter(1, histo->GetMean());    // Mean
0040     fit->SetParameter(2, histo->GetStdDev());  // Sigma
0041     histo->Fit(fit,"QRMNE0");
0042     if (verbosity > 0)  std::cout << "Channel b: " << cb << "\t c: " << cc << "\t" << nameGain.Data() << "  mean: " << fit->GetParameter(1) << "\t width: " << fit->GetParameter(2) << std::endl;
0043     mean = fit->GetParameter(1);
0044     meanErr = fit->GetParError(1);
0045     sigma = fit->GetParameter(2);
0046     sigmaErr = fit->GetParError(2);
0047     return kTRUE;
0048   }
0049
0050   //__________________________________________________________________________________________________________
0051   // Fit noise nicely
0052   //__________________________________________________________________________________________________________
0053   Bool_t FitNoiseWithBG (TH1D* histo, TF1* &fit, Double_t &mean, Double_t &meanErr, Double_t &sigma, Double_t &sigmaErr, Int_t cb, Int_t cc, TString baseName, TString nameGain, Int_t verbosity = 0){
0054
0055     if (!histo) return kFALSE;
0056     if (!(histo->GetEntries() > 0)) return kFALSE;
0057
0058     fit = new TF1(Form("%s_B%d_C%02d",baseName.Data(),cb,cc), "gaus", 0, 100);
0059     // preset parames
0060     Double_t largestContent     = 0;
0061     Int_t minBin = histo->GetXaxis()->FindBin(0.0);
0062     Int_t maxBin = histo->GetXaxis()->FindBin(100)+0.0001;
0063     for (Int_t i= minBin; i < maxBin; i++){
0064         if (largestContent < histo->GetBinContent(i)){
0065             largestContent = histo->GetBinContent(i);
0066         }
0067     }
0068
0069     fit->SetParameter(0, 0.8* largestContent); // Amplitude
0070     fit->SetParLimits(0, 0.1, largestContent); // Amplitude
0071     fit->SetParameter(1, 50);    // Mean
0072     fit->SetParLimits(1, 0, 100);    // Mean
0073     fit->SetParameter(2, 20);  // Sigma
0074     fit->SetParLimits(2, 0, 80);  // Sigma
0075     histo->Fit(fit,"QRMNE0");
0076     if (verbosity > 0)  std::cout << "Channel b: " << cb << "\t c: " << cc << "\t" << nameGain.Data() << "  mean: " << fit->GetParameter(1) << "\t width: " << fit->GetParameter(2) << std::endl;
0077
0078     histo->Fit(fit,"QRMNE0","", fit->GetParameter(1)-3*fit->GetParameter(2), fit->GetParameter(1)+0.05*fit->GetParameter(2));
0079     if (verbosity > 0)  std::cout << "2nd try Channel b: " << cb << "\t c: " << cc << "\t" << nameGain.Data() << "  mean: " << fit->GetParameter(1) << "\t width: " << fit->GetParameter(2) << std::endl;
0080
0081     mean = fit->GetParameter(1);
0082     meanErr = fit->GetParError(1);
0083     sigma = fit->GetParameter(2);
0084     sigmaErr = fit->GetParError(2);
0085     return kTRUE;
0086   }
0087
0088   //__________________________________________________________________________________________________________
0089   // Plot & Fit gain Corr
0090   //__________________________________________________________________________________________________________
0091   void FitAndPlotGainCorr (TH2D* hist, TF1* &fit, TString baseNameFit, Double_t minFX, Double_t maxFX, Double_t maxPX, Double_t maxPY,
0092                           Double_t &slope, Double_t &slopeErr, Double_t &offset, Double_t &offsetErr,
0093                           Int_t cb, Int_t cc, Int_t sl, Int_t sc,
0094                           Bool_t bDetailed, TCanvas* canvas, TString baseNameOut, Double_t textSizeRel = 0.04 ){
0095       fit = new TF1(Form("%s_B%d_C%02d",baseNameFit.Data(), cb,cc), "[0]+[1]*x", minFX, maxFX);
0096       hist->Fit(fit,"QRMNE0");
0097
0098       slope       = fit->GetParameter(1);
0099       slopeErr    = fit->GetParError(1);
0100       offset      = fit->GetParameter(0);
0101       offsetErr   = fit->GetParError(0);
0102
0103       if (!bDetailed) return;
0104       canvas->cd();
0105         SetStyleHistoTH2ForGraphs( hist, hist->GetXaxis()->GetTitle(), hist->GetYaxis()->GetTitle(), 0.85*textSizeRel, textSizeRel, 0.85*textSizeRel, textSizeRel,0.9, 1.25);
0106         SetStyleFit(fit , 0, maxPX, 7, 7, kBlack);
0107         hist->GetXaxis()->SetRangeUser(0,maxPX);
0108         hist->GetYaxis()->SetRangeUser(0,maxPY);
0109         hist->Draw("colz");
0110         fit->Draw("same");
0111
0112         TLegend* legend = GetAndSetLegend2( 0.42, 0.15, 0.95, 0.3,textSizeRel, 1, Form("CAEN B %d, C %d, Stack L %d, C%d",cb, cc, sl, sc), 42,0.2);
0113         legend->AddEntry(fit, "f(x) = a#upoint x + b", "l");
0114         legend->AddEntry((TObject*)0, Form("a = %2.2f #pm %2.2f",slope, slopeErr ) , " ");
0115         legend->AddEntry((TObject*)0, Form("b = %2.2f #pm %2.2f",offset, offsetErr ) , " ");
0116         legend->Draw();
0117       canvas->SaveAs(Form("%s_B%d_C%02d.pdf", baseNameOut.Data(), cb,cc));
0118       return;
0119   }
0120
0121
0122   double langaufun(double *x, double *par) {
0123
0124     //Fit parameters:
0125     //par[0]=Width (scale) parameter of Landau density
0126     //par[1]=Most Probable (MP, location) parameter of Landau density
0127     //par[2]=Total area (integral -inf to inf, normalization constant)
0128     //par[3]=Width (sigma) of convoluted Gaussian function
0129     //
0130     //In the Landau distribution (represented by the CERNLIB approximation),
0131     //the maximum is located at x=-0.22278298 with the location parameter=0.
0132     //This shift is corrected within this function, so that the actual
0133     //maximum is identical to the MP parameter.
0134
0135     // Numeric constants
0136     double invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
0137     double mpshift  = -0.22278298;       // Landau maximum location
0138
0139     // Control constants
0140     double np = 100.0;      // number of convolution steps
0141     double sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
0142
0143     // Variables
0144     double xx;
0145     double mpc;
0146     double fland;
0147     double sum = 0.0;
0148     double xlow,xupp;
0149     double step;
0150     double i;
0151
0152
0153     // MP shift correction
0154     mpc = par[1] - mpshift * par[0];
0155
0156     // Range of convolution integral
0157     xlow = x[0] - sc * par[3];
0158     xupp = x[0] + sc * par[3];
0159
0160     step = (xupp-xlow) / np;
0161
0162     // Convolution integral of Landau and Gaussian by sum
0163     for(i=1.0; i<=np/2; i++) {
0164         xx = xlow + (i-.5) * step;
0165         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
0166         sum += fland * TMath::Gaus(x[0],xx,par[3]);
0167
0168         xx = xupp - (i-.5) * step;
0169         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
0170         sum += fland * TMath::Gaus(x[0],xx,par[3]);
0171     }
0172
0173     return (par[2] * step * sum * invsq2pi / par[3]);
0174   }
0175
0176
0177
0178   TF1 *langaufit(TH1D *his, double *fitrange, double *startvalues, double *parlimitslo, double *parlimitshi, double *fitparams, double *fiterrors, double *ChiSqr, int *NDF, Int_t verbosity = 0)
0179   {
0180     // Once again, here are the Landau * Gaussian parameters:
0181     //   par[0]=Width (scale) parameter of Landau density
0182     //   par[1]=Most Probable (MP, location) parameter of Landau density
0183     //   par[2]=Total area (integral -inf to inf, normalization constant)
0184     //   par[3]=Width (sigma) of convoluted Gaussian function
0185     //
0186     // Variables for langaufit call:
0187     //   his             histogram to fit
0188     //   fitrange[2]     lo and hi boundaries of fit range
0189     //   startvalues[4]  reasonable start values for the fit
0190     //   parlimitslo[4]  lower parameter limits
0191     //   parlimitshi[4]  upper parameter limits
0192     //   fitparams[4]    returns the final fit parameters
0193     //   fiterrors[4]    returns the final fit errors
0194     //   ChiSqr          returns the chi square
0195     //   NDF             returns ndf
0196
0197     int i;
0198     char FunName[100];
0199
0200     sprintf(FunName,"Fitfcn_%s",his->GetName());
0201
0202     TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0203     if (ffitold) delete ffitold;
0204
0205     TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
0206     ffit->SetParameters(startvalues);
0207     ffit->SetParNames("Width","MP","Area","GSigma");
0208
0209     for (i=0; i<4; i++) {
0210       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0211     }
0212
0213     TString fitOption = "QRLMNE0";
0214     if (verbosity > 1)
0215       fitOption = "RLMNE0";
0216
0217     his->Fit(FunName,fitOption);   // fit within specified range, use ParLimits, do not plot
0218
0219     ffit->GetParameters(fitparams);    // obtain fit parameters
0220     for (i=0; i<4; i++) {
0221       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
0222     }
0223     ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
0224     NDF[0] = ffit->GetNDF();           // obtain ndf
0225
0226     return (ffit);              // return fit function
0227
0228   }
0229
0230
0231   int langaupro(double *params, double &maxx, double &FWHM) {
0232
0233     // Searches for the location (x value) at the maximum of the
0234     // Landau-Gaussian convolute and its full width at half-maximum.
0235     //
0236     // The search is probably not very efficient, but it's a first try.
0237     double p,x,fy,fxr,fxl;
0238     double step;
0239     double l,lold;
0240     int i = 0;
0241     int MAXCALLS = 10000;
0242
0243     // Search for maximum
0244     p = params[1] - 0.1 * params[0];
0245     step = 0.05 * params[0];
0246     lold = -2.0;
0247     l    = -1.0;
0248
0249     while ( (l != lold) && (i < MAXCALLS) ) {
0250         i++;
0251         lold = l;
0252         x = p + step;
0253         l = langaufun(&x,params);
0254         if (l < lold)
0255           step = -step/10;
0256         p += step;
0257     }
0258
0259     if (i == MAXCALLS)
0260         return (-1);
0261     maxx = x;
0262     fy = l/2;
0263
0264     // Search for right x location of fy
0265     p = maxx + params[0];
0266     step = params[0];
0267     lold = -2.0;
0268     l    = -1e300;
0269     i    = 0;
0270
0271     while ( (l != lold) && (i < MAXCALLS) ) {
0272         i++;
0273         lold = l;
0274         x = p + step;
0275         l = TMath::Abs(langaufun(&x,params) - fy);
0276         if (l > lold)
0277           step = -step/10;
0278         p += step;
0279     }
0280
0281     if (i == MAXCALLS)
0282         return (-2);
0283     fxr = x;
0284
0285     // Search for left x location of fy
0286     p = maxx - 0.5 * params[0];
0287     step = -params[0];
0288     lold = -2.0;
0289     l    = -1e300;
0290     i    = 0;
0291
0292     while ( (l != lold) && (i < MAXCALLS) ) {
0293       i++;
0294       lold = l;
0295       x = p + step;
0296       l = TMath::Abs(langaufun(&x,params) - fy);
0297       if (l > lold)
0298         step = -step/10;
0299       p += step;
0300     }
0301
0302     if (i == MAXCALLS)
0303       return (-3);
0304
0305
0306     fxl = x;
0307     FWHM = fxr - fxl;
0308     return (0);
0309   }
0310
0311
0312 #endif