Self-Filling Histograms

By
Benno List
and
Jenny List, University of Hamburg

How Do I Get the Code?

The code is located in a cvs repository at DESY afs. To check it out, do

cvs -d /afs/desy.de/user/b/blist/.cvs co jbltools/sfh

What are Self-Filling Histograms?

What we call "Self-Filling Histograms" are histograms, derived from root histograms, where one can define during booking not only name, title, and binning of a histogram, but also what is plotted and under which conditions.

Consider this example:

SFH1F *h1 = 
  new SFH1F ("JPsiElasticTT", 
             "#muODS: m_{J/#Psi}: elastic production, (#mu-#mu)", 
             40, 2., 4.,
             *this,
             jpsiMass, 
             jpsiElastic&&jpsiTrackTrack&&jpsiNoCosmic&&jpsi2mu);

Without going into the details yet it is certainly plausible that this should be a histogram with name "JPsiElasticTT" and title "ODS: : elastic production, ", with 40 bins between 2 and 4 GeV, which shows the mass of meson candidates, after cuts have been applied that demand:

• The J/ has been produced elastically.
• The J/ has 2 tracks.
• It's not a cosmic.
• It decays to 2 muons.

The classes contained in the package SFH provide the means to write this type of code.

At the core of our toolkit are, obviously, the self-filling histogram classes, which are called SFH1F, SFH2F, etc, similarly to the corresponding root classes TH1F, TH2F. etc. from which they are derived.

What do we mean by self-filling?

The main point is that a self-filling histogram implements a method SFH1F::Fill(), which in contrast to the method of Int_t TH1F::Fill(Axis_t x) does not take an argument.

When we loop over the entries of a root Tree, the Fill() method has to be called once per entry. If we have a list of all self-filling histograms, which is provided by class HList, this is easily done in a single call, which in turn calls all Fill() methods of the histograms.

How does the histogram know what to fill?

The secret lies in the objects jpsiMass, jpsiElastic etc, which are declared as follows (at the moment we assume there is only one, and exactly one, J/psi candidate per Tree entry):

FloatFun& jpsiMass       = *new JpsiMassFun;

BaseCut&  jpsiElastic    = *new JpsiElasticCut;
BaseCut&  jpsiTrackTrack = *new JpsiTTCut;
BaseCut&  jpsiNoCosmic   = *new JpsiNoCosmicCut;
BaseCut&  jpsi2mu        = *new Jpsi2MuCut;

So, jpsiMass is an object of type JpsiMassFun.

Assuming we have a ntuple declared in an include file NTuple.h, automatically generated by root, and this ntuple has a member jpsimass, we could have this code:

#include "Ntuple.h"
#include "jbltools/sfh/FloatFun.h"

external Ntuple* ntuple; // global pointer to ntuple

class JpsiMassFun : public FloatFun {
  public: 
    virtual float operator() () const {
      return nt->jpsimass;
    }
};

So, JpsiMassFun is a class that implements operator(), which, when called, returns the jpsi mass stored in the ntuple row that is currently in memory.

This class is so light-weight that it can be completely implemented within a header file, without need for a separate implementation file.

Consequently, the base class FloatFun is an abstract base class which does not much more than declare that any FloatFun subclass must implement operator():

class FloatFun {
  public:
    virtual float operator() () const = 0; 
  protected:  
    virtual ~FloatFun() {}; 
};

(Why the destructor is protected is an issue that will be covered later.)

How does the histogram know how to cut?

By now, it should be almost obvious how the cuts are implemented: They are derived from an abstract base class BaseCut:

class BaseCut {
  public:
    virtual bool operator() () const = 0; 
  protected:
    virtual ~BaseCut() {}; 
};

The only difference between FloatFun and BaseCut is that a BaseCut::operator() returns a bool, not a float.

So, to implement one cut we could write:

#include "Ntuple.h"
#include "jbltools/sfh/BaseCut.h"

external Ntuple* ntuple; // global pointer to ntuple

class Jpsi2MuCut : public BaseCut {
  public: 
    virtual bool operator() () const {
      return nt->trackIsMuon[0] && nt->trackIsMuon[1];
    }
};

How can we form logical expressions of cuts?

Actually, this is very simple. Essentially, the code looks like this:

class AndOfTwoCuts: public BaseCut {
  public:
    AndOfTwoCuts (const BaseCut& lhs_, const BaseCut& rhs_)
    : lhs (&lhs_), rhs (&rhs_) {};
    
    virtual bool operator() () const {return (*lhs)() && (*rhs)();}
  protected:
    const BaseCut *lhs;
    const BaseCut *rhs;
};

inline AndOfTwoCuts& operator&& (const BaseCut& lhs_, const BaseCut& rhs_) {
  return *new AndOfTwoCuts (lhs_, rhs_);
}

(The real code, residing in BaseCut.h, is slightly more complicated, but not much.)

We first define a class AndOfTwoCuts that holds pointers to two other BaseCut objects, and returns the AND of their results. Additionally we overload operator&&(), so that it generates a new AndOfTwoCuts object, which we can pass (by reference!) on to a SFH1F. Voilà!

How can I Use these Self-Filling Histograms?

Typically, one creates one's own class that is derived from class EventLoop, let's name it AnalysisLoop.

The purpose of an AnalysisLoop object is to book, fill, and write out a number of histograms. It might look like this:

#include "jbltools/sfh/EventLoop.h"
#include "jbltools/sfh/Binning.h"
#include "jbltools/sfh/SFH1F.h"
#include "PTMissFun.h"
class AnalysisLoop: public EventLoop {
  public:
    // Constructor books histograms
    AnalysisLoop (Ntuple* ntuple) 
    {
      Binning ptmissbinning (50, 0., 100.);
      FloatFun& ptmiss = *new PTMissFun (ntuple);
      // self-filling histogram for missing pt
      h = new SFH1F ("ptmiss", "Missing pt", ptmissbinning, this, ptmiss);        
    
    }
    
    // Destructor does nothing
    virtual ~AnalysisLoop () 
    {}
    
    // output method writes histos to postscript and file
    virtual void output (const char* rootfile = "", 
                         const char* psfile = "") 
    {
      // create canvas 
      TCanvas *canvas = new TCanvas("canvas", "ptmiss", 600, 800);
      // Write plot to psfile
      TPostScript ps (psfile, 111);
      canvas->Clear();
      h->Draw ("E0");                 // Draws the histo
      canvas->Update();
      ps.Close(); 
      // Write histogram to file
      TFile file (rootfile, "RECREATE");
      h->Write();                     // Writes out the histo
      file.Write(); 
      file.Close();
    }
  protected:
    SFH1F *h;                         // The self-filling histogram
};

The corresponding main program opens an ntuple and loops over its entries, calling the (inherited) loop() method of the AnalysisLoop object for each row:

#define Ntuple_cxx
#include "Ntuple.h"                            // The ntuple
#include "AnalysisLoop.h"                      // The AnalysisLoop
int main(int argc, char** argv) {
  Ntuple nt;                                   // open ntuple
  TApplication theApp("main", &argc, argv);    // for graphics
  AnalysisLoop theAnalysisLoop (&nt);          // book histos 
  // Loop over ntuple entries
  for (int i = 0; i < int(nt.fChain->GetEntriesFast()); ++i) {
    if (nt.LoadTree (i) < 0) break;
    nt.fChain->GetEntry(i);
    theAnalysisLoop.loop();                    // fill the histos
  }
  theAnalysisLoop.output("out.root","out.ps"); // store the histos
  return 0;
}

Here we use the headerfile Ntuple.h, generated by TTree::MakeClass, that describes the ntuple structure.

AnalysisLoop is derived from EventLoop, and thereby from SFROList (list of self-filling registered objects). SFROList defines a Fill() method that fills all self-filling objects by calling their respective Fill() method. EventLoop::loop() calls this Fill() method, and therefore fills all self-filling objects that have registered themselves with the loop, in our example the ptmiss-histogram.

---