conv_linear_test_3.c

#include "Riostream.h"

// Custom two exponential pdf
class TwoExpPdf : public RooAbsPdf {
public:
    TwoExpPdf() {};
    TwoExpPdf(const char *name, const char *title,
        RooAbsReal& _t,
        RooAbsReal& _tau1,
        RooAbsReal& _tau2,
        RooAbsReal& _i2
        );
    TwoExpPdf(const TwoExpPdf& other, const char* name = 0);
    virtual TObject* clone(const char* newname) const { return new TwoExpPdf(*this, newname); }
    inline virtual ~TwoExpPdf() { }

    Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ;
    Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ;
        
protected:
    RooRealProxy t;
    RooRealProxy tau1;
    RooRealProxy tau2;
    RooRealProxy i2;

    Double_t evaluate() const;
    Double_t indefiniteIntegral(Double_t y) const;
private:
	ClassDef(TwoExpPdf, 1)
};

TwoExpPdf::TwoExpPdf(const char *name, const char *title,
	RooAbsReal& _t,
	RooAbsReal& _tau1,
	RooAbsReal& _tau2,
	RooAbsReal& _i2) :
	RooAbsPdf(name, title),
	t("t", "t", this, _t),
	tau1("tau1", "tau1", this, _tau1),
	tau2("tau2", "tau2", this, _tau2),
	i2("i2", "i2", this, _i2){
}

TwoExpPdf::TwoExpPdf(const TwoExpPdf& other, const char* name) :
RooAbsPdf(other, name),
t("t", this, other.t),
tau1("tau1", this, other.tau1),
tau2("tau2", this, other.tau2),
i2("i2", this, other.i2){
}

Double_t TwoExpPdf::evaluate() const {
    if (t < 0) return 0.;
    if (t == 0) return 1;
    return (1. - i2)*exp(-t/tau1) + i2*exp(-t/tau2);
}

Double_t TwoExpPdf::indefiniteIntegral(Double_t y) const {
    if (y < 0) return 0.;
    if (y == 0) return -1/tau1 -1/tau2;
    Double_t int1 = -exp(-y/tau1)*tau1;
    Double_t int2 = -exp(-y/tau2)*tau2;
    return (1. - i2)*int1 + i2*int2;
}

Int_t TwoExpPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const {
   if (matchArgs(allVars,analVars,t)) return 1;
   return 0;
}

Double_t TwoExpPdf::analyticalIntegral(Int_t code, const char* rangeName) const {
    assert(code == 1);
    
    if (code==1){
        Double_t x1 = t.min(rangeName);
        Double_t x2 = t.max(rangeName);

        if (x2 <= 0) return 0;
        x1 = TMath::Max(0.,x1);
        Double_t ret = indefiniteIntegral(x2)-indefiniteIntegral(x1);
        return ret;
    }
    else {
        std::cout << "Error in TwoExpPdf::analyticalIntegral" << std::endl;
    }  
    return 0;
}

// Custom exponential PDF
class ExpPdf : public RooAbsPdf {
public:
  ExpPdf() {};
  ExpPdf(const char *name, const char *title, RooAbsReal& _t, RooAbsReal& _tau);
  ExpPdf(const ExpPdf& other, const char* name = 0);
  virtual TObject* clone(const char* newname) const { return new ExpPdf(*this, newname); }
  inline virtual ~ExpPdf() { }

  Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const ;
  Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const ;
  Double_t indefiniteIntegral(Double_t y) const;

protected:
  RooRealProxy t;
  RooRealProxy tau;
  Double_t evaluate() const;
};

ExpPdf::ExpPdf(const char *name, const char *title, RooAbsReal& _t, RooAbsReal& _tau) :
  RooAbsPdf(name, title), t("t", "t", this, _t), tau("tau", "tau", this, _tau) {
}

ExpPdf::ExpPdf(const ExpPdf& other, const char* name) : RooAbsPdf(other, name), t("t", this, other.t), tau("tau", this, other.tau){
}

Double_t ExpPdf::evaluate() const {
  if (t < 0) return 0.;
  return exp(-t/tau);
}

Double_t ExpPdf::indefiniteIntegral(Double_t y) const {
  return -exp(-y/tau)*tau;
}

Int_t ExpPdf::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName) const {
  if (matchArgs(allVars,analVars,t)) return 1;
  return 0;
}

Double_t ExpPdf::analyticalIntegral(Int_t code, const char* rangeName) const {
  if (code == 1){
    Double_t t1 = t.min(rangeName);
    Double_t t2 = t.max(rangeName);
    if (t2 <= 0) return 0;
    t1 = TMath::Max(0.,t1);
    return indefiniteIntegral(t2) - indefiniteIntegral(t1);
  }
  return 0;
}

// Main macro
void conv_linear_test_3() {
  // Constants
  const Int_t READ_BINS = 3000;
  RooConstVar channelWidth = RooFit::RooConst(0.0061741); // nanoseconds in channel (bin)
  RooConstVar fwhm2dispersion = RooFit::RooConst(0.424661); // coefficient to convert fwhm to dispersion

  // Read histogram from file
  std::ifstream in;
  // in.open("./si-spectrum.dat");
  in.open("./si-spectrum.dat");
  TH1F* hist = new TH1F("hist", "hist", READ_BINS, 0, READ_BINS*channelWidth.getVal());
  for (unsigned i=1; i<=READ_BINS; i++){
    int x; in >> x;
    hist->SetBinContent(i, x+1); // The RooFit chi2 fit does not work when the bins have zero entries.
    hist->SetBinError(i, sqrt(x+1));
  }

  // Create observable
  RooRealVar* t = new RooRealVar("t", "time axis", 0, READ_BINS*channelWidth.getVal(), "ns");
  t->setBins(READ_BINS);
  t->setRange("FULL",t->getMin(),t->getMax()) ;

  // Create data hist
  RooDataHist* data = new RooDataHist("data", "data", RooArgSet(*t), RooFit::Import(*hist));

  // Constant background PDF
  RooPolynomial* bg = new RooPolynomial("bg", "y=1", *t, RooArgSet());
  RooRealVar* I_bg = new RooRealVar("I_bg", "background fraction", 0.07);


  // FIRST TEST: RooAddPdf(ExpPdf + ExpPdf) x gauss
  // Resolution function is RooGaussian
  RooRealVar* mean1 = new RooRealVar("mean1", "gauss1 mean", 5.9, 5.0, 7.0, "ns");
  RooRealVar* fwhm1 = new RooRealVar("fwhm1", "gauss1 fwhm", 0.4, 0.2, 0.8, "ns");
  RooFormulaVar* sigma1 = new RooFormulaVar("sigma1", "gauss1 dispersion", "@0*@1", RooArgList(*fwhm1, fwhm2dispersion));
  RooGaussian* gauss1 = new RooGaussian("gauss1", "gauss pdf", *t, *mean1, *sigma1);

  // Model is sum of two custom ExpPdf's. Their sum is convoluted via RooFFTConvPdf.
  RooRealVar* tauSi1 = new RooRealVar("tauSi1", "positron lifetime in Si", 0.219, "ns");
  RooRealVar* tauKa1 = new RooRealVar("tauKa1", "positron lifetime in Kapton", 0.385, "ns");
  RooRealVar* I_Ka1 = new RooRealVar("I_Ka1", "contribution in Kapton", 0.12, 0.05, 1.0);
  ExpPdf* pdfSi1 = new ExpPdf("pdfSi1", "Silicon pdf", *t, *tauSi1);
  ExpPdf* pdfKa1 = new ExpPdf("pdfKa1", "Kapton pdf", *t, *tauKa1);
  RooAddPdf* pdfSiKa = new RooAddPdf("pdfSiKa", "sum pdf", RooArgList(*pdfKa1, *pdfSi1), *I_Ka1);   // first sum
  pdfSiKa->fixAddCoefNormalization(RooArgSet(*t));
  RooFFTConvPdf* convSiKa1 = new RooFFTConvPdf("convSiKa1", "convolution of sum of pdfs", *t, *pdfSiKa, *gauss1); // then convolute

  // Finally we add same constant background to model
  RooAddPdf* convSiKa1Bg = new RooAddPdf("convSiKa1Bg", "(pdfSi + pdfKa) x gauss", RooArgList(*bg, *convSiKa1), *I_bg);

  // SECOND TEST: (ExpPdf x gauss) + (ExpPdf x gauss)
  // Resolution function is RooGaussian with same parameter starting values
  RooRealVar* mean2 = new RooRealVar("mean2", "gauss2 mean", 5.9, 5.0, 7.0, "ns");
  RooRealVar* fwhm2 = new RooRealVar("fwhm2", "gauss2 fwhm", 0.4, 0.2, 0.8, "ns");
  RooFormulaVar* sigma2 = new RooFormulaVar("sigma2", "gauss2 dispersion", "@0*@1", RooArgList(*fwhm2, fwhm2dispersion));
  RooGaussian* gauss2 = new RooGaussian("gauss2", "gauss pdf", *t, *mean2, *sigma2);

  // Model is sum of two custom ExpPdf's independently convoluted via RooFFTConvPdf
  RooRealVar* tauSi2 = new RooRealVar("tauSi2", "positron lifetime in Si", 0.219, "ns");
  RooRealVar* tauKa2 = new RooRealVar("tauKa2", "positron lifetime in Kapton", 0.385, "ns");
  RooRealVar* I_Ka2 = new RooRealVar("I_Ka2", "contribution in Kapton", 0.12, 0.05, 1.0);
  ExpPdf* pdfSi = new ExpPdf("pdfSi", "Silicon pdf", *t, *tauSi2);
  RooFFTConvPdf* convSi = new RooFFTConvPdf("convSi", "convoluted pdf", *t, *pdfSi, *gauss2);
  ExpPdf* pdfKa = new ExpPdf("pdfKa", "Kapton pdf", *t, *tauKa2);
  RooFFTConvPdf* convKa = new RooFFTConvPdf("convKa", "convoluted pdf", *t, *pdfKa, *gauss2);
  RooAddPdf* convSiKa2 = new RooAddPdf("convSiKa2", "sum of convoluted pdfs", RooArgList(*convKa, *convSi), *I_Ka2);

  // Finally we add same constant background to model
  RooAddPdf* convSiKa2Bg = new RooAddPdf("convSiKa2Bg", "(pdfSi x gauss) + (pdfKa x gauss)", RooArgList(*bg, *convSiKa2), *I_bg);

  // THIRD TEST: TwoExpPdf x gauss
  // Resolution function is RooGaussian
  RooRealVar* mean3 = new RooRealVar("mean3", "gauss3 mean", 5.9, 5.0, 7.0, "ns");
  RooRealVar* fwhm3 = new RooRealVar("fwhm3", "gauss3 fwhm", 0.4, 0.2, 0.8, "ns");
  RooFormulaVar* sigma3 = new RooFormulaVar("sigma3", "gauss3 dispersion", "@0*@1", RooArgList(*fwhm3, fwhm2dispersion));
  RooGaussian* gauss3 = new RooGaussian("gauss3", "gauss pdf", *t, *mean3, *sigma3);

  RooRealVar* tauSi3 = new RooRealVar("tauSi3", "positron lifetime in Si", 0.219, "ns");
  RooRealVar* tauKa3 = new RooRealVar("tauKa3", "positron lifetime in Kapton", 0.385, "ns");
  RooRealVar* I_Ka3 = new RooRealVar("I_Ka3", "contribution in Kapton", 0.12, 0.05, 1.0);
  TwoExpPdf* pdfSiKaSum = new TwoExpPdf("pdfSiKaSum", "Silicon pdf", *t, *tauSi3, *tauKa3, *I_Ka3);
  pdfSiKaSum->fixAddCoefNormalization(RooArgSet(*t)); // does not fix the coefficient??
  RooFFTConvPdf* convSiKa3 = new RooFFTConvPdf("convSiKa3", "convolution of sum of pdfs", *t, *pdfSiKaSum, *gauss3); // then convolute

  // Finally we add same constant background to model
  RooAddPdf* convSiKa3Bg = new RooAddPdf("convSiKa3Bg", "(pdfSi + pdfKa) x gauss", RooArgList(*bg, *convSiKa3), *I_bg);


  // FIT DATA (to the identical dataset)
  convSiKa1Bg->fitTo(*data);
  convSiKa2Bg->fitTo(*data);
  convSiKa3Bg->fitTo(*data);

  // PLOT FITS AND REVEAL THE PROBLEM:
  // I_Ka1 = 7.9% is incorrect parameter value vs correct one I_Ka2 = 13.1%
  TCanvas* canvas = new TCanvas("canvas", "canvas", 1280, 1024);
  canvas->Divide(1, 3);

  // Top plot: convolution of two RooDecay's and RooGaussModel
  canvas->cd(1)->SetLogy();
  RooPlot* frame1 = t->frame(RooFit::Title("RooAddPdf(ExpPdf + ExpPdf) x gauss, gives correct I_Ka2 = 0.131 if I call fixAddCoefNormalization()"));
  data->plotOn(frame1, RooFit::MarkerSize(0.2));
  convSiKa1Bg->plotOn(frame1);
  gauss1->plotOn(frame1, RooFit::LineStyle(kDashed));
  convSiKa1Bg->paramOn(frame1, RooFit::Layout(0.58, 0.9, 0.9), RooFit::Format("NEU", RooFit::AutoPrecision(1)));
  frame1->Draw();

  // Middle plot: convolution of custom ExpPdf's with RooGaussian by means of RooFFTConvPdf
  canvas->cd(2)->SetLogy();
  RooPlot* frame2 = t->frame(RooFit::Title("(ExpPdf x gauss) + (ExpPdf x gauss), gives correct I_Ka2 = 0.131"));
  data->plotOn(frame2, RooFit::MarkerSize(0.2));
  convSiKa2Bg->plotOn(frame2);
  gauss2->plotOn(frame2, RooFit::LineStyle(kDashed));
  convSiKa2Bg->paramOn(frame2, RooFit::Layout(0.58, 0.9, 0.9), RooFit::Format("NEU", RooFit::AutoPrecision(1)));
  frame2->Draw();

  // Bottom plot: convolution of custom ExpPdf's with RooGaussian by means of RooFFTConvPdf
  canvas->cd(3)->SetLogy();
  RooPlot* frame3 = t->frame(RooFit::Title("TwoExpPdf x gauss, gives incorrect I_Ka1 = 0.079 ???"));
  data->plotOn(frame3, RooFit::MarkerSize(0.2));
  convSiKa3Bg->plotOn(frame3);
  gauss3->plotOn(frame3, RooFit::LineStyle(kDashed));
  convSiKa3Bg->paramOn(frame3, RooFit::Layout(0.58, 0.9, 0.9), RooFit::Format("NEU", RooFit::AutoPrecision(1)));
  frame3->Draw();

}