ROOT   6.10/09 Reference Guide
testVectorIO.cxx
Go to the documentation of this file.
1
2 //
3 // Cint macro to test I/O of mathcore Lorentz Vectors in a Tree and compare with a
4 // TLorentzVector. A ROOT tree is written and read in both using either a XYZTVector or /// a TLorentzVector.
5 //
6 // To execute the macro type in:
7 //
8 // root[0]: .x mathcoreVectorIO.C
9
10
11 #ifdef USE_REFLEX
12 #include "Cintex/Cintex.h"
13 #include "Reflex/Reflex.h"
14 #endif
15
16 #include "TRandom3.h"
17 #include "TStopwatch.h"
18 #include "TSystem.h"
19 #include "TFile.h"
20 #include "TTree.h"
21 #include "TH1D.h"
22 #include "TCanvas.h"
23
24 #include <iostream>
25
26 #include "TLorentzVector.h"
27
28 #include "Math/Vector4D.h"
29 #include "Math/Vector3D.h"
30 #include "Math/Point3D.h"
31
32 #include "Track.h"
33
34
35 #define DEBUG
36
37
39
40 //#define USE_REFLEX
41
42 #define DEFVECTOR4D(TYPE) \
43 typedef TYPE AVector4D; \
44 const std::string vector4d_type = #TYPE ;
45
46 #define DEFVECTOR3D(TYPE) \
47 typedef TYPE AVector3D; \
48 const std::string vector3d_type = #TYPE ;
49
50 #define DEFPOINT3D(TYPE) \
51 typedef TYPE APoint3D; \
52 const std::string point3d_type = #TYPE ;
53
54
55
56 //const double tol = 1.0E-16;
57 const double tol = 1.0E-6; // or doublr 32 or float
58
60 //DEFVECTOR4D(ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<Double32_t> >);
61
63
65
66
67 // DEFVECTOR4D(ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<Double32_t> >)
68
69 //using namespace ROOT::Math;
70
71 template<class Vector>
72 inline double getMag2(const Vector & v) {
73  return v.mag2();
74 }
75
76 // inline double getMag2(const VecTrackD & v) {
77 // // print the read points
78 // std::cout << "VecTRackD " << std::endl;
79 // for (VecTrackD::It itr = v.begin() ; itr != v.end(); ++itr)
80 // std::cout << (*itr).Pos() << std::endl;
81
82 // return v.mag2();
83 // }
84
85
86 inline double getMag2(const TVector3 & v) {
87  return v.Mag2();
88 }
89
90 inline double getMag2(const TLorentzVector & v) {
91  return v.Mag2();
92 }
93
94 template<class U>
95 inline void setValues(ROOT::Math::DisplacementVector3D<U> & v, const double x[] ) {
96  v.SetXYZ(x[0],x[1],x[2]);
97 }
98
99 template<class U>
100 inline void setValues(ROOT::Math::PositionVector3D<U> & v, const double x[]) {
101  v.SetXYZ(x[0],x[1],x[2]);
102 }
103
104 template<class U>
105 inline void setValues(ROOT::Math::LorentzVector<U> & v, const double x[]) {
106  v.SetXYZT(x[0],x[1],x[2],x[3]);
107 }
108 // specialization for T -classes
109 inline void setValues(TVector3 & v, const double x[]) {
110  v.SetXYZ(x[0],x[1],x[2]);
111 }
112 inline void setValues(TLorentzVector & v, const double x[]) {
113  v.SetXYZT(x[0],x[1],x[2],x[3]);
114 }
115
116
117 template <class Vector>
118 double testDummy(int n) {
119
120  TRandom3 R(111);
121
123
124  Vector v1;
125  double s = 0;
126  double p[4];
127
128  timer.Start();
129  for (int i = 0; i < n; ++i) {
130  p[0] = R.Gaus(0,10);
131  p[1] = R.Gaus(0,10);
132  p[2] = R.Gaus(0,10);
133  p[3] = R.Gaus(100,10);
134  setValues(v1,p);
135  s += getMag2(v1);
136  }
137
138  timer.Stop();
139
140  double sav = std::sqrt(s/double(n));
141
142  std::cout << " Time for Random gen " << timer.RealTime() << " " << timer.CpuTime() << std::endl;
143  int pr = std::cout.precision(18); std::cout << "Average : " << sav << std::endl; std::cout.precision(pr);
144
145  return sav;
146 }
147
148 //----------------------------------------------------------------
149 /// writing
150 //----------------------------------------------------------------
151
152 template<class Vector>
153 double write(int n, const std::string & file_name, const std::string & vector_type, int compress = 0) {
154
156
157  TRandom3 R(111);
158
159  std::cout << "writing a tree with " << vector_type << std::endl;
160
161  std::string fname = file_name + ".root";
162  TFile f1(fname.c_str(),"RECREATE","",compress);
163
164  // create tree
165  std::string tree_name="Tree with" + vector_type;
166  TTree t1("t1",tree_name.c_str());
167
168  Vector *v1 = new Vector();
169  std::cout << "typeID written : " << typeid(*v1).name() << std::endl;
170
171  t1.Branch("Vector branch",vector_type.c_str(),&v1);
172
173  timer.Start();
174  double p[4];
175  double s = 0;
176  for (int i = 0; i < n; ++i) {
177  p[0] = R.Gaus(0,10);
178  p[1] = R.Gaus(0,10);
179  p[2] = R.Gaus(0,10);
180  p[3] = R.Gaus(100,10);
181  //CylindricalEta4D<double> & c = v1->Coordinates();
182  //c.SetValues(Px,pY,pZ,E);
183  setValues(*v1,p);
184  t1.Fill();
185  s += getMag2(*v1);
186  }
187
188  f1.Write();
189  timer.Stop();
190
191  double sav = std::sqrt(s/double(n));
192
193
194 #ifdef DEBUG
195  t1.Print();
196  std::cout << " Time for Writing " << file_name << "\t: " << timer.RealTime() << " " << timer.CpuTime() << std::endl;
197  int pr = std::cout.precision(18); std::cout << "Average : " << sav << std::endl; std::cout.precision(pr);
198 #endif
199
200  return sav;
201 }
202
203
204 //----------------------------------------------------------------
206 //----------------------------------------------------------------
207
208 template<class Vector>
209 double read(const std::string & file_name) {
210
212
213  std::string fname = file_name + ".root";
214
215  TFile f1(fname.c_str());
216  if (f1.IsZombie() ) {
217  std::cout << " Error opening file " << file_name << std::endl;
218  return -1;
219  }
220
221  //TFile f1("mathcoreVectorIO_D32.root");
222
223  // create tree
224  TTree *t1 = (TTree*)f1.Get("t1");
225
226  Vector *v1 = 0;
227
228  std::cout << "reading typeID : " << typeid(*v1).name() << std::endl;
229
231
232  timer.Start();
233  int n = (int) t1->GetEntries();
234  std::cout << " Tree Entries " << n << std::endl;
235  double s=0;
236  for (int i = 0; i < n; ++i) {
237  t1->GetEntry(i);
238  s += getMag2(*v1);
239  }
240
241
242  timer.Stop();
243
244  double sav = std::sqrt(s/double(n));
245
246 #ifdef DEBUG
247  std::cout << " Time for Reading " << file_name << "\t: " << timer.RealTime() << " " << timer.CpuTime() << std::endl;
248  int pr = std::cout.precision(18); std::cout << "Average : " << sav << std::endl; std::cout.precision(pr);
249 #endif
250
251  return sav;
252 }
253
254 template<class TrackType>
255 double writeTrack(int n, const std::string & file_name, int compress = 0) {
256
257  std::cout << "\n";
258  std::cout << " Test writing .." << file_name << " .....\n";
259  std::cout << "**************************************************\n";
260
261  TRandom3 R(111);
262
264
265  //std::cout << "writing a tree with " << vector_type << std::endl;
266
267  std::string fname = file_name + ".root";
268  TFile f1(fname.c_str(),"RECREATE","",compress);
269
270  // create tree
271  std::string tree_name="Tree with TrackD";
272  TTree t1("t1",tree_name.c_str());
273
274  TrackType *track = new TrackType();
275  std::cout << "typeID written : " << typeid(*track).name() << std::endl;
276
277
278  //t1.Branch("Vector branch",&track,16000,0);
279  // default split model
280  t1.Branch("Vector branch",&track);
281
282  timer.Start();
283  double s = 0;
286  typedef typename TrackType::VectorType V;
287  typedef typename TrackType::PointType P;
288
289  for (int i = 0; i < n; ++i) {
290  q.SetXYZT( R.Gaus(0,10),
291  R.Gaus(0,10),
292  R.Gaus(0,10),
293  R.Gaus(100,10) );
294  p.SetXYZ( q.X(), q.Y(), q.Z() );
295
296  track->Set( V(q), P(p) );
297
298  t1.Fill();
299  s += getMag2( *track );
300  }
301
302  f1.Write();
303  timer.Stop();
304
305  double sav = std::sqrt(s/double(n));
306
307
308 #ifdef DEBUG
309  t1.Print();
310  std::cout << " Time for Writing " << file_name << "\t: " << timer.RealTime() << " " << timer.CpuTime() << std::endl;
311  int pr = std::cout.precision(18); std::cout << "Average : " << sav << std::endl; std::cout.precision(pr);
312 #endif
313
314  return sav;
315
316 }
317
318
319
320 int testResult(double w1, double r1, const std::string & type) {
321
322  int iret = 0;
323  std::cout << w1 << " r " << r1 << std::endl;
324  // do like this to avoid nan's
325  if (!( fabs(w1-r1) < tol)) {
326  std::cout << "\nERROR: Differeces found when reading " << std::endl;
327  int pr = std::cout.precision(18); std::cout << w1 << " != " << r1 << std::endl; std::cout.precision(pr);
328  iret = -1;
329  }
330
331  std::cout << "\n*********************************************************************************************\n";
332  std::cout << "Test :\t " << type << "\t\t";
333  if (iret ==0)
334  std::cout << "OK" << std::endl;
335  else
336  std::cout << "FAILED" << std::endl;
337  std::cout << "********************************************************************************************\n\n";
338
339  return iret;
340 }
341
342
343
344 int testVectorIO(bool readOnly = false) {
345
346  int iret = 0;
347
348 // #ifdef __CINT__
351 // using namespace ROOT::Math;
352 // #endif
353
354 #ifdef USE_REFLEX
355
356 // #ifdef __CINT__
359 // #endif
360
361  std::cout << "Use Reflex dictionary " << std::endl;
362
365
366  ROOT::Cintex::Cintex::SetDebug(1);
367  ROOT::Cintex::Cintex::Enable();
368
369
370
374  if (iret |= 0) {
375  std::cerr <<"Failing to Load Reflex dictionaries " << std::endl;
376  return iret;
377  }
378
379 #endif // endif on using reflex
380
381
382  int nEvents = 100000;
383  //int nEvents = 100;
384
385  double w1, r1 = 0;
386  std::string fname;
387
388  testDummy<AVector4D>(nEvents);
389
390  fname = "lorentzvector";
392  w1 = 98.995527276968474;
393  fname += "_prev";
394  }
395  else
396  w1 = write<AVector4D>(nEvents,fname,vector4d_type);
397
399  iret |= testResult(w1,r1,vector4d_type);
400
401  fname = "displacementvector";
403  w1 = 17.3281570633214095;
404  fname += "_prev";
405  }
406  else
407  w1 = write<AVector3D>(nEvents,fname,vector3d_type);
408
410  iret |= testResult(w1,r1,vector3d_type);
411
412  fname = "positionvector";
414  fname += "_prev";
415  else
416  w1 = write<APoint3D>(nEvents,fname,point3d_type);
417
419  iret |= testResult(w1,r1,point3d_type);
420
421
422  // test TrackD
423  fname = "track";
426
427  //nEvents = 10000;
428
430  fname += "_prev";
431  }
432  else
433  w1 = writeTrack<TrackD>( nEvents,fname);
434
436  iret |= testResult(w1,r1,"TrackD");
437
438  // test TrackD32
439
440  fname = "trackD32";
443
445  fname += "_prev";
446  }
447  else
448  w1 = writeTrack<TrackD32>( nEvents,fname);
449
451  iret |= testResult(w1,r1,"TrackD32");
452
453
454  // test vector of tracks
455  fname = "vectrack";
456
457  nEvents = 10000;
458
460  fname += "_prev";
461  }
462  else
463  w1 = writeTrack<VecTrackD>( nEvents,fname);
464
466  iret |= testResult(w1,r1,"VecTrackD");
467
468  // test ClusterD (cotaining a vector of points)
469  fname = "cluster";
470
471  nEvents = 10000;
472
474  fname += "_prev";
475  }
476  else
477  w1 = writeTrack<ClusterD>( nEvents,fname);
478
480  iret |= testResult(w1,r1,"ClusterD");
481
482  return iret;
483 }
484
485 int main(int argc, char ** ) {
486
488  if (argc > 1) readOnly = true;
489
491  if (iret != 0)
492  std::cerr << "testVectorIO:\t FAILED ! " << std::endl;
493  else
494  std::cerr << "testVectorIO:\t OK ! " << std::endl;
495
496  return iret;
497
498 }
499
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TObject.cxx:778
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
Definition: LorentzVector.h:48
virtual void Print(Option_t *option="") const
Dump this text with its attributes.
Definition: TText.cxx:781
Random number generator class based on M.
Definition: TRandom3.h:27
double testDummy(int n)
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
DisplacementVector3D< CoordSystem, Tag > & SetXYZ(Scalar a, Scalar b, Scalar c)
set the values of the vector from the cartesian components (x,y,z) (if the vector is held in polar or...
double write(int n, const std::string &file_name, const std::string &vector_type, int compress=0)
writing
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format...
Definition: TFile.h:46
Class describing a generic position vector (point) in 3 dimensions.
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Definition: TTree.cxx:5321
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Definition: TSystem.cxx:1825
Double_t Mag2() const
Definition: TVector3.h:339
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
double writeTrack(int n, const std::string &file_name, int compress=0)
LorentzVector< CoordSystem > & SetXYZT(Scalar xx, Scalar yy, Scalar zz, Scalar tt)
set the values of the vector from the cartesian components (x,y,z,t) (if the vector is held in anothe...
TLatex * t1
Definition: textangle.C:20
double sqrt(double)
TStopwatch timer
Definition: pirndm.C:37
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
Change branch address, dealing with clone trees properly.
Definition: TTree.cxx:7873
void SetXYZ(Double_t x, Double_t y, Double_t z)
Definition: TVector3.h:227
#define DEFPOINT3D(TYPE)
TVector3 is a general three vector class, which can be used for the description of different vectors ...
Definition: TVector3.h:22
const int nEvents
Definition: testRooFit.cxx:42
Class describing a generic displacement vector in 3 dimensions.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const double tol
TLorentzVector is a general four-vector class, which can be used either for the description of positi...
R__EXTERN TSystem * gSystem
Definition: TSystem.h:539
void SetXYZT(Double_t x, Double_t y, Double_t z, Double_t t)
SVector< double, 2 > v
Definition: Dict.h:5
int testResult(double w1, double r1, const std::string &type)
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
static double P[]
Double_t Mag2() const
Bool_t IsZombie() const
Definition: TObject.h:122
int type
Definition: TGX11.cxx:120
virtual Long64_t GetEntries() const
Definition: TTree.h:381
#define DEFVECTOR4D(TYPE)
PositionVector3D< CoordSystem, Tag > & SetXYZ(Scalar a, Scalar b, Scalar c)
set the values of the vector from the cartesian components (x,y,z) (if the vector is held in polar or...
int main(int argc, char **)
#define DEFVECTOR3D(TYPE)
TF1 * f1
Definition: legend1.C:11
void setValues(ROOT::Math::DisplacementVector3D< U > &v, const double x[])
A TTree object has a header with a name and a title.
Definition: TTree.h:78
double getMag2(const Vector &v)
float * q
Definition: THbookFile.cxx:87
const Int_t n
Definition: legend1.C:16
TRandom3 R
a TMatrixD.
Definition: testIO.cxx:28
Stopwatch class.
Definition: TStopwatch.h:28