Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-17 01:47:48

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 "RooChi2Var.h"
0026 #include "RooMinimizer.h"
0027 #include "RooBreitWigner.h"
0028 #include "RooFFTConvPdf.h"
0029 
0030 /**
0031  * This macro allows to use RooFit to perform a fit on a TH1 histogram. <br>
0032  * The currently implemented functions are: <br>
0033  * - signal: <br>
0034  * -- gaussian <br>
0035  * -- double gaussian <br>
0036  * -- voigtian <br>
0037  * - background: <br>
0038  * -- exponential <br>
0039  * It is possible to select any combination of signal and background. <br>
0040  * The fit() method receives the TH1, two strings specifying the signal and background type
0041  * and the min and max x of the histogram over which to perform the fit. <br>
0042  * The variables of the fit must be initialized separately. For example when doing a gaussian+exponential
0043  * fit the initMean, initSigma, initFSig and initExpCoeff methods must be used to create and initialize
0044  * the corresponding variables. <br>
0045  * The methods names after the variables return the fit result.
0046  */
0047 
0048 namespace {
0049   typedef std::pair<RooRealVar, RooDataHist*> rooPair;
0050 }
0051 
0052 class FitWithRooFit {
0053 public:
0054   FitWithRooFit()
0055       : useChi2_(false),
0056         mean_(0),
0057         mean2_(0),
0058         mean3_(0),
0059         sigma_(0),
0060         sigma2_(0),
0061         sigma3_(0),
0062         gamma_(0),
0063         gaussFrac_(0),
0064         gaussFrac2_(0),
0065         expCoeffa0_(0),
0066         expCoeffa1_(0),
0067         expCoeffa2_(0),
0068         fsig_(0),
0069         a0_(0),
0070         a1_(0),
0071         a2_(0),
0072         a3_(0),
0073         a4_(0),
0074         a5_(0),
0075         a6_(0),
0076         alpha_(0),
0077         n_(0),
0078         fGCB_(0) {}
0079 
0080   // Import TH1 histogram into a RooDataHist
0081   rooPair importTH1(TH1* histo, const double& inputXmin, const double& inputXmax) {
0082     double xMin = inputXmin;
0083     double xMax = inputXmax;
0084     if ((xMin == xMax) && xMin == 0) {
0085       xMin = histo->GetXaxis()->GetXmin();
0086       xMax = histo->GetXaxis()->GetXmax();
0087     }
0088     // Declare observable x
0089     RooRealVar x("x", "x", xMin, xMax);
0090     // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
0091     return (std::make_pair(x, new RooDataHist("dh", "dh", x, RooFit::Import(*histo))));
0092   }
0093 
0094   // Plot and fit a RooDataHist fitting signal and background
0095   void fit(TH1* histo,
0096            const TString signalType,
0097            const TString backgroundType,
0098            const double& xMin = 0.,
0099            const double& xMax = 0.,
0100            bool sumW2Error = false) {
0101     reinitializeParameters();
0102 
0103     rooPair imported = importTH1(histo, xMin, xMax);
0104     RooRealVar x(imported.first);
0105     RooDataHist* dh = imported.second;
0106 
0107     // Make plot of binned dataset showing Poisson error bars (RooFit default)
0108     RooPlot* frame = x.frame(RooFit::Title("Imported TH1 with Poisson error bars"));
0109     frame->SetName(TString(histo->GetName()) + "_frame");
0110     dh->plotOn(frame);
0111 
0112     // Build the composite model
0113     RooAbsPdf* model = buildModel(&x, signalType, backgroundType);
0114 
0115     RooChi2Var chi2("chi2", "chi2", *model, *dh, RooFit::DataError(RooAbsData::SumW2));
0116 
0117     // Fit the composite model
0118     // -----------------------
0119     // Fit with likelihood
0120     if (!useChi2_) {
0121       if (sumW2Error)
0122         model->fitTo(*dh, RooFit::Save(), RooFit::SumW2Error(kTRUE));
0123       else
0124         model->fitTo(*dh);
0125     }
0126     // Fit with chi^2
0127     else {
0128       std::cout << "FITTING WITH CHI^2" << std::endl;
0129       RooMinimizer m(chi2);
0130       m.migrad();
0131       m.hesse();
0132       // RooFitResult* r_chi2_wgt = m.save();
0133     }
0134     model->plotOn(frame);
0135     model->plotOn(frame, RooFit::Components(backgroundType), RooFit::LineStyle(kDotted), RooFit::LineColor(kRed));
0136     model->paramOn(frame, RooFit::Label("fit result"), RooFit::Format("NEU", RooFit::AutoPrecision(2)));
0137 
0138     // TODO: fix next lines to get the prob(chi2) (ndof should be dynamically set according to the choosen pdf)
0139     // double chi2 = xframe.chiSquare("model","data",ndof);
0140     // double ndoff = xframeGetNbinsX();
0141     // double chi2prob = TMath::Prob(chi2,ndoff);
0142 
0143     // P l o t   a n d   f i t   a   R o o D a t a H i s t   w i t h   i n t e r n a l   e r r o r s
0144     // ---------------------------------------------------------------------------------------------
0145 
0146     // If histogram has custom error (i.e. its contents is does not originate from a Poisson process
0147     // but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
0148     // (same error bars as shown by ROOT)
0149     RooPlot* frame2 = x.frame(RooFit::Title("Imported TH1 with internal errors"));
0150     dh->plotOn(frame2, RooFit::DataError(RooAbsData::SumW2));
0151     model->plotOn(frame2);
0152     model->plotOn(frame2, RooFit::Components(backgroundType), RooFit::LineColor(kRed));
0153     model->paramOn(frame2, RooFit::Label("fit result"), RooFit::Format("NEU", RooFit::AutoPrecision(2)));
0154 
0155     // Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used
0156     // in a maximum likelihood fit
0157     //
0158     // A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition
0159     // of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly
0160     // fitted with a chi^2 fit (see rf602_chi2fit.C)
0161 
0162     // Draw all frames on a canvas
0163     // if( canvas == 0 ) {
0164     //   canvas = new TCanvas("rf102_dataimport","rf102_dataimport",800,800);
0165     //   canvas->cd(1);
0166     // }
0167     // canvas->Divide(2,1);
0168     // canvas->cd(1) ; frame->Draw();
0169     // canvas->cd(2) ; frame2->Draw();
0170 
0171     frame2->Draw();
0172   }
0173 
0174   // Initialization methods for all the parameters
0175   void initMean(const double& value,
0176                 const double& min,
0177                 const double& max,
0178                 const TString& name = "mean",
0179                 const TString& title = "mean") {
0180     if (mean_ != 0)
0181       delete mean_;
0182     mean_ = new RooRealVar(name, title, value, min, max);
0183     initVal_mean = value;
0184   }
0185   void initMean2(const double& value,
0186                  const double& min,
0187                  const double& max,
0188                  const TString& name = "mean2",
0189                  const TString& title = "mean2") {
0190     if (mean2_ != 0)
0191       delete mean2_;
0192     mean2_ = new RooRealVar(name, title, value, min, max);
0193     initVal_mean2 = value;
0194   }
0195   void initMean3(const double& value,
0196                  const double& min,
0197                  const double& max,
0198                  const TString& name = "mean3",
0199                  const TString& title = "mean3") {
0200     if (mean3_ != 0)
0201       delete mean3_;
0202     mean3_ = new RooRealVar(name, title, value, min, max);
0203     initVal_mean3 = value;
0204   }
0205   void initSigma(const double& value,
0206                  const double& min,
0207                  const double& max,
0208                  const TString& name = "sigma",
0209                  const TString& title = "sigma") {
0210     if (sigma_ != 0)
0211       delete sigma_;
0212     sigma_ = new RooRealVar(name, title, value, min, max);
0213     initVal_sigma = value;
0214   }
0215   void initSigma2(const double& value,
0216                   const double& min,
0217                   const double& max,
0218                   const TString& name = "sigma2",
0219                   const TString& title = "sigma2") {
0220     if (sigma2_ != 0)
0221       delete sigma2_;
0222     sigma2_ = new RooRealVar(name, title, value, min, max);
0223     initVal_sigma2 = value;
0224   }
0225   void initSigma3(const double& value,
0226                   const double& min,
0227                   const double& max,
0228                   const TString& name = "sigma3",
0229                   const TString& title = "sigma3") {
0230     if (sigma3_ != 0)
0231       delete sigma3_;
0232     sigma3_ = new RooRealVar(name, title, value, min, max);
0233     initVal_sigma3 = value;
0234   }
0235   void initGamma(const double& value,
0236                  const double& min,
0237                  const double& max,
0238                  const TString& name = "gamma",
0239                  const TString& title = "gamma") {
0240     if (gamma_ != 0)
0241       delete gamma_;
0242     gamma_ = new RooRealVar(name, title, value, min, max);
0243     initVal_gamma = value;
0244   }
0245   void initGaussFrac(const double& value,
0246                      const double& min,
0247                      const double& max,
0248                      const TString& name = "GaussFrac",
0249                      const TString& title = "GaussFrac") {
0250     if (gaussFrac_ != 0)
0251       delete gaussFrac_;
0252     gaussFrac_ = new RooRealVar(name, title, value, min, max);
0253     initVal_gaussFrac = value;
0254   }
0255   void initGaussFrac2(const double& value,
0256                       const double& min,
0257                       const double& max,
0258                       const TString& name = "GaussFrac2",
0259                       const TString& title = "GaussFrac2") {
0260     if (gaussFrac2_ != 0)
0261       delete gaussFrac2_;
0262     gaussFrac2_ = new RooRealVar(name, title, value, min, max);
0263     initVal_gaussFrac2 = value;
0264   }
0265   void initExpCoeffA0(const double& value,
0266                       const double& min,
0267                       const double& max,
0268                       const TString& name = "expCoeffa0",
0269                       const TString& title = "expCoeffa0") {
0270     if (expCoeffa0_ != 0)
0271       delete expCoeffa0_;
0272     expCoeffa0_ = new RooRealVar(name, title, value, min, max);
0273     initVal_expCoeffa0 = value;
0274   }
0275   void initExpCoeffA1(const double& value,
0276                       const double& min,
0277                       const double& max,
0278                       const TString& name = "expCoeffa1",
0279                       const TString& title = "expCoeffa1") {
0280     if (expCoeffa1_ != 0)
0281       delete expCoeffa1_;
0282     expCoeffa1_ = new RooRealVar(name, title, value, min, max);
0283     initVal_expCoeffa1 = value;
0284   }
0285   void initExpCoeffA2(const double& value,
0286                       const double& min,
0287                       const double& max,
0288                       const TString& name = "expCoeffa2",
0289                       const TString& title = "expCoeffa2") {
0290     if (expCoeffa2_ != 0)
0291       delete expCoeffa2_;
0292     expCoeffa2_ = new RooRealVar(name, title, value, min, max);
0293     initVal_expCoeffa2 = value;
0294   }
0295   void initFsig(const double& value,
0296                 const double& min,
0297                 const double& max,
0298                 const TString& name = "fsig",
0299                 const TString& title = "signal fraction") {
0300     if (fsig_ != 0)
0301       delete fsig_;
0302     fsig_ = new RooRealVar(name, title, value, min, max);
0303     initVal_fsig = value;
0304   }
0305   void initA0(const double& value,
0306               const double& min,
0307               const double& max,
0308               const TString& name = "a0",
0309               const TString& title = "a0") {
0310     if (a0_ != 0)
0311       delete a0_;
0312     a0_ = new RooRealVar(name, title, value, min, max);
0313     initVal_a0 = value;
0314   }
0315   void initA1(const double& value,
0316               const double& min,
0317               const double& max,
0318               const TString& name = "a1",
0319               const TString& title = "a1") {
0320     if (a1_ != 0)
0321       delete a1_;
0322     a1_ = new RooRealVar(name, title, value, min, max);
0323     initVal_a1 = value;
0324   }
0325   void initA2(const double& value,
0326               const double& min,
0327               const double& max,
0328               const TString& name = "a2",
0329               const TString& title = "a2") {
0330     if (a2_ != 0)
0331       delete a2_;
0332     a2_ = new RooRealVar(name, title, value, min, max);
0333     initVal_a2 = value;
0334   }
0335   void initA3(const double& value,
0336               const double& min,
0337               const double& max,
0338               const TString& name = "a3",
0339               const TString& title = "a3") {
0340     if (a3_ != 0)
0341       delete a3_;
0342     a3_ = new RooRealVar(name, title, value, min, max);
0343     initVal_a3 = value;
0344   }
0345   void initA4(const double& value,
0346               const double& min,
0347               const double& max,
0348               const TString& name = "a4",
0349               const TString& title = "a4") {
0350     if (a4_ != 0)
0351       delete a4_;
0352     a4_ = new RooRealVar(name, title, value, min, max);
0353     initVal_a4 = value;
0354   }
0355   void initA5(const double& value,
0356               const double& min,
0357               const double& max,
0358               const TString& name = "a5",
0359               const TString& title = "a5") {
0360     if (a5_ != 0)
0361       delete a5_;
0362     a5_ = new RooRealVar(name, title, value, min, max);
0363     initVal_a5 = value;
0364   }
0365   void initA6(const double& value,
0366               const double& min,
0367               const double& max,
0368               const TString& name = "a6",
0369               const TString& title = "a6") {
0370     if (a6_ != 0)
0371       delete a6_;
0372     a6_ = new RooRealVar(name, title, value, min, max);
0373     initVal_a6 = value;
0374   }
0375   void initAlpha(const double& value,
0376                  const double& min,
0377                  const double& max,
0378                  const TString& name = "alpha",
0379                  const TString& title = "alpha") {
0380     if (alpha_ != 0)
0381       delete alpha_;
0382     alpha_ = new RooRealVar(name, title, value, min, max);
0383     initVal_alpha = value;
0384   }
0385   void initN(const double& value,
0386              const double& min,
0387              const double& max,
0388              const TString& name = "n",
0389              const TString& title = "n") {
0390     if (n_ != 0)
0391       delete n_;
0392     n_ = new RooRealVar(name, title, value, min, max);
0393     initVal_n = value;
0394   }
0395   void initFGCB(const double& value,
0396                 const double& min,
0397                 const double& max,
0398                 const TString& name = "fGCB",
0399                 const TString& title = "fGCB") {
0400     if (fGCB_ != 0)
0401       delete fGCB_;
0402     fGCB_ = new RooRealVar(name, title, value, min, max);
0403     initVal_fGCB = value;
0404   }
0405 
0406   void reinitializeParameters() {
0407     if (mean_ != 0)
0408       mean_->setVal(initVal_mean);
0409     if (mean2_ != 0)
0410       mean2_->setVal(initVal_mean2);
0411     if (mean3_ != 0)
0412       mean3_->setVal(initVal_mean3);
0413     if (sigma_ != 0)
0414       sigma_->setVal(initVal_sigma);
0415     if (sigma2_ != 0)
0416       sigma2_->setVal(initVal_sigma2);
0417     if (sigma3_ != 0)
0418       sigma3_->setVal(initVal_sigma3);
0419     if (gamma_ != 0)
0420       gamma_->setVal(initVal_gamma);
0421     if (gaussFrac_ != 0)
0422       gaussFrac_->setVal(initVal_gaussFrac);
0423     if (gaussFrac2_ != 0)
0424       gaussFrac2_->setVal(initVal_gaussFrac2);
0425     if (expCoeffa0_ != 0)
0426       expCoeffa0_->setVal(initVal_expCoeffa0);
0427     if (expCoeffa1_ != 0)
0428       expCoeffa1_->setVal(initVal_expCoeffa1);
0429     if (expCoeffa2_ != 0)
0430       expCoeffa2_->setVal(initVal_expCoeffa2);
0431     if (fsig_ != 0)
0432       fsig_->setVal(initVal_fsig);
0433     if (a0_ != 0)
0434       a0_->setVal(initVal_a0);
0435     if (a1_ != 0)
0436       a1_->setVal(initVal_a1);
0437     if (a2_ != 0)
0438       a2_->setVal(initVal_a2);
0439     if (a3_ != 0)
0440       a3_->setVal(initVal_a3);
0441     if (a4_ != 0)
0442       a4_->setVal(initVal_a4);
0443     if (a5_ != 0)
0444       a5_->setVal(initVal_a5);
0445     if (a6_ != 0)
0446       a6_->setVal(initVal_a6);
0447     if (alpha_ != 0)
0448       alpha_->setVal(initVal_alpha);
0449     if (n_ != 0)
0450       n_->setVal(initVal_n);
0451     if (fGCB_ != 0)
0452       fGCB_->setVal(initVal_fGCB);
0453   }
0454 
0455   inline RooRealVar* mean() { return mean_; }
0456   inline RooRealVar* mean2() { return mean2_; }
0457   inline RooRealVar* mean3() { return mean3_; }
0458   inline RooRealVar* sigma() { return sigma_; }
0459   inline RooRealVar* sigma2() { return sigma2_; }
0460   inline RooRealVar* sigma3() { return sigma3_; }
0461   inline RooRealVar* gamma() { return gamma_; }
0462   inline RooRealVar* gaussFrac() { return gaussFrac_; }
0463   inline RooRealVar* gaussFrac2() { return gaussFrac2_; }
0464   inline RooRealVar* expCoeffa0() { return expCoeffa0_; }
0465   inline RooRealVar* expCoeffa1() { return expCoeffa1_; }
0466   inline RooRealVar* expCoeffa2() { return expCoeffa2_; }
0467   inline RooRealVar* fsig() { return fsig_; }
0468   inline RooRealVar* a0() { return a0_; }
0469   inline RooRealVar* a1() { return a1_; }
0470   inline RooRealVar* a2() { return a2_; }
0471   inline RooRealVar* a3() { return a3_; }
0472   inline RooRealVar* a4() { return a4_; }
0473   inline RooRealVar* a5() { return a5_; }
0474   inline RooRealVar* a6() { return a6_; }
0475   inline RooRealVar* alpha() { return alpha_; }
0476   inline RooRealVar* n() { return n_; }
0477   inline RooRealVar* fGCB() { return fGCB_; }
0478 
0479   /// Build the model for the specified signal type
0480   RooAbsPdf* buildSignalModel(RooRealVar* x, const TString& signalType) {
0481     RooAbsPdf* signal = 0;
0482     if (signalType == "gaussian") {
0483       // Fit a Gaussian p.d.f to the data
0484       if ((mean_ == 0) || (sigma_ == 0)) {
0485         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and sigma"
0486                   << std::endl;
0487         exit(1);
0488       }
0489       signal = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma_);
0490     } else if (signalType == "doubleGaussian") {
0491       // Fit with double gaussian
0492       if ((mean_ == 0) || (sigma_ == 0) || (sigma2_ == 0)) {
0493         std::cout
0494             << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and sigma2"
0495             << std::endl;
0496         exit(1);
0497       }
0498       RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_);
0499       RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean_, *sigma2_);
0500       signal = new RooAddModel("doubleGaussian", "double gaussian", RooArgList(*gaussModel, *gaussModel2), *gaussFrac_);
0501     } else if (signalType == "tripleGaussian") {
0502       // Fit with triple gaussian
0503       if ((mean_ == 0) || (mean2_ == 0) || (mean3_ == 0) || (sigma_ == 0) || (sigma2_ == 0) || (sigma3_ == 0)) {
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       // Fit a Breit-Wigner
0518       if ((mean_ == 0) || (gamma_ == 0)) {
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       // Fit a relativistic Breit-Wigner
0526       if ((mean_ == 0) || (gamma_ == 0)) {
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       // Fit a Voigtian
0537       if ((mean_ == 0) || (sigma_ == 0) || (gamma_ == 0)) {
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       // Fit a CrystalBall
0546       if ((mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
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       // Fit a Breit Wigner convoluted with a CrystalBall
0555       if ((mean_ == 0) || (mean2_ == 0) || (sigma_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
0556         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0557                      "sigma, gamma, alpha and n"
0558                   << std::endl;
0559         exit(1);
0560       }
0561       RooAbsPdf* bw = new RooBreitWigner("breiWigner", "breitWigner", *x, *mean_, *gamma_);
0562       RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0563       signal = new RooFFTConvPdf("breitWignerTimesCB", "breitWignerTimesCB", *x, *bw, *cb);
0564     } else if (signalType == "relBreitWignerTimesCB") {
0565       // Fit a relativistic Breit Wigner convoluted with a CrystalBall
0566       if ((mean_ == 0) || (mean2_ == 0) || (sigma_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
0567         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0568                      "sigma, gamma, alpha and n"
0569                   << std::endl;
0570         exit(1);
0571       }
0572       RooGenericPdf* bw = new RooGenericPdf("Relativistic Breit-Wigner",
0573                                             "RBW",
0574                                             "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0575                                             RooArgList(*x, *mean_, *gamma_));
0576       RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0577       signal = new RooFFTConvPdf("relBreitWignerTimesCB", "relBreitWignerTimesCB", *x, *bw, *cb);
0578     } else if (signalType == "gaussianPlusCrystalBall") {
0579       // Fit a Gaussian + CrystalBall with the same mean
0580       if ((mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0) || (sigma2_ == 0) || (fGCB_ == 0)) {
0581         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0582                      "sigma2, alpha, n and fGCB"
0583                   << std::endl;
0584         exit(1);
0585       }
0586       RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0587       RooAbsPdf* tempGaussian = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma2_);
0588 
0589       signal = new RooAddPdf(
0590           "gaussianPlusCrystalBall", "gaussianPlusCrystalBall", RooArgList(*tempCB, *tempGaussian), *fGCB_);
0591     } else if (signalType == "voigtianPlusCrystalBall") {
0592       // Fit a Voigtian + CrystalBall with the same mean
0593       if ((mean_ == 0) || (sigma_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0) || (sigma2_ == 0) ||
0594           (fGCB_ == 0)) {
0595         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0596                      "sigma, sigma2, alpha, n and fGCB"
0597                   << std::endl;
0598         exit(1);
0599       }
0600       RooAbsPdf* tempVoigt = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0601       RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0602 
0603       signal =
0604           new RooAddPdf("voigtianPlusCrystalBall", "voigtianPlusCrystalBall", RooArgList(*tempCB, *tempVoigt), *fGCB_);
0605     } else if (signalType == "breitWignerPlusCrystalBall") {
0606       // Fit a Breit-Wigner + CrystalBall with the same mean
0607       if ((mean_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0) || (sigma2_ == 0) || (fGCB_ == 0)) {
0608         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0609                      "sigma, alpha, n and fGCB"
0610                   << std::endl;
0611         exit(1);
0612       }
0613       RooAbsPdf* tempBW = new RooBreitWigner("breitWign", "breitWign", *x, *mean_, *gamma_);
0614       RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0615 
0616       signal = new RooAddPdf(
0617           "breitWignerPlusCrystalBall", "breitWignerPlusCrystalBall", RooArgList(*tempCB, *tempBW), *fGCB_);
0618     }
0619 
0620     else if (signalType != "") {
0621       std::cout << "Unknown signal function: " << signalType << ". Signal will not be in the model" << std::endl;
0622     }
0623     return signal;
0624   }
0625 
0626   /// Build the model for the specified background type
0627   RooAbsPdf* buildBackgroundModel(RooRealVar* x, const TString& backgroundType) {
0628     RooAbsPdf* background = 0;
0629     if (backgroundType == "exponential") {
0630       // Add an exponential for the background
0631       if ((expCoeffa1_ == 0) || (fsig_ == 0)) {
0632         std::cout
0633             << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeffa1 and fsig"
0634             << std::endl;
0635         exit(1);
0636       }
0637       background = new RooExponential("exponential", "exponential", *x, *expCoeffa1_);
0638     }
0639 
0640     if (backgroundType == "exponentialpol") {
0641       // Add an exponential for the background
0642       if ((expCoeffa0_ == 0) || (expCoeffa1_ == 0) || (expCoeffa2_ == 0) || (fsig_ == 0)) {
0643         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeff and fsig"
0644                   << std::endl;
0645         exit(1);
0646       }
0647       background = new RooGenericPdf("exponential",
0648                                      "exponential",
0649                                      "TMath::Exp(@1+@2*@0+@3*@0*@0)",
0650                                      RooArgList(*x, *expCoeffa0_, *expCoeffa1_, *expCoeffa2_));
0651     }
0652 
0653     else if (backgroundType == "chebychev0") {
0654       // Add a linear background
0655       if (a0_ == 0) {
0656         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0" << std::endl;
0657         exit(1);
0658       }
0659       background = new RooChebychev("chebychev0", "chebychev0", *x, *a0_);
0660     } else if (backgroundType == "chebychev1") {
0661       // Add a 2nd order chebychev polynomial background
0662       if ((a0_ == 0) || (a1_ == 0)) {
0663         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0 and a1"
0664                   << std::endl;
0665         exit(1);
0666       }
0667       background = new RooChebychev("chebychev1", "chebychev1", *x, RooArgList(*a0_, *a1_));
0668     } else if (backgroundType == "chebychev3") {
0669       // Add a 3rd order chebychev polynomial background
0670       if ((a0_ == 0) || (a1_ == 0) || (a2_ == 0) || (a3_ == 0)) {
0671         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2 and a3"
0672                   << std::endl;
0673         exit(1);
0674       }
0675       background = new RooChebychev("3rdOrderPol", "3rdOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_));
0676     }
0677 
0678     else if (backgroundType == "chebychev6") {
0679       // Add a 6th order chebychev polynomial background
0680       if ((a0_ == 0) || (a1_ == 0) || (a2_ == 0) || (a3_ == 0) || (a4_ == 0) || (a5_ == 0) || (a6_ == 0)) {
0681         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2, a3, "
0682                      "a4, a5 and a6"
0683                   << std::endl;
0684         exit(1);
0685       }
0686       background =
0687           new RooChebychev("6thOrderPol", "6thOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_, *a4_, *a5_, *a6_));
0688     }
0689 
0690     return background;
0691   }
0692 
0693   /// Build the model to fit
0694   RooAbsPdf* buildModel(RooRealVar* x, const TString& signalType, const TString& backgroundType) {
0695     RooAbsPdf* model = 0;
0696 
0697     RooAbsPdf* signal = buildSignalModel(x, signalType);
0698     RooAbsPdf* background = buildBackgroundModel(x, backgroundType);
0699 
0700     if ((signal != 0) && (background != 0)) {
0701       // Combine signal and background pdfs
0702       std::cout << "Building model with signal and backgound" << std::endl;
0703       model = new RooAddPdf("model", "model", RooArgList(*signal, *background), *fsig_);
0704     } else if (signal != 0) {
0705       std::cout << "Building model with signal" << std::endl;
0706       model = signal;
0707     } else if (background != 0) {
0708       std::cout << "Building model with backgound" << std::endl;
0709       model = background;
0710     } else {
0711       std::cout << "Nothing to fit" << std::endl;
0712       exit(0);
0713     }
0714     return model;
0715   }
0716 
0717   bool useChi2_;
0718 
0719 protected:
0720   // Declare all variables
0721   RooRealVar* mean_;
0722   RooRealVar* mean2_;
0723   RooRealVar* mean3_;
0724   RooRealVar* sigma_;
0725   RooRealVar* sigma2_;
0726   RooRealVar* sigma3_;
0727   RooRealVar* gamma_;
0728   RooRealVar* gaussFrac_;
0729   RooRealVar* gaussFrac2_;
0730   RooRealVar* expCoeffa0_;
0731   RooRealVar* expCoeffa1_;
0732   RooRealVar* expCoeffa2_;
0733   RooRealVar* fsig_;
0734   RooRealVar* a0_;
0735   RooRealVar* a1_;
0736   RooRealVar* a2_;
0737   RooRealVar* a3_;
0738   RooRealVar* a4_;
0739   RooRealVar* a5_;
0740   RooRealVar* a6_;
0741   RooRealVar* alpha_;
0742   RooRealVar* n_;
0743   RooRealVar* fGCB_;
0744 
0745   // Initial values
0746   double initVal_mean;
0747   double initVal_mean2;
0748   double initVal_mean3;
0749   double initVal_sigma;
0750   double initVal_sigma2;
0751   double initVal_sigma3;
0752   double initVal_gamma;
0753   double initVal_gaussFrac;
0754   double initVal_gaussFrac2;
0755   double initVal_expCoeffa0;
0756   double initVal_expCoeffa1;
0757   double initVal_expCoeffa2;
0758   double initVal_fsig;
0759   double initVal_a0;
0760   double initVal_a1;
0761   double initVal_a2;
0762   double initVal_a3;
0763   double initVal_a4;
0764   double initVal_a5;
0765   double initVal_a6;
0766   double initVal_alpha;
0767   double initVal_n;
0768   double initVal_fGCB;
0769 };
0770 
0771 #endif