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:/// 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

• Using TTree::Draw: This is very convenient and efficient for small tasks. A TTree::Draw call produces one histogram at the time. The histogram is automatically generated. The selection expression may be specified in the command line.
• Using the TTreeViewer: This is a graphical interface to TTree::Draw with the same functionality.
• Using the code generated by TTree::MakeClass: In this case, the user creates an instance of the analysis class. They have the control over the event loop and he can generate an unlimited number of histograms.
• Using the code generated by TTree::MakeSelector. Like for the code generated by TTree::MakeClass, the user can do complex analysis. However, they cannot control the event loop. The event loop is controlled by TTree::Process called by the user. This solution is illustrated by the current code. The advantage of this method is that it can be run in a parallel environment using PROOF (the Parallel Root Facility).

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.

• Begin(): Called everytime a loop on the tree starts. A convenient place to create your histograms.
• Notify(): This function is called at the first entry of a new Tree in a chain.
• Process(): Called to analyze each entry.
• Terminate(): Called at the end of a loop on a TTree. A convenient place to draw/fit your histograms.

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

Root > TProof *p = TProof::Open("")

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:

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 fdm5(Double_t *xx, Double_t *par)
{
   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 fdm2(Double_t *xx, Double_t *par)
{
   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")) {
      fillList = kTRUE;
      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;
      }
   }
   if (fillList) 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) {
         // 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
      }
   }
}


void h1analysis::SlaveBegin(TTree *tree)
{
// 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")) {
      fillList = kTRUE;
      // 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
            fillList = kFALSE;
      }
   }
   if (fillList) Info("SlaveBegin", "creating an entry-list");
}


Bool_t h1analysis::Process(Long64_t entry)
{
// entry is the entry number in the current Tree
// Selection function to select D* and D0.

   fProcessed++;
   //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_ntracks->GetEntry(entry);
      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
   fStatus++;

   return kTRUE;
}



void h1analysis::SlaveTerminate()
{
   // nothing to be done
}


void h1analysis::Terminate()
{
// 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
   gStyle->SetOptFit();
   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 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);
   hdmd->Fit("f5","lr");

   //create the canvas for tau d0
   gStyle->SetOptFit(0);
   gStyle->SetOptStat(1100);
   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.