Histogram Functors for RDataFrame

RENATO_QUAGLIANI · May 20, 2020, 7:26am

Dear experts,
I am trying to understand why this code doesn’t compile properly. I am failing to understand why this is not working.

Basically i wrote a functor for Reading in a histogram given a input value and return a value , such that i can run MT the weight adding on my tuples using RDataFrame.

// #include <TH1.h>
#include <TH1D.h>
// #include <TH1D.h>
class TH1DHistoAdder {
public:
    TH1DHistoAdder( const TH1D& h , TString _finalBranchName, bool _interpolate ) : m_histo(h) , m_finalBranchName(_finalBranchName), m_interpolate(_interpolate) {}
    TH1DHistoAdder( const TH1DHistoAdder&& other) : m_histo(other.m_histo), m_finalBranchName(other.m_finalBranchName), m_interpolate(other.m_interpolate){}
    TH1DHistoAdder( const TH1DHistoAdder& other)  : m_histo(other.m_histo), m_finalBranchName(other.m_finalBranchName), m_interpolate(other.m_interpolate){}
    double operator()(double _variable) const{ 
      double _val = 1;
      double _var = _variable;
      int _bin = m_histo.FindFixBin(_var);
      if (m_histo.IsBinUnderflow(_bin) || m_histo.IsBinOverflow(_bin)) {
          if (_var <= m_histo.GetXaxis()->GetXmin()) _var = m_histo.GetXaxis()->GetXmin() + m_histo.GetXaxis()->GetBinWidth(1) / 100.;
          if (_var >= m_histo.GetXaxis()->GetXmax()) _var = m_histo.GetXaxis()->GetXmax() - m_histo.GetXaxis()->GetBinWidth(m_histo.GetNbinsX()) / 100.;
          _bin = m_histo.FindFixBin(_var);
      }
      if( !m_interpolate){
          _val = m_histo.GetBinContent(_bin);
      }else{
          _val = m_histo.Interpolate( _var);
      }
      return _val;
    }
    const char* branchName(){
        return m_finalBranchName.Data();
    }
private:
    TH1D    m_histo;
    TString m_finalBranchName;
    bool m_interpolate;
};

This is my class, basically given a input histogram i want to attach a weight in “interpolation” mode or bare mode.

All compiles without Interpolate call, When i try to compile this code i get:

error: no matching member function for call to 'Interpolate'
          _val = m_histo.Interpolate( _var);
                 ~~~~~~~~^~~~~~~~~~~
/cvmfs/sft.cern.ch/lcg/releases/ROOT/6.18.04-c767d/x86_64-centos7-gcc8-opt/include/TH1.h:325:21: note: candidate function not viable: no known conversion from 'const TH1D' to 'TH1' for object argument
   virtual Double_t Interpolate(Double_t x);
                    ^
/cvmfs/sft.cern.ch/lcg/releases/ROOT/6.18.04-c767d/x86_64-centos7-gcc8-opt/include/TH1.h:326:21: note: candidate function not viable: requires 2 arguments, but 1 was provided
   virtual Double_t Interpolate(Double_t x, Double_t y);
                    ^
/cvmfs/sft.cern.ch/lcg/releases/ROOT/6.18.04-c767d/x86_64-centos7-gcc8-opt/include/TH1.h:327:21: note: candidate function not viable: requires 3 arguments, but 1 was provided

Thanks for any help.

Renato

Please read tips for efficient and successful posting and posting code

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

RENATO_QUAGLIANI · May 20, 2020, 7:33am

If i try to change the TH1D to TH1. Things improves but maybe my ROOT versions is not yet including the constification of methods?

/cvmfs/sft.cern.ch/lcg/releases/ROOT/6.18.04-c767d/x86_64-centos7-gcc8-opt/include/TH1.h:325:21: note: candidate function not viable: 'this' argument has type 'const TH1', but method is not marked const
   virtual Double_t Interpolate(Double_t x);

Even if i see from the doxygen that Interpolate is marked const

eguiraud · May 20, 2020, 7:41am

Hi Renato,
the problem is indeed that you are calling a non-const method on a const object. If I understand correctly you are on ROOT v6.18: in that case, the documentation for Interpolate is at https://root.cern.ch/doc/v618/classTH1.html#ad3880f90543c7e60bf3180dc1a196bcc

Cheers,
Enrico

RENATO_QUAGLIANI · May 20, 2020, 7:42am

Hi @eguiraud, is there any way out i can call Interpolate inside the operator() const method ?

I was thinking about making an extension of TH1D which re-implement Interpolate as const . Do you know a safe way to “extend” an existing root class, such as TH1D just adding a new method?
I tried something like

class TH1DExtension : public TH1D {
    public :
    Double_t Interpolate(Double_t x) const{
        if (fBuffer) ((TH1*)this)->BufferEmpty();
        Int_t xbin = fXaxis.FindFixBin(x);
        Double_t x0,x1,y0,y1;
    
        if(x<=GetBinCenter(1)) {
        return RetrieveBinContent(1);
        } else if(x>=GetBinCenter(GetNbinsX())) {
        return RetrieveBinContent(GetNbinsX());
        } else {
        if(x<=GetBinCenter(xbin)) {
            y0 = RetrieveBinContent(xbin-1);
            x0 = GetBinCenter(xbin-1);
            y1 = RetrieveBinContent(xbin);
            x1 = GetBinCenter(xbin);
        } else {
            y0 = RetrieveBinContent(xbin);
            x0 = GetBinCenter(xbin);
            y1 = RetrieveBinContent(xbin+1);
            x1 = GetBinCenter(xbin+1);
        }
        return y0 + (x-x0)*((y1-y0)/(x1-x0));
        }
    }
};

but didn’t work…

RENATO_QUAGLIANI · May 20, 2020, 7:46am

I can also try to make a custom Interpolate1D( const TH1D & histo, double value)
reimplementing this :

Double_t TH1::Interpolate(Double_t x)
 {
    if (fBuffer) ((TH1*)this)->BufferEmpty();
 
    Int_t xbin = FindBin(x);
    Double_t x0,x1,y0,y1;
 
    if(x<=GetBinCenter(1)) {
       return RetrieveBinContent(1);
    } else if(x>=GetBinCenter(GetNbinsX())) {
       return RetrieveBinContent(GetNbinsX());
    } else {
       if(x<=GetBinCenter(xbin)) {
          y0 = RetrieveBinContent(xbin-1);
          x0 = GetBinCenter(xbin-1);
          y1 = RetrieveBinContent(xbin);
          x1 = GetBinCenter(xbin);
       } else {
          y0 = RetrieveBinContent(xbin);
          x0 = GetBinCenter(xbin);
          y1 = RetrieveBinContent(xbin+1);
          x1 = GetBinCenter(xbin+1);
       }
       return y0 + (x-x0)*((y1-y0)/(x1-x0));
    }
 }

I think all the methods here will work fine, but i am not 100% sure about the fBuffer and BufferEmpty() calls how i can access those information since this is using a protected class member.
Unfortunately change ROOT version is not a solution for me.

eguiraud · May 20, 2020, 7:55am

Hi Renato,
about post number 4: you can’t change the constness of a method by overriding it: https://godbolt.org/z/Z8ZjVv

about post number 5: you can access fBuffer via TH1::GetBuffer(). About BufferEmpty(), I’m not sure.

The simplest solution would be to remove the const qualifier from your operator() since, in v6.18 at least, it’s calling non-const methods.

Cheers,
Enrico

RENATO_QUAGLIANI · May 20, 2020, 7:59am

Ok, i remove the const of the operator, nevertheless, i am not sure if doing this will prevent MT to work properly. I will ensure that when this is used i run single thread and i think i will be fine.
Thanks
Renato

eguiraud · May 20, 2020, 8:19am

Yes, if you call that method from multiple threads concurrently, you are in trouble.
I don’t know how TH1DHistoAdder is used in your application, but in RDataFrame typically we duplicate objects we work on, so that each thread has its own copy. That pattern would probably work for you too.

RENATO_QUAGLIANI · May 20, 2020, 8:21am

well, the approach i am following is to just make a

vector<TH1DHistoAdder>

housing all the histograms i need to read to attach a branch, then i call a Define in loop for that and finally a snapshot with all weights definitions.
I am moving to this approach to enable in my analysis workflow 2 properties :

do the weight-definitions on the fly and work on temporary tuples
do the weight-definitions and persistency.

RENATO_QUAGLIANI · May 20, 2020, 8:21am

So maybe for the persistency i am fine running single-threaded, for the working on the fly apporach i may will need to upgrade ROOT at some point i believe.

eguiraud · May 20, 2020, 8:38am

DefineSlot might be useful to make sure you don’t get concurrent calls to a method of the same non-thread-safe object, by providing a different instance of the object to each thread. That might be more memory-consuming but better-performing than just putting a lock to protect the thread-unsafe calls.

Cheers,
Enrico

RENATO_QUAGLIANI · May 20, 2020, 9:22am

Thanks, i think for the moment i will just ensure i run with DisableImplicitMT() the functions adding the weights and persisting the tuples to disk. Once i am fine, i will probably just upgrade ROOT version.

eguiraud · May 20, 2020, 9:32am

Cool! You might want to wait for ROOT v6.22 to update the ROOT version – it’s coming out soon, and it will contain the best RDataFrame so far

system · June 3, 2020, 9:32am

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.