/// \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: [https://root.cern.ch/files/h1/](https://root.cern.ch/files/h1/) /// 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 "TBranch.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(nullptr); 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 = nullptr; } 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(nullptr); //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(fOutput->FindObject("hdmd")); h2 = dynamic_cast(fOutput->FindObject("h2")); if (hdmd == nullptr || h2 == nullptr) { 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); 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) { if (!elist) elist = dynamic_cast(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); }