h1analysis.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_tree
/// \notebook -header -nodraw
/// 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:
/// ~~~{.cpp}
///        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 every time 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
/// ~~~{.cpp}
///  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:
/// ~~~{.cpp}
///  {
///    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
/// ~~~{.cpp}
///  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)
/// ~~~{.cpp}
///  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
/// ~~~{.cpp}
///  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
/// ~~~{.cpp}
///  Root > chain.Process("h1analysis.C+","useList")
/// ~~~
///
/// The same analysis can be run on PROOF. For a quick try start a PROOF-Lite
/// session
/// ~~~{.cpp}
///  Root > TProof *p = TProof::Open("")
/// ~~~
/// 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:
/// ~~~{.cpp}
///  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
/// ~~~{.cpp}
///  Root > TFile f("elist.root")
///  Root > TEntryList *elist = (TEntryList *) f.Get("elist")
/// ~~~
///  set it on the chain:
/// ~~~{.cpp}
///  Root > chain.SetEntryList(elist)
/// ~~~
/// call Process as in step B. Of course this works also for local processing.
///
/// \macro_code
///
/// \author Rene Brun

#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;
      }
      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
      }
   }
}


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");
   if (option.Contains("useList")) useList  = kTRUE;
}


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 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
   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);
}