ROOT  6.06/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 
38 //#define READONLY
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 //----------------------------------------------------------------
205 /// reading
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 
230  t1->SetBranchAddress("Vector branch",&v1);
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__
349 // gSystem->Load("libMathCore");
350 // gSystem->Load("libPhysics");
351 // using namespace ROOT::Math;
352 // #endif
353 
354 #ifdef USE_REFLEX
355 
356 // #ifdef __CINT__
357 // gSystem->Load("libReflex");
358 // gSystem->Load("libCintex");
359 // #endif
360 
361  std::cout << "Use Reflex dictionary " << std::endl;
362 
363  gSystem->Load("libReflex");
364  gSystem->Load("libCintex");
365 
366  ROOT::Cintex::Cintex::SetDebug(1);
367  ROOT::Cintex::Cintex::Enable();
368 
369 
370 
371  //iret |= gSystem->Load("libSmatrixRflx");
372  //iret |= gSystem->Load("libMathAddRflx");
373  //iret |= gSystem->Load("libMathRflx");
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";
391  if (readOnly) {
392  w1 = 98.995527276968474;
393  fname += "_prev";
394  }
395  else
396  w1 = write<AVector4D>(nEvents,fname,vector4d_type);
397 
398  r1 = read<AVector4D>(fname);
399  iret |= testResult(w1,r1,vector4d_type);
400 
401  fname = "displacementvector";
402  if (readOnly) {
403  w1 = 17.3281570633214095;
404  fname += "_prev";
405  }
406  else
407  w1 = write<AVector3D>(nEvents,fname,vector3d_type);
408 
409  r1 = read<AVector3D>(fname);
410  iret |= testResult(w1,r1,vector3d_type);
411 
412  fname = "positionvector";
413  if (readOnly)
414  fname += "_prev";
415  else
416  w1 = write<APoint3D>(nEvents,fname,point3d_type);
417 
418  r1 = read<APoint3D>(fname);
419  iret |= testResult(w1,r1,point3d_type);
420 
421 
422  // test TrackD
423  fname = "track";
424  // load track dictionary
425  gSystem->Load("libTrackDict");
426 
427  //nEvents = 10000;
428 
429  if (readOnly) {
430  fname += "_prev";
431  }
432  else
433  w1 = writeTrack<TrackD>( nEvents,fname);
434 
435  r1 = read<TrackD>(fname);
436  iret |= testResult(w1,r1,"TrackD");
437 
438  // test TrackD32
439 
440  fname = "trackD32";
441  // load track dictionary
442  // gSystem->Load("libTrackDict");
443 
444  if (readOnly) {
445  fname += "_prev";
446  }
447  else
448  w1 = writeTrack<TrackD32>( nEvents,fname);
449 
450  r1 = read<TrackD32>(fname);
451  iret |= testResult(w1,r1,"TrackD32");
452 
453 
454  // test vector of tracks
455  fname = "vectrack";
456 
457  nEvents = 10000;
458 
459  if (readOnly) {
460  fname += "_prev";
461  }
462  else
463  w1 = writeTrack<VecTrackD>( nEvents,fname);
464 
465  r1 = read<VecTrackD>(fname);
466  iret |= testResult(w1,r1,"VecTrackD");
467 
468  // test ClusterD (cotaining a vector of points)
469  fname = "cluster";
470 
471  nEvents = 10000;
472 
473  if (readOnly) {
474  fname += "_prev";
475  }
476  else
477  w1 = writeTrack<ClusterD>( nEvents,fname);
478 
479  r1 = read<ClusterD>(fname);
480  iret |= testResult(w1,r1,"ClusterD");
481 
482  return iret;
483 }
484 
485 int main(int argc, char ** ) {
486 
487  bool readOnly = false;
488  if (argc > 1) readOnly = true;
489 
490  int iret = testVectorIO(readOnly);
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:823
double read(const std::string &file_name)
reading
Class describing a generic LorentzVector in the 4D space-time, using the specified coordinate system ...
Definition: LorentzVector.h:54
Random number generator class based on M.
Definition: TRandom3.h:29
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:108
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
const Double_t * v1
Definition: TArcBall.cxx:33
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:45
Class describing a generic position vector (point) in 3 dimensions.
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all branches of entry and return total number of bytes read.
Definition: TTree.cxx:5163
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1818
Bool_t IsZombie() const
Definition: TObject.h:141
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
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:75
Double_t x[n]
Definition: legend1.C:17
virtual void Print(Option_t *option="") const
Dump this text with its attributes.
Definition: TText.cxx:776
virtual Int_t SetBranchAddress(const char *bname, void *add, TBranch **ptr=0)
Change branch address, dealing with clone trees properly.
Definition: TTree.cxx:7529
void SetXYZ(Double_t x, Double_t y, Double_t z)
Definition: TVector3.h:231
Double_t Mag2() const
#define DEFPOINT3D(TYPE)
int testVectorIO(bool readOnly=false)
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
R__EXTERN TSystem * gSystem
Definition: TSystem.h:549
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
Double_t Mag2() const
Definition: TVector3.h:298
MyComplex< T > P(MyComplex< T > z, T c_real, T c_imag)
[MyComplex]
Definition: mandel.cpp:155
int type
Definition: TGX11.cxx:120
#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[])
virtual Long64_t GetEntries() const
Definition: TTree.h:382
A TTree object has a header with a name and a title.
Definition: TTree.h:94
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:30