[RooFit] How to do Mellin Convoluton

Hello,

Sorry to ask this question again (And I think it is better to move here).
So, I followed Wouter’s instructions in
http://root.cern.ch/phpBB2/viewtopic.php?t=7539&highlight=

To be simple, I just used a single Gaussian for the K factor, but it doesn’t work. I expect that after doing the Mellin convolution, the k is ntegrated out, and the pdf (versus t) should be something like exponential convoluted with a Gaussian. But what I obtained is a flat line. Did I do something wrong?

[code]void tesMellinCov(){
using namespace RooFit;

RooRealVar *t = new RooRealVar( “t”, “proper time”, -0.5 , 10.0 , “ps” );
RooRealVar *st = new RooRealVar(“st”, “proper time resolution”, 0.05, 0.10, “ps” );
RooRealVar *k = new RooRealVar( “k”, “k factor”, 0.3 , 3.0 );

RooRealVar *Tau = new RooRealVar(“Tau”, “lifetime”, 1.5, 0.1, 3.0, “ps”);

RooRealVar *TauResScale = new RooRealVar(“TauResScale”,“resolutions scale”,1.185, 0.5, 3);
RooRealVar *TauResBias = new RooRealVar(“TauResBias”,“bias”,0.,-5,5);
RooRealVar *TauResBiasScale = new RooRealVar(“TauResBiasScale”,“bias scale”,1.0,0.5,3);
RooResolutionModel *TauRes = new RooGaussModel(“TauRes”,“resolution model”,
*t,*TauResBias,*st,*TauResBiasScale,*TauResScale);
// Exp(t|tau) fConv Gau(st) — fConv: Fourier convolution
RooAbsPdf *TauPdf = new RooDecay(“BcTauPdf”, “proper time pdf”, *t, *Tau, *TauRes,RooDecay::SingleSided);

TauResBias->setConstant(kTRUE);
TauResBiasScale->setConstant(kTRUE);

// Take into account K factor
RooFormulaVar PseudoTau = new RooFormulaVar(“PseudoTau”, “pseudo lifetime”, "kTau", RooArgSet(*k,*Tau));
RooAbsPdf *PseudoTauPdf = new RooDecay(“PseudoTauPdf”, “proper time pdf”, *t, *PseudoTau, *TauRes, RooDecay::SingleSided);

// K factor distribution, one Gaussian
RooRealVar *gauMean = new RooRealVar( “gauMean” , “mean” , 1.0, 0.8, 1.2 );
RooRealVar *gauSigma = new RooRealVar( “gauSigma”, “sigma”, 0.2, 0.1, 0.3 );
RooAbsPdf *kPdf = new RooGaussian(“kPdf”, “k distribution”, *k, *gauMean, *gauSigma) ;

// [ Exp(t|k*tau) fConv Gau(st)] mConv Gau(k) — mConv: Mellin convolution, k integrated out
k->setBins(10000,“cache”) ;
RooAbsPdf *PseudoPdf = new RooFFTConvPdf(“PseudoPdf”, “PseudoPdf”, *k, *PseudoTauPdf, *kPdf );

RooPlot *tframe = t->frame();
PseudoPdf->plotOn(tframe);
tframe->Draw();
}[/code]

Thank you very much.

Cheers, Jibo

Hi,

I can reproduce your problem. I’m looking into it now.

Wouter

Hi Jibo,

Sorry for the long wait, this was a rather complicated problem to debug, but I’ve found the problem now:

Inside RooFFTConvPdf a bug was introduced recently that causes improper handling of internal multi-dimensional cache histograms (i.e. when you cache F(x,y) (X) G(x) you have a 2D cache for all values of x and y). Due to this bug the internal copies of F and G were no longer connected to the instance of ‘y’ inside the cache histogram so cache histogram was filled with the same distribution in each ‘slice’ of y. A plot vs y (your ‘t’) will then result in a flat line. With my fix I now get what seems to be a sensible looking plot from your test macro.

I will commit the fix later today to SVN so that it will be included in ROOT 5.22 (due next week).

Wouter

Dear Wouter,

Thank you so much for your kind helps.

Now the above macro can give the plot correctly. But when I tried to use the PseudoPdf to fit some data. It used more than 80% memory of my laptop (4GB) or our cluster (16GB) at once. Then it eventually used up all the memory. I guess there is still something wrong. Would you be so kind to have a look at it. I hope that I didn’t misunderstand something.

And for the distribution of the k factor, I also tried one Gaussian. I saw the same problem. And it seems I can not use RooKeysPdf here (maybe RooKeysPdf is too complicated?).

I’m using the head of svn (today). Hereafter is the marco. The required data files are in the attached tarball.

Thanks a lot.

Cheers, Jibo

void tesMellinCov(){
  using namespace RooFit;
 
  RooRealVar *t  = new RooRealVar( "t", "proper time", -1.0 , 5.0 , "ps" );
  RooRealVar *st = new RooRealVar("st", "proper time resolution",  0.,  0.08, "ps" );
  RooRealVar *k  = new RooRealVar( "k", "k factor",  0.3 ,  1.8  );

  // observables 
  RooArgSet* observables = new RooArgSet(*t,*st);
  
  // read in data
  RooDataSet *tData = RooDataSet::read("lifetimeData.txt", RooArgSet( *observables ) );
  RooDataSet *kData = RooDataSet::read("kFactorData.txt" , RooArgSet( *k ) );

  // check the reading in
  TCanvas *MyCan = new TCanvas();

  RooPlot *tframe = t->frame();
  tData->plotOn(tframe);
  tframe->Draw();
  MyCan->Print("tData.eps"); 

  RooPlot *kframe = k->frame();
  kData->plotOn(kframe);
  kframe->Draw();
  MyCan->Print("kData.eps");

  //
  RooRealVar *Tau = new RooRealVar("Tau", "lifetime", 1.5, 0.1, 3.0, "ps");
   
  RooRealVar *TauResScale     = new RooRealVar("TauResScale","resolutions scale",1.185, 0.5, 3);
  RooRealVar *TauResBias      = new RooRealVar("TauResBias","bias",0.,-5,5);
  RooRealVar *TauResBiasScale = new RooRealVar("TauResBiasScale","bias scale",1.0,0.5,3);
  RooResolutionModel *TauRes = new RooGaussModel("TauRes","resolution model",
                     *t,*TauResBias,*st,*TauResBiasScale,*TauResScale);
  // Exp(t|tau) fConv Gau(st)   --- fConv: Fourier convolution
  RooAbsPdf  *TauPdf = new RooDecay("BcTauPdf", "proper time pdf", *t, *Tau, *TauRes,RooDecay::SingleSided);

  TauResBias->setConstant(kTRUE);
  TauResBiasScale->setConstant(kTRUE);

  // Take into account K factor
  RooFormulaVar *PseudoTau = new RooFormulaVar("PseudoTau", "pseudo lifetime", "k*Tau", RooArgSet(*k,*Tau)); 
  RooAbsPdf  *PseudoTauPdf = new RooDecay("PseudoTauPdf", "proper time pdf", *t, *PseudoTau, *TauRes, RooDecay::SingleSided);

  // K factor distribution, one Gaussian
  // RooRealVar *gauMean  = new RooRealVar( "gauMean" , "mean" , 1.0, 0.8, 1.2 );
  // RooRealVar *gauSigma = new RooRealVar( "gauSigma", "sigma", 0.2, 0.1, 0.3 );
  // RooAbsPdf *kPdf = new RooGaussian("kPdf", "k distribution", *k, *gauMean, *gauSigma) ;

  RooDataHist *kHist = new RooDataHist("kHist","kHist",*k,*kData);
  RooHistPdf  *kPdf  = new RooHistPdf("kPdf","kPdf",*k,*kHist);

  // RooKeysPdf doesn't work
  //RooAbsPdf *kPdf = new RooKeysPdf( "kPdf", "kPdf", *k, *kdata ); 

  // [ Exp(t|k*tau) fConv Gau(st)] mConv Gau(k)  --- mConv: Mellin convolution, k integrated out
  k->setBins(10000,"cache") ;
  RooAbsPdf *PseudoPdf = new RooFFTConvPdf("PseudoPdf", "PseudoPdf", *k, *PseudoTauPdf, *kPdf );

  RooPlot *tframe = t->frame();
  PseudoPdf->plotOn(tframe);
  tframe->Draw();

  PseudoPdf->fitTo(*tData);
}

tesMellinCov.tar.gz (10.1 KB)

Hi Jibo,

Thanks for the update. It looks like there still is a memory leak somewhere.
I will try out your example macro and debug it.

Wouter

Hi, Wouter,

I’m just wondering whether you have got time to look at this problem. Just now I tried with the head of ROOT, RooFit v2.98, the memory leak is still there.

BTW: There is no difference for RooFit between running
svn update
and running
svn co root.cern.ch/svn/root/branches/dev/roostats,
right? You know, I haven’t seen the updates of RooFit until today (it was updated to v2.98 today). And I saw that you suggested to use the latter of the above to check out the dev/roostats branch of ROOT. I doubted that “svn update” doesn’t always give me the head of RooFit. :slight_smile:

Thanks a lot.

Jibo

Hi Jibo,

I think I fixed all the problems (my simplified example now runs fine).

Concerning the updates: you only pick up updates of the branch that
you originally checked out. I’ve been regularly committed fixes to
the dev/roostats branch as I work. Since the ROOT team is about
to build a new release, I have also updated the SVN trunk with all
changes that were made in the dev/roostats branch in the past month.

Let me know if you have more problems. Then I’ll look into it again.

Wouter

Hi, Wouter,

Thanks a lot for your reply.

Would you be so kind to do a test with the tar ball I posted before (tesMellinCov.tar.gz)? You can find it just in this thread. You know, I tired it again using the dev/roostats branch on my laptop and on our cluster. On my laptop, it consumed about 90% memory of the total 4GB for a long time. Then I switched to the cluster (with 16 GB memory). It used about 45% memory for some time. Then the memory usage was increasing. I found that it was keeping complaining that:

[quote][#1] INFO:Caching – RooAbsCachedPdf::getCache(PseudoPdf) creating new cache 0xd4b3e0 with pdf PseudoTauPdf_CONV_kPdf_CACHE_Obs[t,st,k] for nset (t,st,blindState) with code 0
[#0] ERROR:InputArguments – RooDataHist::weight(PseudoPdf_PseudoTauPdf_CONV_kPdf_CACHEHIST_Obs[t,st,k]_BufFrac0.1_BufStrat0) interpolation in 3 dimensions not yet implemented
[#0] ERROR:InputArguments – RooDataHist::weight(PseudoPdf_PseudoTauPdf_CONV_kPdf_CACHEHIST_Obs[t,st,k]_BufFrac0.1_BufStrat0) interpolation in 3 dimensions not yet implemented
[#1] INFO:NumericIntegration – RooRealIntegral::init(PseudoTauPdf_CONV_kPdf_CACHE_Obs[t,st,k]_Int[st,t]) using numeric integrator RooAdaptiveIntegratorND to calculate Int(t,st)[/quote]

Then I killed it before the memory was used up. Did I miss something developed recently? And even before the failure, such a large memory usage is a problem, you know, I only included 10% of the total data in that tar ball.

And for the svn issue, I am a little confused. I thought that the SVN trunk would always include all the updates of all the branches, as it is said that it is equivalent to the cvs head. But it seems it is not. So for RooFit, I should also check out the dev/roostats branch frequently so that we could benefit from the changes you have made? But in this case we would miss the changes made in the other branches and we need to have two version of ROOT (trunk and dev/roostats branch)? Did I misunderstand something?

Thanks again.

Jibo

Hi,

I will run your test again.

The new RooFit will surely fix your memory consumption problem (which is due to an imprudent choice of default caching strategy in earlier versions of RooFFTConvPdf that use large amount of memory with the p.d.f has many dimensions), but I would like to make sure there are no other problems (this type of convolutions is not in my automated regression test suite yet).

Concerning the SVN branches: I think there is some conceptual confusion here. The point of branches is that changes in them are not automatically swept into the trunk (or into each other). The dev/roostats branch allows me to commit development code without disturbing other ROOT developers and vice versa. When the time for a release is close I make sure that the code in the branch works also well in the trunk and only then I synchronize the trunk. NB: If you have a dev/roostats branch checked out, it is sufficient to to do an ‘svn up’ to pick up the latest changes.

[ NB: A new root release 5.23.04 is out meanwhile with the latest RooFit code ]

Wouter

Hi, Wouter,

I’m just wondering whether you have got time to look at this. Is is possible for you to include the fix into 5.24 too. :smiley:

Thanks a lot.

Cheers, Jibo

Hi Jibo,

The type of FFT convolution you are doing now (convolution in a non-observable) works fine in RooFit, so 5.24 should work for you. (It did not in the last 5.23 series release, but problems in FFTconv have been fixed meanwhile)

Your particular example is numerically rather complicated and may require some tuning to get it to work properly, due to the relatively high dimensionality of the problem (t,st,k)
but there are no fundamental issues as far as I can tell now. E.g. if I drop ‘st’ as an observable for a moment, your example runs relatively quickly.

Wouter