From fd30f2e535c10a642f822726a331d1114071a8ab Mon Sep 17 00:00:00 2001 From: Vratislav Chudoba Date: Tue, 7 Feb 2017 17:24:56 +0300 Subject: [PATCH] New macros added. --- FittingDemo.cxx | 102 +++++++++++++++++++++++ fitting.cxx | 122 +++++++++++++++++++++++++++ fittingBetter.cxx | 205 ++++++++++++++++++++++++++++++++++++++++++++++ makefile | 0 opentrees.cxx | 142 ++++++++++++++++++++++++++++++++ 5 files changed, 571 insertions(+) create mode 100644 FittingDemo.cxx create mode 100644 fitting.cxx create mode 100644 fittingBetter.cxx create mode 100644 makefile create mode 100644 opentrees.cxx diff --git a/FittingDemo.cxx b/FittingDemo.cxx new file mode 100644 index 0000000..915a6da --- /dev/null +++ b/FittingDemo.cxx @@ -0,0 +1,102 @@ +//Example for fitting signal/background. +// This example can be executed with: +// root > .x FittingDemo.C (using the CINT interpreter) +// root > .x FittingDemo.C+ (using the native complier via ACLIC) +//Author: Rene Brun + +#include "TH1.h" +#include "TMath.h" +#include "TF1.h" +#include "TLegend.h" +#include "TCanvas.h" + +// Quadratic background function +Double_t background(Double_t *x, Double_t *par) { + return par[0] + par[1]*x[0] + par[2]*x[0]*x[0]; +} + + +// Lorenzian Peak function +Double_t lorentzianPeak(Double_t *x, Double_t *par) { + return (0.5*par[0]*par[1]/TMath::Pi()) / + TMath::Max( 1.e-10,(x[0]-par[2])*(x[0]-par[2]) + + .25*par[1]*par[1]); +} + +// Sum of background and peak function +Double_t fitFunction(Double_t *x, Double_t *par) { + return background(x,par) + lorentzianPeak(x,&par[3]); +} + +void FittingDemo() { + //Bevington Exercise by Peter Malzacher, modified by Rene Brun + + const int nBins = 60; + + Double_t data[nBins] = { 6, 1,10,12, 6,13,23,22,15,21, + 23,26,36,25,27,35,40,44,66,81, + 75,57,48,45,46,41,35,36,53,32, + 40,37,38,31,36,44,42,37,32,32, + 43,44,35,33,33,39,29,41,32,44, + 26,39,29,35,32,21,21,15,25,15}; + TCanvas *c1 = new TCanvas("c1","Fitting Demo",10,10,700,500); + c1->SetFillColor(33); + c1->SetFrameFillColor(41); + c1->SetGrid(); + + TH1F *histo = new TH1F("histo", + "Lorentzian Peak on Quadratic Background",60,0,3); + histo->SetMarkerStyle(21); + histo->SetMarkerSize(0.8); + histo->SetStats(0); + + for(int i=0; i < nBins; i++) histo->SetBinContent(i+1,data[i]); + + // create a TF1 with the range from 0 to 3 and 6 parameters + TF1 *fitFcn = new TF1("fitFcn",fitFunction,0,3,6); + fitFcn->SetNpx(500); + fitFcn->SetLineWidth(4); + fitFcn->SetLineColor(kMagenta); + + // first try without starting values for the parameters + // This defaults to 1 for each param. + // this results in an ok fit for the polynomial function + // however the non-linear part (lorenzian) does not + // respond well. + fitFcn->SetParameters(1,1,1,1,1,1); + histo->Fit("fitFcn","0"); + + // second try: set start values for some parameters + fitFcn->SetParameter(4,0.2); // width + fitFcn->SetParameter(5,1); // peak + + histo->Fit("fitFcn","V+","ep"); + + // improve the picture: + TF1 *backFcn = new TF1("backFcn",background,0,3,3); + backFcn->SetLineColor(kRed); + TF1 *signalFcn = new TF1("signalFcn",lorentzianPeak,0,3,3); + signalFcn->SetLineColor(kBlue); + signalFcn->SetNpx(500); + Double_t par[6]; + + // writes the fit results into the par array + fitFcn->GetParameters(par); + + backFcn->SetParameters(par); + backFcn->Draw("same"); + + signalFcn->SetParameters(&par[3]); + signalFcn->Draw("same"); + + // draw the legend + TLegend *legend=new TLegend(0.6,0.65,0.88,0.85); + legend->SetTextFont(72); + legend->SetTextSize(0.04); + legend->AddEntry(histo,"Data","lpe"); + legend->AddEntry(backFcn,"Background fit","l"); + legend->AddEntry(signalFcn,"Signal fit","l"); + legend->AddEntry(fitFcn,"Global Fit","l"); + legend->Draw(); + +} diff --git a/fitting.cxx b/fitting.cxx new file mode 100644 index 0000000..ec7e004 --- /dev/null +++ b/fitting.cxx @@ -0,0 +1,122 @@ +#include "TFile.h" +#include "TTree.h" +#include "TCanvas.h" +#include "TF1.h" +#include "TH1F.h" +#include "TStyle.h" +#include "TLegend.h" +#include "TPaveStats.h" + +void fitting() +{ + + const Int_t bins = 200; + +// TFile *fr = new TFile("simeventsNew.root"); + TFile *fr = new TFile("simeventsNew1.root"); + + TTree *tr = (TTree*)fr->Get("simevents"); + + TCanvas *c1 = new TCanvas("name1", "title of canvas 1"); + c1->Divide(2,2); + + c1->cd(1); + tr->Draw("E", "E>0"); + + c1->cd(2); + TH1F *hE = new TH1F("hE", "Interesting region", bins, 0.7, 1.3); + tr->Draw("E>>hE", "E>0.7"); + + + c1->cd(3); + TH1F *hEfitSep = new TH1F("hEfitSep", "Fit by two functions", bins, 0.7, 1.3); + hE->Draw(); + tr->Draw("E >> hEfitSep", "", ""); + + TF1 *g1 = new TF1("m1","gaus",0.9,1.05); + TF1 *g2 = new TF1("m2","gaus",1.08,1.15); + // The total is the sum of the three, each has 3 parameters + TF1 *total = new TF1("mstotal","gaus(0)+gaus(3)",0.9,1.3); + + Double_t par[6]; + + // Fit each function and add it to the list of functions + hEfitSep->Fit(g1,"R"); + hEfitSep->Fit(g2,"R+"); + + + c1->cd(4); + TH1F *hEtot = new TH1F("hEtot", "Fit by multiple function", bins, 0.7, 1.3); + hEtot->Draw(); + tr->Draw("E >> hEtot", "", ""); + + // Get the parameters from the fit + g1->GetParameters(&par[0]); + g2->GetParameters(&par[3]); + + // Use the parameters on the sum + total->SetParameters(par); + hEtot->Fit(total,"R+"); + + TCanvas *c2 = new TCanvas("name2", "title of canvas 2"); + c2->Divide(2,2); + + c2->cd(1); + TH1F *hEDecomposition = new TH1F(*hE); + hEDecomposition->Draw(); + + Double_t parTotal[6]; + total->GetParameters(&parTotal[0]); + + TF1 *gLow = new TF1("mLow","gaus",0.8,1.2); + TF1 *gHigh = new TF1("mHigh","gaus",1.0,1.3); + gHigh->SetLineColor(kGreen); + + gLow->SetParameter(0, parTotal[0]); + gLow->SetParameter(1, parTotal[1]); + gLow->SetParameter(2, parTotal[2]); + + gHigh->SetParameter(0, parTotal[3]); + gHigh->SetParameter(1, parTotal[4]); + gHigh->SetParameter(2, parTotal[5]); + + gLow->Draw("same"); + gHigh->Draw("same"); + + c2->cd(2); + hEtot->Draw(); + gHigh->Draw("same"); + + c2->cd(3); + TH1F *hInput = new TH1F(*hEDecomposition); + hInput->SetName("hInput"); + tr->Draw("Egamma >> hInput", ""); + + c2->cd(4); +// TH1F *hEtotZoom = new TH1F(*hEtot); + TH1F *hEtotZoom = (TH1F*)hEtot->Clone(); + hEtotZoom->SetName("hEtotZoom"); +// hEtotZoom-> + hEtotZoom->GetXaxis()->SetRangeUser(1.04, 1.18); + hEtotZoom->Draw(); +// total->Draw("same"); + gHigh->Draw("same"); + + Double_t maxSingle = gHigh->GetMaximumX(1.08, 1.14); + Double_t maxTotal = total->GetMaximumX(1.08, 1.14); + +// cout << maxSingle << endl; + + TLine *l1 = new TLine(maxSingle, 100., maxSingle, 800.); + l1->SetLineWidth(2); + l1->SetLineColor(gHigh->GetLineColor()); + l1->Draw(""); + + TLine *l2 = new TLine(maxTotal, 100., maxTotal, 800.); + l2->SetLineWidth(2); + l2->SetLineColor(total->GetLineColor()); + l2->Draw(""); + + + +} diff --git a/fittingBetter.cxx b/fittingBetter.cxx new file mode 100644 index 0000000..2ea03cd --- /dev/null +++ b/fittingBetter.cxx @@ -0,0 +1,205 @@ +#include "TFile.h" +#include "TTree.h" +#include "TCanvas.h" +#include "TF1.h" +#include "TH1F.h" +#include "TStyle.h" +#include "TLegend.h" +#include "TPaveStats.h" + +void fittingBetter() +{ + +// gStyle->ToggleEditor(); + + TFile *fr = new TFile("simeventsNew.root"); + + TTree *tr = (TTree*)fr->Get("simevents"); + + TCanvas *c1 = new TCanvas("name1", "title of canvas 1"); + c1->Divide(2,2); + + c1->cd(1); + tr->Draw("E", "E>0"); + + c1->cd(2); + TH1F *hE = new TH1F("hE", "Interesting region", 200, 0.7, 1.3); + tr->Draw("E>>hE", "E>0.7"); +// tr1->SetLineColor(kRed); +// tr1->Draw("E", "E>0", "same"); +// tr2->SetLineColor(kBlue); +// tr2->Draw("E", "E>0", "same"); + + + c1->cd(3); + // fr.cd(); + TH1F *hEfitSep = new TH1F("hEfitSep", "Fit by two functions", 200, 0.7, 1.3); + hE->Draw(); + tr->Draw("E >> hEfitSep", "", ""); + + TF1 *g1 = new TF1("m1","gaus",0.9,1.05); + TF1 *g2 = new TF1("m2","gaus",1.08,1.15); + g2->SetLineColor(kBlue); + // The total is the sum of the three, each has 3 parameters + TF1 *total = new TF1("mstotal","gaus(0)+gaus(3)",0.8,1.3); + + TLegend *legend = new TLegend(0.6,0.65,0.88,0.85); + legend->SetTextFont(72); + legend->SetTextSize(0.04); + legend->AddEntry(hEfitSep,"Data","lpe"); + legend->AddEntry(g1,"Broad peak","l"); + legend->AddEntry(g2,"Narrow peak","l"); + + legend->SetX1(0.15); + legend->SetX2(0.35); + + + legend->Draw(); + + cout << legend->GetX1NDC() << endl; + cout << legend->GetX1() << endl; + cout << legend->GetX2() << endl; + + Double_t par[6]; + + // Fit each function and add it to the list of functions + hEfitSep->Fit(g1,"R"); + hEfitSep->Fit(g2,"R+"); + + + c1->cd(4); + TH1F *hEtot = new TH1F("hEtot", "Fit by multiple function", 200, 0.7, 1.3); + hEtot->Draw(); + tr->Draw("E >> hEtot", "", ""); + // Get the parameters from the fit + g1->GetParameters(&par[0]); + g2->GetParameters(&par[3]); + + // Use the parameters on the sum + total->SetParameters(par); + hEtot->Fit(total,"R+"); +// gStyle->SetOptFit(); + gStyle->SetOptStat("ne"); + gStyle->SetOptFit(); +// hEtot->GetObjectStat(); + + TPaveStats *cstats = (TPaveStats*)gPad->GetPrimitive("stats"); + cout << cstats->GetX1NDC() << endl; + cout << cstats->GetX1() << endl; + + cout << cstats->GetX2NDC() << endl; + cout << cstats->GetX2() << endl; + + cout << cstats->GetY1() << endl; + cout << cstats->GetY2() << endl; + + cstats->SetY1(3000); + + c1->Update(); + + gStyle->ToggleEditor(); + + + return; + TCanvas *c2 = new TCanvas("name2", "title of canvas 2"); + c2->Divide(2,2); + + c2->cd(1); + tr->Draw("E", "E>0.95"); + TH1F *h = (TH1F*)gPad->GetPrimitive("htemp"); + + c2->cd(2); +// fr.cd(); + TH1F *hE = new TH1F("hE", "histogram title", 100, 0.95, 1.4); + hE->Draw(); + tr->Draw("E >> hE", "", ""); + + TF1 *g1 = new TF1("m1","gaus",1.,1.15); + TF1 *g2 = new TF1("m2","gaus",1.2,1.3); + // The total is the sum of the three, each has 3 parameters + TF1 *total = new TF1("mstotal","gaus(0)+gaus(3)",1.,1.3); + + + Double_t par[6]; + + // Fit each function and add it to the list of functions + hE->Fit(g1,"R"); + hE->Fit(g2,"R+"); + + + c2->cd(3); + TH1F *hEtot = new TH1F("hEtot", "histogram title", 100, 0.95, 1.4); + hEtot->Draw(); + tr->Draw("E >> hEtot", "", ""); + // Get the parameters from the fit + g1->GetParameters(&par[0]); + g2->GetParameters(&par[3]); + + // Use the parameters on the sum + total->SetParameters(par); + hEtot->Fit(total,"R+"); + + TCanvas *c3 = new TCanvas("name3", "title of canvas 3"); + c3->Divide(2,2); + + c3->cd(1); +// tr1->Draw("E", "E>0.95"); +// TH1F *h = (TH1F*)gPad->GetPrimitive("htemp"); + + c3->cd(2); + // fr.cd(); + TH1F *hE1 = new TH1F("hE1", "histogram title", 100, 0.95, 1.4); + hE1->Draw(); +// tr1->Draw("E >> hE1", "", ""); + + TF1 *g11 = new TF1("m11","gaus",1.,1.15); + TF1 *g21 = new TF1("m21","gaus",1.2,1.3); + // The total is the sum of the three, each has 3 parameters + TF1 *total1 = new TF1("mstotal1","gaus(0)+gaus(3)",1.,1.3); + + + Double_t par1[6]; + + // Fit each function and add it to the list of functions + hE1->Fit(g11,"R"); + hE1->Fit(g21,"R+"); + + + c3->cd(3); + TH1F *hEtot1 = new TH1F("hEtot1", "histogram title", 100, 0.95, 1.4); + hEtot1->Draw(); +// tr1->Draw("E >> hEtot1", "", ""); + // Get the parameters from the fit + g11->GetParameters(&par1[0]); + g21->GetParameters(&par1[3]); + + // Use the parameters on the sum + total1->SetParameters(par1); + hEtot1->Fit(total1,"R+"); + + + c3->cd(4); + TH1F *hEtot11 = new TH1F("hEtot11", "histogram title", 100, 0.95, 1.4); +// tr1->Draw("E >> hEtot11", "", ""); + + Double_t par12[6]; + total1->GetParameters(&par12[0]); + TF1 *gLow = new TF1("mLow","gaus",0.95,1.2); + TF1 *gHigh = new TF1("mHigh","gaus",1.0,1.4); + gHigh->SetLineColor(kGreen); + + gLow->SetParameter(0, par12[0]); + gLow->SetParameter(1, par12[1]); + gLow->SetParameter(2, par12[2]); + + gHigh->SetParameter(0, par12[3]); + gHigh->SetParameter(1, par12[4]); + gHigh->SetParameter(2, par12[5]); + + gLow->Draw("same"); + gHigh->Draw("same"); + + + + +} diff --git a/makefile b/makefile new file mode 100644 index 0000000..e69de29 diff --git a/opentrees.cxx b/opentrees.cxx new file mode 100644 index 0000000..f4fac37 --- /dev/null +++ b/opentrees.cxx @@ -0,0 +1,142 @@ +#include "TFile.h" +#include "TTree.h" +#include "TCanvas.h" +#include "TF1.h" +#include "TH1F.h" + +void opentrees() +{ + + + +// TFile *fr1 = new TFile("simevents-1.root"); + TFile *fr1 = new TFile("simeventsNew.root"); + TFile *fr2 = new TFile("simevents-2.root"); + TFile *fr = new TFile("simevents.root"); + + TTree *tr1 = (TTree*)fr1->Get("simevents"); + TTree *tr2 = (TTree*)fr2->Get("simevents"); + TTree *tr = (TTree*)fr->Get("simevents"); + + TCanvas *c1 = new TCanvas("name1", "title of canvas"); + c1->Divide(2,2); + + c1->cd(1); + tr->Draw("E", "E>0"); + + c1->cd(2); + tr1->Draw("E", "E>0"); + + c1->cd(3); + tr2->Draw("E", "E>0"); + + + c1->cd(4); + tr->Draw("E", "E>0.95"); + tr1->SetLineColor(kRed); + tr1->Draw("E", "E>0", "same"); + tr2->SetLineColor(kBlue); + tr2->Draw("E", "E>0", "same"); + + TCanvas *c2 = new TCanvas("name2", "title of canvas"); + c2->Divide(2,2); + + c2->cd(1); + tr->Draw("E", "E>0.95"); + TH1F *h = (TH1F*)gPad->GetPrimitive("htemp"); + + c2->cd(2); +// fr.cd(); + TH1F *hE = new TH1F("hE", "histogram title", 100, 0.95, 1.4); + hE->Draw(); + tr->Draw("E >> hE", "", ""); + + TF1 *g1 = new TF1("m1","gaus",1.,1.15); + TF1 *g2 = new TF1("m2","gaus",1.2,1.3); + // The total is the sum of the three, each has 3 parameters + TF1 *total = new TF1("mstotal","gaus(0)+gaus(3)",1.,1.3); + + + Double_t par[6]; + + // Fit each function and add it to the list of functions + hE->Fit(g1,"R"); + hE->Fit(g2,"R+"); + + + c2->cd(3); + TH1F *hEtot = new TH1F("hEtot", "histogram title", 100, 0.95, 1.4); + hEtot->Draw(); + tr->Draw("E >> hEtot", "", ""); + // Get the parameters from the fit + g1->GetParameters(&par[0]); + g2->GetParameters(&par[3]); + + // Use the parameters on the sum + total->SetParameters(par); + hEtot->Fit(total,"R+"); + + TCanvas *c3 = new TCanvas("name3", "title of canvas 3"); + c3->Divide(2,2); + + c3->cd(1); + tr1->Draw("E", "E>0.95"); +// TH1F *h = (TH1F*)gPad->GetPrimitive("htemp"); + + c3->cd(2); + // fr.cd(); + TH1F *hE1 = new TH1F("hE1", "histogram title", 100, 0.95, 1.4); + hE1->Draw(); + tr1->Draw("E >> hE1", "", ""); + + TF1 *g11 = new TF1("m11","gaus",1.,1.15); + TF1 *g21 = new TF1("m21","gaus",1.2,1.3); + // The total is the sum of the three, each has 3 parameters + TF1 *total1 = new TF1("mstotal1","gaus(0)+gaus(3)",1.,1.3); + + + Double_t par1[6]; + + // Fit each function and add it to the list of functions + hE1->Fit(g11,"R"); + hE1->Fit(g21,"R+"); + + + c3->cd(3); + TH1F *hEtot1 = new TH1F("hEtot1", "histogram title", 100, 0.95, 1.4); + hEtot1->Draw(); + tr1->Draw("E >> hEtot1", "", ""); + // Get the parameters from the fit + g11->GetParameters(&par1[0]); + g21->GetParameters(&par1[3]); + + // Use the parameters on the sum + total1->SetParameters(par1); + hEtot1->Fit(total1,"R+"); + + + c3->cd(4); + TH1F *hEtot11 = new TH1F("hEtot11", "histogram title", 100, 0.95, 1.4); + tr1->Draw("E >> hEtot11", "", ""); + + Double_t par12[6]; + total1->GetParameters(&par12[0]); + TF1 *gLow = new TF1("mLow","gaus",0.95,1.2); + TF1 *gHigh = new TF1("mHigh","gaus",1.0,1.4); + gHigh->SetLineColor(kGreen); + + gLow->SetParameter(0, par12[0]); + gLow->SetParameter(1, par12[1]); + gLow->SetParameter(2, par12[2]); + + gHigh->SetParameter(0, par12[3]); + gHigh->SetParameter(1, par12[4]); + gHigh->SetParameter(2, par12[5]); + + gLow->Draw("same"); + gHigh->Draw("same"); + + + + +} -- 2.18.1