|
||||||
|
|
This tutorial is based on ROOT 5.26 on a standalone system (i.e.
one not connected to the BaBar
disk system), but root versions as old as ROOT 4.04/02 should also
work.
You can also run it on a yakut machine using the analysis-51
(AKA 24.3.6) release using the modified instruction
located at
the end of some sections. In this case, you start root 5.14 by
typing bbrroot in the workdir
directory. Remember to srtpath first (from the release
directory). In general, I find it tidier to make a new release to
analyze the ntuples, instead of using the same one I used to make the
ntuples in the first place.
A slight familiarity with ROOT is assumed in this tutorial. For example, you should first have completed the WorkBook sections Root1 or Root2 or the FNAL root class.
This WorkBook section does not cover fitting of histograms.
[~]: cat .rootrc
# Path used to find macros.
Unix.*.Root.MacroPath: .:./macros:$(HOME)/Work/macros:/sw/share/root/macros:
[~]:
rootlogin.C is run whenever you start root. I
put mine in the macro directory defined in the previous step, but you
could put it in your analysis directory if you wanted to specialize it
for a particular analysis. Here is mine:
[~]: cat Work/macros/rootlogon.C
{
printf("\nWelcome to rootlogon.C \n\n"); <-- print a message
gROOT->SetStyle("Plain"); <-- plain histogram style
gStyle->SetOptStat(1111111); <-- expanded stats box
gStyle->SetPadTickX(1); <-- tic marks on all axes
gStyle->SetPadTickY(1); <--
}
The macro that is run when starting ROOT is called RooLogon.C,
and is also located in workdir. It
contains a
reasonable default set of commands, but you can edit it if you would
like something different.
Create a directory "data" of your analysis directory and
place the two files in it. While you are at it, you can also copy these
ntuples, which are made from generic B+ MC: SP-1235-Jpsitoll-Run1-R18c-1.root,
SP-1235-Jpsitoll-Run1-R18c-2.root,
SP-1235-Jpsitoll-Run1-R18c-3.root
. We will use these later.
Now use MakeSelector to create a
C++ program
(histAnalysis.h and histAnalysis.C) to
analyze the data:
[~/Work/tutorials/root3]: root
*******************************************
* *
* W E L C O M E to R O O T *
* *
* Version 5.26/00 14 December 2009 *
* *
* You are welcome to visit our Web site *
* http://root.cern.ch *
* *
*******************************************
ROOT 5.26/00 (trunk@31882, Dec 14 2009, 20:18:36 on macosx64)
CINT/ROOT C/C++ Interpreter version 5.17.00, Dec 21, 2008
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
Welcome to rootlogon.C <-- note the message from rootlogin.C
For approved plots use: gROOT->SetStyle("BABAR");
root [0] TFile f("data/SP-989-Run1-1.root") <-- either file will do
root [1] .ls
TFile** data/SP-989-Run1-1.root Created for you by RooTupleManager
TFile* data/SP-989-Run1-1.root Created for you by RooTupleManager
KEY: TTree ntp1;1 Analysis Ntuple
KEY: TTree ntp2;1 myNtuple
root [2] ntp2->MakeSelector("histAnalysis")
Info in <TTreePlayer::MakeClass>: Files: histAnalysis.h and histAnalysis.C generated from TTree: ntp2
(Int_t)0
root [3] .q
This is what the original files look like:
By default, the block sizes are set to the largest value needed for the particular file used with MakeSelector. For example, although we allowed for up to 50 Jpsi candidates when we made the ntuples (see Analysis.tcl in the root3 tutorial), only 1 is allowed for in histAnalysis.h: Float_t JpsiChi2[1]; //[nJpsi]. You will need to go through and edit all of these to the correct values. The resulting file is histAnalysis-blocksize.h.histAnalysis.C as
generated. If you want, you can delete them all to tidy up the file;
you can always check the histAnalysis-orig.C link above.
There are three steps to creating a histogram. First, it
must be
declared prior to the Begin function.
TH1F *hnJpsi, *hmassJpsi;
Next, define the histograms in the Begin
function:
//..Histogram definitions
hnJpsi = new TH1F("hnJpsi", "Number of Jpsi", 20, -0.5, 19.5);
hmassJpsi = new TH1F("hmassJpsi", "Raw Jpsi Mass", 60, 2.9, 3.3);
Process//..Just read the full event
fChain->GetTree()->GetEntry(entry);
//..Fill some simple histograms
hnJpsi->Fill(nJpsi);
for(Int_t i=0; i<nJpsi; i++) {
hmassJpsi->Fill(JpsiUncMass[i]);
}
root [0] TFile f("data/SP-989-Run1-1.root");
root [1] ntp2->Process("histAnalysis.C+");
To run only 100 events, you would say
ntp2->Process("histAnalysis.C+","",100,0);
The following commands create a canvas and display the
two
histograms on it. You can save the canvas using the File menu if you
wish. I normally use eps output format, since it is more convenient for
inserting into documents, but gif is useful for web pages.
root [2] TCanvas MyC("MyC");
root [3] MyC.Divide(2,1);
root [4] MyC.cd(1);
root [5] hnJpsi->Draw()
root [6] MyC.cd(2);
root [7] hmassJpsi->Draw()
Here is a screen capture of the root session:
And the resulting histograms: This version of the code is
TChain chain("ntp2");
chain.Add("data/SP-989-Run1-1.root");
chain.Add("data/SP-989-Run1-2.root");
chain.Process("histAnalysis.C+");
We will use this later.
If your ntuple is in a root subdirectory, the syntax is
(for
example) TChain chain("TauMicroFilter/ntpl"); The
remaining
commands are the same.
a. After filling the histograms, write them to a file:
Begin):
TFile
*fHistFile; //..Output file for histograms
fHistFile = new TFile("histAnalysis.root","RECREATE");
//..Write out histograms to file
fHistFile->cd();
fHistFile->Write();
fHistFile->Close();
cout << "Output file written" << endl;
hnJpsi->Write();
b. running the job now creates the output file:
[~/Work/tutorials/root3]: ls -l histAnalysis.root
-rw-r--r-- 1 hearty staff 4111 Jun 19 20:26 histAnalysis.root
c. Read the file into root to access it. Here is a sample session.
root [0] TFile f("histAnalysis.root");
root [1] .ls
TFile** histAnalysis.root
TFile* histAnalysis.root
KEY: TH1F hnJpsi;1 Number of Jpsi
KEY: TH1F hmassJpsi;1 Raw Jpsi Mass
root [2] TCanvas MyC("MyC");
root [3] MyC.Divide(2,1);
root [4] MyC.cd(1)
(class TVirtualPad*)0x3cea400
root [5] hnJpsi->Draw()
root [6] MyC.cd(2)
(class TVirtualPad*)0x3a04a00
root [7] hmassJpsi->Draw()
d. This version of the code is
The resulting histogram file is histAnalysis-2.root.[~/Work/tutorials/root3]: cat runAnalysis.C
#include "TChain.h"
void runAnalysis () {
TChain chain("ntp2");
#include "SP-989-Run1-Chain.C"
chain.Process("histAnalysis.C+");
}
Note that I include the files to be chained in an include file:
[~/Work/tutorials/root3]: cat SP-989-Run1-Chain.C
chain.Add("data/SP-989-Run1-1.root");
chain.Add("data/SP-989-Run1-2.root");
This is a convenient way to deal with multiple data
types - just
create a similar include file for each (on-peak data, off-peak data,
BB MC, continuum MC and so on). List them all in your macro and
comment out the ones to skip by starting the line with "//". Now we can
compile and execute the macro to run our analysis code on the full
chain. (Actually, since we haven't changed histAnalysis.C, root
doesn't recompile it).
root [0] .x runAnalysis.C+
Info in <TUnixSystem::ACLiC>: creating shared library /Users/hearty/Work/tutorials/root3/./runAnalysis_C.so
histograms written to histAnalysis.root
Here is the resulting histogram file: histAnalysis-3.root.
Note that the include file only works when compiling.
Here are versions that work in interpreted mode, which you run by typing
.x
runAnalysis-i.C. They lack
the flexibility of modifying datasets by adding or subtracting include
files. Note that the analysis code itself is still compiled in this
example. This macro can also be used in batch mode (next step):
[~/Work/tutorials/root3]: root -b -q runAnalysis.C+ >& runAnalysis.log &The "-b" runs in batch mode, "-q" exits root after running the macro. Here is the log: runAnalysis-3.log.
[1] 7544
[~/Work/tutorials/root3]:
[1] Done root -b -q runAnalysis.C+ >& runAnalysis.log
You can skip the log file if you like: root -b -q
runAnalysis-i.C+
b. You can now view the histograms by reading histAnalysis.root. It is convenient to leave an interactive root session running for this purpose. Just close and reopen the file when you remake it:
root [0] TFile f("histAnalysis.root");
root [1] hnJpsi->Draw();
<TCanvas::MakeDefCanvas>: created default TCanvas with name c1
root [2] f.Close();
root [3] hnJpsi->Draw();
Error: Symbol hnJpsi is not defined in current scope FILE:(tmpfile) LINE:1
Error: Failed to evaluate hnJpsi->Draw()Possible candidates are...
filename line:size busy function type and name
*** Interpreter error recovered ***
root [4] TFile f("histAnalysis.root");
root [5] hnJpsi->Draw();
It is likely that your data, and perhaps the ntuples you produce from your data, are located at a remote computing center. A convenient way to operate under these circumstances is to run root in batch mode at the remote site, and to run it interactively on your local machine. You then need to copy only the small histogram file histAnalysis.root over the network. This is very much faster than runnng an X window over the network. If you use AFS, the copying happens seamlessly.
bbrroot -b -q runAnalysis.C+ >& runAnalysis.log &or
bbrroot -b -q runAnalysis.C+
[~/Work/tutorials/root3]: cat runAnalysis.CWe can use TLorentzVector to create a 4-vector of the center of mass, from which we can get the center of mass energy. First, #include "TLorentzVector.h" in histAnalysis.C. Then, in Process:
#include "TChain.h"
#include "TSystem.h"
void runAnalysis () {
gSystem->Load("libPhysics.so");
TChain chain("ntp2");
#include "SP-989-Run1Chain.C"
chain.Process("histAnalysis.C+");
}
TLorentzVector mom4Ups(eePx, eePy, eePz, eeE);See the "Physics Vectors" chapter of the root users manual for more details on this class. Here is this version of the file: runAnalysis-CheckMC.C.
Float_t sqrts = mom4Ups.M();
a. declare the function in histAnalysis.h, right after the Terminate declaration:
virtual void Terminate();
Int_t CheckMC(Int_t mclund[200], Int_t mothidx[200], Int_t daulen[200],
Int_t dauidx[200], Int_t mclen, Float_t sqrts);
ClassDef(histAnalysis,0);
(Note that the second line above is broken only for formatting purposes.) Here is the full file: histAnalysis-CheckMC.h
b. Write the function. I prefer to put it in a separate file, CheckMC.C, but some people like to have all of their code in one file. Here are the first few lines:
// Value of CheckMC tells what type of events this is:
// 0 = on peak data, 1 = off peak, 2 = Jpsi K+, 3 = Other B-->Jpsi
// 4 = Other BB, 5 = continuum
Int_t histAnalysis::CheckMC(Int_t mcLund[200], Int_t mothIdx[200], Int_t dauLen[200],
Int_t dauIdx[200], Int_t mcLen, Float_t sqrts) {
(Note that the line above is broken only for formatting purposes.) Note the "histAnalysis::CheckMC". Here is the full file: CheckMC.C
c. Include the new file into your code with #include
"CheckMC.C":
//-----------------------------------------------------------------------------
#include "CheckMC.C"
d. Lets add a new histogram, hdtype, and fill it with the value returned by CheckMC.
//..Categorize the event
TLorentzVector mom4Ups(eePx, eePy, eePz, eeE);
Float_t sqrts = mom4Ups.M();
Int_t dtype = histAnalysis::CheckMC(mcLund, mothIdx, dauLen, dauIdx, mcLen,
sqrts);
hdtype->Fill(dtype);
Note that you include the class histAnalysis when using
the function.
If you check the histogram, note that every entry is 2, as expected.
For something different, try the generic B+B- ntuples (SP-1235)
instead.
(You will need to create a chain out of them, and modify your
runAnalysis.C). In fact, most
of these are also signal events - it is often the case that when you
run on generic BB MC, you will need to identify the signal events mixed
among the other BB events using code similar to CheckMC.C.
Here is the full code: histAnalysis-CheckMC.C
and the new chain file: SP-1235-Jpsitoll-Run1-Chain.C
In this example, we will use a single class, MyF.
Other classes can be similarly created as
needed.
a. Create the header file MyF.h. This should be put in your Macros directory defined in .rootrc.
[~/Work/tutorials/root3]: cat ~/Work/macros/MyF.h
#include "TLorentzVector.h"
class MyF {
public:
//--------------------------------------------------------------------------
//..FourMomentum returns four vector from p, theta, phi and mass
static TLorentzVector FourMomentum(Float_t mom, Float_t theta, Float_t phi,
Float_t mass);
};
This first function makes a TLorentzVector from the p,theta, phi and mass of the particle. (There is no such standard constructor). The key features:
[~/Work/tutorials/root3]: cat ~/Work/macros/MyF.C
#include "MyF.h"
//--------------------------------------------------------------------------
//..FourMomentum returns four vector from p, Cos(theta), phi and mass
TLorentzVector MyF::FourMomentum(Float_t mom, Float_t costheta, Float_t phi, Float_t mass)
{
TLorentzVector temp(1.0,1.0,1.0,1.0);
Float_t energy = sqrt(mom*mom + mass*mass);
temp.SetRho(mom);
temp.SetTheta(acos(costheta));
temp.SetPhi(phi);
temp.SetE(energy);
return temp;
}
You could also put the code into a file FourMomentum.C
and just include
it in MyF.C using #include "FourMomentum.C"
c. Compile the code and load the library to make the code accessible.
To use the code, we need to compile and load it before running histAnalysis. We do this in runAnalysis.C, right after loading libPhysics.so:
gSystem->Load("libPhysics.so");
gROOT->LoadMacro("MyF.C+");
Recall that the ".C+" extension means that MyF.C is compiled if
changed; otherwise the existing library is loaded. You will also need
to include the corresponding header file for gROOT: #include "TROOT.h"
d. Create a soft link to your Macros directory:
[~/Work/tutorials/root3]: ln -s ~/Work/macros Macros
e. Use the function in your analysis.
First, you need to include the header: #include
"Macros/MyF.h"
Here is a snippet of code to calculate the four vector
and corresponding missing mass
for every B
candidate. Actually, a proper calculation of this quantity using a
kinematic fit is available in the ntuple as BpostFitMmiss[ib]. You
might find it interesting to compare this quantity to Mmiss:
//..Missing mass
for(Int_t ib=0; ib<nb; ib++) {
TLorentzVector mom4B =
MyF::FourMomentum(Bp3[ib], Bcosth[ib], Bphi[ib], BMass[ib]);
Float_t Mmiss = (mom4Ups - mom4B).M();
hmyBmiss->Fill(Mmiss);
hntpBmiss->Fill(BpostFitMmiss[ib]);
hCompMiss->Fill(BpostFitMmiss[ib],Mmiss);
}
Remember to declare and
create the histograms. Here are the complete files.
a. Declare the relevant histograms to be arrays of length 6:
TH1F *hnJpsi[6], *hmassJpsi[6], *hdtype, *hmyBmiss[6], *hntpBmiss[6];
TH2F *hCompMiss[6];
TFile *fHistFile;
Note that hdtype is still only a single histogram.
b. Before defining the histograms in the Begin method of
histAnalysis, create a vector of names and use them to customize each
histogram:
//..Some quantities to handle multiple data types
TString htit[6] = {"_on","_off","_sig","_Jpsi","_BB","_cont"};
//..Define histograms
hdtype = new TH1F("hdtype", "Event Type", 10, -1.5, 8.5);
for( Int_t ih=0; ih<6; ih++ ) {
hnJpsi[ih] = new TH1F("hnJpsi"+htit[ih], "Number of Jpsi"+htit[ih],
20, -0.5, 19.5);
hnJpsi[ih]->Sumw2();
hmassJpsi[ih] = new TH1F("hmassJpsi"+htit[ih], "Raw Jpsi Mass"+htit[ih],
60, 2.9, 3.3);
hmassJpsi[ih]->Sumw2();
hmyBmiss[ih] = new TH1F("hmyBmiss"+htit[ih], "My Missing Mass"+htit[ih],
50, 5.20, 5.30);
hmyBmiss[ih]->Sumw2();
hntpBmiss[ih] = new TH1F("hntpBmiss"+htit[ih],
"Ntp Missing Mass"+htit[ih], 50, 5.20, 5.30);
hntpBmiss[ih]->Sumw2();
hCompMiss[ih] = new TH2F("hCompMiss"+htit[ih],
"My Mmiss vs Ntp Mmiss"+htit[ih],50, 5.20,
5.30, 50, 5.20, 5.30);
hCompMiss[ih]->Sumw2();
}
c. Now create an array of weights corresponding to each data type. The weight = on-peak-luminosity/sample-luminosity. Put this before the Begin function since it will probably be needed in more than one function of histAnalysis. (Of course, you will need to calculate these for your own analysis!)
Float_t hwt[6] = {1., 8.55, 1., 0.241, 0.283, 1.098};
d. Now fill the histograms according to the event type and with the appropriate weight. Note that dtype enters both in the histogram identifier and in the weight.
hnJpsi[dtype]->Fill(nJpsi,hwt[dtype]);
and so forth for hmassJpsi and the others.
e. To test this out, use the B+B- MC (SP
mode 1235). Run the job using root -q -b
runAnalysis.C+ as normal.
f. Take a look at the resulting histograms in your interactive root, starting with the histogram of the event type, hdtype->Draw(). Note that most of the events are actually B+ --> Jpsi K+ signal events. This is not surprising, since we are requiring a reconstructed Jpsi K+ to write the ntuple. Here is a gif version of the plot (or eps). Take a look at the Jpsi masses for each data type as well:
root [4] TCanvas MyC("MyC");
root [5] MyC.Divide(3,2)
root [6] MyC.cd(1)
root [7] hmassJpsi_on->Draw()
root [8] MyC.cd(2)
root [9] hmassJpsi_off->Draw()
root [10] MyC.cd(3)
root [11] hmassJpsi_sig->Draw()
root [12] MyC.cd(4)
root [13] hmassJpsi_Jpsi->Draw()
root [14] MyC.cd(5)
root [15] hmassJpsi_BB->Draw()
root [16] MyC.cd(6)
root [17] hmassJpsi_cont->Draw()
root [18] Info in <TCanvas::Print>: GIF file /Users/hearty/Work/tutorials/root26/Jpsi-mass-SP1235.gif
has been created
Here is resulting plot.
Note that the on-peak, off-peak and continuum plots are
empty,
since we didn't run on those data, and that the Signal sample gives a
nice Jpsi peaks, while the "Other BB" sample does
not.
Also note that for the signal plot "Entries" (which is
the number of
times that Fill was called for this histogram) is equal to "underflow"
+ "overflow" + "integral", which are the sum of weights below, above
and in the plot respectively. This is in contrast to the BB plot,
where the two values differ by the ratio of 0.283 - the BB luminosity
weight.
When doing an analysis, you often find yourself making
the same set of plots over and over. In this case, it is worth writing
a macro to make it. To execute it, type
.x plotMass.C in an interactive root
session.
This is a much easier way to produce the same plot!
Here are the final versions of various files:
RooLogon.C,
located in the workdir
package.
rootlogon.Crootlogon.C
(see "Some basics"), which defines a style
called BABAR. Put it before your normal commands, to be sure that
you
end up in
your normal style. Here is my resulting rootlogon.C.
You then invoke this style by the command
gROOT->SetStyle("BABAR");
workdir/RooAlias.C. Add
this to your macros directory: BABARSmartLabel.C
Note that there is no header file; it is assumed that
this macro
will not be used in compiled macros. You will need to load the macro
either using .L BABARSmartLabel.C interactively, or gROOT->LoadMacro("BABARSmartLabel.C");
in
the
histogramming macro. (You don't need to specify the path to
BABARSmartLabel.C if it's in your macros directory as defined in Some Basics). After loading the macro, the command
BABARSmartLabel(-1,-1,-1,"-1",-1);will put "BaBar" in the upper right-hand corner;
BABARSmartLabel(-1,-1,-1,"~1",-1,-1);will add "Preliminary" underneath. Check the comments in the code for the full usage.
In this macro, we will make a nice version of the Jpsi
mass plot
produced above. Here is the root file containing the produced
histograms, in case you don't have your own: histAnalysis-final.root.
Here is the macro: DrawsighmassJpsi.C
This macro also includes a couple of lines with commands for drawing lines on figures in ROOT.
A key point for manipulating histograms in a macro is that under many circumstances (in particular, if you want to work with more than one in the macro), you need to explicitly create a pointer to each histogram:
TFile f("histAnalysis.root");
TH1* sighmassJpsi = (TH1*) f.Get("sighmassJpsi");
In addition to invoking the BABAR style, this macro does some typical
modifications, like adjusting the marker size and moving the axis
labels. Note that the BABAR style turns off the default titles, so you
will need to add titles yourself. It then writes the plot to an eps
file.
[~/Work/tutorials/root3]: root
*******************************************
* *
* W E L C O M E to R O O T *
* *
* Version 5.26/00 14 December 2009 *
* *
* You are welcome to visit our Web site *
* http://root.cern.ch *
* *
*******************************************
ROOT 5.26/00 (trunk@31882, Dec 14 2009, 20:18:36 on macosx64)
CINT/ROOT C/C++ Interpreter version 5.17.00, Dec 21, 2008
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
Welcome to rootlogon.C
For approved plots use: gROOT->SetStyle("BABAR");
root [0] .x DrawsighmassJpsi.C
Info in <TCanvas::Print>: eps file JpsiMass.eps has been created
root [1] .q
[~/Work/tutorials/root3]:
If you make changes to DrawsighmassJpsi.C, you don't need to quit and restart root; just rerun the macro. Here is the resulting plot: JpsiMass.eps
RooLogon.C
or to your Macros directory.
However, you do need to delete the following command in DrawsighJpsi.C,
which loads BABARSmartLabel:
gROOT->LoadMacro("BABARSmartLabel.C");
Here is the resulting macro: DrawsighmassJpsi-yakut.C
MakeSelector, added histograms, code to
deal with multiple MC types, and external functions, and developed a
macro to load the required libraries and execute the code. These tools
are adequate for a cut-and-count type analysis.
The next step in sophistication would be fitting. Basic fits are discussed in the ROOT II tutorial. The package RooFit provides many convenient and advanced tools for this purpose.
Last modified: 29-Jun-2006
Last significant update: 27-Jun-2006
Updated to analysis-31: 27-Jun-2006
Updated to analysis-51 and root 5.26: 29-Jan-2010
Send comments to Workbook Team.