test/Macros/BackgroundFit.C

0001 #include "TF1.h"
0002 #include <cmath>
0003 #include <iostream>
0004 #include "TFile.h"
0005 #include "TH1F.h"
0006
0007
0008 /**
0009  * This macro can be used to fit the background and signal peak of the Z. <br>
0010  * It contains several different functions that can be activated by uncommenting
0011  * the BackgroundFit function.
0012  */
0013
0014 // e^-(t^2)/(y^2+(x-t)^2)
0015 // x[0] = t
0016 // par[0] = y
0017 // par[1] = x
0018 Double_t integralPart(Double_t *x, Double_t *par)
0019 {
0020   std::cout << "x[0] = " << x[0] << std::endl;
0021   std::cout << "par[0] = " << par[0] << std::endl;
0022   std::cout << "par[1] = " << par[1] << std::endl;
0023   // std::cout << "value = " << exp(-(x[0]-par[1])*(x[0]-par[1])/(par[0]*par[0]))/(par[0]*par[0] + (par[1] - x[0])*(par[1] - x[0])) << std::endl;
0024   return exp(-x[0]*x[0])/(par[0]*par[0] + (par[1] - x[0])*(par[1] - x[0]));
0025 }
0026
0027 // CrystalBall fit
0028 // ---------------
0029 // see for example http://en.wikipedia.org/wiki/Crystal_Ball_function
0030 // par[0] = x0
0031 // par[1] = sigma
0032 // par[2] = A
0033 // par[3] = B
0034 // par[4] = alpha
0035 // par[5] = N
0036 // par[6] = n
0037 Double_t crystalBall(Double_t *x, Double_t *par)
0038 {
0039   double x_0 = par[0];
0040   double sigma = par[1];
0041   double n = par[2];
0042   double alpha = par[3];
0043   double N = par[4];
0044
0045   double A = pow(n/fabs(alpha), n)*exp(-(alpha*alpha)/2);
0046   double B = n/fabs(alpha) - fabs(alpha);
0047
0048   if((x[0]-x_0)/sigma > -alpha) {
0049     return N*exp(-(x[0]-x_0)*(x[0]-x_0)/(2*sigma*sigma));
0050   }
0051   else {
0052     return N*A*pow( (B - (x[0]-x_0)/sigma), -n );
0053   }
0054 }
0055 TF1 * crystalBallFit()
0056 {
0057   TF1 * functionToFit = new TF1("functionToFit", crystalBall, 42, 160, 5);
0058   functionToFit->SetParameter(0, 90);
0059   functionToFit->SetParameter(1, 10);
0060   functionToFit->SetParameter(2, 1.);
0061   functionToFit->SetParameter(3, 1.);
0062   functionToFit->SetParameter(4, 1.);
0063   return functionToFit;
0064 }
0065
0066 // Lorentzian function and fit
0067 // ---------------------------
0068 TF1 * lorentzian = new TF1( "lorentzian", "[2]/((x-[0])*(x-[0])+(([1]/2)*([1]/2)))", 0, 1000);
0069 TF1 * lorentzianFit()
0070 {
0071   TF1 * functionToFit = new TF1("functionToFit", lorentzian, 42, 160, 3);
0072   functionToFit->SetParameter(0, 90);
0073   functionToFit->SetParameter(1, 10);
0074   return functionToFit;
0075 }
0076
0077 // Power law fit
0078 // -------------
0079 TF1 * powerLaw = new TF1( "powerLaw", "[0] + [1]*pow(x, [2])", 0, 1000 );
0080 TF1 * powerLawFit()
0081 {
0082   TF1 * functionToFit = new TF1("functionToFit", powerLaw, 42, 160, 3);
0083   return functionToFit;
0084 }
0085
0086 // Lorentzian + and power law
0087 TF1 * lorentzianAndPowerLaw()
0088 {
0089   TF1 * functionToFit = new TF1("functionToFit", "[2]/((x-[0])*(x-[0])+(([1]/2)*([1]/2))) + ([3] + [4]*pow(x, [5]))", 42, 160);
0090   functionToFit->SetParameter(0, 90);
0091   functionToFit->SetParameter(1, 10);
0092   functionToFit->SetParameter(2, 0.33);
0093   functionToFit->SetParameter(3, -0.0220578);
0094   functionToFit->SetParameter(4, 0.0357716);
0095   functionToFit->SetParameter(5, -0.0962408);
0096   return functionToFit;
0097 }
0098
0099 // Exponential fit
0100 // ---------------
0101 TF1 * exponential = new TF1 ("exponential", "[0]*exp([1]*x)", 0, 1000 );
0102 TF1 * exponentialFit()
0103 {
0104   TF1 * functionToFit = new TF1("functionToFit", exponential, 42, 160, 2);
0105   return functionToFit;
0106 }
0107
0108 // Lorentzian + exponential fit
0109 // ----------------------------
0110 TF1 * lorenzianAndExponentialFit()
0111 {
0112   TF1 * functionToFit = new TF1("functionToFit", "[2]/((x-[0])*(x-[0])+(([1]/2)*([1]/2))) + [3]*exp([4]*x)", 42, 160);
0113   functionToFit->SetParameter(0, 90);
0114   functionToFit->SetParameter(1, 10);
0115   functionToFit->SetParameter(2, 0.33);
0116   functionToFit->SetParameter(3, 0.0115272);
0117   functionToFit->SetParameter(4, -0.0294229);
0118   return functionToFit;
0119 }
0120
0121
0122
0123 void BackgroundFit() {
0124   TFile* inputFile = new TFile("0_MuScleFit.root", "READ");
0125   TH1F* histo = (TH1F*)inputFile->FindObjectAny("hRecBestRes_Mass");
0126   histo->Rebin(30);
0127   histo->Scale(1/histo->GetEntries());
0128
0129   // TF1 * functionToFit = lorentzianFit();
0130   // TF1 * functionToFit = crystalBallFit();
0131   // TF1 * functionToFit = powerLawFit();
0132   // TF1 * functionToFit = lorentzianAndPowerLaw();
0133   // TF1 * functionToFit = exponentialFit();
0134   TF1 * functionToFit = lorenzianAndExponentialFit();
0135
0136   histo->Fit(functionToFit, "M", "", 42, 160);
0137
0138 }