#ifndef __CINT__
#include <iostream>
//#include <sstream>

#include <TString.h>
#include <TFile.h>
#include <TROOT.h>
#include <TStopwatch.h>
#include <TVector3.h>
#include <TGeoManager.h>
#include <TSystem.h>
#include <TPython.h>
#include <TGeant3.h>

#include "FairRun.h"
#include "FairRunAna.h"
#include "FairField.h"
#include "FairRuntimeDb.h"
#include "FairParAsciiFileIo.h"
#include "FairParGenericSet.h"
#include "FairParRootFileIo.h"

#include "PndFieldAdaptor.h"

#include "GFTGeoMaterialInterface.h"
#include "GFMaterialEffects.h"
#include "GFFieldManager.h"

#include "TpcGas.h"
#include "FopiDataReaderTask.h"
#include "TpcDataReaderTask.h"
#include "FopiRpcInitTask.h"
#include "FopiTrackInitTask.h"
#include "FopiBarInitTask.h"
#include "FopiTrackMergingTask.h"
#include "TpcTrackInitTask.h"
#include "FopiRpcBarMergingTask.h"
#include "TpcCdcMatchingTask.h"
#include "CdcCorrectionTask.h"
#include "TpcAlignmentTask.h"
#include "TpcRiemannTrackingTask.h"
#include "TpcClusterFinderTask.h"
#include "TpcDigiPar.h"
#include "TpcDigiAmpCorrectionTask.h"
#include "TpcPSATask.h"
#include "TpcDigiMapper.h"
#include "TpcAlignmentManager.h"
#include "FOPIField.h"

#include "KalmanTask.h"
#endif
// -----------------------------------------------------------------------------
void runRecoFOPI_batch_matching_t0(int runNumber,
				   TString jobname,
				   TString outpath="$RECOOUT",
				   bool forceDriftVel = false,
				   unsigned int deepPersist = 0,  //persist Digis, clusters
				   unsigned int smoothing = 0,
				   unsigned int gaincorr = 1,
				   unsigned int nEvents = 0,
				   unsigned int nEvStart = 0,
				   TString decodedInPath="/nfs/hicran/data/tpc/fopi/2011/decoded/",
				   //TString fopiRecoPath="/nfs/nas/data/panda/tpc/SIM/fcusanno/FopiFiles/", 
				   TString fopiRecoPath="/nfs/mds/data/tpc/fopi/2011/fopi_reco_sverre/rev176/",
				   TString SGE_IDENTIFIER="",
				   TString fopiRecoSuffix="",
				   double t0corr=0,
				   bool flightTimeCorr = true,
				   double wirePlaneShift = 0.18
				   ) 
{
  
  // ========================================================================
  // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
  // Load basic libraries in rootlogon
  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  rootlogon();

  Int_t iVerbose = 1;

  TStopwatch timer;
  timer.Start();

  unsigned int nFopiEvents = 0;
  
  gROOT->SetBatch(); //to ensure running in batch

  TString basedir = gSystem->Getenv("VMCWORKDIR");
  FairRunAna* fRun = new FairRunAna();
  
  if(smoothing)
    jobname += "_smoothed.reco.root";
  else
    jobname += ".reco.root";
  
  TString outFile = outpath+"/";
  outFile += jobname;

  TString zCorrFile = outpath+"/analysis/zShifts.root";
  //std::cout<<zCorrFile.Data()<<std::endl;
  
    
    
  fRun->SetOutputFile(outFile);

  std::stringstream ss;
  ss<<runNumber;
  TString runNr(ss.str().c_str());  

  //--- Parameter database ---
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo();
  //  TString parFile = "/tpc/FOPI/par/tpc.run"+runNr+"_davidDrift.par";
  TString parFile = "/tpc/FOPI/par/tpc.run"+runNr+".par";
  parFile = basedir+parFile;
  std::cout << "parFile: " << parFile << std::endl;
  parInput1->open(parFile.Data(),"in");
  rtdb->setFirstInput(parInput1);

  TpcDigiPar* par = (TpcDigiPar*) rtdb->getContainer("TpcDigiPar");
  par->setInputVersion(fRun->GetRunId(),1);
  par->setChanged(kTRUE);
  par->init();

  //root file parameters
  TString copy = jobname;
  TString parfile = outFile.ReplaceAll("reco.root","param.root");
  FairParRootFileIo* parOut=new FairParRootFileIo(kTRUE);
  parOut->open(parfile.Data());
  rtdb->setOutput(parOut);
 
   
  //TpcAlignmentManager::init(par->getAlignmentFile());
  TString alignFile = basedir;
  alignFile += "/tpc/FOPI/par/alignment_carbon_manual_phi_2015_03_25.txt";
  TpcAlignmentManager::init(alignFile);  
  TpcAlignmentManager::getInstance();
  TpcDigiMapper::getInstance()->init(par);
  TpcDigiMapper::getInstance()->forceManualDriftVel(forceDriftVel);
  double vdrift=TpcDigiMapper::getInstance()->getGas()->VDrift(0,0);
  std::cout<<vdrift<<std::endl;
  
  //################################################################################
  //################################################################################ 
  // build a run specific GEOMETRY
  // run the FopiGeom.C in a separate root session
  //
  // taroff: distance between target and center of magnet
  // targetpos: distance between real target position and nominal FOPI target
  // position
  // In experiment S339 the distance between the nominal FOPI target position
  // and the center of the magnet is 41 cm.
  //################################################################################
  Double_t taroff = -61.;  //FIXED from 65->61 on 30.7.2014 F.B
  Double_t targetpos = taroff+41.;
  // possible materials: copper (default), lead, carbon
  TString targetmat = TString("carbon");
  TString alignmentFileName =alignFile;//par->getAlignmentFile();
  TString geoFile;

  //IF we run parallel, each job has to write his OWN geometry file to avoid conflicts :
  if(SGE_IDENTIFIER.Length()>0)
    geoFile.Form("/tmp/FopiGeom_s339%5.1f_%s_%s.root",
		 targetpos,targetmat.Data(),SGE_IDENTIFIER.Data());
  else 
    geoFile.Form("$FOPI2ROOT/fopigeometry/FopiGeom_s339%5.1f_%s.root",targetpos,targetmat.Data());
  
  TString teststring = gSystem->GetFromPipe(TString::Format("if [ -f %s/macro/geometry/FopiGeoms339.root ]; then echo yes; else echo no; fi", outpath.Data()));
  bool rebuildGeom=true;
  if (teststring=="yes"){
    std::cout<<"yes"<<std::endl;
    rebuildGeom=false;
    geoFile=outpath+"/macro/geometry/FopiGeoms339.root";
  }
    
  TString cmd;
  cmd.Form(".! root -b -q '$FOPI2ROOT/fopigeometry/FopiGeom.C(%f,\"%s\",\"%s\",\"%s\",kFALSE)'",targetpos,targetmat.Data(),alignmentFileName.Data(),geoFile.Data());
  std::cout<<std::endl<<std::endl;
  std::cout<<"################################################################################"<<std::endl;
  std::cout<<"################################################################################"<<std::endl;
  if(rebuildGeom){
    gROOT->ProcessLine(cmd);
  }
  std::cout<<"################################################################################"<<std::endl;
  std::cout<<"################################################################################"<<std::endl;
  std::cout<<std::endl<<std::endl;
  //################################################################################


  // add geometry file to fRun
  // work around a bug in FAIRsoft... I hate these clowns
  unsigned int ntries = 0;
  while(fRun->GetGeoFile() == NULL) {
    std::cout<<"\n   trying to set geometry file... try #"<<ntries<<std::endl;
    fRun->SetGeomFile(geoFile);
    ntries++;
  }

   
  // FOPI magnetic field
  //For pion beam data, FOPIField(0.616) has to be used according to CDC/RPC analysis!
  //if the constructor is called with no argument, 0.6T is used
  FOPIField *fMagField = new FOPIField(0.616); 
  fMagField->SetTargetOffset(taroff);

  //const field:
  //PndConstField *fMagField=new PndConstField();
  //fMagField->SetField(0., 0., 6.16 ); // values are in kG
  //fMagField->SetFieldRegion(-500, 500,-500, 500, -2000, 2000); // values are in cm

  fRun->SetField(fMagField);

  //initialize GF field
  GFFieldManager::getInstance()->init(new PndFieldAdaptor(fRun->GetField()));
  GFMaterialEffects::getInstance()->init(new GFTGeoMaterialInterface());
  //GFMaterialEffects::getInstance()->setNoEffects();
  
  //parse input files:
  TString pythonScript=basedir+"/macro/tpc/FOPI/parseFiles.py";
  TPython::LoadMacro(pythonScript);
  PyFileParser fp;
  if(decodedInPath.Last('/') != decodedInPath.Length()-1)
    decodedInPath+="/";
  fp.setFiles(TString(decodedInPath+"*"+runNr+"*").Data());
  int nFiles = (int)fp.getNum();
  std::cout<<"FOUND "<<nFiles<<" decoded files"<<std::endl;
  TChain* decodedTestChain = new TChain("tpcEvent");
  for(int i=0; i<nFiles; i++) {
    decodedTestChain->AddFile((char*)fp.getFile(i));
    std::cout<<"Add file "<<(char*)fp.getFile(i)<<" to decoded chain"<<std::endl;
  }
  if(nFiles>1)
    Fatal("MACRO:","More than 1 decoded file found in input path");

  if(fopiRecoPath.Last('/') != fopiRecoPath.Length()-1)
    fopiRecoPath+="/";
  fp.setFiles(TString(fopiRecoPath+"FOPI_*"+runNr+fopiRecoSuffix+"*.root").Data());
  nFiles = (int)fp.getNum();
  TChain* recoTestChain = new TChain("FOPITree");
  for(int i=0; i<nFiles; i++) {
    recoTestChain->AddFile((char*)fp.getFile(i));
    std::cout<<"Add file "<<(char*)fp.getFile(i)<<" to reco chain"<<std::endl;
  }
  if(nFiles>1)
    Fatal("MACRO:","More than 1 reco file found in input path");

 
  if(nEvents==0)
    nEvents=decodedTestChain->GetEntries();
  std::cout<<"Found "<<decodedTestChain->GetEntries()<<" events in TPC input data file"<<std::endl;
  std::cout<<"Found "<<recoTestChain->GetEntries()<<" events in fopi input data file"<<std::endl;


  TString pythonScript2=basedir+"/macro/tpc/FOPI/getCdcPars.py";
  TPython::LoadMacro(pythonScript2);
  PyGetCdcPars cdcParLoader;
  TString bla = basedir+"/tpc/FOPI/par/cdc_par_AE.par";
  cdcParLoader.loadInputFile(bla.Data());

  double cdcT0=cdcParLoader.getT0(runNumber);
  double cdcVDrift=cdcParLoader.getVDrift(runNumber);
  double cdcLorentzAngle=cdcParLoader.getLorentzAngle(runNumber);

  std::cout<<cdcLorentzAngle<<std::endl;
  std::cout<<cdcT0<<std::endl;
  std::cout<<cdcVDrift<<std::endl;


  //delete decodedTestChain;
  //delete recoTestChain;
  
  bool DP = false;
  if(deepPersist == 1)
    DP = true;
    
  //--------------------SET UP TASKS ------------------------------

    
  TpcDataReaderTask* read = new TpcDataReaderTask();
  read->SetPersistence(kFALSE);
  read->SetDatafile(decodedTestChain->GetFile()->GetName());
  read->SetSampleBranchName("TpcSample");
  read->SetStartEvent(nEvStart);
  //read->SetCutSmallPad();
  //read->SetMinSamples(1000);
  fRun->AddTask(read);
    
  FopiDataReaderTask* readFopi = new FopiDataReaderTask();
  readFopi->SetPersistence(false);
  readFopi->SetCdcHitPersistence(false);
  readFopi->SetCdcTrackPersistence(false);
  readFopi->SetRpcTrackPersistence(false);
  //readFopi->SetBarrelOn();
  //not (re)implemented yet, but it is only a double-check
  //readFopi->SetRunNr(runNr);
  readFopi->SetInputFile(recoTestChain->GetFile()->GetName());
  //readFopi->SetNevent(nFopiEvents);
  readFopi->SetDisableHeliTracks(true);
  readFopi->SetBarrelOn();
  fRun->AddTask(readFopi);

  //CORRECT FOR CDC SHIFTS
  CdcCorrectionTask* corTask=new CdcCorrectionTask();
  corTask->SetCdcVdr(cdcVDrift);
  corTask->SetCdcLorentz(cdcLorentzAngle);
  corTask->SetCdcT0(cdcT0);
  corTask->SetEnableFlightTimeCorr(flightTimeCorr);
  corTask->SetCorrection(wirePlaneShift); //default is 0.2
  //corTask->SetCorrectionFile("testShifts.root");
  corTask->SetCorrectionFile(zCorrFile);
  corTask->SetCdcT0Corr(t0corr);
  fRun->AddTask(corTask);


  //Conversion from fopi2root RpcTrack objects to GF strip hits
  FopiRpcInitTask* rpcInit = new FopiRpcInitTask();
  rpcInit->SetPersistence(true);
  rpcInit->SetUseConstantOffset(true);
  rpcInit->SetConstantErrorZ(2.5);
  rpcInit->SetConstantErrorPhi(0.2);
  rpcInit->SetOutBranchName("RpcPixHit");
  fRun->AddTask(rpcInit);

  //Conversion from fopi2root BarTrack objects to GF strip hits
  FopiBarInitTask* barInit = new FopiBarInitTask();
  barInit->SetPersistence(true);
  barInit->SetConstantErrorZ(8.);
  barInit->SetConstantErrorPhi(3./TMath::Sqrt(12.));
  barInit->SetOutBranchName("BarPixHit");
  //barInit->SetConstantErrorPhi(100.); //basically turning phi info off
  fRun->AddTask(barInit);
  
  
  TpcPSATask* tpsa = new  TpcPSATask();
  tpsa->SetPersistence(DP);
  tpsa->SetSamplePersistence(DP);
  tpsa->SetSampleBranchName("TpcSample"); // Input of PSA
  //tpsa->UseNoiseMaxSampleCut(decodedTestChain->GetFile()->GetName(),7.);
  fRun->AddTask(tpsa);

  if(gaincorr==1){
    TString gainfile = basedir;
    TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask();
    tpcCalib->SetVerbose(kFALSE);
    tpcCalib->SetPersistence(DP);
    fRun->AddTask(tpcCalib);
  }
  
  //bool SimpleClustering = true;
  TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask();
  tpcCF->SetPersistence(DP); // keep Digis refs in clusters
  tpcCF->SetDigiPersistence(DP);
  //tpcCF->SetMode(2); // individual timeslice
  tpcCF->timeslice(6); //in samples
  //tpcCF->SetSingleDigiClusterAmpCut(10);
  //tpcCF->SetClusterAmpCut(5); // cut on mean digi amplitude
  tpcCF->SetErrorPars(-1,10.);
  tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple
  fRun->AddTask(tpcCF);

  //find TpcRiemannTracks in the TPC alone
  TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask();
  tpcSPR->SetClusterBranchName("TpcCluster");
  tpcSPR->SetSortingParameters(true,3,42.78);
  tpcSPR->SetPersistence(DP);
  //tpcSPR->SetVerbose(1);
  fRun->AddTask(tpcSPR);

  //build GFTracks from TpcRiemannTracks
  TpcTrackInitTask* trackInit=new TpcTrackInitTask();
  trackInit->SetPersistence(DP);
  trackInit->SetAllRecoHitPersistence(true);
  trackInit->SetTrackOutBranchName("TpcTrackPreFit");
  trackInit->SetClusterBranchName("TpcCluster");
  trackInit->SetRecoHitOutBranchName("TpcSPHit");
  TVector3 IP(0,0,42.78);
  trackInit->SetInteractionPoint(IP);
  trackInit->SetPDG(211);
  trackInit->SetSmoothing(smoothing);
  trackInit->SetWeightedPlaneConstruction(true);
  trackInit->SetCutCov(true);
  trackInit->SetVerbose(0);
  fRun->AddTask(trackInit);
  

  TpcAlignmentTask* align = new TpcAlignmentTask();
  align->SetRecoHitBranchName("TpcSPHit");
  fRun->AddTask(align);

  //TPC standalone fitting. I guess we don't really need to do that...
  KalmanTask* kalman = new KalmanTask();
  //kalman->SetUseDAF();
  kalman->SetPersistence(true);
  //kalman->SetTpcClusterBranchName("TpcCluster");
  kalman->SetTpcClusterBranchName("TpcSPHit");
  kalman->SetTrackBranchName("TpcTrackPreFit");
  kalman->SetOutBranchName("TpcTrackPostFit");
  kalman->SetNumIterations(1); // number of fitting iterations (back and forth)
  fRun->AddTask(kalman);
  
  //build GFTracks for CDC hits only. NO BARREL/RPC
  FopiTrackInitTask* trackCdcInit=new FopiTrackInitTask();
  trackCdcInit->SetPersistence(true);                   //affects hits
  trackCdcInit->SetAllRecoHitPersistence(true);
  trackCdcInit->SetExtrapToZ(true);
  //trackCdcInit->SetTrackPersistence(true);
  trackCdcInit->SetHitBranchName("CdcHit");
  trackCdcInit->SetTrackBranchName("CdcTrack");
  trackCdcInit->SetOutBranchName("CdcTrackPreTempFit");  //still missing BAR/RPC
  trackCdcInit->SetRecoHitOutBranchName("CdcSpacepoint");
  trackCdcInit->SetTrackPersistence(DP);
  //trackCdcInit->SetPiKaonThreshold(0.32);
  //trackCdcInit->SetMinCDCHits(5);
  trackCdcInit->SetSmoothing(smoothing);
  
  //trackCdcInit->SetUsePseudoSPHits();
  trackCdcInit->SetWeightedPlaneConstruction(true);
  trackCdcInit->SetCutCov(true);
  
  //in cm:
  trackCdcInit->SetConstantPointErrors(3.5e-2, 3.5e-2, 10.);
  
  fRun->AddTask(trackCdcInit);
  
  
  //refit with added BAR / RPC hits
  KalmanTask* KalmanCdc1 = new KalmanTask();
  //KalmanCdc1->SetUseDAF();
  KalmanCdc1->SetPersistence(false);
  KalmanCdc1->SetTpcClusterBranchName("TpcSPHit");
  KalmanCdc1->SetFopiCDCBranchName("CdcSpacepoint");
  KalmanCdc1->SetTrackBranchName("CdcTrackPreTempFit");
  KalmanCdc1->SetOutBranchName("CdcTrackPostTempFit");
  KalmanCdc1->SetNumIterations(1); // number of fitting iterations
  fRun->AddTask(KalmanCdc1);
  
   
  //match TPC and CDC tracks (no GFTracks!)
  TpcCdcMatchingTask* tpcCdcMatch = new TpcCdcMatchingTask();
  tpcCdcMatch->SetPersistence(false);
  tpcCdcMatch->SetTpcTrackBranchName("TpcTrackPostFit");
  tpcCdcMatch->SetCdcTrackBranchName("CdcTrack");
  tpcCdcMatch->SetCdcHitBranchName("CdcHit");
  //  tpcCdcMatch->SetWriteFailedMatches(); 
  fRun->AddTask(tpcCdcMatch);  

  
  
  //merge CDC and TPC
  FopiTrackMergingTask* merge = new FopiTrackMergingTask();
  merge->SetPersistence();
  merge->SetTrackPersistence(false); //prefit tracks
  merge->SetTpcTrackBranchName("TpcTrackPostFit");
  merge->AddTupleBranchName("TpcCdcMatchingPair");
  //merge->SetTrackOutBranchName("FopiTrackPreFit");     //TEST
  merge->SetTrackOutBranchName("FopiTrackPreTempFit");
  merge->SetCdcGFTrackBranchName("CdcTrackPreTempFit");
  fRun->AddTask(merge);

  //obtain fits for TPC-CDC tracks. Needed for RPC/BAR matching decision
  KalmanTask* KalmanComb = new KalmanTask();
  KalmanComb->SetUseDAF();
  KalmanComb->SetPersistence(false);
  KalmanComb->SetTpcClusterBranchName("TpcSPHit");
  KalmanComb->SetFopiCDCBranchName("CdcSpacepoint");
  KalmanComb->SetTrackBranchName("FopiTrackPreTempFit");
  KalmanComb->SetOutBranchName("FopiTrackPostTempFit");
  //KalmanComb->SetNumIterations(1); // number of fitting iterations
  fRun->AddTask(KalmanComb);

  //RPC / Barrel merging goes here. RPC residual check is performed 
  //creates both TPC-CDC + BAR/RPC and CDC + BAR/RPC
  FopiRpcBarMergingTask* addBarRpc = new FopiRpcBarMergingTask();
  addBarRpc->SetPersistence(true);
  addBarRpc->SetExtrapToZ(true);
  addBarRpc->SetSmoothing(smoothing);
  addBarRpc->SetTrackBranchName("FopiTrackPreTempFit");
  addBarRpc->SetCdcGFTrackBranchName("CdcTrackPostTempFit");
  addBarRpc->SetTrackOutBranchName("FopiTrackPreFit");
  addBarRpc->SetCdcTrackOutBranchName("CdcTrackPreFit");
  addBarRpc->SetDissolveFOPIRpcMatching(false); //DISSOLVE FOPI's CDC-RPC matching if residual too large
  //addBarRpc->SetExtrapResMeanZ(2.82e-1);
  //addBarRpc->SetExtrapResSigZ(1.40);
  //addBarRpc->SetExtrapWindow(2.);  
  fRun->AddTask(addBarRpc);
  

  //refit with added BAR / RPC hits
  KalmanTask* KalmanCdc = new KalmanTask();
  KalmanCdc->SetUseDAF();
  KalmanCdc->SetPersistence(kTRUE);
  KalmanCdc->SetTpcClusterBranchName("TpcSPHit");
  KalmanCdc->SetFopiCDCBranchName("CdcSpacepoint");
  KalmanCdc->SetTrackBranchName("CdcTrackPreFit");
  KalmanCdc->SetOutBranchName("CdcTrackPostFit");
  KalmanCdc->SetNumIterations(2); // number of fitting iterations
  fRun->AddTask(KalmanCdc);

  //refit with added BAR / RPC hits
  KalmanTask* KalmanTpcRpc = new KalmanTask();
  KalmanTpcRpc->SetUseDAF();
  KalmanTpcRpc->SetPersistence(kTRUE);
  KalmanTpcRpc->SetTpcClusterBranchName("TpcSPHit");
  KalmanTpcRpc->SetFopiCDCBranchName("CdcSpacepoint");
  KalmanTpcRpc->SetTrackBranchName("TpcRpcTrackPreFit");
  KalmanTpcRpc->SetOutBranchName("TpcRpcTrackPostFit");
  KalmanTpcRpc->SetNumIterations(2); // number of fitting iterations
  fRun->AddTask(KalmanTpcRpc);
  
  
  
  // Chi2CleanupTask* cleanupCDC =  new Chi2CleanupTask();
  // cleanupCDC->SetTrackBranchName("CdcTrackPostFit");
  // //cleanupCDC->AddDiffCutWindow(-1.e6,-10.);
  // cleanupCDC->AddDiffCutWindow(10.,1.e6);
  // cleanupCDC->SetHitWindow(0.93,1.);
  // //fRun->AddTask(cleanupCDC);
  

  // AddTrackRepTask* addReps = new AddTrackRepTask();
  // addReps->SetTrackBranchName("FopiTrackPreFit");
  // addReps->AddHypothesis("proton");
  // addReps->AddHypothesis("kaon");
  // //  fRun->AddTask(addReps);
  

  // //refit with added BAR / RPC hits
  // KalmanTask* KalmanComb2 = new KalmanTask();
  // KalmanComb2->SetUseDAF();
  // KalmanComb2->SetPersistence(kFALSE);
  // KalmanComb2->SetTpcClusterBranchName("TpcSPHit");
  // KalmanComb2->SetFopiCDCBranchName("CdcSpacepoint");
  // KalmanComb2->SetTrackBranchName("FopiTrackPreFit");
  // KalmanComb2->SetOutBranchName("FopiTrackPostFit");
  // //KalmanComb2->SetNumIterations(2); // number of fitting iterations
  // KalmanComb2->AddIgnorePDG(-2212);  //ignore antiprotons in the fit
  // //  fRun->AddTask(KalmanComb2);

  

  // Chi2CleanupTask* cleanupCOMB =  new Chi2CleanupTask();
  // cleanupCOMB->SetTrackBranchName("FopiTrackPostFit");
  // //cleanupCOMB->AddDiffCutWindow(-1.e6,-10.);
  // cleanupCOMB->AddDiffCutWindow(10.,1.e6);
  // cleanupCOMB->SetHitWindow(0.95,1.);
  // //fRun->AddTask(cleanupCOMB);


  // // dE/dx task
  // TpcdEdxTask* dEdx = new TpcdEdxTask();
  // dEdx->SetTrackBranchName("FopiTrackPostFit");
  // dEdx->SetClusterBranchName("TpcSPHit");
  // //dEdx->SetIdealdEdx();
  // dEdx->SetPersistence();
  // //dEdx->SetVerbose(1);
  // dEdx->SetUseFirstDigiPos();
  // dEdx->SetIgnoreEdgeDigis(5.8,14.5);
  // dEdx->SetDXgrid(0.5);
  // //  fRun->AddTask(dEdx);

  //  //the actual calculator:
  // TpcCdc2DMatchedResCalc* cdcRefRes = new TpcCdc2DMatchedResCalc();
  // cdcRefRes->addBranchName("FopiTrackTuples");
  // cdcRefRes->addBranchName("TpcSPHit");
  // cdcRefRes->setTpcTrackBranchName("TpcTrackPostFit");
  // cdcRefRes->setCdcTrackBranchName("CdcTrack");
  // //the task:
  // TpcRefTrackResidualTask* res = new TpcRefTrackResidualTask();
  // res->SetResCalculator(cdcRefRes);
  // res->SetOutBranch("TpcCdcRefResiduals");
  // res->SetPersistence();   
  // //fRun->AddTask(res);

  // TpcQATask* qa = new TpcQATask();
  // qa->SetTpcCdcTrackBranchName("FopiTrackPostFit");
  // qa->SetTpcDecodedFile(decodedTestChain->GetFile()->GetName());
  // //  fRun->AddTask(qa);
  
  
  // -----   Intialise and run   --------------------------------------------
  fRun->Init();
  std::cout<<"\n****** gGeoManager->Print() ***************** \n"<<std::endl;
  gGeoManager->Print();
  std::cout<<"\n*********************************************\n\n"<<std::endl;
  
  rtdb->saveOutput();
  rtdb->print();

  std::cout<<"start:"<<nEvStart<<" end:"<<nEvStart+nEvents<<std::endl;
  fRun->Run(nEvStart,nEvStart+nEvents); 

  // close pdf-file
  // TO BE USED WITH TpcHelDisplayTask TASK
  //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);
  
  std::cout<<"OutputFile: "<<outFile<<std::endl;
  
  //testFile.Close();
}