fluka_gfortran7_hit/helium/langaus.C

 /// \ingroup tutorial_fit
 /// \notebook
 /// Convoluted Landau and Gaussian Fitting Function
 ///         (using ROOT's Landau and Gauss functions)
 ///
 ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
 ///
 ///  to execute this example, do:
 ///
 /// ~~~{.cpp}
 ///  root > .x langaus.C
 /// ~~~
 ///
 /// or
 ///
 /// ~~~{.cpp}
 ///  root > .x langaus.C++
 /// ~~~
 ///
 /// \macro_image
 /// \macro_output
 /// \macro_code
 ///
 /// \authors H.Pernegger, Markus Friedl
 
 #include "TH1.h"
 #include "TF1.h"
 #include "TROOT.h"
 #include "TStyle.h"
 #include "TMath.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <iostream>
#include <fstream>
#include <iomanip> 
#include <cstdlib>
#include <cmath>
#include <string.h>
#include <TLorentzVector.h>
#include <vector>
#include <TApplication.h>
#include <TAxis.h>
#include <stdio.h>
#include <stdlib.h>
#include <TFile.h>
#include <TTree.h>
 
 Double_t langaufun(Double_t *x, Double_t *par) {
 
    //Fit parameters:
    //par[0]=Width (scale) parameter of Landau density
    //par[1]=Most Probable (MP, location) parameter of Landau density
    //par[2]=Total area (integral -inf to inf, normalization constant)
    //par[3]=Width (sigma) of convoluted Gaussian function
    //
    //In the Landau distribution (represented by the CERNLIB approximation),
    //the maximum is located at x=-0.22278298 with the location parameter=0.
    //This shift is corrected within this function, so that the actual
    //maximum is identical to the MP parameter.
 
       // Numeric constants
       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
       Double_t mpshift  = -0.22278298;       // Landau maximum location
 
       // Control constants
       Double_t np = 200.0;      // number of convolution steps
       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
 
       // Variables
       Double_t xx;
       Double_t mpc;
       Double_t fland;
       Double_t sum = 0.0;
       Double_t xlow,xupp;
       Double_t step;
       Double_t i;
 
 
       // MP shift correction
       mpc = par[1] - mpshift * par[0];
 
       // Range of convolution integral
       xlow = x[0] - sc * par[3];
       xupp = x[0] + 2*sc * par[3];
 
       step = (xupp-xlow) / np;
 
       // Convolution integral of Landau and Gaussian by sum
       for(i=1.0; i<=np/2; i++) {
          xx = xlow + (i-.5) * step;
          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
          sum += fland * TMath::Gaus(x[0],xx,par[3]);
 
          xx = xupp - (i-.5) * step;
          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
          sum += fland * TMath::Gaus(x[0],xx,par[3]);
       }
 
       return (par[2] * step * sum * invsq2pi / par[3]);
 }
 
 
 
 TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
 {
    // Once again, here are the Landau * Gaussian parameters:
    //   par[0]=Width (scale) parameter of Landau density
    //   par[1]=Most Probable (MP, location) parameter of Landau density
    //   par[2]=Total area (integral -inf to inf, normalization constant)
    //   par[3]=Width (sigma) of convoluted Gaussian function
    //
    // Variables for langaufit call:
    //   his             histogram to fit
    //   fitrange[2]     lo and hi boundaries of fit range
    //   startvalues[4]  reasonable start values for the fit
    //   parlimitslo[4]  lower parameter limits
    //   parlimitshi[4]  upper parameter limits
    //   fitparams[4]    returns the final fit parameters
    //   fiterrors[4]    returns the final fit errors
    //   ChiSqr          returns the chi square
    //   NDF             returns ndf
 
    Int_t i;
    Char_t FunName[100];
 
    sprintf(FunName,"Fitfcn_%s",his->GetName());
 
    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
    if (ffitold) delete ffitold;
 
    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
    ffit->SetParameters(startvalues);
    ffit->SetParNames("Width","MP","Area","GSigma");
 
    for (i=0; i<4; i++) {
       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
    }
 
    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
 
    ffit->GetParameters(fitparams);    // obtain fit parameters
    for (i=0; i<4; i++) {
       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
    }
    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
    NDF[0] = ffit->GetNDF();           // obtain ndf
 
    return (ffit);              // return fit function
 
 }
 
 
 Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
 
    // Seaches for the location (x value) at the maximum of the
    // Landau-Gaussian convolute and its full width at half-maximum.
    //
    // The search is probably not very efficient, but it's a first try.
 
    Double_t p,x,fy,fxr,fxl;
    Double_t step;
    Double_t l,lold;
    Int_t i = 0;
    Int_t MAXCALLS = 10000;
 
 
    // Search for maximum
 
    p = params[1] - 0.1 * params[0];
    step = 0.05 * params[0];
    lold = -2.0;
    l    = -1.0;
 
 
    while ( (l != lold) && (i < MAXCALLS) ) {
       i++;
 
       lold = l;
       x = p + step;
       l = langaufun(&x,params);
 
       if (l < lold)
          step = -step/10;
 
       p += step;
    }
 
    if (i == MAXCALLS)
       return (-1);
 
    maxx = x;
 
    fy = l/2;
 
 
    // Search for right x location of fy
 
    p = maxx + params[0];
    step = params[0];
    lold = -2.0;
    l    = -1e300;
    i    = 0;
 
 
    while ( (l != lold) && (i < MAXCALLS) ) {
       i++;
 
       lold = l;
       x = p + step;
       l = TMath::Abs(langaufun(&x,params) - fy);
 
       if (l > lold)
          step = -step/10;
 
       p += step;
    }
 
    if (i == MAXCALLS)
       return (-2);
 
    fxr = x;
 
 
    // Search for left x location of fy
 
    p = maxx - 0.5 * params[0];
    step = -params[0];
    lold = -2.0;
    l    = -1e300;
    i    = 0;
 
    while ( (l != lold) && (i < MAXCALLS) ) {
       i++;
 
       lold = l;
       x = p + step;
       l = TMath::Abs(langaufun(&x,params) - fy);
 
       if (l > lold)
          step = -step/10;
 
       p += step;
    }
 
    if (i == MAXCALLS)
       return (-3);
 
 
    fxl = x;
 
    FWHM = fxr - fxl;
    return (0);
 }
 
void langaus() {
    // Fill Histogram
   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
   */
   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
   TH1D * hSNR;

    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
   
    
    // hSNR = h_beamSignal_b0[13];
    // Fitting SNR histo
    printf("Fitting...\n");
    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
    
    // Setting fit range and start values
    Double_t fr[2];
    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
    Double_t chisqr;
    Int_t    ndf;
    TF1 *fitsnr;
    Double_t SNRPeak, SNRFWHM;
    char rootfilename[50] = "";
    char saveplotname[50] = "";
    char fout_mpv_name[50] = "";
    Double_t norm;
    
    ofstream myfile, fout_mpv;
    int j = 1;

    myfile.open ("MPVcorrection_helium0mm05.txt");
    sprintf(fout_mpv_name, "jobs%i/plots/mpv.txt",j);
     fout_mpv.open (fout_mpv_name);
    
    for (int i = 0; i<26;i++){
      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
      TFile *rootFile = new TFile(rootfilename,"OPEN");

      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
      norm = hSNR->GetEntries();
      hSNR->Scale(1/norm);
      
      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
      fr[1]=hSNR->GetMean() + 3*hSNR->GetRMS();
 
      pllo[0]=0.0006; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01;
      plhi[0]=0.01; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.05;
      sv[0]= hSNR->GetRMS()/10.;;
      //sv[1]=2.1;
      sv[2]=0.01;
      //   sv[3]=0.03;
      sv[3] = hSNR->GetRMS()/2.;
      sv[1] = hSNR->GetMean();
    
      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);

      
      langaupro(fp,SNRPeak,SNRFWHM);
      graph_x[i] = beta_helium[i];
      graph_y[i] = fp[1];
      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
      fout_mpv << graph_x[i] << " " << graph_y[i] << " " << graph_yerr[i] << endl;
      printf("Fitting done\nPlotting results...\n");
            myfile << i << " " << graph_y[i] << endl;

      // Global style settings
      gStyle->SetOptStat(1111);
      gStyle->SetOptFit(111);
      //gStyle->SetLabelSize(0.03,"x");
      //gStyle->SetLabelSize(0.03,"y");
      
      // hSNR->GetXaxis()->SetRange(0,70);
      hSNR->Draw();
      fitsnr->Draw("lsame");
      c1->Update();
      sprintf(saveplotname, "jobs%i/plots/runjob%i001_h_spratio.pdf",j,i);
      c1->SaveAs(saveplotname);
      //  c1->WaitPrimitive();
      hSNR->Delete();
      rootFile->Close();
    }
     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
    c2->cd();
    graph_1->Draw("A*");
    c2->SaveAs("MPVcorrection_helium0mm05.pdf");
    myfile.close();
    fout_mpv.close();

    
 }