File indexing completed on 2023-03-17 11:14:51
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
0041
0042
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
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
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
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
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
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
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
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
0225
0226 class ResidualFit
0227 {
0228 public:
0229 ResidualFit()
0230 {
0231 parNum_ = 6;
0232 }
0233 double operator() (double *x, double *p)
0234 {
0235
0236
0237
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
0240 if( x[0] < 90.56 ) return p[6] + p[7]*x[0] + p[8]*x[0]*x[0];
0241
0242 if( x[0] < 91.76 ) return p[0] + p[1]*x[0] + p[2]*x[0]*x[0];
0243
0244
0245
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
0275
0276
0277
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
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
0334
0335
0336 functionToFit = expRelativisticBWintPhotFit(ss.str());
0337
0338
0339
0340
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
0367
0368 canvas->cd(2+2*iCanvas);
0369 TH1F * residuals = (TH1F*)histo->Clone();
0370 residuals->Add(functionHisto, -1);
0371 residuals->SetName("Residuals");
0372
0373
0374
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
0393
0394 TF1 * residuals1 = iterateFitter(histo, 0, 0, canvas);
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406 }