Logo ROOT  
Reference Guide
mathcoreVectorCollection.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_math
3/// \notebook -js
4/// Example showing how to write and read a std vector of ROOT::Math LorentzVector in a ROOT tree.
5/// In the write() function a variable number of track Vectors is generated
6/// according to a Poisson distribution with random momentum uniformly distributed
7/// in phi and eta.
8/// In the read() the vectors are read back and the content analysed and
9/// some information such as number of tracks per event or the track pt
10/// distributions are displayed in a canvas.
11///
12/// To execute the macro type in:
13///
14/// ~~~{.cpp}
15/// root[0]: .x mathcoreVectorCollection.C
16/// ~~~
17///
18/// \macro_image
19/// \macro_output
20/// \macro_code
21///
22/// \author Andras Zsenei
23
24#include "TRandom.h"
25#include "TStopwatch.h"
26#include "TSystem.h"
27#include "TFile.h"
28#include "TTree.h"
29#include "TH1D.h"
30#include "TCanvas.h"
31#include "TMath.h"
32
33#include <iostream>
34
35// CINT does not understand some files included by LorentzVector
36#include "Math/Vector3D.h"
37#include "Math/Vector4D.h"
38
39using namespace ROOT::Math;
40
41double write(int n) {
42 TRandom R;
43 TStopwatch timer;
44
45 TFile f1("mathcoreLV.root","RECREATE");
46
47 // create tree
48 TTree t1("t1","Tree with new LorentzVector");
49
50 std::vector<ROOT::Math::XYZTVector> tracks;
51 std::vector<ROOT::Math::XYZTVector> * pTracks = &tracks;
52 t1.Branch("tracks","std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > >",&pTracks);
53
54 double M = 0.13957; // set pi+ mass
55
56 timer.Start();
57 double sum = 0;
58 for (int i = 0; i < n; ++i) {
59 int nPart = R.Poisson(5);
60 pTracks->clear();
61 pTracks->reserve(nPart);
62 for (int j = 0; j < nPart; ++j) {
63 double px = R.Gaus(0,10);
64 double py = R.Gaus(0,10);
65 double pt = sqrt(px*px +py*py);
66 double eta = R.Uniform(-3,3);
67 double phi = R.Uniform(0.0 , 2*TMath::Pi() );
68 RhoEtaPhiVector vcyl( pt, eta, phi);
69 // set energy
70 double E = sqrt( vcyl.R()*vcyl.R() + M*M);
71 XYZTVector q( vcyl.X(), vcyl.Y(), vcyl.Z(), E);
72 // fill track vector
73 pTracks->push_back(q);
74 // evaluate sum of components to check
75 sum += q.x()+q.y()+q.z()+q.t();
76 }
77 t1.Fill();
78 }
79
80 f1.Write();
81 timer.Stop();
82 std::cout << " Time for new Vector " << timer.RealTime() << " " << timer.CpuTime() << std::endl;
83
84 t1.Print();
85 return sum;
86}
87
88double read() {
89 TRandom R;
90 TStopwatch timer;
91
92 TH1D * h1 = new TH1D("h1","total event energy ",100,0,1000.);
93 TH1D * h2 = new TH1D("h2","Number of track per event",21,-0.5,20.5);
94 TH1D * h3 = new TH1D("h3","Track Energy",100,0,200);
95 TH1D * h4 = new TH1D("h4","Track Pt",100,0,100);
96 TH1D * h5 = new TH1D("h5","Track Eta",100,-5,5);
97 TH1D * h6 = new TH1D("h6","Track Cos(theta)",100,-1,1);
98
99 TFile f1("mathcoreLV.root");
100
101 // create tree
102 TTree *t1 = (TTree*)f1.Get("t1");
103
104 std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > * pTracks = 0;
105 t1->SetBranchAddress("tracks",&pTracks);
106
107 timer.Start();
108 int n = (int) t1->GetEntries();
109 std::cout << " Tree Entries " << n << std::endl;
110 double sum=0;
111 for (int i = 0; i < n; ++i) {
112 t1->GetEntry(i);
113 int ntrk = pTracks->size();
114 h3->Fill(ntrk);
116 for (int j = 0; j < ntrk; ++j) {
117 XYZTVector v = (*pTracks)[j];
118 q += v;
119 h3->Fill(v.E());
120 h4->Fill(v.Pt());
121 h5->Fill(v.Eta());
122 h6->Fill(cos(v.Theta()));
123 sum += v.x() + v.y() + v.z() + v.t();
124 }
125 h1->Fill(q.E() );
126 h2->Fill(ntrk);
127 }
128
129 timer.Stop();
130 std::cout << " Time for new Vector " << timer.RealTime() << " " << timer.CpuTime() << std::endl;
131
132 TCanvas *c1 = new TCanvas("c1","demo of Trees",10,10,600,800);
133 c1->Divide(2,3);
134
135 c1->cd(1);
136 h1->Draw();
137 c1->cd(2);
138 h2->Draw();
139 c1->cd(3);
140 h3->Draw();
141 c1->cd(3);
142 h3->Draw();
143 c1->cd(4);
144 h4->Draw();
145 c1->cd(5);
146 h5->Draw();
147 c1->cd(6);
148 h6->Draw();
149
150 return sum;
151}
152
153int mathcoreVectorCollection() {
154
155 int nEvents = 10000;
156 double s1 = write(nEvents);
157 double s2 = read();
158
159 if (fabs(s1-s2) > s1*1.E-15 ) {
160 std::cout << "ERROR: Found difference in Vector when reading ( " << s1 << " != " << s2 << " diff = " << fabs(s1-s2) << " ) " << std::endl;
161 return -1;
162 }
163 return 0;
164}
165
166int main() {
167 return mathcoreVectorCollection();
168}
169
#define s1(x)
Definition: RSha256.hxx:91
#define R(a, b, c, d, e, f, g, h, i)
Definition: RSha256.hxx:110
float * q
Definition: THbookFile.cxx:89
int main(int argc, char *argv[])
Definition: cef_main.cxx:54
Class describing a generic displacement vector in 3 dimensions.
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
Definition: LorentzVector.h:59
The Canvas class.
Definition: TCanvas.h:23
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:620
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3338
void Draw(Option_t *option="") override
Draw this histogram with options.
Definition: TH1.cxx:3060
virtual Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:875
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
Stopwatch class.
Definition: TStopwatch.h:28
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
A TTree represents a columnar dataset.
Definition: TTree.h:79
TPaveText * pt
RVec< PromoteType< T > > cos(const RVec< T > &v)
Definition: RVec.hxx:1798
return c1
Definition: legend1.C:41
const Int_t n
Definition: legend1.C:16
TH1F * h1
Definition: legend1.C:5
TF1 * f1
Definition: legend1.C:11
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
constexpr Double_t E()
Base of natural log: .
Definition: TMath.h:93
constexpr Double_t Pi()
Definition: TMath.h:37
auto * t1
Definition: textangle.C:20
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345
void tracks()
Definition: tracks.C:49