epic-lfhcal-tbana/OldStructure/LanGauData.C

0001 //This code fits the noise subtracted data to a Landau-Gauss fit
0002 //This creates a per channel (1-64) TGraph of the MPV values and writes it to a root file
0003 #include <TCanvas.h>
0004 #include <TH1F.h>
0005 #include <TF1.h>
0006 #include <TRandom3.h>
0007 #include <TMath.h>
0008
0009 const int gMaxChannels = 64;
0010 const int gMaxBoard    = 1;
0011
0012 //void CheckFitStatus(TF1* fitFunction, TH1* histogram, int numpara)
0013 //{
0014 //    TF1 *result = (TF1*)histogram->Fit(fitFunction, "S", ""); // "S" for silent fit
0015 //    // Check the fit status
0016 ////    int fitStatus = result->Status();
0017 //    if (fitStatus == 0)
0018 //    {
0019 //        for(int i=0; i<numpara; i++ )
0020 //        {
0021 //            // Fit was successful
0022 //            fitval = fitFunction->GetParameter(i);
0023 //        }
0024 //    }
0025 //    else {
0026 //        // Fit failed
0027 //        std::cout << "Fit failed. Fit status: " << fitStatus << std::endl;
0028 //
0029 //        // Handle the failure or take appropriate action
0030 //    }
0031 //}
0032
0033 void FitWithParameterLimitExpansion(int maxIterations, float increment, TF1* fitFunction, TH1* histogram) {
0034
0035     histogram->Fit(fitFunction, "Q");
0036
0037     int numFitParameters = fitFunction->GetNpar();
0038
0039     for(int i=0; i < numFitParameters; i++)
0040     {
0041         // Access the fit parameter 'Mean' and its limits
0042         double fitval = fitFunction->GetParameter(i);
0043         double fitvalerr = fitFunction->GetParError(i);
0044         double fitvalLowerLimit, fitvalUpperLimit;
0045         fitFunction->GetParLimits(i, fitvalLowerLimit, fitvalUpperLimit);
0046
0047         // Perform an iterative fit with parameter limit expansion
0048         for (int j = 0; j < maxIterations; j++) {
0049             // Check if the fit parameter is at its upper limit
0050             if (fitval == fitvalUpperLimit) {
0051                 // Increase the upper limit by a small increment (you can adjust this)
0052                 fitvalUpperLimit += 0.1;
0053
0054                 // Set the new parameter limits
0055                 fitFunction->SetParLimits(i, fitvalLowerLimit, fitvalUpperLimit);
0056
0057                 // Perform the fit again
0058                 histogram->Fit(fitFunction, "Q");
0059
0060                 // Update the fit parameter and chisq
0061                 fitval = fitFunction->GetParameter(i);
0062                 //                chisq = fitsnr->GetChisquare();
0063                 //                ndf = fitsnr->GetNDF();
0064             } else {
0065                 // If the parameter is not at its limit, exit the loop
0066                 break;
0067             }
0068         }
0069         if (fitval == fitvalUpperLimit) {
0070             printf("Fit parameteris still at its upper limit after iterations.\n");
0071         } else {
0072             printf("Fit parameter is within its limits.\n");
0073         }
0074
0075 //        CheckFitStatus(fitFunction, histogram, numFitParameters);
0076
0077         // Draw the histogram
0078         histogram->Draw();
0079
0080         //    c1->Update();
0081     }
0082 }
0083
0084
0085 double langaufun(double *x, double *par) {
0086
0087    //Fit parameters:
0088    //par[0]=Width (scale) parameter of Landau density
0089    //par[1]=Most Probable (MP, location) parameter of Landau density
0090    //par[2]=Total area (integral -inf to inf, normalization constant)
0091    //par[3]=Width (sigma) of convoluted Gaussian function
0092    //
0093    //In the Landau distribution (represented by the CERNLIB approximation),
0094    //the maximum is located at x=-0.22278298 with the location parameter=0.
0095    //This shift is corrected within this function, so that the actual
0096    //maximum is identical to the MP parameter.
0097
0098       // Numeric constants
0099       double invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
0100       double mpshift  = -0.22278298;       // Landau maximum location
0101
0102       // Control constants
0103       double np = 100.0;      // number of convolution steps
0104       double sc =   5.0;      // convolution extends to +-sc Gaussian sigmas
0105
0106       // Variables
0107       double xx;
0108       double mpc;
0109       double fland;
0110       double sum = 0.0;
0111       double xlow,xupp;
0112       double step;
0113       double i;
0114
0115
0116       // MP shift correction
0117       mpc = par[1] - mpshift * par[0];
0118
0119       // Range of convolution integral
0120       xlow = x[0] - sc * par[3];
0121       xupp = x[0] + sc * par[3];
0122
0123       step = (xupp-xlow) / np;
0124
0125       // Convolution integral of Landau and Gaussian by sum
0126       for(i=1.0; i<=np/2; i++) {
0127          xx = xlow + (i-.5) * step;
0128          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
0129          sum += fland * TMath::Gaus(x[0],xx,par[3]);
0130
0131          xx = xupp - (i-.5) * step;
0132          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
0133          sum += fland * TMath::Gaus(x[0],xx,par[3]);
0134       }
0135
0136       return (par[2] * step * sum * invsq2pi / par[3]);
0137 }
0138
0139
0140
0141 TF1 *langaufit(TH1F *his, double *fitrange, double *startvalues, double *parlimitslo, double *parlimitshi, double *fitparams, double *fiterrors, double *ChiSqr, int *NDF)
0142 {
0143    // Once again, here are the Landau * Gaussian parameters:
0144    //   par[0]=Width (scale) parameter of Landau density
0145    //   par[1]=Most Probable (MP, location) parameter of Landau density
0146    //   par[2]=Total area (integral -inf to inf, normalization constant)
0147    //   par[3]=Width (sigma) of convoluted Gaussian function
0148    //
0149    // Variables for langaufit call:
0150    //   his             histogram to fit
0151    //   fitrange[2]     lo and hi boundaries of fit range
0152    //   startvalues[4]  reasonable start values for the fit
0153    //   parlimitslo[4]  lower parameter limits
0154    //   parlimitshi[4]  upper parameter limits
0155    //   fitparams[4]    returns the final fit parameters
0156    //   fiterrors[4]    returns the final fit errors
0157    //   ChiSqr          returns the chi square
0158    //   NDF             returns ndf
0159
0160    int i;
0161    char FunName[100];
0162
0163    sprintf(FunName,"Fitfcn_%s",his->GetName());
0164
0165    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
0166    if (ffitold) delete ffitold;
0167
0168    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
0169    ffit->SetParameters(startvalues);
0170    ffit->SetParNames("Width","MP","Area","GSigma");
0171
0172    for (i=0; i<4; i++) {
0173       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
0174    }
0175
0176    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
0177
0178    ffit->GetParameters(fitparams);    // obtain fit parameters
0179    for (i=0; i<4; i++) {
0180       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
0181    }
0182    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
0183    NDF[0] = ffit->GetNDF();           // obtain ndf
0184
0185    return (ffit);              // return fit function
0186
0187 }
0188
0189
0190 int langaupro(double *params, double &maxx, double &FWHM) {
0191
0192    // Searches for the location (x value) at the maximum of the
0193    // Landau-Gaussian convolute and its full width at half-maximum.
0194    //
0195    // The search is probably not very efficient, but it's a first try.
0196
0197    double p,x,fy,fxr,fxl;
0198    double step;
0199    double l,lold;
0200    int i = 0;
0201    int MAXCALLS = 10000;
0202
0203
0204    // Search for maximum
0205
0206    p = params[1] - 0.1 * params[0];
0207    step = 0.05 * params[0];
0208    lold = -2.0;
0209    l    = -1.0;
0210
0211
0212    while ( (l != lold) && (i < MAXCALLS) ) {
0213       i++;
0214
0215       lold = l;
0216       x = p + step;
0217       l = langaufun(&x,params);
0218
0219       if (l < lold)
0220          step = -step/10;
0221
0222       p += step;
0223    }
0224
0225    if (i == MAXCALLS)
0226       return (-1);
0227
0228    maxx = x;
0229
0230    fy = l/2;
0231
0232
0233    // Search for right x location of fy
0234
0235    p = maxx + params[0];
0236    step = params[0];
0237    lold = -2.0;
0238    l    = -1e300;
0239    i    = 0;
0240
0241
0242    while ( (l != lold) && (i < MAXCALLS) ) {
0243       i++;
0244
0245       lold = l;
0246       x = p + step;
0247       l = TMath::Abs(langaufun(&x,params) - fy);
0248
0249       if (l > lold)
0250          step = -step/10;
0251
0252       p += step;
0253    }
0254
0255    if (i == MAXCALLS)
0256       return (-2);
0257
0258    fxr = x;
0259
0260
0261    // Search for left x location of fy
0262
0263    p = maxx - 0.5 * params[0];
0264    step = -params[0];
0265    lold = -2.0;
0266    l    = -1e300;
0267    i    = 0;
0268
0269    while ( (l != lold) && (i < MAXCALLS) ) {
0270       i++;
0271
0272       lold = l;
0273       x = p + step;
0274       l = TMath::Abs(langaufun(&x,params) - fy);
0275
0276       if (l > lold)
0277          step = -step/10;
0278
0279       p += step;
0280    }
0281
0282    if (i == MAXCALLS)
0283       return (-3);
0284
0285
0286    fxl = x;
0287
0288    FWHM = fxr - fxl;
0289    return (0);
0290 }
0291
0292 void LanGauData(const char* datafile, int runnumber, TString outputDir    = "LanGauFitData/")
0293 {
0294
0295
0296     outputDir = Form("%s/Run%05d",outputDir.Data(), runnumber );
0297     gSystem->Exec("mkdir -p "+outputDir);
0298     // Create a ROOT canvas to display the histogram and fit results
0299     TCanvas *c1 = new TCanvas("c1", "LanGau Fit Example", 800, 600);
0300
0301     //Get the data file with subtracted data
0302     TFile* fdatsub = new TFile(datafile,"READ");
0303
0304     // Create a 2D array of histograms
0305 //    TH2F *histHGLG[gMaxBoard][gMaxChannels];
0306     TH1F* histLG[gMaxBoard][gMaxChannels];
0307     TH1F* histHG[gMaxBoard][gMaxChannels];
0308     TH1F* newhistHG[gMaxBoard][gMaxChannels];
0309
0310     //Get the noise subtracted files at the right voltage
0311     for (Int_t j = 0; j < gMaxBoard; j++){
0312         for (Int_t i = 0; i < gMaxChannels; i++){
0313             //noise histogram
0314             newhistHG[j][i] = (TH1F*)fdatsub->Get(Form("sh_HG_B%d_C%02d",j,i));
0315             newhistHG[j][i]->RebinX(5);
0316             histHG[j][i] = (TH1F*)newhistHG[j][i]->Clone();
0317             for(int k=1; k < 701; k++)
0318             {
0319                 histHG[j][i]->SetBinContent(k,newhistHG[j][i]->GetBinContent(k+10));
0320             }
0321             histLG[j][i] = (TH1F*)fdatsub->Get(Form("sh_LG_B%d_C%02d",j,i));
0322
0323         }
0324     }
0325
0326     //Now we want to fit these with Landau-Gauss fit and store the fit values in a TGraph with the errors representing the fit errors.
0327     //The TGraph will have 64 values corresponding to the channels on one digitizer board
0328     //The parameters of the LG fit are
0329         //   par[0]=Width (scale) parameter of Landau density
0330         //   par[1]=Most Probable (MP, location) parameter of Landau density
0331         //   par[2]=Total area (integral -inf to inf, normalization constant)
0332         //   par[3]=Width (sigma) of convoluted Gaussian function
0333     //Additionally, we want to keep track of the chisq and ndf for the fits
0334     TGraphAsymmErrors* mpv_gr = new TGraphAsymmErrors();;
0335     TGraphAsymmErrors* scale_gr = new TGraphAsymmErrors();;
0336     TGraphAsymmErrors* area_gr = new TGraphAsymmErrors();;
0337     TGraphAsymmErrors* sigma_gr = new TGraphAsymmErrors();;
0338     TGraph* chisq_gr = new TGraph();
0339     TGraph* ndf_gr = new TGraph();
0340
0341     mpv_gr->GetYaxis()->SetTitle("MPV");
0342     mpv_gr->SetName("MPV");
0343     scale_gr->GetYaxis()->SetTitle("Scale");
0344     scale_gr->SetName("Scale");
0345     area_gr->GetYaxis()->SetTitle("Area");
0346     area_gr->SetName("Area");
0347     sigma_gr->GetYaxis()->SetTitle("Sigma");
0348     sigma_gr->SetName("Sigma");
0349
0350     mpv_gr->GetXaxis()->SetTitle("Channel");
0351     scale_gr->GetXaxis()->SetTitle("Channel");
0352     area_gr->GetXaxis()->SetTitle("Channel");
0353     sigma_gr->GetXaxis()->SetTitle("Channel");
0354
0355     TH1D* mpv_hist = new TH1D("mpv_hist","MPV fit values",500,0,500);
0356     TH1* scale_hist = new TH1D("scale_hist","Scale fit values",500,0,500);
0357     TH1* area_hist = new TH1D("area_hist","Area fit values",500,0,1000);
0358     TH1* sigma_hist = new TH1D("sigma_hist","Sigma fit values",100,0,100);
0359     TH1* chisq_hist = new TH1D("chisq_hist","ChiSQ",100,0,100);
0360     TH1* ndf_hist=new TH1D("ndf_hist","NDF",100,0,100);
0361 //
0362     //Defining fit parameter names
0363     Double_t mpv;
0364     Double_t mpv_err;
0365     float_t chisq;
0366     Int_t ndf;
0367     float_t area;
0368     float_t area_err;
0369     float_t sigma;
0370     float_t sigma_err;
0371     float_t scale;
0372     float_t scale_err;
0373
0374     //Defining the x array of channels
0375 //    Double_t x[gMaxChannels];
0376 //    for(int i =0; i < 64; i++)
0377 //    {
0378 //      x[i]=i;
0379 //    }
0380 //
0381     //Now loop over histograms in the file and run the fitter
0382      for (Int_t j = 0; j < gMaxBoard; j++)
0383     {
0384         for (Int_t i = 0; i < gMaxChannels; i++)
0385         {
0386
0387             double fr[2];
0388             double sv[4], pllo[4], plhi[4], fp[4], fpe[4];
0389             double chisqr;
0390             int    ndf;
0391             int binmax = histHG[j][i]->GetMaximumBin(); double x = histHG[j][i]->GetXaxis()->GetBinCenter(binmax);
0392
0393             //Set low and high parameters for the fit
0394             fr[0]=0.15*histHG[j][i]->GetMean();
0395             fr[1]=2*histHG[j][i]->GetMean();
0396
0397             pllo[0] = 0;
0398             pllo[1] = x - histHG[j][i]->GetMean();//MPV
0399             pllo[2] = 0.9*(histHG[j][i]->Integral());
0400             pllo[3] = 2;
0401
0402             plhi[0] = 0.5*histHG[j][i]->GetStdDev(1);
0403             plhi[1] = x + histHG[j][i]->GetMean();//MPV
0404             plhi[2] = 1000*(histHG[j][i]->Integral());
0405             plhi[3] = 18;
0406
0407             //Fix start values
0408             sv[0]=0.25*histHG[j][i]->GetStdDev(1) ; sv[1]=x; sv[2]=1.5*(histHG[j][i]->Integral()); sv[3]=(pllo[3]+plhi[3])/2;
0409
0410             //Perform the fit first time
0411             TF1 *fitsnr = langaufit(histHG[j][i],fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
0412
0413             scale = fp[0]; scale_err = fpe[0];
0414             mpv = fp[1]; mpv_err = fpe[1];
0415             area = fp[2]; area_err = fpe[2];
0416             sigma = fp[3]; sigma_err = fpe[3];
0417
0418             FitWithParameterLimitExpansion(30, 0.1, fitsnr, histHG[j][i]);
0419
0420
0421
0422             chisq = fitsnr->GetChisquare();
0423             ndf = fitsnr->GetNDF();
0424
0425             //Add values to the graph
0426             int pointIndex = mpv_gr->GetN(); // Get the current number of points
0427             mpv_gr->SetPoint(pointIndex, i, mpv);
0428             mpv_gr->SetPointError(pointIndex, 0, 0, mpv_err, mpv_err);
0429
0430             scale_gr->SetPoint(pointIndex, i, scale);
0431             scale_gr->SetPointError(pointIndex, 0, 0, scale_err, scale_err);
0432
0433             area_gr->SetPoint(pointIndex, i, area);
0434             area_gr->SetPointError(pointIndex, 0, 0, area_err, area_err);
0435
0436             sigma_gr->SetPoint(pointIndex, i, sigma);
0437             sigma_gr->SetPointError(pointIndex, 0, 0, sigma_err, sigma_err);
0438
0439             //Fill histogram with fit values
0440             mpv_hist->Fill(mpv);
0441             scale_hist->Fill(scale);
0442             area_hist->Fill(area);
0443             sigma_hist->Fill(sigma);
0444             chisq_hist->Fill(chisq);
0445             ndf_hist->Fill(ndf);
0446
0447         }
0448     }
0449
0450 //histHG[0][1]->Draw();
0451 ////            fitsnr->Draw("lsame");
0452 //            gStyle->SetOptStat(1111);
0453 //            gStyle->SetOptFit(111);
0454 //            gStyle->SetLabelSize(0.03,"x");
0455 //            gStyle->SetLabelSize(0.03,"y");
0456
0457     TFile* fileOutput = new TFile(Form("%s/%d.root","LanGauFitData/",runnumber),"RECREATE");
0458     for (Int_t j = 0; j < gMaxBoard; j++)
0459     {
0460             mpv_hist->Write();
0461             scale_hist->Write();
0462             area_hist->Write();
0463             sigma_hist->Write();
0464             chisq_hist->Write();
0465             ndf_hist->Write();
0466
0467             //Write TGraphs to output file
0468             mpv_gr->Write();
0469             scale_gr->Write();
0470             area_gr->Write();
0471             sigma_gr->Write();
0472
0473     }
0474
0475     histHG[0][24]->Draw();
0476 }
0477
0478
0479
0480
0481