2D interpolation using moment morphing

Roofit and RooStats
1

Dear all,

I would like to perform a 2d interpolation using a moment morphing technique discussed in this paper: https://arxiv.org/pdf/1410.7388.pdf
There is minimal 1D example code in the paper that I wanted to extend to a 2D, however the 1D example code itself doesn’t work for me and I guess a bit of changes is required that I am having trouble to do:
here is the example code:

// This example builds two normal distributions and uses moment morphing to
// interpolate between the templates using RooFit.
using namespace RooFit;

// Create a persistable container for RooFit projects, allowing to use a simplified
// scripting language to build the p.d.f.s needed in this example
RooWorkspace w("w", 1);

// Build two normal distributions, corresponding to different values in the morph
 // parameter space. They share the same observable, but have otherwise different
 // moments, i.e. mean and width.
 w.factory("RooGaussian::gaussian1(obs[0,400],50,15)");
 w.factory("RooGaussian::gaussian2(obs,300,40)");

 // Build a RooMomentMorph p.d.f. which interpolates between the normal distributions
 // created before. The interpolation is parametrized by the parameter alpha and the
 // reference templates map to alpha=0 and alpha=1 respectively.
 w.factory("RooMomentMorph::morphpdf(alpha[0,1],obs,{gaussian1,gaussian2},{0,1})");

 // Set the morphing parameter alpha explicitly to 0.5.
 w::alpha->setVal(0.5);

 // Create a frame to draw the p.d.f. from before and show the input templates as
 // solid blue curves and the moment morph p.d.f. at alpha=0.5 in dashed red.
 RooPlot* frame = w::obs->frame();
 w::gaussian1->plotOn(frame, LineColor(kBlue), LineStyle(kSolid));
 w::gaussian2->plotOn(frame, LineColor(kBlue), LineStyle(kSolid));
 w::morphpdf->plotOn(frame, LineColor(kRed), LineStyle(kDashed));
 frame->Draw();

Actually the trouble comes from this line:

w::alpha->setVal(0.5);

Does any have any idea how to fix that?

Best regards,

Diallo.

2

Hi ,

You should not use w::alpha to access workspace elements.
The correct code is:

w.var("alpha")->setVal(0.5);

and similar later for w::obs. For the pdf you do:

w.pdf("gaussian1")->plotOn(....)

Lorenzo

3

Hi Lorenzo,

Thanks for your answer after the changes you suggested it works but these errors come up which are related to the plotting.

[#1] INFO:ObjectHandling -- RooWorkspace::exportToCint(w) INFO: references to all objects in this workspace will be created in CINT in 'namespace w'
[#0] ERROR:Plotting -- RooGaussian::gaussian1:plotOn: WARNING: variable is not an explicit dependent: alpha
[#0] ERROR:Plotting -- RooGaussian::gaussian1:plotOn: WARNING: variable is not an explicit dependent: alpha
[#0] ERROR:Plotting -- RooGaussian::gaussian1:createPlotProjection: "alpha" is not a dependent and will be ignored.
[#0] ERROR:Plotting -- RooGaussian::gaussian2:plotOn: WARNING: variable is not an explicit dependent: alpha
[#0] ERROR:Plotting -- RooGaussian::gaussian2:plotOn: WARNING: variable is not an explicit dependent: alpha
[#0] ERROR:Plotting -- RooGaussian::gaussian2:createPlotProjection: "alpha" is not a dependent and will be ignored.

The code has been changed like this after your suggestions:

// This example builds two normal distributions and uses moment morphing to
// interpolate between the templates using RooFit.
using namespace RooFit;

// Create a persistable container for RooFit projects, allowing to use a simplified
// scripting language to build the p.d.f.s needed in this example
RooWorkspace w("w", 1);

// Build two normal distributions, corresponding to different values in the morph
 // parameter space. They share the same observable, but have otherwise different
 // moments, i.e. mean and width.
 w.factory("RooGaussian::gaussian1(obs[0,400],50,15)");
 w.factory("RooGaussian::gaussian2(obs,300,40)");

 // Build a RooMomentMorph p.d.f. which interpolates between the normal distributions
 // created before. The interpolation is parametrized by the parameter alpha and the
 // reference templates map to alpha=0 and alpha=1 respectively.
 w.factory("RooMomentMorph::morphpdf(alpha[0,1],obs,{gaussian1,gaussian2},{0,1})");

 // Set the morphing parameter alpha explicitly to 0.5.
 //w::alpha->setVal(0.5);
 w.var("alpha")->setVal(0.5); 

 // Create a frame to draw the p.d.f. from before and show the input templates as
 // solid blue curves and the moment morph p.d.f. at alpha=0.5 in dashed red.
 RooPlot* frame = w.var("alpha")->frame();
 w.pdf("gaussian1")->plotOn(frame, LineColor(kBlue), LineStyle(kSolid));
 w.pdf("gaussian2")->plotOn(frame, LineColor(kBlue), LineStyle(kSolid));
 w.pdf("morphpdf")->plotOn(frame, LineColor(kRed), LineStyle(kDashed));
 frame->Draw();

Best regards,

Diallo.

4

Hi Lorenzo,

I used the obs variable for the frame definition and it works. So I changed this:

w.var("alpha")->frame();

to this:

w.var("obs")->frame();
5

Hi Lorenzo,

So I moved to 2D with these changes:

 //first mass point
 w.factory("RooGaussian::gaussian1(obs1[0,400],50,15)"); //mll
 w.factory("RooGaussian::gaussian2(obs2[0,400],300,40)"); //m4l

 //second mass point
 w.factory("RooGaussian::gaussian3(obs1,250,25)"); //mll
 w.factory("RooGaussian::gaussian4(obs2,100,20)"); //m4l

 //make a 2D of the mass points
 w.factory("PROD::gaussian12(gaussian1,gaussian2)"); 
 w.factory("PROD::gaussian34(gaussian3,gaussian4)");

 TH1 *hh12_pdf = w.pdf("gaussian12")->createHistogram("obs1,obs2", 50, 50);
 hh12_pdf->SetLineColor(kBlue);
 
 TH1 *hh34_pdf = w.pdf("gaussian34")->createHistogram("obs1,obs2", 50, 50);
 hh34_pdf->SetLineColor(kBlue);

 w.factory("RooMomentMorphND::morphpdfND({alpha1[0,1]}, {alpha2[0,1]},{obs, obs1},{gaussian12,gaussian34},{0,1},{0,1})");
 w.var("alpha1")->setVal(0.5);
w.var("alpha2")->setVal(0.5);
 RooPlot* frame;
 frame = (w.var("obs1"),w.var("obs2"))->frame(); //this doesn't seem right
 w.pdf("gaussian12")->plotOn(frame, LineColor(kRed), LineStyle(kSolid));
 w.pdf("gaussian34")->plotOn(frame, LineColor(kBlue), LineStyle(kSolid));
 w.pdf("morphpdfND")->plotOn(frame, LineColor(kRed), LineStyle(kDashed));
 frame->Draw();

 hh12_pdf->Draw("surf4");
 hh34_pdf->Draw("surf4 same");

The code actually only works if I comment out these lines:

w.factory("RooMomentMorphND::morphpdfND({alpha1[0,1]}, {alpha2[0,1]},{obs, obs1},{gaussian12,gaussian34},{0,1},{0,1})");
w.pdf("morphpdfND")->plotOn(frame, LineColor(kRed), LineStyle(kDashed));
 w.var("alpha1")->setVal(0.5);
w.var("alpha2")->setVal(0.5);

That I need the most. Could you advise please how to fix this?

Best regards,

Diallo.

6

I did it actually, it should be constructed like this:

 TVectorD paramVec = TVectorD(2);
 paramVec[0]=0;
 paramVec[1]=1;
 /*paramVec[2]=2;
 paramVec[3]=3;
 paramVec[4]=4;
 paramVec[5]=5;
 paramVec[6]=6;*/
 RooRealVar alpha("alpha","alpha",0,1);
 RooAbsPdf* pdf1 = w.pdf("gaussian12");
 RooAbsPdf* pdf2 = w.pdf("gaussian34");
 //mu.setVal(25);
 RooArgList pdfs;
 pdfs.add(*pdf1);
 pdfs.add(*pdf2);
 RooArgList varlist;
 varlist.add(*w.var("obs1"));
 varlist.add(*w.var("obs2"));

 RooMomentMorph morphpdfND("morphpdfND", "test morph",alpha ,varlist, pdfs, paramVec,RooMomentMorph::Linear);
7

Yes, if you have issues building a RooFit class using the WooWorkspace::factory method you can always call the C++ constructor externally and if you need import later the pdf using w.import(pdf)

Lorenzo

8

Hi Lorenzo,

I am quite advanced in my project of processing a 2d interpolation using RooMomentMorphND in root. Unfortunately the documentation of this is very poor (or maybe it’s me who could’t find it), so I really start from scratch. So what I am doing is as following:
I have 33 h1 histograms and 33 h2 histograms which are 1D histograms. I fit them with a double sided crystall function and save the parameters fits in text file. So my code is reading the text file and make pdfs from the fit parameters. The observable in h1 is mll and the one in h2 is mS. I use PROD:: to construct the pairs (mll, mS) in a 2d plan, wich will give me 33 pairs of (mll, mS). So now I want to populate the 2D plane by interpolating using RooMomentMorphND. The results I got so far are promising as you can in see the file I attached:

Screen Shot 2021-02-04 at 9.02.03 PM
Screen Shot 2021-02-04 at 9.02.03 PM2880×1800 511 KB

The interpolated points are the ones plotted with the option cont1 which form the vertical line along the Y axis at mS ~ 175 GeV.
Here is my code:

#include "RooWorkspace.h"
#include <RooRealVar.h>
#include "RooDataHist.h"
#include "RooHistPdf.h"
#include "RooAbsArg.h"
#include "RooMomentMorphND.h"
#include "RooFit.h"
#include "TRandom3.h"
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooConstVar.h"
#include "RooPolynomial.h"
#include "RooIntegralMorph.h"
#include "RooNLLVar.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "RooPlot.h"
#include "TH1.h"
#include "RooMyPdf.cxx"
using namespace RooFit;



void morphingBackup()
{
  

  
 TString file_input[]={"./FitParameterM4l_2e2m.txt","./FitParameterAvgMll_2e2m.txt"};
   TString tmps;
   ifstream fin_mS(file_input[0]);
   ifstream fin_mZd(file_input[1]);
   int  Npoints = 33;
   double Param_mS[7][Npoints];
   double Param_mZd[7][Npoints];
   double PmS[33];
   double PmZd[33];
   fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;fin_mS>>tmps;//read first line.
   
   for (Int_t i=0; i<Npoints; i++)
     {
       
       
       fin_mS>>tmps; PmS[i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; PmZd[i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[0][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[1][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[2][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[3][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[4][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[5][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mS>>tmps; Param_mS[6][i]=tmps.IsFloat()?tmps.Atof():0;
       
     }
 
   fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;fin_mZd>>tmps;//read first line.
   for (Int_t i=0; i<Npoints; i++)
     {
       
       
       fin_mZd>>tmps; PmS[i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; PmZd[i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[0][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[1][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[2][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[3][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[4][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[5][i]=tmps.IsFloat()?tmps.Atof():0;
       fin_mZd>>tmps; Param_mZd[6][i]=tmps.IsFloat()?tmps.Atof():0;
       
     }

   /* for (int i =0; i< 33; i++)
     {
       cout << "mS = " << Param_mZd[6][i] << endl;
     }*/

   RooRealVar mll("mll","mll",0,300) ;
 


 RooRealVar mS("mS","mS", 120, 900) ;



 RooArgList pdfs;

  TVectorD paramVec = TVectorD(33);

RooWorkspace w("w", 2);
 RooAbsPdf *pdf;

 TH1 *hmZdmS[33];
 TH1 *hmZdmS_morph[10];

 //RooBinning dim1(1, 0.0, 1.0);
 //RooBinning dim2(1, 0.0, 1.0);
 //RooMomentMorphND::Grid referenceGrid(dim1, dim2);
 //referenceGrid.addPdf(*gauss1, 0, 0);

 for( int i=0; i<33; ++i)
   {
     w.factory(TString::Format("RooMyPdf::myPdf_mZd_%d(mll[0,300], %f, %f,%f, %f,%f,%f, %f)",i, Param_mZd[0][i],Param_mZd[1][i],Param_mZd[2][i],Param_mZd[3][i],Param_mZd[4][i],Param_mZd[5][i],Param_mZd[6][i]));
     w.factory(TString::Format("RooMyPdf::myPdf_mS_%d(mS[120,900], %f, %f,%f, %f,%f,%f, %f)",i, Param_mS[0][i],Param_mS[1][i],Param_mS[2][i],Param_mS[3][i],Param_mS[4][i],Param_mS[5][i],Param_mS[6][i]));
     w.factory(TString::Format("PROD::myPdf_mZdmS_%d(myPdf_mZd_%d,myPdf_mS_%d)", i,i,i));
     
     pdf = w.pdf(TString::Format("myPdf_mZdmS_%d",i));
     
     pdfs.add(*pdf);
     //referenceGrid.addPdf(*pdf, 0, 0);
     paramVec[i]=i;
     w.Print() ;

     hmZdmS[i] = pdf->createHistogram(TString::Format("hmZdmS_%d",i),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
   }
 pdfs.Print();
 RooPlot* mZdFrame = w.var("mll")->frame();
 RooPlot* mSFrame = w.var("mS")->frame();
 RooRealVar alpha("alpha","alpha",0,1);
 RooArgList varlist;
 varlist.add(mS);
 varlist.add(mll);
 RooArgList parlist;
 RooRealVar alpha1("alpha1","alpha1",0,1);
 //alpha1=*(w.var("alphaL"));
 //parlist.add(*w.var("alphaL"));
 TFile f("outputHist.root","recreate");
 
 alpha.setVal(0.1);
 RooMomentMorphND morphpdfND_0("morphpdfND_0","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[0] = morphpdfND_0.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[0]->Write();
 
 alpha.setVal(0.2);
 RooMomentMorphND morphpdfND_1("morphpdfND_1","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[1] = morphpdfND_1.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[1]->Write();
 
 alpha.setVal(0.3);
 RooMomentMorphND morphpdfND_2("morphpdfND_2","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[2] = morphpdfND_2.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[2]->Write();

 alpha.setVal(0.4);
 RooMomentMorphND morphpdfND_3("morphpdfND_3","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[3] = morphpdfND_3.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[3]->Write();

 alpha.setVal(0.5);
 RooMomentMorphND morphpdfND_4("morphpdfND_4","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[4] = morphpdfND_4.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[4]->Write();

 alpha.setVal(0.6);
 RooMomentMorphND morphpdfND_5("morphpdfND_5","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[5] = morphpdfND_5.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[5]->Write();

 alpha.setVal(0.7);
 RooMomentMorphND morphpdfND_6("morphpdfND_6","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[6] = morphpdfND_6.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[6]->Write();

 alpha.setVal(0.8);
 RooMomentMorphND morphpdfND_7("morphpdfND_7","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[7] = morphpdfND_7.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[7]->Write();

 alpha.setVal(0.9);
 RooMomentMorphND morphpdfND_8("morphpdfND_8","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);
 hmZdmS_morph[8] = morphpdfND_8.createHistogram(TString::Format("hmZdmS_morph_alpha%d",0),mS, Binning(460,120,1500), YVar(mll,Binning(100)));
 hmZdmS_morph[8]->Write();
 alpha.setVal(0.9);



 hmZdmS[0]->Draw("cont");
 for(int i =1; i<33; i++)
   {
     hmZdmS[i]->Draw("cont same");
     hmZdmS[i]->Write();
   }
  for(int i =0; i<8; i++)
   {
     hmZdmS_morph[i]->Draw("cont1 same");
   }
 //hmZdmS_morph[1]->Draw("cont1 same");
  f.Close();

}

I varied alpha from 0.1 to 0.9 hoping it will covers the plane but it’s not the case. It’s like alpha is varying in 1 direction along the Y axis. So my question is how should I change alpha to covers the plane? I think I should define a referenceGrid but I don’t know how should I change this line

 RooMomentMorphND morphpdfND_8("morphpdfND_8","2D morph",alpha ,varlist, pdfs, paramVec,RooMomentMorphND::Linear);

To include a referenceGrid. Could you advise please?

Tell me if you want to run my code I can give you the other files that are included.

Best regards,

Diallo.

9

Hi Diallo,

You are right, there is no documentation of the class and this is not good. We need to improve this. There is only the archive paper you have linked before.
Unfortunately I am not the author of the class and I am not expert on how using it.
I would need to have a look at the code and understand how to do what you are asking.
If you can post all the needed code and files to run, it would be good

Best

Lorenzo

10

Hi Lorenzo,

Ok I posted them.FitParameterAvgMll_2e2m.txt (1.2 KB) FitParameterM4l_2e2m.txt (1.3 KB) RooMyPdf.cxx (2.5 KB) RooMyPdf.h (1.3 KB)

To run, you just need to do: root -l morphingBackup.C
Unfortunately, it takes around 1h30 to run.

11

Hi Lorenzo,

Did you get a chance to run the code?

Best,

Diallo.

12

Hi Diallo,

Apologies, not yet. I have been scared by the long running time. I see however you did not attach the macro, morphingBackup.C. Is it the same included in the text above ?

Lorenzo

13

Hi @moneta indeed the macro was there. Not attached but in the text of a post. Here it is:
morphingBackup.C (7.0 KB)

14

Hi,
The code can be run, but too slow. RooFit is doing an integral of an integral, something one should avoid. Since you have an analytical expression for the PDF, I would first add an analytical integral.

Lorenzo

15

Hi Lorenzo and Couet,

Thanks for your reply. Actually I get in touch with Stefan Gadatsch one of the authors who wrote the paper I attached before. As I was suspected 2 morph parameters have to be defined as in in this example he gave me:

#include "RooAbsArg.h"
#include "RooConstVar.h"
#include "RooDataHist.h"
#include "RooDataSet.h"
#include "RooFit.h"
#include "RooGaussian.h"
#include "RooHistPdf.h"
#include "RooIntegralMorph.h"
#include "RooMomentMorphND.h"
#include "RooMyPdf.cxx"
#include "RooNLLVar.h"
#include "RooPlot.h"
#include "RooPolynomial.h"
#include "RooRealVar.h"
#include "RooWorkspace.h"
#include "TAxis.h"
#include "TCanvas.h"
#include "TH1.h"
#include "TRandom3.h"
#include <RooRealVar.h>

using namespace RooFit;



void morphing() {

  TString file_input[] = {"./FitParameterM4l_2e2m.txt",
                          "./FitParameterAvgMll_2e2m.txt"};
  TString tmps;
  ifstream fin_mS(file_input[0]);
  ifstream fin_mZd(file_input[1]);
  int Npoints = 33;
  double Param_mS[7][Npoints];
  double Param_mZd[7][Npoints];
  double PmS[33];
  double PmZd[33];
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps;
  fin_mS >> tmps; // read first line.

  for (Int_t i = 0; i < Npoints; i++) {

    fin_mS >> tmps;
    PmS[i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    PmZd[i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[0][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[1][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[2][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[3][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[4][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[5][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mS >> tmps;
    Param_mS[6][i] = tmps.IsFloat() ? tmps.Atof() : 0;
  }

  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps;
  fin_mZd >> tmps; // read first line.
  for (Int_t i = 0; i < Npoints; i++) {

    fin_mZd >> tmps;
    PmS[i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    PmZd[i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[0][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[1][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[2][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[3][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[4][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[5][i] = tmps.IsFloat() ? tmps.Atof() : 0;
    fin_mZd >> tmps;
    Param_mZd[6][i] = tmps.IsFloat() ? tmps.Atof() : 0;
  }

  RooRealVar mll("mll", "mll", 0, 300);
  RooRealVar mS("mS", "mS", 120, 900);


  RooArgList pdfs;

  TVectorD paramVec = TVectorD(33);

  RooWorkspace w("w", 2);
  RooAbsPdf *pdf;

  TH1 *hmZdmS[33];
  TH1 *hmZdmS_morph[10];

  RooBinning dim1(1, 0.0, 1.0);
  RooBinning dim2(1, 0.0, 1.0);
  RooMomentMorphND::Grid referenceGrid(dim1, dim2);

  for (int i = 0; i < 33; ++i) {
    w.factory(TString::Format(
        "RooMyPdf::myPdf_mZd_%d(mll[0,300], %f, %f,%f, %f,%f,%f, %f)", i,
        Param_mZd[0][i], Param_mZd[1][i], Param_mZd[2][i], Param_mZd[3][i],
        Param_mZd[4][i], Param_mZd[5][i], Param_mZd[6][i]));
    w.factory(TString::Format(
        "RooMyPdf::myPdf_mS_%d(mS[120,900], %f, %f,%f, %f,%f,%f, %f)", i,
        Param_mS[0][i], Param_mS[1][i], Param_mS[2][i], Param_mS[3][i],
        Param_mS[4][i], Param_mS[5][i], Param_mS[6][i]));
    w.factory(TString::Format("PROD::myPdf_mZdmS_%d(myPdf_mZd_%d,myPdf_mS_%d)",
                              i, i, i));

    pdf = w.pdf(TString::Format("myPdf_mZdmS_%d", i));

    pdfs.add(*pdf);
    paramVec[i] = i;
    w.Print();

    hmZdmS[i] =
        pdf->createHistogram(TString::Format("hmZdmS_%d", i), mS,
                             Binning(460, 120, 1500), YVar(mll, Binning(100)));
  }
  
  pdfs.Print();
  
  RooPlot *mZdFrame = w.var("mll")->frame();
  RooPlot *mSFrame = w.var("mS")->frame();
  
  RooRealVar alpha1("alpha1", "alpha1", 0.5, 0, 1);
  RooRealVar alpha2("alpha2", "alpha2", 0.5, 0, 1);

  RooArgList varlist;
  varlist.add(mS);
  varlist.add(mll);
  referenceGrid.addPdf(*w.pdf("myPdf_mZdmS_0"), 0, 0);
  referenceGrid.addPdf(*w.pdf("myPdf_mZdmS_1"), 1, 0);
  referenceGrid.addPdf(*w.pdf("myPdf_mZdmS_2"), 0, 1);
  referenceGrid.addPdf(*w.pdf("myPdf_mZdmS_3"), 1, 1);
  TFile f("outputHist.root", "recreate");

  RooMomentMorphND morphpdfND("morphpdfND", "2D morph",
                              RooArgList(alpha1, alpha2), RooArgList(mS, mll),
                              referenceGrid, RooMomentMorphND::Linear);

  alpha1.setVal(0.0);
  alpha2.setVal(0.0);
  hmZdmS_morph[0] = morphpdfND.createHistogram(
      TString::Format("hmZdmS_morph_alpha%d", 0), mS, Binning(460, 120, 1500),
      YVar(mll, Binning(100)));
  hmZdmS_morph[0]->Write();

  alpha1.setVal(1.0);
  alpha2.setVal(1.0);
  hmZdmS_morph[1] = morphpdfND.createHistogram(
      TString::Format("hmZdmS_morph_alpha%d", 0), mS, Binning(460, 120, 1500),
      YVar(mll, Binning(100)));
  hmZdmS_morph[1]->Write();

  hmZdmS[0]->Draw("cont");
  for (int i = 1; i < 33; i++) {
    hmZdmS[i]->Draw("cont same");
    hmZdmS[i]->Write();
  }
  for (int i = 0; i < 2; i++) {
    hmZdmS_morph[i]->Draw("cont1 same");
  }
  f.Close();
}

It’s already a good start and I am trying to develop it further. I agree with you Lorenzo that the code is very slow I guess by adding an analytical integral you mean you calculating the integral of the Double Sided crystal ball from it’s function, right?
For that I calculate it as coded it like this:

double DoubleSidedCrystalballIntegral(double *x, double *par)
 { 
   double alpha_l = par[0]; 
   double alpha_h = par[1]; 
   double n_l     = par[2]; 
   double n_h     = par[3]; 
   double mean	= par[4]; 
   double sigma	= par[5];
   double N	= par[6];
   double t = (x[0]-mean)/sigma;
   double Int;
   double ratio_alphaN_l = n_l/alpha_l;
   double ratio_alphaN_h = n_h/alpha_h;
   
   double firstNum_l = (alpha_l*(ratio_alphaN_l-alpha_l-t))/n_l;
   double secondNum_l = alpha_l*(mean-x[0]) -sigma*alpha_l*alpha_l + n_l*sigma;
   //double secondNum_l = -alpha_l*alpha_l*sigma + alpha_l*(mean-x[0]) + n_h*sigma;
   
   double firstNum_h = (alpha_h*(ratio_alphaN_h-alpha_h+t))/n_h;
   double secondNum_h = alpha_h*(mean-x[0]) +sigma*alpha_h*alpha_h - n_h*sigma;
   //double secondNum_h = -alpha_h*alpha_h*sigma + alpha_h*(x[0]-mean) - n_h*sigma;
   
   double Denom_l=(n_l-1)*alpha_l;
   double Denom_h=(n_h-1)*alpha_h;
   
   if (-alpha_l <= t && alpha_h >= t)
     {
       Int = sqrt(M_PI*2)*sigma;
       //Int = sqrt(M_PI/2)*sigma*erf(-t/sqrt(2));
     }
   else if (t < -alpha_l)
     {
       
       Int = (exp((-alpha_l*alpha_l/2))*pow(firstNum_l, -n_l)*secondNum_l)/Denom_l;

     }
   else if (t > alpha_h)
     {
       Int = (exp((-alpha_h*alpha_h/2))*pow(firstNum_h, -n_h)*secondNum_h)/Denom_h;
       
     }
   return N*Int;
 }

The DSCB function can be found in RooMyPdf.cxx file that I coded based on this paper, page 10: https://dash.harvard.edu/bitstream/handle/1/29362185/1606.03833.pdf?sequence=1
I get the analytical integral from https://www.wolframalpha.com/ where you can enter some functions as input and get its integral as output.
But I didn’t get how could I add it to get the code faster.

Best regards,

Diallo.

16

Hi,

Good you contact the author and found a solution. To integrate the analytical integration in the PDF you have to implement 2 functions :

  • RooMyPdf::getAnalyticalntegral that advertises the computation of the analytical integral
  • RooMyPdf::analyticalIntegral that performs the actual computation.

You can see as example what is done for some of the pdf’s like the RooGaussian.
See ROOT: RooGaussian Class Reference and ROOT: RooGaussian Class Reference.

Lorenzo

17

Hi,

I recreated a skeleton of my PDF using RooClassFactory via the command:

RooClassFactory::makePdf("RooMyPdf", "x,alphaL,alphaH,nL,nH,mean,sigma,N", "", "", kTRUE, kFALSE, "");

This created skeletons of the 2 functions needed for the integral inside the .cxx file. I implemented my pdf function (which is a Double sided crystall ball) and its analytical integral. I attached the file containing the DSCB function and its integral. DSCB .pdf (223.5 KB)
You can see that one term of the DSCB function is a Gaussian, another term is a simple cirstall ball function and their analytical integrals exist already in root that can be found for instance in these links respectively: ROOT: roofit/roofit/src/RooGaussian.cxx Source File and here:
ROOT: roofit/roofit/src/RooSDSCBShape.cxx Source File. Since the DSCB doesn’t exist in root I was trying to implement it following these Gaussian and the CB. But this is not straightforward for me and what I did so far doesn’t seem to converge.
I am attaching again my files: morphing.C (5.2 KB) RooMyPdf.cxx (4.8 KB) RooMyPdf.h (1.5 KB) FitParameterAvgMll_2e2m.txt (1.2 KB) FitParameterM4l_2e2m.txt (1.3 KB)

Best

18

Hi,

Do you have any idea about this issue? I am stuck here, so your help would be really appreciated.

Best

19

Hi,
The problem is in the integral calculation. If I just call RooMyPdf.getVal(x), i.e. the normalised pdf value, I get a NaN. We need to debug your integral code

Best,

Lorenzo

20

Hi,

I am pretty sure the integration calculation is correct, but I am afraid the implementation should be more complex that what I did. I think it should be similar to the CB one already available in root, but I found it a bit complex to adapt it to the DSCB case.

Best regards,

Diallo.