#include "TVector3.h"
#include <iostream>
using namespace std;

void print(TVector3 a) {
  cout << a.X() << " " << a.Y() << " " << a.Z() << endl;
}


void checkspu(){

  // ========================================================================
  // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
  Int_t iVerbose = 1;

  // Number of events to process
  Int_t nEvents = 5; // process all events in input file

  // Output file
  TString outFile = "testgenfit.root";
  
  TString inFile  = "../testrun.root";
  TString parFile = "../testparams.root"; 
  
  // ----  Load libraries   -------------------------------------------------
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  rootlogon();
  // ------------------------------------------------------------------------
  
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------

  // -----   Reco run   -------------------------------------------
  FairRunAna *fRun= new FairRunAna();
  fRun->SetInputFile(inFile);
  fRun->SetOutputFile(outFile);

  // ----- Prepare GEANE --------------------------------------------
  // this will load Geant3 and execute setup macros to initialize geometry:
  FairGeane *Geane = new FairGeane();
  fRun->AddTask(Geane);
  // ------------------------------------------------------------------------

  // -----  Parameter database   --------------------------------------------
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParRootFileIo* parInput1 = new FairParRootFileIo();
  parInput1->open(parFile.Data());
  rtdb->setFirstInput(parInput1);

  //-------------------------------------------------------------------------
  // -----   Intialise and run   --------------------------------------------
  fRun->Init();
  rtdb->print();
  Geane->SetField(fRun->GetField());
  fRun->Run(0, nEvents);

  //----------------------------
  // start parameters
  TVector3 StartPos(0.,0.,0.);
  TVector3 StartPosErr(0.01,0.01,0.01);
  TVector3 StartMom(1.,-0.5,0.);
  StartMom.SetMag(1.);
  TVector3 StartMomErr(0.01,0.01,0.01);
  // particle
  Double_t  Charge= -1;
  Int_t PDGCode= 13;
  // start plane
  TVector3 so(0.,0.,0.);
  TVector3 sdj(0.,1.,0.);
  TVector3 sdk(0.,0.,1.);
  // end plane
  TVector3 eo(10.,0.,0.);
  TVector3 edj(0.,0.,1.);
  TVector3 edk(0.,1.,0.);
  
  // ---------------------- GEANE -------------------------
  FairTrackParP *start = new FairTrackParP(StartPos, StartMom, 
					   StartPosErr, StartMomErr, 
					   Charge, so, sdj, sdk);

  FairGeanePro *pro = new FairGeanePro();
  pro->PropagateFromPlane(sdj, sdk);
  pro->PropagateToPlane(eo, edj, edk);

  FairTrackParP *end = new FairTrackParP();
  pro->Propagate(start, end, PDGCode);

  cout << "print each momentum 3 times! " << endl;
  cout << endl;

  cout << "from GEANE alone --------------------------" << endl;
  cout << "start spu = " << start->GetSPU() 
       << "\nstart momentum " << endl;
  print(start->GetMomentum());
  print(start->GetMomentum());
  print(start->GetMomentum());
  cout << endl;
  cout << "end spu = " << end->GetSPU()
       << "\nend momentum " << endl;
  print(end->GetMomentum());
  print(end->GetMomentum());
  print(end->GetMomentum());
  cout << endl;
  
  // ---------------------- GENFIT -------------------------
  DetPlane splane(so, sdj, sdk);

  GeaneTrackRep *srep = new GeaneTrackRep(pro,
					  splane,
					  StartMom,
					  StartPosErr,
					  StartMomErr,
					  Charge,
					  PDGCode);
  DetPlane eplane(eo, edj, edk);

 cout << "from GeaneTrackRep ------------------------" << endl;
  cout << "GeaneTrackRep::getMom(ref plane)" << endl;
  print(((TVector3) srep->getMom(splane)));
  print(((TVector3) srep->getMom(splane)));
  print(((TVector3) srep->getMom(splane)));
  cout << endl;
  cout << "GeaneTrackRep::getMom(end plane)" << endl;
  print(((TVector3) srep->getMom(eplane)));
  print(((TVector3) srep->getMom(eplane)));
  print(((TVector3) srep->getMom(eplane)));
  
  cout << endl;
  cout << "======================================================" << endl;
  cout << "the problem is in SPU, which changes from the first call 
to the second and to the third, when it should stay the 
same until we decide to update the track representation.
In fact, let' s call once again the getMom(end plane)
and print the spu:" << endl;
  cout << "1st time starting spu is " << srep->getSPU() << " (WRONG) and momentum: " << endl;
  print(((TVector3) srep->getMom(eplane)));
  cout << "2nd time starting spu is " << srep->getSPU() << " (RIGHT) and momentum: " << endl;
  print(((TVector3) srep->getMom(eplane)));
  cout << "3rd time starting spu is " << srep->getSPU() << " (WRONG) and momentum: " << endl;
  print(((TVector3) srep->getMom(eplane)));
  cout << "and so on, alternating wrong and right values. :(" << endl;

  // ------------------------------------------------------------------------
  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished succesfully." << endl;
  cout << "Output file is "    << outFile << endl;
  cout << "Parameter file is " << parFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
  cout << endl;
  // ------------------------------------------------------------------------
}