test/Macros/ProbsFitter.C

0001 #include "TFile.h"
0002 #include "TCanvas.h"
0003 #include "TH1F.h"
0004 #include "TF1.h"
0005 #include <cmath>
0006 #include <sstream>
0007 #include "TROOT.h"
0008 #include "TStyle.h"
0009 #include "TMath.h"
0010
0011 Double_t crystalBall(Double_t *x, Double_t *par)
0012 {
0013   double x_0 = par[0];
0014   double sigma = par[1];
0015   double n = par[2];
0016   double alpha = par[3];
0017   double N = par[4];
0018
0019   double A = pow(n/fabs(alpha), n)*exp(-(alpha*alpha)/2);
0020   double B = n/fabs(alpha) - fabs(alpha);
0021
0022   if((x[0]-x_0)/sigma > -alpha) {
0023     return N*exp(-(x[0]-x_0)*(x[0]-x_0)/(2*sigma*sigma));
0024   }
0025   else {
0026     return N*A*pow( (B - (x[0]-x_0)/sigma), -n );
0027   }
0028 }
0029 TF1 * crystalBallFit()
0030 {
0031   TF1 * functionToFit = new TF1("functionToFit", crystalBall, 60, 120, 5);
0032   functionToFit->SetParameter(0, 90);
0033   functionToFit->SetParameter(1, 10);
0034   functionToFit->SetParameter(2, 1.);
0035   functionToFit->SetParameter(3, 1.);
0036   functionToFit->SetParameter(4, 1.);
0037   return functionToFit;
0038 }
0039
0040 // Double_t reversedCrystalBall(Double_t *x, Double_t *par)
0041 // {
0042 //   double x_0 = par[0];
0043 //   double sigma = par[1];
0044 //   double n = par[2];
0045 //   double alpha = par[3];
0046 //   double N = par[4];
0047
0048 //   double A = pow(n/fabs(alpha), n)*exp(-(alpha*alpha)/2);
0049 //   double B = n/fabs(alpha) - fabs(alpha);
0050
0051 //   if((x[0]-x_0)/sigma < alpha) {
0052 //     return N*exp(-(x[0]-x_0)*(x[0]-x_0)/(2*sigma*sigma));
0053 //   }
0054 //   else {
0055 //     return N*A*pow( (B + (x[0]-x_0)/sigma), -n );
0056 //   }
0057 // }
0058 // TF1 * reversedCrystalBallFit()
0059 // {
0060 //   TF1 * functionToFit = new TF1("functionToFit", reversedCrystalBall, 60, 120, 5);
0061 //   functionToFit->SetParameter(0, 90);
0062 //   functionToFit->SetParameter(1, 10);
0063 //   functionToFit->SetParameter(2, 1.);
0064 //   functionToFit->SetParameter(3, 1.);
0065 //   functionToFit->SetParameter(4, 1.);
0066 //   return functionToFit;
0067 // }
0068
0069 // Lorentzian function and fit (non relativistic Breit-Wigner)
0070 // -----------------------------------------------------------
0071 TF1 * lorentzian = new TF1( "lorentzian", "[2]/((x-[0])*(x-[0])+(([1]/2)*([1]/2)))", 0, 1000);
0072 TF1 * lorentzianFit()
0073 {
0074   TF1 * functionToFit = new TF1("functionToFit", lorentzian, 60, 120, 3);
0075   functionToFit->SetParameter(0, 90);
0076   functionToFit->SetParameter(1, 10);
0077   return functionToFit;
0078 }
0079
0080 // Relativistic Breit-Wigner
0081 // -------------------------
0082 TF1 * relativisticBW = new TF1 ("relativisticBW", "[2]*pow([0]*[1],2)/(pow(x*x-[1]*[1],2)+pow([0]*[1],2))", 0., 1000.);
0083 TF1 * relativisticBWFit(const std::string & index)
0084 {
0085   TF1 * functionToFit = new TF1(("functionToFit"+index).c_str(), relativisticBW, 60, 120, 3);
0086   functionToFit->SetParameter(1, 90);
0087   functionToFit->SetParameter(0, 2);
0088   return functionToFit;
0089 }
0090
0091 // Relativistic Breit-Wigner with Z/gamma interference term
0092 // --------------------------------------------------------
0093 class RelativisticBWwithZGammaInterference
0094 {
0095  public:
0096   RelativisticBWwithZGammaInterference()
0097   {
0098     parNum_ = 4;
0099     twoOverPi_ = 2./TMath::Pi();
0100   }
0101   double operator() (double *x, double *p)
0102   {
0103     double squaredMassDiff = pow((x[0]*x[0] - p[1]*p[1]), 2);
0104     double denominator = squaredMassDiff + pow(x[0], 4)*pow(p[0]/p[1], 2);
0105     return p[2]*( p[3]*twoOverPi_*pow(p[1]*p[0], 2)/denominator + (1-p[3])*p[1]*squaredMassDiff/denominator );
0106   }
0107   int parNum() const { return parNum_; }
0108  protected:
0109   int parNum_;
0110   double twoOverPi_;
0111 };
0112
0113 TF1 * relativisticBWintFit(const std::string & index)
0114 {
0115   RelativisticBWwithZGammaInterference * fobj = new RelativisticBWwithZGammaInterference;
0116   TF1 * functionToFit = new TF1(("functionToFit"+index).c_str(), fobj, 60, 120, fobj->parNum(), "RelativisticBWwithZGammaInterference");
0117   functionToFit->SetParameter(0, 2.);
0118   functionToFit->SetParameter(1, 90.);
0119   functionToFit->SetParameter(2, 1.);
0120   functionToFit->SetParameter(3, 1.);
0121
0122   functionToFit->SetParLimits(3, 0., 1.);
0123
0124   return functionToFit;
0125 }
0126
0127
0128
0129
0130 // Relativistic Breit-Wigner with Z/gamma interference term and photon propagator
0131 // ------------------------------------------------------------------------------
0132 class RelativisticBWwithZGammaInterferenceAndPhotonPropagator
0133 {
0134  public:
0135   RelativisticBWwithZGammaInterferenceAndPhotonPropagator()
0136   {
0137     parNum_ = 5;
0138     twoOverPi_ = 2./TMath::Pi();
0139   }
0140   double operator() (double *x, double *p)
0141   {
0142
0143     // if( p[3]+p[4] > 1 ) return -10000.;
0144
0145     double squaredMassDiff = pow((x[0]*x[0] - p[1]*p[1]), 2);
0146     double denominator = squaredMassDiff + pow(x[0], 4)*pow(p[0]/p[1], 2);
0147     return p[2]*( p[3]*twoOverPi_*pow(p[1]*p[0], 2)/denominator + (1-p[3]-p[4])*p[1]*squaredMassDiff/denominator + p[4]/(x[0]*x[0]));
0148   }
0149   int parNum() const { return parNum_; }
0150  protected:
0151   int parNum_;
0152   double twoOverPi_;
0153 };
0154
0155 TF1 * relativisticBWintPhotFit(const std::string & index)
0156 {
0157   RelativisticBWwithZGammaInterferenceAndPhotonPropagator * fobj = new RelativisticBWwithZGammaInterferenceAndPhotonPropagator;
0158   TF1 * functionToFit = new TF1(("functionToFit"+index).c_str(), fobj, 60, 120, fobj->parNum(), "RelativisticBWwithZGammaInterferenceAndPhotonPropagator");
0159   functionToFit->SetParameter(0, 2.);
0160   functionToFit->SetParameter(1, 90.);
0161   functionToFit->SetParameter(2, 1.);
0162   functionToFit->SetParameter(3, 1.);
0163   functionToFit->SetParameter(4, 0.);
0164
0165   functionToFit->SetParLimits(3, 0., 1.);
0166   functionToFit->SetParLimits(4, 0., 1.);
0167
0168   return functionToFit;
0169 }
0170
0171
0172
0173
0174
0175
0176 // Product between an exponential term and the relativistic Breit-Wigner with Z/gamma interference term and photon propagator
0177 // --------------------------------------------------------------------------------------------------------------------------
0178 class ExpRelativisticBWwithZGammaInterferenceAndPhotonPropagator
0179 {
0180  public:
0181   ExpRelativisticBWwithZGammaInterferenceAndPhotonPropagator()
0182   {
0183     parNum_ = 6;
0184     twoOverPi_ = 2./TMath::Pi();
0185   }
0186   double operator() (double *x, double *p)
0187   {
0188
0189     // if( p[3]+p[4] > 1 ) return -10000.;
0190
0191     double squaredMassDiff = pow((x[0]*x[0] - p[1]*p[1]), 2);
0192     double denominator = squaredMassDiff + pow(x[0], 4)*pow(p[0]/p[1], 2);
0193     return p[2]*exp(-p[5]*x[0])*( p[3]*twoOverPi_*pow(p[1]*p[0], 2)/denominator + (1-p[3]-p[4])*p[1]*squaredMassDiff/denominator + p[4]/(x[0]*x[0]));
0194   }
0195   int parNum() const { return parNum_; }
0196  protected:
0197   int parNum_;
0198   double twoOverPi_;
0199 };
0200
0201 TF1 * expRelativisticBWintPhotFit(const std::string & index)
0202 {
0203   ExpRelativisticBWwithZGammaInterferenceAndPhotonPropagator * fobj = new ExpRelativisticBWwithZGammaInterferenceAndPhotonPropagator;
0204   TF1 * functionToFit = new TF1(("functionToFit"+index).c_str(), fobj, 60, 120, fobj->parNum(), "ExpRelativisticBWwithZGammaInterferenceAndPhotonPropagator");
0205   functionToFit->SetParameter(0, 2.);
0206   functionToFit->SetParameter(1, 90.);
0207   functionToFit->SetParameter(2, 1.);
0208   functionToFit->SetParameter(3, 1.);
0209   functionToFit->SetParameter(4, 0.);
0210   functionToFit->SetParameter(5, 0.);
0211
0212   functionToFit->SetParLimits(3, 0., 1.);
0213   functionToFit->SetParLimits(4, 0., 1.);
0214
0215   return functionToFit;
0216 }
0217
0218
0219
0220
0221
0222
0223
0224 // Residual fit
0225 // ------------
0226 class ResidualFit
0227 {
0228  public:
0229   ResidualFit()
0230   {
0231     parNum_ = 6;
0232   }
0233   double operator() (double *x, double *p)
0234   {
0235 //     return p[0] + p[1]*x[0] + p[2]*x[0]*x[0] + p[3]/x[0] + p[4]/(x[0]*x[0]) + p[5]/(pow(x[0],3)) + p[6]*pow(x[0],p[7]);
0236
0237     // if( x[0] < 90.25 ) return (p[0] + p[1]/x[0] + p[2]*x[0] + p[3]*x[0]*x[0] + p[4]*pow(x[0],3))/(p[5] + p[6]/x[0] + p[7]*x[0] + p[8]*x[0]*x[0] + p[9]*pow(x[0],3));
0238     if( x[0] < 89.5 ) return p[0] + p[1]*x[0] + p[2]*x[0]*x[0] + p[3]/x[0] + p[4]/(x[0]*x[0]) + p[5]/(pow(x[0],3));
0239     // if( x[0] < 92 ) return p[6] + p[7]/x[0];
0240     if( x[0] < 90.56 ) return p[6] + p[7]*x[0] + p[8]*x[0]*x[0];
0241     // if( x[0] < 94.5 ) return p[8] + p[9]/x[0];
0242     if( x[0] < 91.76 ) return p[0] + p[1]*x[0] + p[2]*x[0]*x[0];
0243 //     if( x[0] < 94 ) return p[8] + p[9]*x[0] + p[14]*x[0]*x[0];
0244 //     return p[10] + p[11]*x[0] + p[12]*x[0]*x[0] + p[15]*pow(x[0],3) + p[16]*exp(p[17]*x[0]);
0245 //     // else return p[10] + p[11]/x[0] + p[12]*x[0] + p[13]*x[0] + p[14]*x[0] + p[15]*exp(p[16]*x[0]);
0246
0247     return p[3] + p[4]/x[0] + p[5]/pow(x[0], 2);
0248
0249   }
0250   int parNum() const { return parNum_; }
0251  protected:
0252   int parNum_;
0253 };
0254
0255 TF1 * residualFit(const std::string & index)
0256 {
0257   ResidualFit * fobj = new ResidualFit;
0258   TF1 * functionToFit = new TF1(("f"+index).c_str(), fobj, 60, 120, fobj->parNum(), "ResidualFit");
0259
0260   double pars[] = { -2.30972e-02, 7.94253e-01, 7.40701e-05, 2.79889e-06,
0261                     -2.04844e-08, -5.78370e-43, 6.84145e-02, -6.21316e+00,
0262                     -6.90792e-03, 6.43131e-01, -2.03049e-03, 3.70415e-05,
0263                     -1.69014e-07 };
0264   functionToFit->SetParameters(pars);
0265
0266   return functionToFit;
0267 }
0268
0269 TF1 * residualFit(double pars[], const std::string & index)
0270 {
0271   ResidualFit * fobj = new ResidualFit;
0272   TF1 * functionToFit = new TF1(("fp"+index).c_str(), fobj, 60, 120, fobj->parNum(), "ResidualFit");
0273
0274 //   double pars[] = { -2.30972e-02, 7.94253e-01, 7.40701e-05, 2.79889e-06,
0275 //                     -2.04844e-08, -5.78370e-43, 6.84145e-02, -6.21316e+00,
0276 //                     -6.90792e-03, 6.43131e-01, -2.03049e-03, 3.70415e-05,
0277 //                     -1.69014e-07 };
0278   functionToFit->SetParameters(pars);
0279
0280   return functionToFit;
0281 }
0282
0283
0284
0285
0286
0287 class CombinedFit
0288 {
0289  public:
0290   CombinedFit(double pars[], const std::string & index)
0291   {
0292     residuals_ = residualFit(pars, index);
0293
0294     parNum_ = 3;
0295   }
0296   double operator() (double *x, double *p)
0297   {
0298     return p[2]/((x[0]-p[0])*(x[0]-p[0])+((p[1]/2)*(p[1]/2))) + residuals_->Eval(x[0]);
0299     // return p[2]/((x[0]-p[0])*(x[0]-p[0])+((p[1]/2)*(p[1]/2)));
0300   }
0301   int parNum() const { return parNum_; }
0302  protected:
0303   int parNum_;
0304   TF1 * residuals_;
0305 };
0306
0307 TF1 * combinedFit(double pars[], const std::string & index)
0308 {
0309   CombinedFit * combined = new CombinedFit(pars, index);
0310   TF1 * functionToFit = new TF1(("functionToFit"+index).c_str(), combined, 60, 120, combined->parNum(), "functionToFit");
0311
0312
0313   functionToFit->SetParameter(0, 90);
0314   functionToFit->SetParameter(1, 10);
0315   functionToFit->SetParameter(2, 1.);
0316
0317   return functionToFit;
0318 }
0319
0320
0321
0322 TF1 * iterateFitter(TH1F * histo, int i, TF1 * previousResiduals = 0, TCanvas * inputCanvas = 0)
0323 {
0324   TCanvas * canvas = inputCanvas;
0325   std::stringstream ss;
0326   int iCanvas = i;
0327   ss << i;
0328
0329   int nBins = histo->GetNbinsX();
0330
0331   TF1 * functionToFit;
0332   if( previousResiduals == 0 ) {
0333     // functionToFit = relativisticBWFit(ss.str());
0334     // functionToFit = relativisticBWintFit(ss.str());
0335     // functionToFit = relativisticBWintPhotFit(ss.str());
0336     functionToFit = expRelativisticBWintPhotFit(ss.str());
0337     // functionToFit = crystalBallFit();
0338     // functionToFit = reversedCrystalBallFit();
0339     // functionToFit = lorentzianFit();
0340     // functionToFit = relativisticBWFit();
0341   }
0342   else {
0343     functionToFit = combinedFit(previousResiduals->GetParameters(), ss.str());
0344   }
0345   histo->Fit(functionToFit, "MN", "", 60, 120);
0346
0347
0348   double xMin = histo->GetXaxis()->GetXmin();
0349   double xMax = histo->GetXaxis()->GetXmax();
0350   double step = (xMax-xMin)/(double(nBins));
0351   TH1F * functionHisto = new TH1F(("functionHisto"+ss.str()).c_str(), "functionHisto", nBins, xMin, xMax);
0352   for( int i=0; i<nBins; ++i ) {
0353     functionHisto->SetBinContent( i+1, functionToFit->Eval(xMin + (i+0.5)*step) );
0354   }
0355
0356   if( canvas == 0 ) {
0357     canvas = new TCanvas(("canvasResiduals"+ss.str()).c_str(), ("canvasResiduals"+ss.str()).c_str(), 1000, 800);
0358     canvas->Divide(2,1);
0359     canvas->Draw();
0360     iCanvas = 0;
0361   }
0362   canvas->cd(1+2*iCanvas);
0363   histo->Draw();
0364   functionToFit->SetLineColor(kRed);
0365   functionToFit->Draw("same");
0366   // functionHisto->Draw("same");
0367   // functionHisto->SetLineColor(kGreen);
0368   canvas->cd(2+2*iCanvas);
0369   TH1F * residuals = (TH1F*)histo->Clone();
0370   residuals->Add(functionHisto, -1);
0371   residuals->SetName("Residuals");
0372
0373   // TF1 * residualFit = new TF1("residualFit", "-[0] + [1]*x+sqrt( ([1]-1)*([1]-1)*x*x + [0]*[0] )", 0., 1000. );
0374   // TF1 * residualFit = new TF1("residualFit", "[0]*(x - [1])/([2]*x*x + [3]*x + [4])", 0., 1000. );
0375   TF1 * residualFitFunction = residualFit(ss.str());
0376   residuals->Fit(residualFitFunction, "ME", "", 90.56, 120);
0377
0378   residuals->Draw();
0379
0380   return residualFitFunction;
0381 }
0382
0383 void ProbsFitter()
0384 {
0385   TFile * inputFile = new TFile("Sherpa_nocuts.root", "READ");
0386   TH1F * histo = (TH1F*)inputFile->FindObjectAny("HistAllZMass");
0387   histo->Scale(1/histo->GetEntries());
0388   histo->Rebin(6);
0389   TCanvas * canvas = new TCanvas("canvas", "canvas", 1000, 800);
0390   canvas->Divide(2,3);
0391
0392   // gStyle->SetOptFit(1);
0393
0394   TF1 * residuals1 = iterateFitter(histo, 0, 0, canvas);
0395 //   TF1 * residuals2 = iterateFitter(histo, 1, residuals1, canvas);
0396 //   iterateFitter(histo, 2, residuals2, canvas);
0397
0398 //   canvas->cd(3);
0399 //   TH1F * hclone = (TH1F*)histo->Clone();
0400 //   TF1 * combinedFunctionToFit = combinedFit(residualFitFunction->GetParameters());
0401 //   hclone->Fit(combinedFunctionToFit, "MN", "", 60, 120);
0402 //   hclone->Draw();
0403 //   combinedFunctionToFit->Draw("same");
0404 //   combinedFunctionToFit->SetLineColor(kRed);
0405
0406 }