/// \file
/// \ingroup tutorial_math
/// \notebook -nodraw
/// Macro illustrating  I/O with Lorentz Vectors of floats
/// The dictionary for LorentzVector of float is not in the libMathCore, therefore
/// is generated when parsed the file with CINT.
///
///  To run  this macro you must do
///
/// ~~~{.cpp}
/// root[0] .L  mathcoreVectorFloatIO.C+
/// root[1] runIt();
/// ~~~
///
/// \macro_code
///
/// \author Lorenzo Moneta

#include "TRandom.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"


// Now the dictionary contains the vector's with float types
// No need to force dictionary generation
// You need to run ACLIC with old ROOT version
// and uncomment these lines below 
// #ifdef __MAKECINT__
// #pragma link C++ class ROOT::Math::PxPyPzE4D<float>+;
// #pragma link C++ class ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<float> >+;
// #pragma link C++ typedef ROOT::Math::XYZTVectorF;
// #endif

using namespace ROOT::Math;

void write(int n) {

   TRandom R;
   TStopwatch timer;

   TFile f1("mathcoreVectorIO_F.root","RECREATE");

   // create tree
   TTree t1("t1","Tree with new Float LorentzVector");

   XYZTVectorF *v1 = new XYZTVectorF();
   t1.Branch("LV branch","ROOT::Math::XYZTVectorF",&v1);

   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 Float Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;
   t1.Print();
}

void read() {

   TRandom R;
   TStopwatch timer;

   TFile f1("mathcoreVectorIO_F.root");

   // create tree
   TTree *t1 = (TTree*)f1.Get("t1");

   XYZTVectorF *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->E();
   }

   timer.Stop();
   std::cout << " Time for new Float Vector " << timer.RealTime() << "  " << timer.CpuTime() << std::endl;
   std::cout << " E average" << n<< "  " << etot << "  " << etot/double(n) << endl;
}

void runIt() {

#if defined(__CINT__) && !defined(__MAKECINT__)
   gSystem->Load("libMathCore");
   gSystem->Load("libPhysics");
   using namespace ROOT::Math ;

   cout << "This tutorial can run only using ACliC, you must run it by doing: " << endl;
   cout << "\t  .L tutorials/math/mathcoreVectorFloatIO.C+" << endl;
   cout << "\t  runIt()" << endl;
#endif
   int nEvents = 100000;
   write(nEvents);
   read();
}

void mathcoreVectorFloatIO() {
#if defined(__CINT__) && !defined(__MAKECINT__)
   gSystem->Load("libMathCore");
   gSystem->Load("libPhysics");
   using namespace ROOT::Math ;

   cout << "This tutorial can run only using ACliC, you must run it by doing: " << endl;
   cout << "\t  .L tutorials/math/mathcoreVectorFloatIO.C+" << endl;
   cout << "\t  runIt()" << endl;

#endif
}