Skip to content

Commit

Permalink
Added custom Double-Crystal-Ball pdf.
Browse files Browse the repository at this point in the history
 - added the required files to allow storage in a ROOT file
 - created workspace is now stored in ROOT file
 - fit with this function is not yet properly converging
  • Loading branch information
gregor-mittag committed Nov 24, 2015
1 parent c3d69d5 commit 1d79bf4
Show file tree
Hide file tree
Showing 7 changed files with 250 additions and 23 deletions.
26 changes: 16 additions & 10 deletions interface/FitContainer.h
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,8 @@
#include <string>
#include <memory>
#include "TH1.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
Expand All @@ -24,33 +26,37 @@ namespace analysis {
virtual ~FitContainer();

inline static const std::vector<std::string>& availableModels() {
return availableModels_;
};

return availableModels_; };
void setModel(const std::string& type, const std::string& model);
void setModel(const std::string& type, const std::string& model,
const std::vector<ParamModifier>& modifiers);
std::unique_ptr<RooFitResult> backgroundOnlyFit();

private:
std::string getOutputPath_(const std::string& subdirectory = "");
bool checkType_(const std::string& type);
double getMaxPosition_(const RooDataHist& data);
int getNonZeroBins_(const RooDataHist& data);
bool applyModifiers_(RooAbsPdf& pdf,
const std::vector<ParamModifier>& modifiers);

// methods to set the fit model
static std::unique_ptr<RooArgList>
getCoefficients_(const int numCoeffs, const std::string& name);
void setNovosibirsk_(const std::string& type);
void setCrystalBall_(const std::string& type);
void setNovoEffProd_(const std::string& type);
void setExpEffProd_(const std::string& type);
void setDoubleCB_(const std::string& type);
void setBernstein_(const std::string& type, const int numCoeffs);
void setChebychev_(const std::string& type, const int numCoeffs);
void setBernEffProd_(const std::string& type, const int numCoeffs);
static const std::vector<std::string> availableModels_;

// internal methods
std::string getOutputPath_(const std::string& subdirectory = "");
void prepareCanvas_(TCanvas& raw);
void prepareFrame_(RooPlot& raw);
bool checkType_(const std::string& type);
double getMaxPosition_(const RooDataHist& data);
int getNonZeroBins_(const RooDataHist& data);
bool applyModifiers_(RooAbsPdf& pdf,
const std::vector<ParamModifier>& modifiers);

// data member
static const int defaultNumberOfCoefficients_;
std::string plotDir_;
RooWorkspace workspace_;
Expand Down
47 changes: 47 additions & 0 deletions interface/RooDoubleCB.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
#ifndef Analysis_BackgroundModel_RooDoubleCB_h
#define Analysis_BackgroundModel_RooDoubleCB_h 1

#include "RooAbsPdf.h"
#include "RooRealProxy.h"
#include "RooAbsReal.h"

namespace analysis {
namespace backgroundmodel {

class RooDoubleCB : public RooAbsPdf {
public:
inline RooDoubleCB() {};
RooDoubleCB(const char* name, const char* title,
RooAbsReal& x,
RooAbsReal& mean,
RooAbsReal& width,
RooAbsReal& alpha1,
RooAbsReal& n1,
RooAbsReal& alpha2,
RooAbsReal& n2
);
RooDoubleCB(const RooDoubleCB& other, const char* name = 0) ;
virtual TObject* clone(const char* newname) const;
inline virtual ~RooDoubleCB() {};
int getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars,
const char* name = 0) const;
double analyticalIntegral(int code, const char* rangeName = 0) const;

protected:
double evaluate() const;

RooRealProxy x_;
RooRealProxy mean_;
RooRealProxy width_;
RooRealProxy alpha1_;
RooRealProxy n1_;
RooRealProxy alpha2_;
RooRealProxy n2_;

ClassDef(RooDoubleCB,1);
};

}
}

#endif // Analysis_BackgroundModel_RooDoubleCB_h
62 changes: 49 additions & 13 deletions src/FitContainer.cc
Original file line number Diff line number Diff line change
Expand Up @@ -3,11 +3,8 @@
#include <boost/algorithm/string.hpp>
#include <boost/algorithm/string/join.hpp>
#include "TSystem.h"
#include "TCanvas.h"
#include "TLatex.h"
#include "RooPlot.h"
#include "RooArgList.h"
#include "RooRealVar.h"
#include "RooFormulaVar.h"
#include "RooAbsPdf.h"
#include "RooEffProd.h"
Expand All @@ -17,37 +14,42 @@
#include "RooChebychev.h"
#include "RooNovosibirsk.h"
#include "RooCBShape.h"
#include "Analysis/BackgroundModel/interface/RooDoubleCB.h"
#include "Analysis/BackgroundModel/interface/FitContainer.h"
#include "Analysis/BackgroundModel/interface/Tools.h"


using namespace analysis::backgroundmodel;


FitContainer::FitContainer(TH1& data, TH1& signal, TH1& background)
: plotDir_(getOutputPath_("plots")), workspace_(RooWorkspace("workspace")),
mbb_("mbb", "m_{bb}",
data.GetXaxis()->GetXmin(), data.GetXaxis()->GetXmax(), "GeV"),
data_("data_container", "", mbb_, RooFit::Import(data)),
signal_("signal_container", "", mbb_, RooFit::Import(signal)),
background_("background_container", "", mbb_, RooFit::Import(background)) {
FitContainer::FitContainer(TH1& data, TH1& signal, TH1& background) :
plotDir_(getOutputPath_("plots")),
workspace_(RooWorkspace("workspace")),
mbb_("mbb", "m_{bb}",
data.GetXaxis()->GetXmin(), data.GetXaxis()->GetXmax(), "GeV"),
data_("data_container", "", mbb_, RooFit::Import(data)),
signal_("signal_container", "", mbb_, RooFit::Import(signal)),
background_("background_container", "", mbb_, RooFit::Import(background)) {

// some preliminary test code
std::unique_ptr<RooPlot> frame(mbb_.frame());
data_.plotOn(frame.get());
TCanvas canvas("canvas", "", 600, 600);
canvas.cd();
prepareCanvas_(canvas);
prepareFrame_(*frame);
frame->Draw();
canvas.SaveAs((plotDir_+"input_data.pdf").c_str());
}


FitContainer::FitContainer(const DataContainer& container)
: FitContainer(*(container.data()), *(container.bbH()), *(container.background())) {
FitContainer::FitContainer(const DataContainer& container) :
FitContainer(*(container.data()), *(container.bbH()), *(container.background())) {
}


FitContainer::~FitContainer() {
workspace_.writeToFile((getOutputPath_("workspace")+"workspace.root").c_str());
}


Expand Down Expand Up @@ -83,6 +85,7 @@ void FitContainer::setModel(const std::string& type, const std::string& name,
else if (nameSplitted[0] == "crystalball") setCrystalBall_(type);
else if (nameSplitted[0] == "novoeffprod") setNovoEffProd_(type);
else if (nameSplitted[0] == "expeffprod") setExpEffProd_(type);
else if (nameSplitted[0] == "doublecb") setDoubleCB_(type);
else if (nameSplitted[0] == "bernstein") setBernstein_(type, numCoeffs);
else if (nameSplitted[0] == "chebychev") setChebychev_(type, numCoeffs);
else if (nameSplitted[0] == "berneffprod") setBernEffProd_(type, numCoeffs);
Expand All @@ -109,11 +112,13 @@ std::unique_ptr<RooFitResult> FitContainer::backgroundOnlyFit() {
std::cout << "\nfloating parameters (final):" << std::endl;
fitResult->floatParsFinal().Print("v");

std::unique_ptr<RooPlot> frame(mbb_.frame(RooFit::Title(" ")));
std::unique_ptr<RooPlot> frame(mbb_.frame());
data_.plotOn(frame.get(), RooFit::Name("data_curve"));
bkg.plotOn(frame.get(), RooFit::Name("background_curve"));
TCanvas canvas("canvas", "", 600, 600);
canvas.cd();
prepareCanvas_(canvas);
prepareFrame_(*frame);
frame->Draw();

int nPars = fitResult->floatParsFinal().getSize();
Expand Down Expand Up @@ -199,6 +204,24 @@ void FitContainer::setExpEffProd_(const std::string& type) {
}


void FitContainer::setDoubleCB_(const std::string& type) {
double meanStart = (mbb_.getMin() + mbb_.getMax()) / 2.0;
if (type == "signal") meanStart = getMaxPosition_(signal_);
else if (type == "background") meanStart = getMaxPosition_(background_);
RooRealVar mean("mean", "", meanStart, mbb_.getMin(), mbb_.getMax(), "GeV");
RooRealVar width("width", "", 35.0, 5.0, 100.0, "GeV");
RooRealVar alpha1("alpha1", "", -1.0, 0.0);
RooRealVar n1("n1", "", 1.0);
RooRealVar alpha2("alpha2", "", -1.0, 0.0);
RooRealVar n2("n2", "", 1.0);
n1.setConstant(false);
n2.setConstant(false);
RooDoubleCB doubleCB(type.c_str(), (type+"_doublecb").c_str(),
mbb_, mean, width, alpha1, n1, alpha2, n2);
workspace_.import(doubleCB);
}


void FitContainer::setBernstein_(const std::string& type, const int numCoeffs) {
std::unique_ptr<RooArgList> coeffs(getCoefficients_(numCoeffs, "bernstein"));
RooBernstein bern(type.c_str(), (type+"_bernstein").c_str(), mbb_, *coeffs);
Expand Down Expand Up @@ -240,6 +263,18 @@ std::string FitContainer::getOutputPath_(const std::string& subdirectory) {
}


void FitContainer::prepareCanvas_(TCanvas& raw) {
raw.SetLeftMargin(0.15);
raw.SetRightMargin(0.05);
}


void FitContainer::prepareFrame_(RooPlot& raw) {
raw.GetYaxis()->SetTitleOffset(2);
raw.SetTitle("");
}


bool FitContainer::checkType_(const std::string& type) {
// update this list if needed
const std::vector<std::string> allowedTypes = {"signal", "background"};
Expand Down Expand Up @@ -318,6 +353,7 @@ const std::vector<std::string> FitContainer::availableModels_ =
"crystalball",
"novoeffprod",
"expeffprod",
"doublecb",
"bernstein",
"chebychev",
"berneffprod"};
Expand Down
133 changes: 133 additions & 0 deletions src/RooDoubleCB.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
#include <iostream>
#include <cmath>
#include "TMath.h"

#include "RooRealVar.h"
#include "RooRealConstant.h"
#include "Analysis/BackgroundModel/interface/RooDoubleCB.h"


using namespace analysis::backgroundmodel;


ClassImp(RooDoubleCB)


RooDoubleCB::RooDoubleCB(const char* name, const char* title,
RooAbsReal& x,
RooAbsReal& mean,
RooAbsReal& width,
RooAbsReal& alpha1,
RooAbsReal& n1,
RooAbsReal& alpha2,
RooAbsReal& n2
) :
RooAbsPdf(name, title),
x_("x", "x", this, x),
mean_("mean", "mean", this, mean),
width_("width", "width", this, width),
alpha1_("alpha1", "alpha1", this, alpha1),
n1_("n1", "n1", this, n1),
alpha2_("alpha2", "alpha2", this, alpha2),
n2_("n2", "n2", this, n2) {
}


RooDoubleCB::RooDoubleCB(const RooDoubleCB& other, const char* name) :
RooAbsPdf(other, name),
x_("x", this, other.x_),
mean_("mean", this, other.mean_),
width_("width", this, other.width_),
alpha1_("alpha1", this, other.alpha1_),
n1_("n1", this, other.n1_),
alpha2_("alpha2", this, other.alpha2_),
n2_("n2", this, other.n2_) {
}


TObject* RooDoubleCB::clone(const char* newname) const {
return new RooDoubleCB(*this, newname);
}


double RooDoubleCB::evaluate() const {
double t = (x_ - mean_)/width_;
if (t > -alpha1_ && t < alpha2_) {
return std::exp(-0.5*t*t);
} else if (t < -alpha1_) {
double A1 = std::pow(n1_/fabs(alpha1_), n1_) * std::exp(-alpha1_*alpha1_/2);
double B1 = n1_/std::fabs(alpha1_) - std::fabs(alpha1_);
return A1 * std::pow(B1-t, -n1_);
} else if (t > alpha2_) {
double A2 = std::pow(n2_/fabs(alpha2_), n2_) * std::exp(-alpha2_*alpha2_/2);
double B2 = n2_/std::fabs(alpha2_)-std::fabs(alpha2_);
return A2 * std::pow(B2+t, -n2_);
} else {
std::cerr << "ERROR evaluating range..." << std::endl;
return 99;
}
}


int RooDoubleCB::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars,
const char* range) const {
if (matchArgs(allVars, analVars, x_)) return 1;
return 0;
}


double RooDoubleCB::analyticalIntegral(int code,
const char* rangeName) const {
assert(code == 1) ;

double central = 0;
double left = 0;
double right = 0;

static const double root2 = std::sqrt(2) ;
static const double rootPiBy2 = std::sqrt(std::atan2(0.0, -1.0) / 2.0);
double xscale = root2 * width_;

//compute gaussian contribution
double centralLow = std::max(x_.min(rangeName), mean_ - alpha1_ * width_);
double centralHigh = std::min(x_.max(rangeName), mean_ + alpha2_ * width_);
if (centralLow < centralHigh) // is the gaussian part in range?
central = rootPiBy2 * width_ * (TMath::Erf((centralHigh - mean_)/xscale) -
TMath::Erf((centralLow - mean_)/xscale));

//compute left tail;
double A1 = (std::pow(n1_/std::fabs(alpha1_), n1_) *
std::exp(-alpha1_*alpha1_/2));
double B1 = n1_/std::fabs(alpha1_) - std::fabs(alpha1_);

double leftLow = x_.min(rangeName);
double leftHigh = std::min(x_.max(rangeName), mean_ - alpha1_*width_);
if (leftLow < leftHigh) { //is the left tail in range?
if (std::fabs(n1_ - 1.0) > 1.e-5)
left = A1/(-n1_+1.0) * width_ *
(std::pow(B1-(leftLow-mean_)/width_,-n1_+1.) -
std::pow(B1-(leftHigh-mean_)/width_,-n1_+1.));
else
left = A1 * width_ * (std::log(B1-(leftLow-mean_)/width_) -
std::log(B1-(leftHigh-mean_)/width_));
}

//compute right tail;
double A2 = (std::pow(n2_/std::fabs(alpha2_), n2_) *
std::exp(-alpha2_*alpha2_/2));
double B2 = n2_/std::fabs(alpha2_) - std::fabs(alpha2_);

double rightLow = std::max(x_.min(rangeName), mean_ + alpha2_*width_);
double rightHigh = x_.max(rangeName);
if (rightLow < rightHigh) { //is the right tail in range?
if (std::fabs(n2_-1.0)>1.e-5)
right = A2/(-n2_+1.0) * width_ *
(std::pow(B2+(rightHigh-mean_)/width_, -n2_+1.) -
std::pow(B2 + (rightLow-mean_)/width_, -n2_+1.));
else
right = A2 * width_ * (std::log(B2+(rightHigh-mean_)/width_) -
std::log(B2+(rightLow-mean_)/width_));
}

return left + central + right;
}
1 change: 1 addition & 0 deletions src/classes.h
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
#include "Analysis/BackgroundModel/interface/RooDoubleCB.h"
3 changes: 3 additions & 0 deletions src/classes_def.xml
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
<lcgdict>
<class name="analysis::backgroundmodel::RooDoubleCB" />
</lcgdict>
1 change: 1 addition & 0 deletions test/.gitignore
Original file line number Diff line number Diff line change
@@ -1 +1,2 @@
plots
workspace

0 comments on commit 1d79bf4

Please sign in to comment.