Carlos Escobar

To create histograms and ntuples in your Athena algorithm it is probably easiest to copy code from an example. You can find an example in the Athena examples in the ATLAS repository offline/Control/AthenaExamples/AthExHistNtup and also the example analysis job described in the page: Running an example job at NIKHEF. Below is an explanation of what is needed to include histogramming and ntuples in your Athena algorithms.

It is best to declare the pointers to you histograms as data members of you algorithm class. First, however, you need to predeclare the histogram classes you will use. There are only two, 1-D and 2-D histograms.

// Predeclare histogram classes that you use.
class IHistogram1D;
class IHistogram2D;

class MyAlgorithm : public Algorithm { ...

Then add private data members for all the histograms that you will declare.

class MyAlgorithm : public Algorithm {

  ...

  private:
    IHistogram1D* m_hist1;  // A 1-D histogram
    IHistogram2D* m_hist2;  // A 2-D histogram

}

What's needed in the implementation file (*.cxx)

You need to add the following include statements.

#include "GaudiKernel/IHistogramSvc.h"
#include "AIDA/IHistogram1D.h"
#include "AIDA/IHistogram2D.h"

Booking your histograms

Booking your histograms is normally done in the initialize() method of your algorithm. There are a number of alternative methods available for booking. Two that are common in examples are below and these should be sufficient for most users. There are also methods for variable binning. See the file GaudiKernel/IHistogramSvc.h for all the possibilities. In the two methods illustrated here, you can either specify the histogram id in the directory string or as a separate integer parameter. For Root there is no restriction on what you use for the histogram ID but to keep your code independent of persistency technology it is recommended to use a valid HBOOK identifier (ie an integer). The histSvc() function is a member of the Algorithm base class and so you automatically have access to this function.

Booking 1-D histograms:
m_hist1 = histoSvc()->book("/stat/mydirectory", 100, "Px of track",10, 0.,10.);
or
m_hist1 = histoSvc()->book("/stat/mydirectory/100", "Px of track", 10, 0.,10.);

The first parameter is the directory name of where to store the histogram in the data store. This directory name must start with "/stat". This is a special area in the data store that doesn't get deleted after each event (as is necessary with histograms). The next parameter is the histogram ID (or you can specify it in the directory name as in the second example). The next parameter is a descriptive name for the histogram. The parameters after that are number of bins, lower histogram edge and upper upper histogram edge. For 2-D histograms use

Booking 2-D histograms:
m_hist2 = histoSvc()->book("/stat/mydirectory", 110, "Px vs Py", 10, 0., 10., 10, 0., 10.);
or
m_hist2 = histoSvc()->book("/stat/mydirectory/110", "Px vs Py", 10, 0., 10., 10, 0., 10.);

For 2-D histograms you must specify the number of X bins, lower edge on X axis, upper edge on X axis, number of Y bins, lower edge on Y axis, and upper edge on Y axis.

Filling your histograms

You normally fill your histograms in the execute() method of your algorithm. To fill just specify the value and optionally the weight. For example:

Filling 1-D and 2-D histograms:

double px = ...;
double py = ...;
m_hist1->fill(px);     // 1-D histogram. Weight will be 1
m_hist2->fill(px, py); // 2-D histogram. Weight will be 1

or

double weight = ...;
m_hist1->fill(px, weight);      // 1-D histogram.
m_hist2->fill(px, py, weight);  // 2-D histogram.

Accessing your histograms from other algorithms

If you need to access you histograms in another algorithm (for example if you fill your histogram in one algorithm and do some manipulation on that histogram in another algorithm) it is possible to get the pointer to the histogram from the data store. This is done as follows:

SmartDataPtr<IHistogram1D> hist1D( histoSvc(), "/stat/mydirectory/100" );
if( 0 != hist1D ) {  
  histoSvc()->print( hist1D );  // Print the found histogram  
}

Other histogram manipulation

Although you are most likely to do any further manipulation in PAW or Root it is possible to do some manipulation in Athena. You can do things such as getting bin contents, extracting the rms, etc. You have available all the operations defined in the AIDA interface. See the AIDA web pages and follow the links to the header files for more details.

What's needed in the jobOptions file.

You must tell Athena what persistency method you wish for the histogram files. The choices are Root or HBOOK.

For HBOOK:

// Select the appropriate shared library
ApplicationMgr.DLLs += { "HbookCnv" };

// Select HBOOK or ROOT persistency 
ApplicationMgr.HistogramPersistency = "HBOOK";

// Specify the output file name
HistogramPersistencySvc.OutputFile  = "histo.hbook";

For Root:

// Select the appropriate shared library
ApplicationMgr.DLLs += { "RootHistCnv" };

// Select HBOOK or ROOT persistency 
ApplicationMgr.HistogramPersistency = "ROOT";

// Specify the output file name
HistogramPersistencySvc.OutputFile  = "histo.root";

Tip: Use the extension root for Root files so that the Root browser recognizes it as such.

Ntuples in ATHENA

What's needed in the header file (*.h)

You need the following include file

#include "GaudiKernel/NTuple.h"

You then need to declare all the ntuple variables (tags) as data members in your header file. You can have single values (Item), arrays (Array) and 2-d arrays (Matrix). You have available the following types: float, long, bool and double. However, booleans don't appear to be supported when using Root. If you include booleans they won't appear in the Root output file. Also, it is probably best not to use double as they require special attention in HBOOK and are not readable by PAW. It is recommended to use float instead. NB. There is no int, use long instead. Some examples of how to declare ntuple tags:

private:

  NTuple::Item<float>  m_higgsMass;
  NTuple::Item<bool>   m_flag;
  NTuple::Item<long>   m_nTracks;  
  NTuple::Item<long>   m_number;  
  NTuple::Array<float> m_pt;
  NTuple::Array<long>  m_index;  
  NTuple::Matrix<long> m_hits;

What's needed in the implementation file (*.cxx)

You need to add the following include statements.

#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"

Booking your ntuples

In the initialize() method you need to book your ntuple and declare the tags to the ntuple. This is done with the following commands.

    // Access the output file
    NTupleFilePtr file1(ntupleSvc(), "/NTUPLES/FILE1");

    NTuplePtr nt(ntupleSvc(), "/NTUPLES/FILE1/MyDir/10");
    if ( !nt )    {    // Check if already booked
      nt = ntupleSvc()->book ("/NTUPLES/FILE1/MyDir/10", CLID_ColumnWiseTuple, "Some Description");
      if ( nt )    {
        log << MSG::DEBUG << "booked ntuple " << endreq;

        status = nt->addItem ("PARTICLE/HiggsMass"  , m_higgsMass ) ;
        status = nt->addItem ("Flag"  ,    m_flag);
        status = nt->addItem ("index2",    m_number);
        status = nt->addItem ("NumTracks", m_nTracks, 0, 5000);
        status = nt->addItem ("Pt" ,       m_nTracks, m_pt) ;
        status = nt->addItem ("index",     m_nTracks, m_index, 0, 5000 );
        status = nt->addItem ("hits",      m_nTracks, m_hits, 2, 0, 8 );

      } else {   // did not manage to book the N tuple....
      }
    }

The top level directory of the n-tuple transient data store is called "/NTUPLES". The next directory layer is connected to the different output streams: e.g. "/NTUPLES/FILE1", where FILE1 is the logical name of the requested output file for a given stream. There can be several output streams connected to the service. After that you can have normal directory structures. Directory names should be less than or equal to 8 characters in order to be compatible with HBOOK. For other variants of the book method see the GaudiKernel/INTupleSvc.h file and the various examples listed above.

The case of the tag names is irrelevant in HBOOK but in Root the case will be preserved. You can have long variable names. The length restriction is probably 32 characters which is the limit in HBOOK.

You can specify block names (like in HBOOK column wise ntuples) by prepending the tag with the block name. eg if you want HiggsMass to be in a BLOCK named PARTICLE then call the tag "PARTICLE/HiggsMass". If no block name is specified then the BLOCK name is AUTO_BLK. These block names are ignored for root as far as I can tell. Block names have an 8 character limit.

WARNING: Tag names should be unique regardless of whether you use block names or not. ie you cannot have both tag names "PARTICLE/NUM" and "TRACK/NUM".

Filling your ntuples

The tags can be treated as normal variables. Just assign whatever is appropriate to these variable in your algorithm's execute() method

  m_higgsMass = ...
  m_nTracks = ...
  for (int i = 0; i < m_nTracks; i++) {
    m_pt = ...
  }

and then commit the ntuple row.

status = ntupleSvc()->writeRecord("/NTUPLES/FILE1/MyDir/10");

What's needed in the jobOptions file.

You must tell Athena the persistency method (HBOOK or ROOT) and the file name to associate with each stream.

For HBook:

NTupleSvc.Output     = { "FILE1 DATAFILE='ntuple1.hbook' TYP='HBOOK' OPT='NEW'" };

For Root:

NTupleSvc.Output     = { "FILE1 DATAFILE='ntuple1.root'  TYP='ROOT'  OPT='NEW'" };

The first parameter (FILE1) is the stream name. It can be any string and must correspond to the stream name you use in directory name above (eg "/NTUPLES/FILE1/MyDir"). As mentioned before you can have multiple streams.

For multiple streams you can have something like:

NTupleSvc.Output     = { "STREAM1 DATAFILE='ntuple1.hbook' TYP='HBOOK' OPT='NEW'",
                         "STREAM2 DATAFILE='ntuple2.hbook' TYP='HBOOK' OPT='NEW'" };

OPT can be

• NEW,CREATE,WRITE for creating a new file (the most common usage) or
• OLD,READ for reading an existing file or
• UPDATE for adding records to an existing file.

Reading your histogram or ntuples files in root

An example of how you can access your histograms and ntuples in root follows: (These examples assume you have created a root ntuple using the HiggsAnal as described in the Running an example job at NIKHEF.) and you have called your histogram file atlfast1.root and your ntuple file atlfast2.root.

The easiest way to look at your files is to use the ROOT browser.

> root

root[0] TBrowser b;

To access the root file and its contents from the Root command line you could do:

> root

For the histogram file:
root[0] TFile *f1 = new TFile("atlfast1.root");
root[1] h100 = f1->Get("simple1D/100");
root[2] h100->Draw();

For the tree (ntuple):
root[3] TFile *f2 = new TFile("atlfast2.root");
root[4] TTree *tree = (TTree *)f2->Get("/Atlfast/2002");
root[5] tree->Draw("H_M");

or write a root macro , for example

In a file named readfile.C

            
// Example of how to read a root file and access the histograms or tree
// produced in an Athena job.

// Make tree a global variable so that we can access it outside this function.
// You can then also access it from the command line. 
TTree *tree;

readHist(TString filename = "atlfast1.root")
{
  TFile *f = new TFile(filename);
  TH1F *h100 = (TH1F *)f->Get("simple1D/100");
  h100->Draw();
}

readTree(TString filename = "atlfast2.root")
{
  TFile *f = new TFile(filename);
  tree = (TTree *)f->Get("/Atlfast/2002");
  tree->Draw("H_M");

}

// An example of how to access the tree/ntuple variables
eventLoop()
{
  float h_m;
  tree->SetBranchAddress("H_M",&h_m);
  for (int i=0; i < tree->GetEntries(); i++){
    tree->GetEntry(i);
    cout << h_m << endl;
  }
}

Then from root (assuming the above file is called readfile.C) you can do the something like the following

> root

root[0] .L readfile.C
root[1] readHist("atlfast1.root");
root[2] readTree("atlfast2.root");
root[3] tree->Draw("H_M","H_PT>10");
root[4] eventLoop();

Further Information

Further information (although not a whole lot more) about using the histogramming and Ntuple services in Athena can be found in the Athena Developers guide which is available from the Athena Web Site.