void cascade_rec(TString subDir = "")
{
  // Macro created 20/09/2006 by S.Spataro
  // It loads a digi file and performs tracking

  // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
  Int_t iVerbose = 0; // just forget about it, for the moment

	// Number of events to process
  Int_t nEvents = 0;  // if 0 all the vents will be processed

  // Parameter file
  TString parFile = subDir + "simparams.root"; // at the moment you do not need it

  // Digitisation file (ascii)
  TString digiFile = "all.par";

  // Output file
  TString outFile = subDir + "reco_complete.root";

  // which propagator should be used if any?
  Bool_t useGeane = kFALSE;
  Bool_t useGeane2 = kTRUE;

  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
    // ------------------------------------------------------------------------

  // -----   Reconstruction run   -------------------------------------------
  FairRunAna *fRun= new FairRunAna();
  fRun->SetInputFile(subDir + "sim_complete.root");
  fRun->AddFriend(subDir + "digi_complete.root");
  fRun->SetOutputFile(outFile);
  fRun->SetWriteRunInfoFile(kFALSE);
  fRun->SetUseFairLinks(kTRUE);
  FairGeane *Geane = new FairGeane();
  fRun->AddTask(Geane);

  // -----  Parameter database   --------------------------------------------
  TString emcDigiFile = gSystem->Getenv("VMCWORKDIR");
  emcDigiFile += "/macro/params/";
  emcDigiFile += digiFile;

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

  FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
  parIo1->open(emcDigiFile.Data(),"in");

  rtdb->setFirstInput(parInput1);
  rtdb->setSecondInput(parIo1);


  // -----------------------------------------------------------------------
  // -----   TRACKING   ----------------------------------------------------
  // -----------------------------------------------------------------------

  TString branchTrackBaseName = "IdealTrack";
  TString branchKalmanBaseName = "IdealGenTrack";
  TString branchTrackBaseNameFwd = "FtsIdealTrack";
  TString branchKalmanBaseNameFwd = "FtsIdealGenTrack";

  std::vector<TString> branchTrackExtensionName;
  // branchTrackExtensionName.push_back(TString(""));
  branchTrackExtensionName.push_back(TString("Filtered"));


  // -----   ideal track finder   ------------------------------------------
  PndMCIdealTrackFinderNewLinks* trackIdeal = new PndMCIdealTrackFinderNewLinks();
  trackIdeal->AddBranchName("MVDHitsPixel");
  trackIdeal->AddBranchName("MVDHitsStrip");
  trackIdeal->AddBranchName("STTHit");
  trackIdeal->AddBranchName("GEMHit");
  trackIdeal->AddBranchName("FTSHit");  // filter out these tracks later
  trackIdeal->SetRelativeMomentumSmearing(0.02);
  trackIdeal->SetVertexSmearing(0.05, 0.05, 0.05);
  fRun->AddTask(trackIdeal);

  PndMCIdealTrackFinderNewLinks* trackIdealFwd = new PndMCIdealTrackFinderNewLinks();
  trackIdealFwd->SetOutputBranchName("FtsIdealTrack");
  trackIdealFwd->AddBranchName("MVDHitsPixel");
  trackIdealFwd->AddBranchName("MVDHitsStrip");
  trackIdealFwd->AddBranchName("GEMHit");
  trackIdealFwd->AddBranchName("FTSHit");
  trackIdealFwd->SetRelativeMomentumSmearing(0.02);
  trackIdealFwd->SetVertexSmearing(0.05, 0.05, 0.05);
  fRun->AddTask(trackIdealFwd);


  // -----   Count the hits for track candidates   --------------------------
  PndTrackCandInfoTask * trackCandInfo = new PndTrackCandInfoTask();
  trackCandInfo->SetTrackInputBranch(branchTrackBaseName + "Cand");
  trackCandInfo->SetTrackOutputBranch(branchTrackBaseName + "CandInfo");
  fRun->AddTask(trackCandInfo);

  PndTrackCandInfoTask * trackCandInfoFwd = new PndTrackCandInfoTask();
  trackCandInfoFwd->SetTrackInputBranch(branchTrackBaseNameFwd + "Cand");
  trackCandInfoFwd->SetTrackOutputBranch(branchTrackBaseNameFwd + "CandInfo");
  fRun->AddTask(trackCandInfoFwd);


  // -----   Track filtering   ----------------------------------------------
  for ( std::vector<TString>::iterator it = branchTrackExtensionName.begin();
        it != branchTrackExtensionName.end();
        ++it ) {
    if (  *it != "" ) {
      PndTrackFilterTask* trackFilter = new PndTrackFilterTask();
      trackFilter->SetTrackInBranchName(branchTrackBaseName);
      trackFilter->SetFilteredTrackOutBranchName(branchTrackBaseName + *it);
      trackFilter->SetTrackCandInfoInBranchName(branchTrackBaseName + "CandInfo");
      trackFilter->SetFilteredTrackCandInfoOutBranchName(branchTrackBaseName + "CandInfo" + *it);
      // trackFilter->SetVerbose();
      trackFilter->AddMinHitsMvdFilter(4);
      trackFilter->AddMinHitsSttFilter(4);
      trackFilter->AddMinHitsGemFilter(4);
      trackFilter->AddMaxHitsFtsFilter(0);  // no FTS hits allowed here
      fRun->AddTask(trackFilter);

      PndTrackFilterTask* trackFilterFwd = new PndTrackFilterTask();
      trackFilterFwd->SetTrackInBranchName(branchTrackBaseNameFwd);
      trackFilterFwd->SetFilteredTrackOutBranchName(branchTrackBaseNameFwd + *it);
      trackFilterFwd->SetTrackCandInfoInBranchName(branchTrackBaseNameFwd + "CandInfo");
      trackFilterFwd->SetFilteredTrackCandInfoOutBranchName(branchTrackBaseNameFwd + "CandInfo" + *it);
      // trackFilterFwd->SetVerbose();
      trackFilterFwd->AddMinHitsFilter(4,0,0,1,0);  // if MVD hits, require >0 FTS hits as well
      trackFilterFwd->AddMinHitsFilter(0,0,4,1,0);  // if GEM hits, require >0 FTS hits as well
      // trackFilterFwd->AddMinHitsGemFilter(4);
      trackFilterFwd->AddMinHitsFtsFilter(4);
      fRun->AddTask(trackFilterFwd);
    }

    // track association and filtering (target spectrometer)
    PndMCTrackAssociator* trackMC = new PndMCTrackAssociator();
    trackMC->SetTrackInBranchName(branchTrackBaseName + *it);
    trackMC->SetTrackOutBranchName(branchTrackBaseName + *it + "ID");
    fRun->AddTask(trackMC);

    if ( useGeane ) {
      // central part
      PndRecoKalmanTask* recoKalman = new PndRecoKalmanTask();
      recoKalman->SetTrackInBranchName(branchTrackBaseName + *it);
      recoKalman->SetTrackInIDBranchName(branchTrackBaseName + *it + "ID");
      recoKalman->SetTrackOutBranchName(branchKalmanBaseName + *it);
      recoKalman->SetBusyCut(50); // CHECK to be tuned
      //recoKalman->SetIdealHyp(kTRUE);
      //recoKalman->SetNumIterations(3);
      // recoKalman->SetPropagateToIP(kFALSE);
      fRun->AddTask(recoKalman);

      PndMCTrackAssociator* trackMC2 = new PndMCTrackAssociator();
      trackMC2->SetTrackInBranchName(branchKalmanBaseNameFwd + *it);
      trackMC2->SetTrackOutBranchName(branchKalmanBaseNameFwd + *it + "ID");
      fRun->AddTask(trackMC2);

      // forward part
      PndRecoKalmanTask* recoKalmanFwd = new PndRecoKalmanTask();
      recoKalmanFwd->SetTrackInBranchName(branchTrackBaseNameFwd + *it);
      recoKalmanFwd->SetTrackInIDBranchName(branchTrackBaseNameFwd + *it + "ID");
      recoKalmanFwd->SetTrackOutBranchName(branchKalmanBaseNameFwd + *it);
      recoKalmanFwd->SetBusyCut(50); // CHECK to be tuned
      //recoKalmanFwd->SetIdealHyp(kTRUE);
      //recoKalmanFwd->SetNumIterations(3);
      // recoKalmanFwd->SetPropagateToIP(kFALSE);
      fRun->AddTask(recoKalmanFwd);

      PndMCTrackAssociator* trackMC2Fwd = new PndMCTrackAssociator();
      trackMC2Fwd->SetTrackInBranchName(branchKalmanBaseNameFwd + *it);
      trackMC2Fwd->SetTrackOutBranchName(branchKalmanBaseNameFwd + *it + "ID");
      fRun->AddTask(trackMC2Fwd);
    }

    if ( useGeane2 ) {
      // central part
      PndRecoKalmanTask2* recoKalman2 = new PndRecoKalmanTask2();
      recoKalman2->SetTrackInBranchName(branchTrackBaseName + *it);
      recoKalman2->SetTrackInIDBranchName(branchTrackBaseName + *it + "ID");
      recoKalman2->SetTrackOutBranchName(branchKalmanBaseName + "2" + *it);
      recoKalman2->SetBusyCut(50); // CHECK to be tuned
      recoKalman2->SetIdealHyp(kTRUE);
      recoKalman2->SetNumIterations(5);
      // recoKalman2->SetPropagateToIP(kFALSE);
      fRun->AddTask(recoKalman2);

      PndMCTrackAssociator* trackMC3 = new PndMCTrackAssociator();
      trackMC3->SetTrackInBranchName(branchKalmanBaseName + "2" + *it);
      trackMC3->SetTrackOutBranchName(branchKalmanBaseName + "2" + *it + "ID");
      fRun->AddTask(trackMC3);

      // forward part
      PndRecoKalmanTask2* recoKalman2Fwd = new PndRecoKalmanTask2();
      recoKalman2Fwd->SetTrackInBranchName(branchTrackBaseNameFwd + *it);
      recoKalman2Fwd->SetTrackInIDBranchName(branchTrackBaseNameFwd + *it + "ID");
      recoKalman2Fwd->SetTrackOutBranchName(branchKalmanBaseNameFwd + "2" + *it);
      recoKalman2Fwd->SetBusyCut(50); // CHECK to be tuned
      recoKalman2Fwd->SetIdealHyp(kTRUE);
      recoKalman2Fwd->SetNumIterations(5);
      // recoKalman2Fwd->SetPropagateToIP(kFALSE);
      fRun->AddTask(recoKalman2Fwd);

      PndMCTrackAssociator* trackMC3Fwd = new PndMCTrackAssociator();
      trackMC3Fwd->SetTrackInBranchName(branchKalmanBaseNameFwd + "2" + *it);
      trackMC3Fwd->SetTrackOutBranchName(branchKalmanBaseNameFwd + "2" + *it + "ID");
      fRun->AddTask(trackMC3Fwd);
    }
  }


  // -----   Intialise and run   --------------------------------------------
  PndEmcMapper::Init(1);
  cout << "fRun->Init()" << endl;
  fRun->Init();

  timer.Start();
  fRun->Run(0,nEvents);
  // ------------------------------------------------------------------------


  // -----   Finish   -------------------------------------------------------
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  cout << endl << endl;
  cout << "Macro finished successfully." << 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;
  // ------------------------------------------------------------------------
  cout << " Test passed" << endl;
  cout << " All ok " << endl;
  //exit(0);
}