ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
h1analysis.C File Reference

Example of analysis class for the H1 data. More...

Go to the source code of this file.

Detailed Description

Example of analysis class for the H1 data.

This file uses 4 large data sets from the H1 collaboration at DESY Hamburg. One can access these data sets (277 MBytes) from the standard Root web site at: ftp:/// The Physics plots below generated by this example cannot be produced when using smaller data sets.

There are several ways to analyze data stored in a Root Tree

A chain of 4 files (originally converted from PAW ntuples) is used to illustrate the various ways to loop on Root data sets. Each data set contains a Root Tree named "h42" The class definition in h1analysis.h has been generated automatically by the Root utility TTree::MakeSelector using one of the files with the following statement:


This produces two files: h1analysis.h and h1analysis.C (skeleton of this file) The h1analysis class is derived from the Root class TSelector.

The following members functions are called by the TTree::Process functions.

To use this file, try the following sessions

Root > gROOT->Time(); /// will show RT & CPU time per command

Case A: Create a TChain with the 4 H1 data files

The chain can be created by executed the short macro h1chain.C below:

TChain chain("h42");
chain.Add("$H1/dstarmb.root"); /// 21330730 bytes 21920 events
chain.Add("$H1/dstarp1a.root"); /// 71464503 bytes 73243 events
chain.Add("$H1/dstarp1b.root"); /// 83827959 bytes 85597 events
chain.Add("$H1/dstarp2.root"); /// 100675234 bytes 103053 events
/// where $H1 is a system symbol pointing to the H1 data directory.

Case B: Loop on all events

Root > chain.Process("h1analysis.C")

Case C: Same as B, but in addition fill the entry list with selected entries.

The entry list is saved to a file "elist.root" by the Terminate function. To see the list of selected events, you can do elist->Print("all"). The selection function has selected 7525 events out of the 283813 events in the chain of files. (2.65 per cent)

Root > chain.Process("h1analysis.C","fillList")

Case D: Process only entries in the entry list

The entry list is read from the file in elist.root generated by step C

Root > chain.Process("h1analysis.C","useList")

Case E: The above steps have been executed via the interpreter.

You can repeat the steps B, C and D using the script compiler by replacing "h1analysis.C" by "h1analysis.C+" or "h1analysis.C++" in a new session (see F).

Case F: Create the chain as in A, then execute

Root > chain.Process("h1analysis.C+","useList")

The same analysis can be run on PROOF. For a quick try start a PROOF-Lite session

create (if mot already done) the chain by executing the 'h1chain.C' macro mentioned above, and then tell ROOT to use PROOF to process the chain:

You can then repeat step B above. Step C can also be executed in PROOF. However, step D cannot be executed in PROOF as in the local session (i.e. just passing option 'useList'): to use the entry list you have to

Case G: Load first in the session the list form the file

Root > TFile f("elist.root")
Root > TEntryList *elist = (TEntryList *) f.Get("elist")

set it on the chain:

call Process as in step B. Of course this works also for local processing.

#include "h1analysis.h"
#include "TH2.h"
#include "TF1.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TPaveStats.h"
#include "TLine.h"
#include "TMath.h"
const Double_t dxbin = (0.17-0.13)/40; // Bin-width
const Double_t sigma = 0.0012;
Double_t x = xx[0];
if (x <= 0.13957) return 0;
Double_t xp3 = (x-par[3])*(x-par[3]);
Double_t res = dxbin*(par[0]*TMath::Power(x-0.13957, par[1])
+ par[2] / 2.5066/par[4]*TMath::Exp(-xp3/2/par[4]/par[4]));
return res;
Double_t x = xx[0];
if (x <= 0.13957) return 0;
Double_t xp3 = (x-0.1454)*(x-0.1454);
Double_t res = dxbin*(par[0]*TMath::Power(x-0.13957, 0.25)
+ par[1] / 2.5066/sigma*TMath::Exp(-xp3/2/sigma/sigma));
return res;
void h1analysis::Begin(TTree * /*tree*/)
// function called before starting the event loop
// -it performs some cleanup
// -it creates histograms
// -it sets some initialisation for the entry list
// This is needed when re-processing the object
//print the option specified in the Process function.
TString option = GetOption();
Info("Begin", "starting h1analysis with process option: %s", option.Data());
//process cases with entry list
delete gDirectory->GetList()->FindObject("elist");
// case when one creates/fills the entry list
if (option.Contains("fillList")) {
elist = new TEntryList("elist", "H1 selection from Cut");
// Add to the input list for processing in PROOF, if needed
if (fInput) {
fInput->Add(new TNamed("fillList",""));
// We send a clone to avoid double deletes when importing the result
// This is needed to avoid warnings from output-to-members mapping
elist = 0;
if (fillList) Info("Begin", "creating an entry-list");
// case when one uses the entry list generated in a previous call
if (option.Contains("useList")) {
if (fInput) {
// Option "useList" not supported in PROOF directly
Warning("Begin", "option 'useList' not supported in PROOF - ignoring");
Warning("Begin", "the entry list must be set on the chain *before* calling Process");
} else {
TFile f("elist.root");
elist = (TEntryList*)f.Get("elist");
if (elist) elist->SetDirectory(0); //otherwise the file destructor will delete elist
// function called before starting the event loop
// -it performs some cleanup
// -it creates histograms
// -it sets some initialisation for the entry list
//initialize the Tree branch addresses
//print the option specified in the Process function.
TString option = GetOption();
"starting h1analysis with process option: %s (tree: %p)", option.Data(), tree);
//create histograms
hdmd = new TH1F("hdmd","dm_d",40,0.13,0.17);
h2 = new TH2F("h2","ptD0 vs dm_d",30,0.135,0.165,30,-3,6);
// Entry list stuff (re-parse option because on PROOF only SlaveBegin is called)
if (option.Contains("fillList")) {
// Get the list
if (fInput) {
if ((elist = (TEntryList *) fInput->FindObject("elist")))
// Need to clone to avoid problems when destroying the selector
if (elist)
if (fillList) Info("SlaveBegin", "creating an entry-list");
// entry is the entry number in the current Tree
// Selection function to select D* and D0.
//in case one entry list is given in input, the selection has already been done.
if (!useList) {
// Read only the necessary branches to select entries.
// return as soon as a bad entry is detected
// to read complete event, call fChain->GetTree()->GetEntry(entry)
b_md0_d->GetEntry(entry); if (TMath::Abs(md0_d-1.8646) >= 0.04) return kFALSE;
b_ptds_d->GetEntry(entry); if (ptds_d <= 2.5) return kFALSE;
b_etads_d->GetEntry(entry); if (TMath::Abs(etads_d) >= 1.5) return kFALSE;
b_ik->GetEntry(entry); ik--; //original ik used f77 convention starting at 1
b_ipi->GetEntry(entry); ipi--;
if (nhitrp[ik]*nhitrp[ipi] <= 1) return kFALSE;
if (rend[ik] -rstart[ik] <= 22) return kFALSE;
if (rend[ipi]-rstart[ipi] <= 22) return kFALSE;
b_nlhk->GetEntry(entry); if (nlhk[ik] <= 0.1) return kFALSE;
b_nlhpi->GetEntry(entry); if (nlhpi[ipi] <= 0.1) return kFALSE;
b_ipis->GetEntry(entry); ipis--; if (nlhpi[ipis] <= 0.1) return kFALSE;
b_njets->GetEntry(entry); if (njets < 1) return kFALSE;
// if option fillList, fill the entry list
if (fillList) elist->Enter(entry);
// to read complete event, call fChain->GetTree()->GetEntry(entry)
// read branches not processed in ProcessCut
b_dm_d->GetEntry(entry); //read branch holding dm_d
b_rpd0_t->GetEntry(entry); //read branch holding rpd0_t
b_ptd0_d->GetEntry(entry); //read branch holding ptd0_d
//fill some histograms
// Count the number of selected events
return kTRUE;
// nothing to be done
// function called at the end of the event loop
hdmd = dynamic_cast<TH1F*>(fOutput->FindObject("hdmd"));
h2 = dynamic_cast<TH2F*>(fOutput->FindObject("h2"));
if (hdmd == 0 || h2 == 0) {
Error("Terminate", "hdmd = %p , h2 = %p", hdmd, h2);
//create the canvas for the h1analysis fit
TCanvas *c1 = new TCanvas("c1","h1analysis analysis",10,10,800,600);
hdmd->GetXaxis()->SetTitle("m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]");
//fit histogram hdmd with function f5 using the loglikelihood option
if (gROOT->GetListOfFunctions()->FindObject("f5"))
delete gROOT->GetFunction("f5");
TF1 *f5 = new TF1("f5",fdm5,0.139,0.17,5);
f5->SetParameters(1000000, .25, 2000, .1454, .001);
//create the canvas for tau d0
TCanvas *c2 = new TCanvas("c2","tauD0",100,100,800,600);
// Project slices of 2-d histogram h2 along X , then fit each slice
// with function f2 and make a histogram for each fit parameter
// Note that the generated histograms are added to the list of objects
// in the current directory.
if (gROOT->GetListOfFunctions()->FindObject("f2"))
delete gROOT->GetFunction("f2");
TF1 *f2 = new TF1("f2",fdm2,0.139,0.17,2);
f2->SetParameters(10000, 10);
// Restrict to three bins in this example
Info("Fit Slices","Restricting fit to two bins only in this example...");
h2->FitSlicesX(f2,10,20,10,"g5 l");
TH1D *h2_1 = (TH1D*)gDirectory->Get("h2_1");
TLine *line = new TLine(0,0,0,c2->GetUymax());
// Have the number of entries on the first histogram (to cross check when running
// with entry lists)
TPaveStats *psdmd = (TPaveStats *)hdmd->GetListOfFunctions()->FindObject("stats");
//save the entry list to a Root file if one was produced
if (fillList) {
if (!elist)
elist = dynamic_cast<TEntryList*>(fOutput->FindObject("elist"));
if (elist) {
Printf("Entry list 'elist' created:");
TFile efile("elist.root","recreate");
} else {
Error("Terminate", "entry list requested but not found in output");
// Notify the amount of processed events
if (!fInput) Info("Terminate", "processed %lld events", fProcessed);
Rene Brun

Definition in file h1analysis.C.