Revision 9a9e5b4cf3f43c391212ffc2c292ca1ffdba8767 authored by Unknown Author on 08 May 2006, 14:01:31 UTC, committed by Unknown Author on 08 May 2006, 14:01:31 UTC
git-svn-id: http://root.cern.ch/svn/root/tags/v5-11-02@14953 27541ba8-7e3a-0410-8455-c3a389f83636
1 parent 28e762f
h1analysisProxy.C
// 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. 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"
// 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.
// SlaveBegin():
//
// Notify(): This function is called at the first entry of a new Tree
// in a chain.
// Process(): called to analyze each entry.
//
// SlaveTerminate():
//
// 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.Process("h1analysis.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.Process("h1analysis.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.Process("h1analysis.C","useList")
//
//==> 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 with ,eg:
//
//==> F- Create the chain as in A, then execute
// Root > chain.Process("h1analysis.C+","useList")
//
// The commands executed with the 4 different methods B,C,D and E
// 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
TEventList *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());
//some cleanup in case this function had already been executed
//delete any previously generated histograms or functions
gDirectory->Delete("hdmd");
gDirectory->Delete("h2*");
delete gROOT->GetFunction("f5");
delete gROOT->GetFunction("f2");
}
//_____________________________________________________________________
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 event 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());
//some cleanup in case this function had already been executed
//delete any previously generated histograms or functions
gDirectory->Delete("hdmd");
gDirectory->Delete("h2*");
delete gROOT->GetFunction("f5");
delete gROOT->GetFunction("f2");
//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 event list
fillList = kFALSE;
useList = kFALSE;
tree->SetEventList(0);
delete gDirectory->GetList()->FindObject("elist");
// case when one creates/fills the event list
if (option.Contains("fillList")) {
fillList = kTRUE;
elist = new TEventList("elist","selection from Cut",5000);
} else elist = 0;
// case when one uses the event list generated in a previous call
if (option.Contains("useList")) {
useList = kTRUE;
TFile f("elist.root");
elist = (TEventList*)f.Get("elist");
if (elist) elist->SetDirectory(0); //otherwise the file destructor will delete elist
tree->SetEventList(elist);
}
}
Double_t h1analysisProxy() {
return 0;
}
//_____________________________________________________________________
Bool_t h1analysisProxy_Process(Int_t entry)
{
// entry is the entry number in the current Tree
// Selection function to select D* and D0.
//in case one event 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 #%d 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 #%d 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(fChain->GetChainEntryNumber(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,0,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();
//save the event list to a Root file if one was produced
if (fillList) {
TFile efile("elist.root","recreate");
elist->Write();
}
}
Computing file changes ...