// -----------------------------------------------------------------------------
void runRecoFOPI_display_matching_Felix(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 = 100,
				      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/rev170/",
				      TString SGE_IDENTIFIER=""
				      ) 
{
  
// ========================================================================
// 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;
    
  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+".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());
  TpcAlignmentManager::getInstance();
  TpcDigiMapper::getInstance()->init(par);
  TpcDigiMapper::getInstance()->forceManualDriftVel(forceDriftVel);
  double vdrift=TpcDigiMapper::getInstance()->getGas()->VDrift();
  
  //################################################################################
  //################################################################################ 
  // 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 = 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 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;
  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);
  fRun->SetField(fMagField);

  //initialize GF field
  GFFieldManager::getInstance()->init(new PndFieldAdaptor(fRun->GetField()));
  GFMaterialEffects::getInstance()->init(new GFTGeoMaterialInterface());
  
  //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+"*").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;

  //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(DP);
  readFopi->SetCdcTrackPersistence(true);
  readFopi->SetRpcTrackPersistence(true);
  //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);
  
  //Conversion from fopi2root RpcTrack objects to GF strip hits
  FopiRpcInitTask* rpcInit = new FopiRpcInitTask();
  rpcInit->SetPersistence(true);
  rpcInit->SetUseConstantOffset(true);
  rpcInit->SetConstantErrorZ(1.5);
  rpcInit->SetConstantErrorPhi(0.2); //basically turning phi info off
  fRun->AddTask(rpcInit);

  //Conversion from fopi2root BarTrack objects to GF strip hits
  FopiBarInitTask* barInit = new FopiBarInitTask();
  barInit->SetPersistence(true);
  barInit->SetConstantErrorZ(10.);
  //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(600,14);
  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->useGeane(); // uses RKTrackrep and GeaneTrackrep
  trackInit->SetSmoothing(smoothing);
  trackInit->SetWeightedPlaneConstruction(true);
  trackInit->SetCutCov(true);
  fRun->AddTask(trackInit);

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


  KalmanTask* kalman = new KalmanTask();
  //kalman->SetUseDAF();
  kalman->SetPersistence(true);
  //kalman->SetTpcClusterBranchName("TpcCluster");
  kalman->SetTpcClusterBranchName("TpcSPHit");
  kalman->SetTrackBranchName("TpcTrackPreFit");
  kalman->SetOutBranchName("TpcTrackPostFit");
  kalman->SetNumIterations(2); // number of fitting iterations (back and forth)
  fRun->AddTask(kalman);
  
  //build GFTracks for CDC hits only 
  FopiTrackInitTask* trackCdcInit=new FopiTrackInitTask();
  trackCdcInit->SetPersistence();  //affects prefit tracks
  trackCdcInit->SetAllRecoHitPersistence(true);
  //trackCdcInit->SetTrackPersistence(true);
  trackCdcInit->SetHitBranchName("CdcHit");
  trackCdcInit->SetTrackBranchName("CdcTrack");
  trackCdcInit->SetOutBranchName("CdcTrackPreFit");
  trackCdcInit->SetTrackPersistence(DP);
  //trackCdcInit->SetPiKaonThreshold(0.32);
  //trackCdcInit->SetMinCDCHits(5);
  trackCdcInit->SetSmoothing(smoothing);
  
  //trackCdcInit->SetUsePseudoSPHits();
  trackCdcInit->SetWeightedPlaneConstruction(true);
  trackCdcInit->SetCutCov(true);

  //do we want to include RPC/BAR hits in the fit?
  trackCdcInit->SetMatchRpc(true);
  trackCdcInit->SetMatchBar(true);
  
  //trackCdcInit->SetRejectRpcMirrorHits(true);

  //in cm:
  trackCdcInit->SetConstantPointErrors(4.5e-2, 4.5e-2, 10.);
  
  fRun->AddTask(trackCdcInit);
  
  KalmanTask* KalmanCdc = new KalmanTask(); 
  //KalmanCdc->SetUseDAF();
  KalmanCdc->SetPersistence(true);
  KalmanCdc->SetTrackBranchName("CdcTrackPreFit");
  KalmanCdc->SetOutBranchName("CdcTrackPostFit");
  KalmanCdc->SetFopiCDCBranchName("CdcSpacepoint");
  KalmanCdc->SetNumIterations(2); // number of fitting iterations
  fRun->AddTask(KalmanCdc);
   
  //match TPC and CDC tracks (no GFTracks!)
  TpcCdcMatchingTask* tpcCdcMatch = new TpcCdcMatchingTask();
  tpcCdcMatch->SetPersistence(kTRUE);
  tpcCdcMatch->SetTpcTrackBranchName("TpcTrackPostFit");
  tpcCdcMatch->SetCdcTrackBranchName("CdcTrack");
  tpcCdcMatch->SetCdcHitBranchName("CdcHit");
  tpcCdcMatch->SetWriteFailedMatches(); 
  fRun->AddTask(tpcCdcMatch);  
  
   
  FopiTrackMergingTask* merge = new FopiTrackMergingTask();
  merge->SetPersistence();
  merge->SetTrackPersistence(false); //prefit tracks
  merge->SetTpcTrackBranchName("TpcTrackPostFit");
  merge->AddTupleBranchName("TpcCdcMatchingPair");
  merge->SetCdcGFTrackBranchName("CdcTrackPostFit");
  fRun->AddTask(merge);

  KalmanTask* KalmanComb = new KalmanTask();
  //KalmanComb->SetUseDAF();
  KalmanComb->SetPersistence(kTRUE);
  KalmanComb->SetTpcClusterBranchName("TpcSPHit");
  KalmanComb->SetFopiCDCBranchName("CdcSpacepoint");
  KalmanComb->SetTrackBranchName("FopiTrackPreFit");
  KalmanComb->SetOutBranchName("FopiTrackPostFit");
  KalmanComb->SetNumIterations(3); // number of fitting iterations
  fRun->AddTask(KalmanComb);

  // 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);

  FopiEventDisplayTask* dt = new FopiEventDisplayTask();
  dt->SetGFTrackBranchName("CdcTrackPostFit");
  fRun->AddTask(dt);
  
  
  // -----   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:"<<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();

  GenfitDisplay::getInstance()->setOptions("DHMPTA");
  GenfitDisplay::getInstance()->open();
}