test/Macros/DoubleGaussianFit.C

0001 #include <TFile.h>
0002 #include <TH1F.h>
0003 #include <TF1.h>
0004 #include <TCanvas.h>
0005 #include <TPaveText.h>
0006 #include <TStyle.h>
0007
0008 #include <sstream>
0009 #include <iostream>
0010 #include <iomanip>
0011
0012 /// Helper class holding a TPaveText for better formatting and predefined options
0013 class PaveText
0014 {
0015  public:
0016   PaveText(const double & textX = 0.7, const double & textY = 0.4 )
0017   {
0018     paveText_ = new TPaveText(textX, textY, textX+0.2, textY+0.17, "NDC");
0019   }
0020   void AddText(const TString & text)
0021   {
0022     paveText_->AddText(text);
0023   }
0024   void Draw(const TString & option)
0025   {
0026     paveText_->SetFillColor(0); // text is black on white
0027     paveText_->SetTextSize(0.03);
0028     paveText_->SetBorderSize(0);
0029     paveText_->SetTextAlign(12);
0030     paveText_->Draw(option);
0031   }
0032   void SetTextColor(const int color)
0033   {
0034     paveText_->SetTextColor(color);
0035   }
0036  protected:
0037   TPaveText * paveText_;
0038 };
0039
0040 /**
0041  * Compute the precision to give to the stream operator so that the passed number
0042  * will be printed with two significant figures.
0043  */
0044 int precision( const double & value )
0045 {
0046   // Counter gives the precision
0047   int precision = 1;
0048   int k=1;
0049   while( int(value*k) == 0 ) {
0050     k*=10;
0051     ++precision;
0052   }
0053   return precision;
0054 }
0055
0056 /// Helper function to extract and format the text for the fitted parameters
0057 void getParameters( const TF1 * func, TString & fit1, TString & fit2, TString & fit3 )
0058 {
0059   std::stringstream a;
0060
0061   double error = func->GetParError(1);
0062   a << std::setprecision(precision(error)) << std::fixed << func->GetParameter(1);
0063   fit1 += a.str() + "+-";
0064   a.str("");
0065   a << error;
0066   fit1 += a.str();
0067   a.str("");
0068
0069   error = func->GetParError(2);
0070
0071   a << std::setprecision(precision(error)) << std::fixed << func->GetParameter(2);
0072   fit2 += a.str() + "+-";
0073   a.str("");
0074   a << func->GetParError(2);
0075   fit2 += a.str();
0076   a.str("");
0077   a << std::setprecision(1) << std::fixed << func->GetChisquare();
0078   fit3 += a.str() + "/";
0079   a.str("");
0080   a << std::setprecision(0) << std::fixed << func->GetNDF();
0081   fit3 += a.str();
0082 }
0083
0084 void DoubleGaussianFit()
0085 {
0086   TFile * file_0 = new TFile("0_MuScleFit.root", "READ");
0087   TH1F * histo_0 = (TH1F*)file_0->Get("hRecBestRes_Mass");
0088   TFile * file_1 = new TFile("3_MuScleFit.root", "READ");
0089   TH1F * histo_1 = (TH1F*)file_1->Get("hRecBestRes_Mass");
0090
0091   histo_0->Rebin(2);
0092   histo_1->Rebin(2);
0093
0094   // TF1 *f1_0 = new TF1("f1", "gaus(0) + gaus(3)", 1., 5.);
0095   // f1_0->SetParameters(136., 3.096916, 0.03, 136, 3., 0.03);
0096   // f1_0->FixParameter(1, 3.096916);
0097
0098   TF1 *f1_0 = new TF1("f1", "gaus(0)", 1., 5.);
0099   f1_0->SetParameters(136., 3.096916, 0.03);
0100
0101   // TF1 *f1_1 = new TF1("f1", "gaus(0) + gaus(3)", 1., 5.);
0102   // f1_1->SetParameters(136., 3.096916, 0.03, 136, 3., 0.03);
0103   // f1_1->FixParameter(1, 3.096916);
0104   // histo_1->SetLineColor(kRed);
0105   // f1_1->SetLineColor(kRed);
0106
0107   TF1 *f1_1 = new TF1("f1", "gaus(0)", 1., 5.);
0108   f1_1->SetParameters(136., 3.096916, 0.03);
0109   histo_1->SetLineColor(kRed);
0110   f1_1->SetLineColor(kRed);
0111
0112   // TF1 *f1 = new TF1("f1", "gaus(0)", 1., 5.);
0113   // f1->SetParameters(136., 3., 0.03);
0114
0115   histo_0->Fit(f1_0);
0116   histo_1->Fit(f1_1);
0117
0118   TCanvas * canvas = new TCanvas("canvas", "canvas", 1000, 800);
0119   canvas->Draw();
0120
0121   histo_0->Draw();
0122   histo_0->GetXaxis()->SetTitle("Mass (GeV)");
0123   histo_0->GetYaxis()->SetTitle("a.u.");
0124   histo_0->SetTitle("");
0125   histo_1->Draw("SAME");
0126
0127
0128
0129   TString fit11("mean = ");
0130   TString fit12("sigma = ");
0131   TString fit13("Chi2/NDF = ");
0132   getParameters(f1_0, fit11, fit12, fit13);
0133
0134   PaveText pt1(0.45, 0.15);
0135   pt1.AddText("before:");
0136   pt1.AddText(fit11);
0137   pt1.AddText(fit12);
0138   pt1.AddText(fit13);
0139   pt1.Draw("same");
0140
0141   TString fit21("mean = ");
0142   TString fit22("sigma = ");
0143   TString fit23("Chi2/NDF = ");
0144   getParameters(f1_1, fit21, fit22, fit23);
0145
0146   PaveText pt2(0.65, 0.15);
0147   pt2.SetTextColor(2);
0148   pt2.AddText("after:");
0149   pt2.AddText(fit21);
0150   pt2.AddText(fit22);
0151   pt2.AddText(fit23);
0152   pt2.Draw("same");
0153
0154   gStyle->SetOptStat(0);
0155
0156 }