Macros/RooFit/FitWithRooFit.cc

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  * This macro allows to use RooFit to perform a fit on a TH1 histogram. <br>
0031  * The currently implemented functions are: <br>
0032  * - signal: <br>
0033  * -- gaussian <br>
0034  * -- double gaussian <br>
0035  * -- voigtian <br>
0036  * - background: <br>
0037  * -- exponential <br>
0038  * It is possible to select any combination of signal and background. <br>
0039  * The fit() method receives the TH1, two strings specifying the signal and background type
0040  * and the min and max x of the histogram over which to perform the fit. <br>
0041  * The variables of the fit must be initialized separately. For example when doing a gaussian+exponential
0042  * fit the initMean, initSigma, initFSig and initExpCoeff methods must be used to create and initialize
0043  * the corresponding variables. <br>
0044  * The methods names after the variables return the fit result.
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_(0),
0056         mean2_(0),
0057         mean3_(0),
0058         sigma_(0),
0059         sigma2_(0),
0060         sigma3_(0),
0061         gamma_(0),
0062         gaussFrac_(0),
0063         gaussFrac2_(0),
0064         expCoeffa0_(0),
0065         expCoeffa1_(0),
0066         expCoeffa2_(0),
0067         fsig_(0),
0068         a0_(0),
0069         a1_(0),
0070         a2_(0),
0071         a3_(0),
0072         a4_(0),
0073         a5_(0),
0074         a6_(0),
0075         alpha_(0),
0076         n_(0),
0077         fGCB_(0) {}
0078
0079   // Import TH1 histogram into a RooDataHist
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     // Declare observable x
0088     RooRealVar x("x", "x", xMin, xMax);
0089     // Create a binned dataset that imports contents of TH1 and associates its contents to observable 'x'
0090     return (std::make_pair(x, new RooDataHist("dh", "dh", x, RooFit::Import(*histo))));
0091   }
0092
0093   // Plot and fit a RooDataHist fitting signal and background
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     // Make plot of binned dataset showing Poisson error bars (RooFit default)
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     // Build the composite model
0112     RooAbsPdf* model = buildModel(&x, signalType, backgroundType);
0113
0114     std::unique_ptr<RooAbsReal> chi2{model->createChi2(*dh, RooFit::DataError(RooAbsData::SumW2))};
0115
0116     // Fit the composite model
0117     // -----------------------
0118     // Fit with likelihood
0119     if (!useChi2_) {
0120       if (sumW2Error)
0121         model->fitTo(*dh, RooFit::SumW2Error(true));
0122       else
0123         model->fitTo(*dh);
0124     }
0125     // Fit with chi^2
0126     else {
0127       std::cout << "FITTING WITH CHI^2" << std::endl;
0128       RooMinimizer m(*chi2);
0129       m.migrad();
0130       m.hesse();
0131       // RooFitResult* r_chi2_wgt = m.save();
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     // TODO: fix next lines to get the prob(chi2) (ndof should be dynamically set according to the choosen pdf)
0138     // double chi2 = xframe.chiSquare("model","data",ndof);
0139     // double ndoff = xframeGetNbinsX();
0140     // double chi2prob = TMath::Prob(chi2,ndoff);
0141
0142     // 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
0143     // ---------------------------------------------------------------------------------------------
0144
0145     // If histogram has custom error (i.e. its contents is does not originate from a Poisson process
0146     // but e.g. is a sum of weighted events) you can data with symmetric 'sum-of-weights' error instead
0147     // (same error bars as shown by ROOT)
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     // Please note that error bars shown (Poisson or SumW2) are for visualization only, the are NOT used
0155     // in a maximum likelihood fit
0156     //
0157     // A (binned) ML fit will ALWAYS assume the Poisson error interpretation of data (the mathematical definition
0158     // of likelihood does not take any external definition of errors). Data with non-unit weights can only be correctly
0159     // fitted with a chi^2 fit (see rf602_chi2fit.C)
0160
0161     // Draw all frames on a canvas
0162     // if( canvas == 0 ) {
0163     //   canvas = new TCanvas("rf102_dataimport","rf102_dataimport",800,800);
0164     //   canvas->cd(1);
0165     // }
0166     // canvas->Divide(2,1);
0167     // canvas->cd(1) ; frame->Draw();
0168     // canvas->cd(2) ; frame2->Draw();
0169
0170     frame2->Draw();
0171   }
0172
0173   // Initialization methods for all the parameters
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
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_ != 0)
0400       delete fGCB_;
0401     fGCB_ = new RooRealVar(name, title, value, min, max);
0402     initVal_fGCB = value;
0403   }
0404
0405   void reinitializeParameters() {
0406     if (mean_ != 0)
0407       mean_->setVal(initVal_mean);
0408     if (mean2_ != 0)
0409       mean2_->setVal(initVal_mean2);
0410     if (mean3_ != 0)
0411       mean3_->setVal(initVal_mean3);
0412     if (sigma_ != 0)
0413       sigma_->setVal(initVal_sigma);
0414     if (sigma2_ != 0)
0415       sigma2_->setVal(initVal_sigma2);
0416     if (sigma3_ != 0)
0417       sigma3_->setVal(initVal_sigma3);
0418     if (gamma_ != 0)
0419       gamma_->setVal(initVal_gamma);
0420     if (gaussFrac_ != 0)
0421       gaussFrac_->setVal(initVal_gaussFrac);
0422     if (gaussFrac2_ != 0)
0423       gaussFrac2_->setVal(initVal_gaussFrac2);
0424     if (expCoeffa0_ != 0)
0425       expCoeffa0_->setVal(initVal_expCoeffa0);
0426     if (expCoeffa1_ != 0)
0427       expCoeffa1_->setVal(initVal_expCoeffa1);
0428     if (expCoeffa2_ != 0)
0429       expCoeffa2_->setVal(initVal_expCoeffa2);
0430     if (fsig_ != 0)
0431       fsig_->setVal(initVal_fsig);
0432     if (a0_ != 0)
0433       a0_->setVal(initVal_a0);
0434     if (a1_ != 0)
0435       a1_->setVal(initVal_a1);
0436     if (a2_ != 0)
0437       a2_->setVal(initVal_a2);
0438     if (a3_ != 0)
0439       a3_->setVal(initVal_a3);
0440     if (a4_ != 0)
0441       a4_->setVal(initVal_a4);
0442     if (a5_ != 0)
0443       a5_->setVal(initVal_a5);
0444     if (a6_ != 0)
0445       a6_->setVal(initVal_a6);
0446     if (alpha_ != 0)
0447       alpha_->setVal(initVal_alpha);
0448     if (n_ != 0)
0449       n_->setVal(initVal_n);
0450     if (fGCB_ != 0)
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   /// Build the model for the specified signal type
0479   RooAbsPdf* buildSignalModel(RooRealVar* x, const TString& signalType) {
0480     RooAbsPdf* signal = 0;
0481     if (signalType == "gaussian") {
0482       // Fit a Gaussian p.d.f to the data
0483       if ((mean_ == 0) || (sigma_ == 0)) {
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       // Fit with double gaussian
0491       if ((mean_ == 0) || (sigma_ == 0) || (sigma2_ == 0)) {
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       // Fit with triple gaussian
0502       if ((mean_ == 0) || (mean2_ == 0) || (mean3_ == 0) || (sigma_ == 0) || (sigma2_ == 0) || (sigma3_ == 0)) {
0503         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0504                      "mean3, sigma, sigma2, sigma3"
0505                   << std::endl;
0506         exit(1);
0507       }
0508       RooGaussModel* gaussModel = new RooGaussModel("gaussModel", "gaussModel", *x, *mean_, *sigma_);
0509       RooGaussModel* gaussModel2 = new RooGaussModel("gaussModel2", "gaussModel2", *x, *mean2_, *sigma2_);
0510       RooGaussModel* gaussModel3 = new RooGaussModel("gaussModel3", "gaussModel3", *x, *mean3_, *sigma3_);
0511       signal = new RooAddModel("tripleGaussian",
0512                                "triple gaussian",
0513                                RooArgList(*gaussModel, *gaussModel2, *gaussModel3),
0514                                RooArgList(*gaussFrac_, *gaussFrac2_));
0515     } else if (signalType == "breitWigner") {
0516       // Fit a Breit-Wigner
0517       if ((mean_ == 0) || (gamma_ == 0)) {
0518         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma"
0519                   << std::endl;
0520         exit(1);
0521       }
0522       signal = new RooBreitWigner("breiWign", "breitWign", *x, *mean_, *gamma_);
0523     } else if (signalType == "relBreitWigner") {
0524       // Fit a relativistic Breit-Wigner
0525       if ((mean_ == 0) || (gamma_ == 0)) {
0526         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean and gamma"
0527                   << std::endl;
0528         exit(1);
0529       }
0530       signal = new RooGenericPdf("Relativistic Breit-Wigner",
0531                                  "RBW",
0532                                  "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0533                                  RooArgList(*x, *mean_, *gamma_));
0534     } else if (signalType == "voigtian") {
0535       // Fit a Voigtian
0536       if ((mean_ == 0) || (sigma_ == 0) || (gamma_ == 0)) {
0537         std::cout
0538             << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma and gamma"
0539             << std::endl;
0540         exit(1);
0541       }
0542       signal = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0543     } else if (signalType == "crystalBall") {
0544       // Fit a CrystalBall
0545       if ((mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
0546         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0547                      "alpha and n"
0548                   << std::endl;
0549         exit(1);
0550       }
0551       signal = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0552     } else if (signalType == "breitWignerTimesCB") {
0553       // Fit a Breit Wigner convoluted with a CrystalBall
0554       if ((mean_ == 0) || (mean2_ == 0) || (sigma_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
0555         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0556                      "sigma, gamma, alpha and n"
0557                   << std::endl;
0558         exit(1);
0559       }
0560       RooAbsPdf* bw = new RooBreitWigner("breiWigner", "breitWigner", *x, *mean_, *gamma_);
0561       RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0562       signal = new RooFFTConvPdf("breitWignerTimesCB", "breitWignerTimesCB", *x, *bw, *cb);
0563     } else if (signalType == "relBreitWignerTimesCB") {
0564       // Fit a relativistic Breit Wigner convoluted with a CrystalBall
0565       if ((mean_ == 0) || (mean2_ == 0) || (sigma_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0)) {
0566         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, mean2, "
0567                      "sigma, gamma, alpha and n"
0568                   << std::endl;
0569         exit(1);
0570       }
0571       RooGenericPdf* bw = new RooGenericPdf("Relativistic Breit-Wigner",
0572                                             "RBW",
0573                                             "@0/(pow(@0*@0 - @1*@1,2) + @2*@2*@0*@0*@0*@0/(@1*@1))",
0574                                             RooArgList(*x, *mean_, *gamma_));
0575       RooAbsPdf* cb = new RooCBShape("crystalBall", "crystalBall", *x, *mean2_, *sigma_, *alpha_, *n_);
0576       signal = new RooFFTConvPdf("relBreitWignerTimesCB", "relBreitWignerTimesCB", *x, *bw, *cb);
0577     } else if (signalType == "gaussianPlusCrystalBall") {
0578       // Fit a Gaussian + CrystalBall with the same mean
0579       if ((mean_ == 0) || (sigma_ == 0) || (alpha_ == 0) || (n_ == 0) || (sigma2_ == 0) || (fGCB_ == 0)) {
0580         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, sigma, "
0581                      "sigma2, alpha, n and fGCB"
0582                   << std::endl;
0583         exit(1);
0584       }
0585       RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma_, *alpha_, *n_);
0586       RooAbsPdf* tempGaussian = new RooGaussian("gauss", "gauss", *x, *mean_, *sigma2_);
0587
0588       signal = new RooAddPdf(
0589           "gaussianPlusCrystalBall", "gaussianPlusCrystalBall", RooArgList(*tempCB, *tempGaussian), *fGCB_);
0590     } else if (signalType == "voigtianPlusCrystalBall") {
0591       // Fit a Voigtian + CrystalBall with the same mean
0592       if ((mean_ == 0) || (sigma_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0) || (sigma2_ == 0) ||
0593           (fGCB_ == 0)) {
0594         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0595                      "sigma, sigma2, alpha, n and fGCB"
0596                   << std::endl;
0597         exit(1);
0598       }
0599       RooAbsPdf* tempVoigt = new RooVoigtian("voigt", "voigt", *x, *mean_, *gamma_, *sigma_);
0600       RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0601
0602       signal =
0603           new RooAddPdf("voigtianPlusCrystalBall", "voigtianPlusCrystalBall", RooArgList(*tempCB, *tempVoigt), *fGCB_);
0604     } else if (signalType == "breitWignerPlusCrystalBall") {
0605       // Fit a Breit-Wigner + CrystalBall with the same mean
0606       if ((mean_ == 0) || (gamma_ == 0) || (alpha_ == 0) || (n_ == 0) || (sigma2_ == 0) || (fGCB_ == 0)) {
0607         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize mean, gamma, "
0608                      "sigma, alpha, n and fGCB"
0609                   << std::endl;
0610         exit(1);
0611       }
0612       RooAbsPdf* tempBW = new RooBreitWigner("breitWign", "breitWign", *x, *mean_, *gamma_);
0613       RooAbsPdf* tempCB = new RooCBShape("crystalBall", "crystalBall", *x, *mean_, *sigma2_, *alpha_, *n_);
0614
0615       signal = new RooAddPdf(
0616           "breitWignerPlusCrystalBall", "breitWignerPlusCrystalBall", RooArgList(*tempCB, *tempBW), *fGCB_);
0617     }
0618
0619     else if (signalType != "") {
0620       std::cout << "Unknown signal function: " << signalType << ". Signal will not be in the model" << std::endl;
0621     }
0622     return signal;
0623   }
0624
0625   /// Build the model for the specified background type
0626   RooAbsPdf* buildBackgroundModel(RooRealVar* x, const TString& backgroundType) {
0627     RooAbsPdf* background = 0;
0628     if (backgroundType == "exponential") {
0629       // Add an exponential for the background
0630       if ((expCoeffa1_ == 0) || (fsig_ == 0)) {
0631         std::cout
0632             << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeffa1 and fsig"
0633             << std::endl;
0634         exit(1);
0635       }
0636       background = new RooExponential("exponential", "exponential", *x, *expCoeffa1_);
0637     }
0638
0639     if (backgroundType == "exponentialpol") {
0640       // Add an exponential for the background
0641       if ((expCoeffa0_ == 0) || (expCoeffa1_ == 0) || (expCoeffa2_ == 0) || (fsig_ == 0)) {
0642         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize expCoeff and fsig"
0643                   << std::endl;
0644         exit(1);
0645       }
0646       background = new RooGenericPdf("exponential",
0647                                      "exponential",
0648                                      "TMath::Exp(@1+@2*@0+@3*@0*@0)",
0649                                      RooArgList(*x, *expCoeffa0_, *expCoeffa1_, *expCoeffa2_));
0650     }
0651
0652     else if (backgroundType == "chebychev0") {
0653       // Add a linear background
0654       if (a0_ == 0) {
0655         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0" << std::endl;
0656         exit(1);
0657       }
0658       background = new RooChebychev("chebychev0", "chebychev0", *x, *a0_);
0659     } else if (backgroundType == "chebychev1") {
0660       // Add a 2nd order chebychev polynomial background
0661       if ((a0_ == 0) || (a1_ == 0)) {
0662         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0 and a1"
0663                   << std::endl;
0664         exit(1);
0665       }
0666       background = new RooChebychev("chebychev1", "chebychev1", *x, RooArgList(*a0_, *a1_));
0667     } else if (backgroundType == "chebychev3") {
0668       // Add a 3rd order chebychev polynomial background
0669       if ((a0_ == 0) || (a1_ == 0) || (a2_ == 0) || (a3_ == 0)) {
0670         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2 and a3"
0671                   << std::endl;
0672         exit(1);
0673       }
0674       background = new RooChebychev("3rdOrderPol", "3rdOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_));
0675     }
0676
0677     else if (backgroundType == "chebychev6") {
0678       // Add a 6th order chebychev polynomial background
0679       if ((a0_ == 0) || (a1_ == 0) || (a2_ == 0) || (a3_ == 0) || (a4_ == 0) || (a5_ == 0) || (a6_ == 0)) {
0680         std::cout << "Error: one or more parameters are not initialized. Please be sure to initialize a0, a1, a2, a3, "
0681                      "a4, a5 and a6"
0682                   << std::endl;
0683         exit(1);
0684       }
0685       background =
0686           new RooChebychev("6thOrderPol", "6thOrderPol", *x, RooArgList(*a0_, *a1_, *a2_, *a3_, *a4_, *a5_, *a6_));
0687     }
0688
0689     return background;
0690   }
0691
0692   /// Build the model to fit
0693   RooAbsPdf* buildModel(RooRealVar* x, const TString& signalType, const TString& backgroundType) {
0694     RooAbsPdf* model = 0;
0695
0696     RooAbsPdf* signal = buildSignalModel(x, signalType);
0697     RooAbsPdf* background = buildBackgroundModel(x, backgroundType);
0698
0699     if ((signal != 0) && (background != 0)) {
0700       // Combine signal and background pdfs
0701       std::cout << "Building model with signal and backgound" << std::endl;
0702       model = new RooAddPdf("model", "model", RooArgList(*signal, *background), *fsig_);
0703     } else if (signal != 0) {
0704       std::cout << "Building model with signal" << std::endl;
0705       model = signal;
0706     } else if (background != 0) {
0707       std::cout << "Building model with backgound" << std::endl;
0708       model = background;
0709     } else {
0710       std::cout << "Nothing to fit" << std::endl;
0711       exit(0);
0712     }
0713     return model;
0714   }
0715
0716   bool useChi2_;
0717
0718 protected:
0719   // Declare all variables
0720   RooRealVar* mean_;
0721   RooRealVar* mean2_;
0722   RooRealVar* mean3_;
0723   RooRealVar* sigma_;
0724   RooRealVar* sigma2_;
0725   RooRealVar* sigma3_;
0726   RooRealVar* gamma_;
0727   RooRealVar* gaussFrac_;
0728   RooRealVar* gaussFrac2_;
0729   RooRealVar* expCoeffa0_;
0730   RooRealVar* expCoeffa1_;
0731   RooRealVar* expCoeffa2_;
0732   RooRealVar* fsig_;
0733   RooRealVar* a0_;
0734   RooRealVar* a1_;
0735   RooRealVar* a2_;
0736   RooRealVar* a3_;
0737   RooRealVar* a4_;
0738   RooRealVar* a5_;
0739   RooRealVar* a6_;
0740   RooRealVar* alpha_;
0741   RooRealVar* n_;
0742   RooRealVar* fGCB_;
0743
0744   // Initial values
0745   double initVal_mean;
0746   double initVal_mean2;
0747   double initVal_mean3;
0748   double initVal_sigma;
0749   double initVal_sigma2;
0750   double initVal_sigma3;
0751   double initVal_gamma;
0752   double initVal_gaussFrac;
0753   double initVal_gaussFrac2;
0754   double initVal_expCoeffa0;
0755   double initVal_expCoeffa1;
0756   double initVal_expCoeffa2;
0757   double initVal_fsig;
0758   double initVal_a0;
0759   double initVal_a1;
0760   double initVal_a2;
0761   double initVal_a3;
0762   double initVal_a4;
0763   double initVal_a5;
0764   double initVal_a6;
0765   double initVal_alpha;
0766   double initVal_n;
0767   double initVal_fGCB;
0768 };
0769
0770 #endif