mathcoreVectorIO.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_math
/// \notebook -nodraw
/// Example of  I/O of a GenVector Lorentz Vectors in a Tree and comparison with legacy TLorentzVector.
///
/// A ROOT tree is written and read in both using either a XYZTVector or a TLorentzVector.
///
///  To execute the macro type in:
///
/// ~~~{.cpp}
/// root[0] .x  mathcoreVectorIO.C
/// ~~~
///
/// \macro_output
/// \macro_code
///
/// \author Lorenzo Moneta
 
 
#include "TRandom2.h"
#include "TStopwatch.h"
#include "TSystem.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1D.h"
#include "TCanvas.h"
 
#include <iostream>
 
#include "TLorentzVector.h"
 
#include "Math/Vector4D.h"
 
using namespace ROOT::Math;
 
void write(int n) {
   TRandom2 R;
   TStopwatch timer;
 
 
   R.SetSeed(1);
   timer.Start();
   double s = 0;
   for (int i = 0; i < n; ++i) {
      s  += R.Gaus(0,10);
      s  += R.Gaus(0,10);
      s  += R.Gaus(0,10);
      s  += R.Gaus(100,10);
   }
 
   timer.Stop();
   std::cout << s/double(n) << std::endl;
   std::cout << " Time for Random gen " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;
 
 
   TFile f1("mathcoreVectorIO_1.root","RECREATE");
 
   // create tree
   TTree t1("t1","Tree with new LorentzVector");
 
   XYZTVector *v1 = new XYZTVector();
   t1.Branch("LV branch","ROOT::Math::XYZTVector",&v1);
 
   R.SetSeed(1);
   timer.Start();
   for (int i = 0; i < n; ++i) {
      double Px = R.Gaus(0,10);
      double Py = R.Gaus(0,10);
      double Pz = R.Gaus(0,10);
      double E  = R.Gaus(100,10);
      v1->SetCoordinates(Px,Py,Pz,E);
      t1.Fill();
   }
 
   f1.Write();
   timer.Stop();
   std::cout << " Time for new Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;
 
   t1.Print();
 
   // create tree with old LV
 
   TFile f2("mathcoreVectorIO_2.root","RECREATE");
   TTree t2("t2","Tree with TLorentzVector");
 
   TLorentzVector * v2 = new TLorentzVector();
   TLorentzVector::Class()->IgnoreTObjectStreamer();
   TVector3::Class()->IgnoreTObjectStreamer();
 
   t2.Branch("TLV branch","TLorentzVector",&v2,16000,2);
 
   R.SetSeed(1);
   timer.Start();
   for (int i = 0; i < n; ++i) {
      double Px = R.Gaus(0,10);
      double Py = R.Gaus(0,10);
      double Pz = R.Gaus(0,10);
      double E  = R.Gaus(100,10);
      v2->SetPxPyPzE(Px,Py,Pz,E);
      t2.Fill();
   }
 
   f2.Write();
   timer.Stop();
   std::cout << " Time for old Vector " << timer.RealTime() << "  " << timer.CpuTime() << endl;
   t2.Print();
}
 
 
void read() {
 
   TRandom R;
   TStopwatch timer;
 
   TFile f1("mathcoreVectorIO_1.root");
 
   // create tree
   TTree *t1 = (TTree*)f1.Get("t1");
 
   XYZTVector *v1 = 0;
   t1->SetBranchAddress("LV branch",&v1);
 
   timer.Start();
   int n = (int) t1->GetEntries();
   std::cout << " Tree Entries " << n << std::endl;
   double etot=0;
   for (int i = 0; i < n; ++i) {
      t1->GetEntry(i);
      etot += v1->Px();
      etot += v1->Py();
      etot += v1->Pz();
      etot += v1->E();
   }
   timer.Stop();
   std::cout << " Time for new Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;
 
   std::cout << " TOT average : n = " << n << "\t " << etot/double(n) << endl;
 
   // create tree with old LV
   TFile f2("mathcoreVectorIO_2.root");
   TTree *t2 = (TTree*)f2.Get("t2");
 
   TLorentzVector * v2 = 0;
   t2->SetBranchAddress("TLV branch",&v2);
 
   timer.Start();
   n = (int) t2->GetEntries();
   std::cout << " Tree Entries " << n << std::endl;
   etot = 0;
   for (int i = 0; i < n; ++i) {
      t2->GetEntry(i);
      etot  += v2->Px();
      etot  += v2->Py();
      etot  += v2->Pz();
      etot  += v2->E();
   }
 
   timer.Stop();
   std::cout << " Time for old Vector " << timer.RealTime() << "  " << timer.CpuTime() << endl;
   std::cout << " TOT average:\t" << etot/double(n) << endl;
}
 
void mathcoreVectorIO() {
   
   int nEvents = 100000;
   write(nEvents);
   read();
}