h1analysisProxy.C File Reference

Example of analysis class for the H1 data using code generated by MakeProxy. More...

Go to the source code of this file.

Detailed Description

Example of analysis class for the H1 data using code generated by MakeProxy.

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"

h1analysProxy.C can be used either via TTree::Draw:

h42->Draw("h1analysisProxy.C");

or it can be used directly with TTree::MakeProxy, for example to generate a shared library. TTree::MakeProxy will generate a TSelector skeleton that include h1analysProxy.C:

h42->MakeProxy("h1sel","h1analysisProxy.C");

This produces one file: h1sel.h which does a #include "h1analysProxy.C" The h1sel class is derived from the Root class TSelector and can then be used as:

h42->Process("h1sel.h+");

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

• h1analysProxy_Begin(): Called everytime a loop on the tree starts. a convenient place to create your histograms.
• h1analysProxy_Notify(): This function is called at the first entry of a new Tree in a chain.
• h1analysProxy_Process(): called to analyze each entry.
• h1analysProxy_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 session

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.Draw("h1analysisProxy.C")

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

The event 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.Draw("h1analysisProxy.C","","fillList")

Case D: Process only entries in the event list

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

Root > chain.Draw("h1analysisProxy.C","","useList")

The commands executed with the 3 different methods B,C and D produce two canvases shown below: begin_html the Dstar plot end_html begin_html the Tau D0 plot end_html

TEntryList *elist;
Bool_t useList, fillList;
TH1F *hdmd;
TH2F *h2;


void h1analysisProxy_Begin(TTree *tree)
{
// function called before starting the event loop
//  -it performs some cleanup
//  -it creates histograms
//  -it sets some initialisation for the event list

   //print the option specified in the Process function.
   TString option = GetOption();
   printf("Starting (begin) h1analysis with process option: %s\n",option.Data());

   //process cases with event list
   fillList = kFALSE;
   useList  = kFALSE;
   if (fChain) fChain->SetEntryList(0);
   delete gDirectory->GetList()->FindObject("elist");

   // case when one creates/fills the event 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",""));
         fInput->Add(elist);
      }
   } else elist = 0;

   // case when one uses the event list generated in a previous call
   if (option.Contains("useList")) {
      useList  = kTRUE;
      if (fInput) {
         tree->SetEntryList(elist);
         TFile f("elist.root");
         elist = (TEntryList*)f.Get("elist");
         if (elist) elist->SetDirectory(0); //otherwise the file destructor will delete elist
      } else {
         // 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");
      }
   }
}


void h1analysisProxy_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();
   printf("Starting (slave) h1analysis with process option: %s\n",option.Data());

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

   //process cases with entry list
   fillList = kFALSE;
   useList  = kFALSE;

   // case when one creates/fills the entry list
   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;
   } else elist = 0;

   // case when one uses the entry list generated in a previous call
   if (option.Contains("useList")) {
      useList  = kTRUE;
      TFile f("elist.root");
      elist = (TEntryList*)f.Get("elist");
      if (elist) elist->SetDirectory(0); //otherwise the file destructor will delete elist
      if (tree) tree->SetEntryList(elist);
      else {
         // 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");
      }
   }

}

Double_t h1analysisProxy() {
   return 0;
}


Bool_t h1analysisProxy_Process(Long64_t entry)
{
// 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) {

      float f1 = md0_d;
      float f2 = md0_d-1.8646;
      bool test = TMath::Abs(md0_d-1.8646) >= 0.04;
      if (gDebug>0) fprintf(stderr,"entry #%lld f1=%f f2=%f test=%d\n",
                            fChain->GetReadEntry(),f1,f2,test);

      if (TMath::Abs(md0_d-1.8646) >= 0.04) return kFALSE;
      if (ptds_d <= 2.5) return kFALSE;
      if (TMath::Abs(etads_d) >= 1.5) return kFALSE;

      int cik = ik-1;    //original ik used f77 convention starting at 1
      int cipi = ipi-1;  //original ipi used f77 convention starting at 1

      f1 = nhitrp[cik];
      f2 = nhitrp[cipi];
      test = nhitrp[cik]*nhitrp[cipi] <= 1;
      if (gDebug>0) fprintf(stderr,"entry #%lld f1=%f f2=%f test=%d\n",
                            fChain->GetReadEntry(),f1,f2,test);

      if (nhitrp[cik]*nhitrp[cipi] <= 1) return kFALSE;
      if (rend[cik] -rstart[cik]  <= 22) return kFALSE;
      if (rend[cipi]-rstart[cipi] <= 22) return kFALSE;
      if (nlhk[cik] <= 0.1)    return kFALSE;
      if (nlhpi[cipi] <= 0.1)  return kFALSE;
      // fix because read-only
      if (nlhpi[ipis-1] <= 0.1) return kFALSE;
      if (njets < 1)          return kFALSE;

   }
   // if option fillList, fill the event list
   if (fillList) elist->Enter(entry);

   //fill some histograms
   hdmd->Fill(dm_d);
   h2->Fill(dm_d,rpd0_t/0.029979*1.8646/ptd0_d);

   return kTRUE;
}



void h1analysisProxy_SlaveTerminate()
{
   // nothing to be done
   printf("Terminate (slave) h1analysis\n");
}


void h1analysisProxy_Terminate()
{
   printf("Terminate (final) h1analysis\n");

   // 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
   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.
   TF1 *f2 = new TF1("f2",fdm2,0.139,0.17,2);
   f2->SetParameters(10000, 10);
   h2->FitSlicesX(f2,0,-1,1,"qln");
   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) {
      elist = dynamic_cast<TEntryList*>(fOutput->FindObject("elist"));
      if (elist) {
         TFile efile("elist.root","recreate");
         elist->Write();
      } else {
         Error("Terminate", "entry list requested but not found in output");
      }
   }
}

Author

Philippe Canal from original h1analysis.C by Rene Brun

Definition in file h1analysisProxy.C.