#include <"GFException">
// -----------------------------------------------------------------------------
void runTPC_Helitron() {

  // start timer
	TStopwatch timer;
  timer.Start();

  unsigned int nEvents = 0;
  unsigned int nFopiEvents = 0;
  
  // Load basic libraries in rootlogon
  TString OutPath="/home/pbuehler/physics/projects/fopi/tpc/results/";
  TString outfname=OutPath+"axibixi.root";
	TString basedir = gSystem->Getenv("VMCWORKDIR");
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");rootlogon();
	
	// define input files
	TString tpcfname=
     "/home/pbuehler/physics/projects/fopi/tpc/data/tpc/decoded/run_3280.lmd_decoded.root";
	TString fopifname=
	  "/home/pbuehler/physics/projects/fopi/tpc/data/fopi_new/FOPI_0.root";

  //extract number of entries in external data tree
  TFile testFile(tpcfname);
  if(nEvents==0) nEvents = ((TTree*)testFile.Get("tpcEvent"))->GetEntries();
  TFile testFile2(fopifname);
  if(nFopiEvents==0) nFopiEvents = ((TTree*)testFile2.Get("FOPITree"))->GetEntries();
  
  std::cout<<"Found "<<nEvents<<" events in  input data file"<<std::endl;
  std::cout<<"Found "<<nFopiEvents<<" events in fopi input data file"<<std::endl;
  
   // create run manager
  FairRunAna* fRun = new FairRunAna();

  // Create alignment manager
  TpcAlignmentManager * align = TpcAlignmentManager::getInstance(
	  "/home/pbuehler/physics/projects/fopi/tpc/fopiroot/dev/tpc/FOPI/par/june11_alignment.txt");

  // define output file
	fRun->SetOutputFile(outfname);

  // have to define digi-parameters for TPC
	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParAsciiFileIo* tpcdigiinp = new FairParAsciiFileIo();
  TString tpcDigiFile = basedir;
  tpcDigiFile += "/tpc/FOPI/par/tpc.run3280.par";
  tpcdigiinp->open(tpcDigiFile.Data(),"in");
  rtdb->setFirstInput(tpcdigiinp);
  rtdb->Print();

  // FOPI geometry
  TString geoFile = basedir;
  geoFile="/home/pbuehler/physics/projects/fopi/tpc/FOPI-dst2root/FopiGeom/FopiGeom.root";
  fRun->SetGeomFile(geoFile);
  
	// FOPI magnetic field
  FOPIField *fMagField = new FOPIField();
  fMagField->SetTargetOffset(-65.);
	fRun->SetField(fMagField);

  // Geane task
  FairGeane *Geane = new FairGeane();
  fRun->AddTask(Geane);
    
  TpcDataReaderTask* tpcread = new TpcDataReaderTask();
  tpcread->SetPersistence(1);
  tpcread->SetDatafile(tpcfname);
  tpcread->SetStartEvent(0);
  fRun->AddTask(tpcread);
 
  FopiDataReaderTask* fopiread = new FopiDataReaderTask();
  fopiread->SetPersistence(1);
  fopiread->SetInputFile(fopifname);
  fopiread->SetDebug(kFALSE);
  fRun->AddTask(fopiread);
 
  // TPC clustering 
  TpcPSATask* tpsa = new  TpcPSATask();
  tpsa->SetSampleBranchName("TpcSample");
  tpsa->SetDigiBranchName("TpcDigi");
  tpsa->SetPersistence();
  fRun->AddTask(tpsa);

  TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask();
  tpcCF->SetDigiBranchName("TpcDigi");
  tpcCF->SetClusterBranchName("TpcCluster");
  tpcCF->SetDigiPersistence(1); // keep Digis refs in clusters
  tpcCF->SetPersistence(1); // keep Clusters
  tpcCF->timeslice(6); //in samples
  tpcCF->SetSingleDigiClusterAmpCut(20);
  tpcCF->SetClusterAmpCut(9.1); // cut on mean digi amplitude
  tpcCF->SetErrorPars(600,400);
  tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple
  fRun->AddTask(tpcCF);
 
  // TPC tracking
  TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask();
  tpcSPR->SetClusterBranchName("TpcCluster");
  tpcSPR->SetPersistence(0);
  tpcSPR->SetVerbose(kTRUE);
  fRun->AddTask(tpcSPR);

  // build GFTracks from TpcRiemannTracks
  TpcTrackInitTask* trackInit=new TpcTrackInitTask();
  trackInit->SetClusterBranchName("TpcCluster");
  trackInit->SetRiemannBranchName("RiemannTrack");
  trackInit->SetOutBranchName("TpcTrackPreFit");
  trackInit->SetPersistence(0);
  trackInit->SetVerbose(kFALSE);
  //trackInit->SetMCPid(); // use ideal particle identification
  trackInit->SetPDG(211);
  trackInit->useGeane(); // uses RKTrackrep and GeaneTrackrep
  trackInit->SetSmoothing(true);
  fRun->AddTask(trackInit);
  
  KalmanTask* kalman =new KalmanTask();
  kalman->SetTrackBranchName("TpcTrackPreFit");
  kalman->SetClusterBranchName("TpcCluster");
  kalman->SetOutBranchName("TpcTrackPostFit");
  kalman->SetPersistence(0);
  //kalman->SetClusterBranchName("TpcCluster_cut");
  kalman->SetNumIterations(3); // number of fitting iterations (back and forth)
  kalman->SetVerbose(kFALSE);
  fRun->AddTask(kalman);

  // match TPC and Helitron tracks
	TpcHelMatchingTask* tpcHelMatch = new TpcHelMatchingTask();
  tpcHelMatch->SetPersistence(kTRUE);
  tpcHelMatch->SetVerbose(kTRUE);
  tpcHelMatch->setNHeliHits(1);
  tpcHelMatch->setWithPlawaHit(0);
  tpcHelMatch->SetPosErr(0.1,0.1,0.1);
  tpcHelMatch->SetMomErr(0.2,0.2,0.2);
  tpcHelMatch->SetHHitRes(0.1,1.,0.1);
  tpcHelMatch->SetdrMax(10.0);
  tpcHelMatch->SetResScale(10.0);
  tpcHelMatch->SetPHitRes(0.1,0.1,0.1);
	fRun->AddTask(tpcHelMatch);

  // fit matched TPC and Helitron tracks
	TpcHelTrackFittingTask* tpcHelTrackFit = new TpcHelTrackFittingTask();
  tpcHelTrackFit->SetPersistence(kTRUE);
  tpcHelTrackFit->SetVerbose(kTRUE);
  tpcHelTrackFit->SetDefRep(0);
	fRun->AddTask(tpcHelTrackFit);

  // display tpc and Helitron/Plawa tracks together
	TpcHelDisplayTask* tpcHelDisp = new TpcHelDisplayTask();
  tpcHelDisp->SetOutPath(OutPath);
  tpcHelDisp->SetDisplayType(1);
  tpcHelDisp->SetVerbose(kFALSE);
  //                     tt,hf,ht, v, t, c, h, p,th     
	tpcHelDisp->SetToShow( 1, 1, 1, 1, 1, 1, 1, 1, 1);
	tpcHelDisp->SetToShow( 1, 0, 1, 1, 0, 1, 0, 1, 1);
	fRun->AddTask(tpcHelDisp);

  // open pdf file
	TCanvas * c1 = new TCanvas();
  c1->Print(OutPath+"TpcHel_output.pdf[");

  // loop over events
  fRun->Init();

  // Initialize the Digimapper:
  TpcDigiPar* tpcpar = FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar");
  TpcDigiMapper::getInstance()->init(tpcpar);
  TpcAlignmentManager::init(tpcpar->getAlignmentFile());
  // rtdb->print();

  fRun->Run(0,1);
  
  // close pdf file
  c1->Print(OutPath+"TpcHel_output.pdf]");

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime(); 
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
  
  exit();
}

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