h1analysisProxy.C

Go to the documentation of this file.

// 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. He has 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, he 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_SlaveBegin():
//
//    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_SlaveTerminate():
//
//    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
//
//==>   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.
// }
//
// Root > .x h1chain.C
//
//==>   B- loop on all events
// Root > chain.Draw("h1analysisProxy.C")
//
//==>   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")
//
//==>   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 <a href="gif/h1analysis_dstar.gif" >the Dstar plot</a> end_html
// begin_html <a href="gif/h1analysis_tau.gif" >the Tau D0 plot</a> end_html
//
//Author; Philippe Canal from original h1analysis.C by Rene Brun

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

---