File indexing completed on 2023-03-17 10:39:57
0001 #ifndef FitWithRooFit_cc
0002 #define FitWithRooFit_cc
0003
0004 #ifndef __CINT__
0005 #include "RooGlobalFunc.h"
0006 #endif
0007 #include "TCanvas.h"
0008 #include "TTree.h"
0009 #include "TH1D.h"
0010 #include "TRandom.h"
0011 #include "RooRealVar.h"
0012 #include "RooDataSet.h"
0013 #include "RooGaussian.h"
0014 #include "RooVoigtian.h"
0015 #include "RooExponential.h"
0016 #include "RooPlot.h"
0017 #include "RooDataHist.h"
0018 #include "RooAddPdf.h"
0019 #include "RooChebychev.h"
0020 #include "RooGenericPdf.h"
0021 #include "RooGaussModel.h"
0022 #include "RooAddModel.h"
0023 #include "RooPolynomial.h"
0024 #include "RooCBShape.h"
0025 #include "RooMinimizer.h"
0026 #include "RooBreitWigner.h"
0027 #include "RooFFTConvPdf.h"
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047 namespace {
0048 typedef std::pair<RooRealVar, RooDataHist*> rooPair;
0049 }
0050
0051 class FitWithRooFit {
0052 public:
0053 FitWithRooFit()
0054 : useChi2_(false),
0055 mean_(nullptr),
0056 mean2_(nullptr),
0057 mean3_(nullptr),
0058 sigma_(nullptr),
0059 sigma2_(nullptr),
0060 sigma3_(nullptr),
0061 gamma_(nullptr),
0062 gaussFrac_(nullptr),
0063 gaussFrac2_(nullptr),
0064 expCoeffa0_(nullptr),
0065 expCoeffa1_(nullptr),
0066 expCoeffa2_(nullptr),
0067 fsig_(nullptr),
0068 a0_(nullptr),
0069 a1_(nullptr),
0070 a2_(nullptr),
0071 a3_(nullptr),
0072 a4_(nullptr),
0073 a5_(nullptr),
0074 a6_(nullptr),
0075 alpha_(nullptr),
0076 n_(nullptr),
0077 fGCB_(nullptr) {}
0078
0079
0080 rooPair importTH1(TH1* histo, const double& inputXmin, const double& inputXmax) {
0081 double xMin = inputXmin;
0082 double xMax = inputXmax;
0083 if ((xMin == xMax) && xMin == 0) {
0084 xMin = histo->GetXaxis()->GetXmin();
0085 xMax = histo->GetXaxis()->GetXmax();
0086 }
0087
0088 RooRealVar x("x", "x", xMin, xMax);
0089
0090 return (std::make_pair(x, new RooDataHist("dh", "dh", x, RooFit::Import(*histo))));
0091 }
0092
0093
0094 void fit(TH1* histo,
0095 const TString signalType,
0096 const TString backgroundType,
0097 const double& xMin = 0.,
0098 const double& xMax = 0.,
0099 bool sumW2Error = false) {
0100 reinitializeParameters();
0101
0102 rooPair imported = importTH1(histo, xMin, xMax);
0103 RooRealVar x(imported.first);
0104 RooDataHist* dh = imported.second;
0105
0106
0107 RooPlot* frame = x.frame(RooFit::Title("Imported TH1 with Poisson error bars"));
0108 frame->SetName(TString(histo->GetName()) + "_frame");
0109 dh->plotOn(frame);
0110
0111
0112 RooAbsPdf* model = buildModel(&x, signalType, backgroundType);
0113
0114 std::unique_ptr<RooAbsReal> chi2{model->createChi2(*dh, RooFit::DataError(RooAbsData::SumW2))};
0115
0116
0117
0118
0119 if (!useChi2_) {
0120 if (sumW2Error)
0121 model->fitTo(*dh, RooFit::Save(), RooFit::SumW2Error(kTRUE));
0122 else
0123 model->fitTo(*dh);
0124 }
0125
0126 else {
0127 std::cout << "FITTING WITH CHI^2" << std::endl;
0128 RooMinimizer m(*chi2);
0129 m.migrad();
0130 m.hesse();
0131
0132 }
0133 model->plotOn(frame);
0134 model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDotted), RooFit::LineColor(kRed));
0135 model->paramOn(frame, RooFit::Label("fit result"), RooFit::Format("NEU", RooFit::AutoPrecision(2)));
0136
0137
0138
0139
0140
0141
0142
0143
0144
0145
0146
0147
0148 RooPlot* frame2 = x.frame(RooFit::Title("Imported TH1 with internal errors"));
0149 dh->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2));
0150 model->plotOn(frame2);
0151 model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineColor(kRed));
0152 model->paramOn(frame2, RooFit::Label("fit result"), RooFit::Format("NEU", RooFit::AutoPrecision(2)));
0153
0154
0155
0156
0157
0158
0159
0160
0161
0162
0163
0164
0165
0166
0167
0168
0169
0170 frame2->Draw();
0171 }
0172
0173
0174 void initMean(const double& value,
0175 const double& min,
0176 const double& max,
0177 const TString& name = "mean",
0178 const TString& title = "mean") {
0179 if (mean_ != nullptr)
0180 delete mean_;
0181 mean_ = new RooRealVar(name, title, value, min, max);
0182 initVal_mean = value;
0183 }
0184 void initMean2(const double& value,
0185 const double& min,
0186 const double& max,
0187 const TString& name = "mean2",
0188 const TString& title = "mean2") {
0189 if (mean2_ != nullptr)
0190 delete mean2_;
0191 mean2_ = new RooRealVar(name, title, value, min, max);
0192 initVal_mean2 = value;
0193 }
0194 void initMean3(const double& value,
0195 const double& min,
0196 const double& max,
0197 const TString& name = "mean3",
0198 const TString& title = "mean3") {
0199 if (mean3_ != nullptr)
0200 delete mean3_;
0201 mean3_ = new RooRealVar(name, title, value, min, max);
0202 initVal_mean3 = value;
0203 }
0204 void initSigma(const double& value,
0205 const double& min,
0206 const double& max,
0207 const TString& name = "sigma",
0208 const TString& title = "sigma") {
0209 if (sigma_ != nullptr)
0210 delete sigma_;
0211 sigma_ = new RooRealVar(name, title, value, min, max);
0212 initVal_sigma = value;
0213 }
0214 void initSigma2(const double& value,
0215 const double& min,
0216 const double& max,
0217 const TString& name = "sigma2",
0218 const TString& title = "sigma2") {
0219 if (sigma2_ != nullptr)
0220 delete sigma2_;
0221 sigma2_ = new RooRealVar(name, title, value, min, max);
0222 initVal_sigma2 = value;
0223 }
0224 void initSigma3(const double& value,
0225 const double& min,
0226 const double& max,
0227 const TString& name = "sigma3",
0228 const TString& title = "sigma3") {
0229 if (sigma3_ != nullptr)
0230 delete sigma3_;
0231 sigma3_ = new RooRealVar(name, title, value, min, max);
0232 initVal_sigma3 = value;
0233 }
0234 void initGamma(const double& value,
0235 const double& min,
0236 const double& max,
0237 const TString& name = "gamma",
0238 const TString& title = "gamma") {
0239 if (gamma_ != nullptr)
0240 delete gamma_;
0241 gamma_ = new RooRealVar(name, title, value, min, max);
0242 initVal_gamma = value;
0243 }
0244 void initGaussFrac(const double& value,
0245 const double& min,
0246 const double& max,
0247 const TString& name = "GaussFrac",
0248 const TString& title = "GaussFrac") {
0249 if (gaussFrac_ != nullptr)
0250 delete gaussFrac_;
0251 gaussFrac_ = new RooRealVar(name, title, value, min, max);
0252 initVal_gaussFrac = value;
0253 }
0254 void initGaussFrac2(const double& value,
0255 const double& min,
0256 const double& max,
0257 const TString& name = "GaussFrac2",
0258 const TString& title = "GaussFrac2") {
0259 if (gaussFrac2_ != nullptr)
0260 delete gaussFrac2_;
0261 gaussFrac2_ = new RooRealVar(name, title, value, min, max);
0262 initVal_gaussFrac2 = value;
0263 }
0264 void initExpCoeffA0(const double& value,
0265 const double& min,
0266 const double& max,
0267 const TString& name = "expCoeffa0",
0268 const TString& title = "expCoeffa0") {
0269 if (expCoeffa0_ != nullptr)
0270 delete expCoeffa0_;
0271 expCoeffa0_ = new RooRealVar(name, title, value, min, max);
0272 initVal_expCoeffa0 = value;
0273 }
0274 void initExpCoeffA1(const double& value,
0275 const double& min,
0276 const double& max,
0277 const TString& name = "expCoeffa1",
0278 const TString& title = "expCoeffa1") {
0279 if (expCoeffa1_ != nullptr)
0280 delete expCoeffa1_;
0281 expCoeffa1_ = new RooRealVar(name, title, value, min, max);
0282 initVal_expCoeffa1 = value;
0283 }
0284 void initExpCoeffA2(const double& value,
0285 const double& min,
0286 const double& max,
0287 const TString& name = "expCoeffa2",
0288 const TString& title = "expCoeffa2") {
0289 if (expCoeffa2_ != nullptr)
0290 delete expCoeffa2_;
0291 expCoeffa2_ = new RooRealVar(name, title, value, min, max);
0292 initVal_expCoeffa2 = value;
0293 }
0294 void initFsig(const double& value,
0295 const double& min,
0296 const double& max,
0297 const TString& name = "fsig",
0298 const TString& title = "signal fraction") {
0299 if (fsig_ != nullptr)
0300 delete fsig_;
0301 fsig_ = new RooRealVar(name, title, value, min, max);
0302 initVal_fsig = value;
0303 }
0304 void initA0(const double& value,
0305 const double& min,
0306 const double& max,
0307 const TString& name = "a0",
0308 const TString& title = "a0") {
0309 if (a0_ != nullptr)
0310 delete a0_;
0311 a0_ = new RooRealVar(name, title, value, min, max);
0312 initVal_a0 = value;
0313 }
0314 void initA1(const double& value,
0315 const double& min,
0316 const double& max,
0317 const TString& name = "a1",
0318 const TString& title = "a1") {
0319 if (a1_ != nullptr)
0320 delete a1_;
0321 a1_ = new RooRealVar(name, title, value, min, max);
0322 initVal_a1 = value;
0323 }
0324 void initA2(const double& value,
0325 const double& min,
0326 const double& max,
0327 const TString& name = "a2",
0328 const TString& title = "a2") {
0329 if (a2_ != nullptr)
0330 delete a2_;
0331 a2_ = new RooRealVar(name, title, value, min, max);
0332 initVal_a2 = value;
0333 }
0334 void initA3(const double& value,
0335 const double& min,
0336 const double& max,
0337 const TString& name = "a3",
0338 const TString& title = "a3") {
0339 if (a3_ != nullptr)
0340 delete a3_;
0341 a3_ = new RooRealVar(name, title, value, min, max);
0342 initVal_a3 = value;
0343 }
0344 void initA4(const double& value,
0345 const double& min,
0346 const double& max,
0347 const TString& name = "a4",
0348 const TString& title = "a4") {
0349 if (a4_ != nullptr)
0350 delete a4_;
0351 a4_ = new RooRealVar(name, title, value, min, max);
0352 initVal_a4 = value;
0353 }
0354 void initA5(const double& value,
0355 const double& min,
0356 const double& max,
0357 const TString& name = "a5",
0358 const TString& title = "a5") {
0359 if (a5_ != nullptr)
0360 delete a5_;
0361 a5_ = new RooRealVar(name, title, value, min, max);
0362 initVal_a5 = value;
0363 }
0364 void initA6(const double& value,
0365 const double& min,
0366 const double& max,
0367 const TString& name = "a6",
0368 const TString& title = "a6") {
0369 if (a6_ != nullptr)
0370 delete a6_;
0371 a6_ = new RooRealVar(name, title, value, min, max);
0372 initVal_a6 = value;
0373 }
0374 void initAlpha(const double& value,
0375 const double& min,
0376 const double& max,
0377 const TString& name = "alpha",
0378 const TString& title = "alpha") {
0379 if (alpha_ != nullptr)
0380 delete alpha_;
0381 alpha_ = new RooRealVar(name, title, value, min, max);
0382 initVal_alpha = value;
0383 }
0384 void initN(const double& value,
0385 const double& min,
0386 const double& max,
0387 const TString& name = "n",
0388 const TString& title = "n") {
0389 if (n_ != nullptr)
0390 delete n_;
0391 n_ = new RooRealVar(name, title, value, min, max);
0392 initVal_n = value;
0393 }
0394 void initFGCB(const double& value,
0395 const double& min,
0396 const double& max,
0397 const TString& name = "fGCB",
0398 const TString& title = "fGCB") {
0399 if (fGCB_ != nullptr)
0400 delete fGCB_;
0401 fGCB_ = new RooRealVar(name, title, value, min, max);
0402 initVal_fGCB = value;
0403 }
0404
0405 void reinitializeParameters() {
0406 if (mean_ != nullptr)
0407 mean_->setVal(initVal_mean);
0408 if (mean2_ != nullptr)
0409 mean2_->setVal(initVal_mean2);
0410 if (mean3_ != nullptr)
0411 mean3_->setVal(initVal_mean3);
0412 if (sigma_ != nullptr)
0413 sigma_->setVal(initVal_sigma);
0414 if (sigma2_ != nullptr)
0415 sigma2_->setVal(initVal_sigma2);
0416 if (sigma3_ != nullptr)
0417 sigma3_->setVal(initVal_sigma3);
0418 if (gamma_ != nullptr)
0419 gamma_->setVal(initVal_gamma);
0420 if (gaussFrac_ != nullptr)
0421 gaussFrac_->setVal(initVal_gaussFrac);
0422 if (gaussFrac2_ != nullptr)
0423 gaussFrac2_->setVal(initVal_gaussFrac2);
0424 if (expCoeffa0_ != nullptr)
0425 expCoeffa0_->setVal(initVal_expCoeffa0);
0426 if (expCoeffa1_ != nullptr)
0427 expCoeffa1_->setVal(initVal_expCoeffa1);
0428 if (expCoeffa2_ != nullptr)
0429 expCoeffa2_->setVal(initVal_expCoeffa2);
0430 if (fsig_ != nullptr)
0431 fsig_->setVal(initVal_fsig);
0432 if (a0_ != nullptr)
0433 a0_->setVal(initVal_a0);
0434 if (a1_ != nullptr)
0435 a1_->setVal(initVal_a1);
0436 if (a2_ != nullptr)
0437 a2_->setVal(initVal_a2);
0438 if (a3_ != nullptr)
0439 a3_->setVal(initVal_a3);
0440 if (a4_ != nullptr)
0441 a4_->setVal(initVal_a4);
0442 if (a5_ != nullptr)
0443 a5_->setVal(initVal_a5);
0444 if (a6_ != nullptr)
0445 a6_->setVal(initVal_a6);
0446 if (alpha_ != nullptr)
0447 alpha_->setVal(initVal_alpha);
0448 if (n_ != nullptr)
0449 n_->setVal(initVal_n);
0450 if (fGCB_ != nullptr)
0451 fGCB_->setVal(initVal_fGCB);
0452 }
0453
0454 inline RooRealVar* mean() { return mean_; }
0455 inline RooRealVar* mean2() { return mean2_; }
0456 inline RooRealVar* mean3() { return mean3_; }
0457 inline RooRealVar* sigma() { return sigma_; }
0458 inline RooRealVar* sigma2() { return sigma2_; }
0459 inline RooRealVar* sigma3() { return sigma3_; }
0460 inline RooRealVar* gamma() { return gamma_; }
0461 inline RooRealVar* gaussFrac() { return gaussFrac_; }
0462 inline RooRealVar* gaussFrac2() { return gaussFrac2_; }
0463 inline RooRealVar* expCoeffa0() { return expCoeffa0_; }
0464 inline RooRealVar* expCoeffa1() { return expCoeffa1_; }
0465 inline RooRealVar* expCoeffa2() { return expCoeffa2_; }
0466 inline RooRealVar* fsig() { return fsig_; }
0467 inline RooRealVar* a0() { return a0_; }
0468 inline RooRealVar* a1() { return a1_; }
0469 inline RooRealVar* a2() { return a2_; }
0470 inline RooRealVar* a3() { return a3_; }
0471 inline RooRealVar* a4() { return a4_; }
0472 inline RooRealVar* a5() { return a5_; }
0473 inline RooRealVar* a6() { return a6_; }
0474 inline RooRealVar* alpha() { return alpha_; }
0475 inline RooRealVar* n() { return n_; }
0476 inline RooRealVar* fGCB() { return fGCB_; }
0477
0478
0479 RooAbsPdf* buildSignalModel(RooRealVar* x, const TString& signalType) {
0480 RooAbsPdf* signal = nullptr;
0481 if (signalType == "gaussian") {
0482
0483 if ((mean_ == nullptr) || (sigma_ == nullptr)) {
0484 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and sigma"
0485 << std::endl;
0486 exit(1);
0487 }
0488 signal = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma_);
0489 } else if (signalType == "doubleGaussian") {
0490
0491 if ((mean_ == nullptr) || (sigma_ == nullptr) || (sigma2_ == nullptr)) {
0492 std::cout
0493 << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and sigma2"
0494 << std::endl;
0495 exit(1);
0496 }
0497 RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_);
0498 RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean_, *sigma2_);
0499 signal = new RooAddModel("doubleGaussian", "double gaussian", RooArgList(*gaussModel, *gaussModel2), *gaussFrac_);
0500 } else if (signalType == "tripleGaussian") {
0501
0502 if ((mean_ == nullptr) || (mean2_ == nullptr) || (mean3_ == nullptr) || (sigma_ == nullptr) ||
0503 (sigma2_ == nullptr) || (sigma3_ == nullptr)) {
0504 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0505 "mean3, sigma, sigma2, sigma3"
0506 << std::endl;
0507 exit(1);
0508 }
0509 RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_);
0510 RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean2_, *sigma2_);
0511 RooGaussModel* gaussModel3 = new RooGaussModel("gaussModel3", "gaussModel3", *x, *mean3_, *sigma3_);
0512 signal = new RooAddModel("tripleGaussian",
0513 "triple gaussian",
0514 RooArgList(*gaussModel, *gaussModel2, *gaussModel3),
0515 RooArgList(*gaussFrac_, *gaussFrac2_));
0516 } else if (signalType == "breitWigner") {
0517
0518 if ((mean_ == nullptr) || (gamma_ == nullptr)) {
0519 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma"
0520 << std::endl;
0521 exit(1);
0522 }
0523 signal = new RooBreitWigner("breiWign", "breitWign", *x, *mean_, *gamma_);
0524 } else if (signalType == "relBreitWigner") {
0525
0526 if ((mean_ == nullptr) || (gamma_ == nullptr)) {
0527 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma"
0528 << std::endl;
0529 exit(1);
0530 }
0531 signal = new RooGenericPdf("Relativistic Breit-Wigner",
0532 "RBW",
0533 "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0534 RooArgList(*x, *mean_, *gamma_));
0535 } else if (signalType == "voigtian") {
0536
0537 if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr)) {
0538 std::cout
0539 << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and gamma"
0540 << std::endl;
0541 exit(1);
0542 }
0543 signal = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0544 } else if (signalType == "crystalBall") {
0545
0546 if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr)) {
0547 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0548 "alpha and n"
0549 << std::endl;
0550 exit(1);
0551 }
0552 signal = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0553 } else if (signalType == "breitWignerTimesCB") {
0554
0555 if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) ||
0556 (alpha_ == nullptr) || (n_ == nullptr)) {
0557 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0558 "sigma, gamma, alpha and n"
0559 << std::endl;
0560 exit(1);
0561 }
0562 RooAbsPdf* bw = new RooBreitWigner("breiWigner", "breitWigner", *x, *mean_, *gamma_);
0563 RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0564 signal = new RooFFTConvPdf("breitWignerTimesCB", "breitWignerTimesCB", *x, *bw, *cb);
0565 } else if (signalType == "relBreitWignerTimesCB") {
0566
0567 if ((mean_ == nullptr) || (mean2_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) ||
0568 (alpha_ == nullptr) || (n_ == nullptr)) {
0569 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0570 "sigma, gamma, alpha and n"
0571 << std::endl;
0572 exit(1);
0573 }
0574 RooGenericPdf* bw = new RooGenericPdf("Relativistic Breit-Wigner",
0575 "RBW",
0576 "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0577 RooArgList(*x, *mean_, *gamma_));
0578 RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0579 signal = new RooFFTConvPdf("relBreitWignerTimesCB", "relBreitWignerTimesCB", *x, *bw, *cb);
0580 } else if (signalType == "gaussianPlusCrystalBall") {
0581
0582 if ((mean_ == nullptr) || (sigma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) ||
0583 (fGCB_ == nullptr)) {
0584 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0585 "sigma2, alpha, n and fGCB"
0586 << std::endl;
0587 exit(1);
0588 }
0589 RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0590 RooAbsPdf* tempGaussian = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma2_);
0591
0592 signal = new RooAddPdf(
0593 "gaussianPlusCrystalBall", "gaussianPlusCrystalBall", RooArgList(*tempCB, *tempGaussian), *fGCB_);
0594 } else if (signalType == "voigtianPlusCrystalBall") {
0595
0596 if ((mean_ == nullptr) || (sigma_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) ||
0597 (sigma2_ == nullptr) || (fGCB_ == nullptr)) {
0598 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0599 "sigma, sigma2, alpha, n and fGCB"
0600 << std::endl;
0601 exit(1);
0602 }
0603 RooAbsPdf* tempVoigt = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0604 RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0605
0606 signal =
0607 new RooAddPdf("voigtianPlusCrystalBall", "voigtianPlusCrystalBall", RooArgList(*tempCB, *tempVoigt), *fGCB_);
0608 } else if (signalType == "breitWignerPlusCrystalBall") {
0609
0610 if ((mean_ == nullptr) || (gamma_ == nullptr) || (alpha_ == nullptr) || (n_ == nullptr) || (sigma2_ == nullptr) ||
0611 (fGCB_ == nullptr)) {
0612 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0613 "sigma, alpha, n and fGCB"
0614 << std::endl;
0615 exit(1);
0616 }
0617 RooAbsPdf* tempBW = new RooBreitWigner("breitWign", "breitWign", *x, *mean_, *gamma_);
0618 RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0619
0620 signal = new RooAddPdf(
0621 "breitWignerPlusCrystalBall", "breitWignerPlusCrystalBall", RooArgList(*tempCB, *tempBW), *fGCB_);
0622 }
0623
0624 else if (signalType != "") {
0625 std::cout << "Unknown signal function: " << signalType << ". Signal will not be in the model" << std::endl;
0626 }
0627 return signal;
0628 }
0629
0630
0631 RooAbsPdf* buildBackgroundModel(RooRealVar* x, const TString& backgroundType) {
0632 RooAbsPdf* background = nullptr;
0633 if (backgroundType == "exponential") {
0634
0635 if ((expCoeffa1_ == nullptr) || (fsig_ == nullptr)) {
0636 std::cout
0637 << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeffa1 and fsig"
0638 << std::endl;
0639 exit(1);
0640 }
0641 background = new RooExponential("exponential", "exponential", *x, *expCoeffa1_);
0642 }
0643
0644 if (backgroundType == "exponentialpol") {
0645
0646 if ((expCoeffa0_ == nullptr) || (expCoeffa1_ == nullptr) || (expCoeffa2_ == nullptr) || (fsig_ == nullptr)) {
0647 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeff and fsig"
0648 << std::endl;
0649 exit(1);
0650 }
0651 background = new RooGenericPdf("exponential",
0652 "exponential",
0653 "TMath::Exp(@1+@2*@0+@3*@0*@0)",
0654 RooArgList(*x, *expCoeffa0_, *expCoeffa1_, *expCoeffa2_));
0655 }
0656
0657 else if (backgroundType == "chebychev0") {
0658
0659 if (a0_ == nullptr) {
0660 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0" << std::endl;
0661 exit(1);
0662 }
0663 background = new RooChebychev("chebychev0", "chebychev0", *x, *a0_);
0664 } else if (backgroundType == "chebychev1") {
0665
0666 if ((a0_ == nullptr) || (a1_ == nullptr)) {
0667 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0 and a1"
0668 << std::endl;
0669 exit(1);
0670 }
0671 background = new RooChebychev("chebychev1", "chebychev1", *x, RooArgList(*a0_, *a1_));
0672 } else if (backgroundType == "chebychev3") {
0673
0674 if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr)) {
0675 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2 and a3"
0676 << std::endl;
0677 exit(1);
0678 }
0679 background = new RooChebychev("3rdOrderPol", "3rdOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_));
0680 }
0681
0682 else if (backgroundType == "chebychev6") {
0683
0684 if ((a0_ == nullptr) || (a1_ == nullptr) || (a2_ == nullptr) || (a3_ == nullptr) || (a4_ == nullptr) ||
0685 (a5_ == nullptr) || (a6_ == nullptr)) {
0686 std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2, a3, "
0687 "a4, a5 and a6"
0688 << std::endl;
0689 exit(1);
0690 }
0691 background =
0692 new RooChebychev("6thOrderPol", "6thOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_, *a4_, *a5_, *a6_));
0693 }
0694
0695 return background;
0696 }
0697
0698
0699 RooAbsPdf* buildModel(RooRealVar* x, const TString& signalType, const TString& backgroundType) {
0700 RooAbsPdf* model = nullptr;
0701
0702 RooAbsPdf* signal = buildSignalModel(x, signalType);
0703 RooAbsPdf* background = buildBackgroundModel(x, backgroundType);
0704
0705 if ((signal != nullptr) && (background != nullptr)) {
0706
0707 std::cout << "Building model with signal and backgound" << std::endl;
0708 model = new RooAddPdf("model", "model", RooArgList(*signal, *background), *fsig_);
0709 } else if (signal != nullptr) {
0710 std::cout << "Building model with signal" << std::endl;
0711 model = signal;
0712 } else if (background != nullptr) {
0713 std::cout << "Building model with backgound" << std::endl;
0714 model = background;
0715 } else {
0716 std::cout << "Nothing to fit" << std::endl;
0717 exit(0);
0718 }
0719 return model;
0720 }
0721
0722 bool useChi2_;
0723
0724 protected:
0725
0726 RooRealVar* mean_;
0727 RooRealVar* mean2_;
0728 RooRealVar* mean3_;
0729 RooRealVar* sigma_;
0730 RooRealVar* sigma2_;
0731 RooRealVar* sigma3_;
0732 RooRealVar* gamma_;
0733 RooRealVar* gaussFrac_;
0734 RooRealVar* gaussFrac2_;
0735 RooRealVar* expCoeffa0_;
0736 RooRealVar* expCoeffa1_;
0737 RooRealVar* expCoeffa2_;
0738 RooRealVar* fsig_;
0739 RooRealVar* a0_;
0740 RooRealVar* a1_;
0741 RooRealVar* a2_;
0742 RooRealVar* a3_;
0743 RooRealVar* a4_;
0744 RooRealVar* a5_;
0745 RooRealVar* a6_;
0746 RooRealVar* alpha_;
0747 RooRealVar* n_;
0748 RooRealVar* fGCB_;
0749
0750
0751 double initVal_mean;
0752 double initVal_mean2;
0753 double initVal_mean3;
0754 double initVal_sigma;
0755 double initVal_sigma2;
0756 double initVal_sigma3;
0757 double initVal_gamma;
0758 double initVal_gaussFrac;
0759 double initVal_gaussFrac2;
0760 double initVal_expCoeffa0;
0761 double initVal_expCoeffa1;
0762 double initVal_expCoeffa2;
0763 double initVal_fsig;
0764 double initVal_a0;
0765 double initVal_a1;
0766 double initVal_a2;
0767 double initVal_a3;
0768 double initVal_a4;
0769 double initVal_a5;
0770 double initVal_a6;
0771 double initVal_alpha;
0772 double initVal_n;
0773 double initVal_fGCB;
0774 };
0775
0776 #endif