Back to home page

Project CMSSW displayed by LXR
 
 

    

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  * 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_(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   // 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::Save(), RooFit::SumW2Error(kTRUE));
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_ != 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   /// Build the model for the specified signal type
0479   RooAbsPdf* buildSignalModel(RooRealVar* x, const TString& signalType) {
0480     RooAbsPdf* signal = nullptr;
0481     if (signalType == "gaussian") {
0482       // Fit a Gaussian p.d.f to the data
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       // Fit with double gaussian
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       // Fit with triple gaussian
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       // Fit a Breit-Wigner
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       // Fit a relativistic Breit-Wigner
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       // Fit a Voigtian
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       // Fit a CrystalBall
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       // Fit a Breit Wigner convoluted with a CrystalBall
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       // Fit a relativistic Breit Wigner convoluted with a CrystalBall
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       // Fit a Gaussian + CrystalBall with the same mean
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       // Fit a Voigtian + CrystalBall with the same mean
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       // Fit a Breit-Wigner + CrystalBall with the same mean
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   /// Build the model for the specified background type
0631   RooAbsPdf* buildBackgroundModel(RooRealVar* x, const TString& backgroundType) {
0632     RooAbsPdf* background = nullptr;
0633     if (backgroundType == "exponential") {
0634       // Add an exponential for the background
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       // Add an exponential for the background
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       // Add a linear background
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       // Add a 2nd order chebychev polynomial background
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       // Add a 3rd order chebychev polynomial background
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       // Add a 6th order chebychev polynomial background
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   /// Build the model to fit
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       // Combine signal and background pdfs
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   // Declare all variables
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   // Initial values
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
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!