// Macro for running Cbm  with Geant3  or Geant4 (M. Al-Turany , D. Bertini)
// Modified 22/06/2005 D.Bertini
{
  gDebug=0;
  // Load basic libraries
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
  basiclibs();

  // Load this example libraries
  gSystem->Load("libGeoBase");
  gSystem->Load("libParBase");
  gSystem->Load("libBase");
  gSystem->Load("libField");
  gSystem->Load("libMCStack");
  gSystem->Load("libPassive");
  gSystem->Load("libTrkBase");  
  gSystem->Load("libGeane");
  gSystem->Load("libStt");
  //gSystem->Load("libPlane");
  gSystem->Load("libGen");
  //gSystem->Load("libGeaneEx");
  gSystem->Load("libgenfit");
  gSystem->Load("librecotasks");

  CbmRunAna *fRun = new CbmRunAna();

  TString base="data/test1";
  TString infile=base+".mc.root";
  TString outfile=base+".g.root";
  fRun->SetInputFile(infile);
  fRun->SetOutputFile(outfile);

  // setup the run for  use with Geane:
  // Geane object will call scripts to load Geant3
  // the object will not be used afterwards!
  CbmGeane *Geane = new CbmGeane(infile);

 // -----  Parameter database   --------------------------------------------

  CbmRuntimeDb* rtdb = fRun->GetRuntimeDb();
  CbmParRootFileIo* parInput1 = new CbmParRootFileIo();

  TString dbfile=base+".param.root";
  parInput1->open(dbfile.Data()); 
  rtdb->setFirstInput(parInput1);

  // ------------------------------------------------------------------------

  fRun->Init();

  rtdb->print();
 
  // Set the field(if any) to Geane
  Geane->SetField(fRun->GetField());
     
  
  // Setup a trackrep which uses geane
  CbmGeanePro* gePro=new CbmGeanePro();
  
  CbmTrackParP dummy;

  TVector3 pos(0.,0.,0.);
  TVector3 mom(0.1,0.0001,1.);
  TVector3 poserr=TVector3(0.1,0.1,0.1);
  TVector3 momerr=0.1*mom;

  DetPlane initplane(TVector3(pos),TVector3(1.,0.,0.),TVector3(0.,1.,0.));
  double q=1;
  
  // The track representation:
  AbsTrackRep* rep=new GeaneTrackRep(gePro,pos,mom,poserr,momerr,q,initplane);
  rep.Print();

  // now exrapolate the parameters to another plane
  DetPlane plane(TVector3(0,0,20),TVector3(1,0,0),TVector3(0,1,0));

  TMatrixT<double> newpar(5,1);
  TMatrixT<double> newcov(5,5);
  TMatrixT<double> jac(5,5); // will be empty!

  rep->extrapolate(plane,newpar,newcov,jac);
  
  newpar.Print();
  newcov.Print();

  // now create a simple recohit:
  // the DemoRecoHit uses the DetPlane (0,0,z)(1,0,0)(0,1,0) -> u=x v=y

  AbsRecoHit* hit=new DemoRecoHit(10,8,0,      // coordinates (u,v,z)
				  0.01,0.01);   // errors (du,dv)

  //hit.Print();

  TMatrixT<double> res=hit.residualVector(rep,rep->getState());
  res->Print();
  // Note that the residual calcualtion already assumes
  // the track to be extrapolated to the right frame!

  
  AbsRecoHit* hit2=new DemoRecoHit(14,12,5,      // coordinates (u,v,z)
				  0.01,0.01);   // errors (du,dv)

  rep->extrapolate(hit2->getDetPlane(rep),newpar,newcov,jac);

  TMatrixT<double> res2=hit.residualVector(rep,newpar);
  std::cout<<"Residual at z=5cm"<<std::endl;
  res2->Print();
  

}