Logo ROOT   6.10/09
Reference Guide
h1analysis.C File Reference

Detailed Description

View in nbviewer Open in SWAN 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:/// root.cern.ch/root/h1analysis 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:

h42->MakeSelector("h1analysis");

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 not already done) the chain by executing the 'h1chain.C' macro mentioned above, and then tell ROOT to use PROOF to process the chain:

Root > chain.SetProof()

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:

Root > chain.SetEntryList(elist)

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
Reset();
//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
if (fChain) fChain->SetEntryList(0);
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
fInput->Add(elist);
// This is needed to avoid warnings from output-to-members mapping
elist = 0;
}
Info("Begin", "creating an entry-list");
}
// case when one uses the entry list generated in a previous call
if (option.Contains("useList")) {
useList = kTRUE;
if (fInput) {
// In PROOF option "useList" is processed in SlaveBegin and we do not need
// to do anything here
} 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
Init(tree);
//print the option specified in the Process function.
TString option = GetOption();
Info("SlaveBegin",
"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);
fOutput->Add(hdmd);
fOutput->Add(h2);
// 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
elist = (TEntryList *) elist->Clone();
if (elist)
fOutput->Add(elist);
else
}
}
if (fillList) Info("SlaveBegin", "creating an entry-list");
if (option.Contains("useList")) useList = kTRUE;
}
{
// 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--;
b_nhitrp->GetEntry(entry);
if (nhitrp[ik]*nhitrp[ipi] <= 1) return kFALSE;
b_rend->GetEntry(entry);
b_rstart->GetEntry(entry);
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
hdmd->Fill(dm_d);
h2->Fill(dm_d,rpd0_t/0.029979*1.8646/ptd0_d);
// 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);
return;
}
//create the canvas for the h1analysis fit
TCanvas *c1 = new TCanvas("c1","h1analysis analysis",10,10,800,600);
c1->SetBottomMargin(0.15);
hdmd->GetXaxis()->SetTitle("m_{K#pi#pi} - m_{K#pi}[GeV/c^{2}]");
hdmd->GetXaxis()->SetTitleOffset(1.4);
//fit histogram hdmd with function f5 using the log-likelihood 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);
hdmd->Fit("f5","lr");
//create the canvas for tau d0
TCanvas *c2 = new TCanvas("c2","tauD0",100,100,800,600);
c2->SetGrid();
c2->SetBottomMargin(0.15);
// 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");
h2_1->GetXaxis()->SetTitle("#tau[ps]");
h2_1->SetMarkerStyle(21);
h2_1->Draw();
c2->Update();
TLine *line = new TLine(0,0,0,c2->GetUymax());
line->Draw();
// Have the number of entries on the first histogram (to cross check when running
// with entry lists)
TPaveStats *psdmd = (TPaveStats *)hdmd->GetListOfFunctions()->FindObject("stats");
psdmd->SetOptStat(1110);
c1->Modified();
//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:");
elist->Print();
TFile efile("elist.root","recreate");
elist->Write();
} 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);
}
Author
Rene Brun

Definition in file h1analysis.C.