void runReco_mb2(unsigned int run, TString outpath, bool clcorr=false, bool reclust=false,
		unsigned int nEvents = 0, unsigned int nEvStart = 0,
		unsigned int bfield=1,unsigned int smoothing = 1,
	        int outnum=-1, unsigned int gaincorr=1) {
// ========================================================================
// Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
  Int_t iVerbose = 1;
  Int_t clustvis = 0;

  TStopwatch timer;
  timer.Start();
  
  // Load basic libraries in rootlogon
  //  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
  //  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
  //  cout<<"$VMCWORKDIR/macro/tpc/FOPI/mberger/startClustVis.C"<<endl;
  //gROOT->LoadMacro("./macro/tpc/FOPI/mberger/startClustVis.C");
  //  rootlogon.C
  //  FairLogger::GetLogger()->SetLogToScreen(true);
  //  FairLogger::GetLogger()->SetLogScreenLevel("logINFO");
  TString filename="";
  if (clustvis==0)
    {
      //      TString filename="/nfs/hicran/data/tpc/fopi/2011/mberger/decoded/";
      filename+="decoded/";
      filename+="runC_";
      filename+=run;
      filename+=".lmd_decoded.root";
    }
  
  TString basedir = gSystem->Getenv("VMCWORKDIR");
  FairRunAna* fRun = new FairRunAna();
  //FairRunSim* fSim = new FairRunSim();

  TString jobdir = outpath; 
  std::string jobname(filename.Data());
  int last = jobname.rfind("/");
  if(last>0)
    jobname = jobname.substr(last+1,jobname.size()+1);  

  //create output file -----------------------------------------------  
  TString outName(jobname);
  if(outName.Contains("repaired."))
    outName.ReplaceAll(".lmd_decoded_repaired","");
  if(outName.Contains("decoded."))
    outName.ReplaceAll(".lmd_decoded","");
  if(smoothing)
    outName.ReplaceAll(".root", "_smoothed.root");
  if(clcorr)
    outName.ReplaceAll(".root","_clcorr.root");
  if (reclust)
    outName.ReplaceAll(".root","_recl.root");
  if(outnum>-1)
    outName.ReplaceAll(".root",Form("_%03i.root",outnum));

  std::cout<<outName<<"\n";
  //  return;

  outName.ReplaceAll(".root",".reco.root");
  TString outFile = outpath+"/";
  outFile += outName; 
    
  // TString PROutFile = outFile;
  // PROutFile.ReplaceAll(".reco.root",".patternReco.root");
  // TFile test(PROutFile, "recreate");
  // if(!test.IsZombie()) { //delete file
  //   gSystem->Setenv("PROUTFILENAME", PROutFile.Data());
  //   gROOT->ProcessLine(".! rm $PROUTFILENAME");
  //   gSystem->Unsetenv("PROUTFILENAME");
  // }
  fRun->SetOutputFile(outFile);
  //DONE create output file -----------------------------------------------  
 
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo();
  TString tpcDigiFile = basedir;
  tpcDigiFile+="/tpc/FOPI/par/tpc.run";
  tpcDigiFile+=run;
  tpcDigiFile +=".par";
  cout<<tpcDigiFile<<endl;
  fstream fin;
  fin.open(tpcDigiFile.Data(),ios::in);
  if( !fin.is_open() )
    {
      cout<<"parfile doesnt exist"<<endl;
      return;
    }
  fin.close();
  parInput1->open(tpcDigiFile.Data(),"in"); 
  rtdb->setFirstInput(parInput1);


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

  Double_t taroff = -65;
  Double_t targetpos = taroff+41.;
  // possible materials: copper (default), lead, carbon
  TString targetmat = TString("carbon");
  TString alignmentFileName = par->getAlignmentFile();
  TString geoFile;
  geoFile.Form("$FOPI2ROOT/fopigeometry/FopiGeom_s339%5.1f.root",targetpos);

  TString cmd;
  cmd.Form(".! root -b -q '$FOPI2ROOT/fopigeometry/FopiGeom.C(%f,\"%s\",\"%s\",\"%s\",kFALSE)'",targetpos,targetmat.Data(),alignmentFileName.Data(),geoFile.Data(),false,false);
  gROOT->ProcessLine(cmd);
  //TString geoFile = basedir;
  fRun->SetGeomFile(geoFile);


  

  FOPIField *fMagField = new FOPIField(0.616);
  fMagField->SetTargetOffset(taroff);
  // PndConstField *fMagField=new PndConstField();
  // if (bfield == 1)
  //   fMagField->SetField(0., 0. , 6. ); // values are in kG
  // if (bfield == 0)
  //   fMagField->SetField(0., 0. , 0. ); // values are in kG
  // // values are in cm
  // fMagField->SetFieldRegion(-50, 50,-50, 50, -2000, 2000);
  fRun->SetField(fMagField);
  
  //extract number of entries in external data tree
  cout<<"filename: "<<filename<<endl;
  TFile testFile(filename);
  TTree* testtree = testFile.Get("tpcEvent");
  unsigned int numEvents = testtree->GetEntries();
  //      unsigned int numEvents = ((TTree*)testFile.Get("cbmsim"))->GetEntries();
  std::cout<<"Found "<<numEvents<<" events in input data file"<<std::endl;
  if(nEvents==0 || nEvents>numEvents) nEvents=numEvents;
  testFile.Close();
  
  //  TpcAlignmentManager * align = TpcAlignmentManager::getInstance("/nfs/hicran/project/panda/SIM/bergerm/fopiroot/tpc/FOPI/par/dummy_alignment.txt");
  //--------------------SET UP TASKS ------------------------------
  TpcEventCounter* evCount =  new TpcEventCounter();
  //evCount->SetnEvts(nEvents);
  //evCount->SetStep(1);
  fRun->AddTask(evCount); 


  TpcDataReaderTask* read = new TpcDataReaderTask();
  read->SetPersistence();
  read->SetDatafile(filename);
  read->SetSampleBranchName("TpcSample");
  read->SetStartEvent(nEvStart);
  //read->SetCutSmallPad();
  //read->SetMinSamples(1000);
  fRun->AddTask(read);
  

  TpcPSATask* tpsa = new  TpcPSATask();
  tpsa->SetPersistence();
  tpsa->SetSamplePersistence();
  tpsa->SetSampleBranchName("TpcSample"); // Input of PSA
  fRun->AddTask(tpsa);

  if (gaincorr==1)
    {
      TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask();
      fRun->AddTask(tpcCalib);
    }

  TpcClusterFinderTask* tpcCF1 = new TpcClusterFinderTask();
  if(clcorr)
    tpcCF1->SetClusterBranchName("TpcPreCluster_raw");
  else
    tpcCF1->SetClusterBranchName("TpcPreCluster");
  //tpcCF1->SetDigiPersistence(); // keep Digis copys in clusters
  //tpcCF1->SetPersistence(); // keep Clusters
  tpcCF1->timeslice(6); //in samples
  tpcCF1->SetErrorPars(600.,400.);
  tpcCF1->SetSimpleClustering(); // use TpcClusterFinderSimple
  fRun->AddTask(tpcCF1);
  
  if (clcorr)
    {
      cout<<"Setting up cluster correction task"<<endl;
      TpcClusterCorrectionTask* tpcCC = new TpcClusterCorrectionTask();
      tpcCC->SetClusterBranchName("TpcPreCluster_raw");
      tpcCC->SetClusterOutBranchName("TpcPreCluster");
      //tpcCC->SetDevMap("datafiles/full_field_-148.93_-312.89_4s3_4.7_merged.txt");
      tpcCC->SetUseDevMap();
      fRun->AddTask(tpcCC);
    }


  cout<<"Setting up Rieman Task"<<endl;
  //find TpcRiemannTracks in the TPC alone
  TpcRiemannTrackingTask* tpcSPR1 = new TpcRiemannTrackingTask();
  tpcSPR1->SetClusterBranchName("TpcPreCluster");
  tpcSPR1->SetRiemannTrackBranchName("RiemannTrack");
  tpcSPR1->SetRiemannHitBranchName("RiemannHit");
  //tpcSPR1->SetPersistence();
  tpcSPR1->SetSortingParameters(true,1,0);
  fRun->AddTask(tpcSPR1);

  if (reclust)
    {
      cout<<"Setting up TrackInitTask"<<endl;
      //build GFTracks from TpcRiemannTracks
      TpcTrackInitTask* trackInit1=new TpcTrackInitTask();
      //trackInit1->SetPersistence();
      trackInit1->SetRiemannBranchName("RiemannTrack");
      trackInit1->SetTrackOutBranchName("TpcPreTrackPreFit");
      trackInit1->SetClusterBranchName("TpcPreCluster");
      trackInit1->SetRecoHitOutBranchName("TpcPreSPHit");
      trackInit1->SetCosmicSorting();
      trackInit1->SetVerbose(0);
      trackInit1->SetPDG(211);
      //trackInit1->SetAlign(kFALSE);
      if (smoothing==1)
	trackInit1->SetSmoothing(true);
      fRun->AddTask(trackInit1);

      cout<<"Setting up KalmanTask 1"<<endl;
      KalmanTask* kalman1 =new KalmanTask();
      kalman1->SetTpcClusterBranchName("TpcPreSPHit");
      kalman1->SetTrackBranchName("TpcPreTrackPreFit");
      kalman1->SetOutBranchName("TrackPrePostFit");
      //kalman1->SetPersistence();
      kalman1->SetNumIterations(1); // number of fitting iterations (back and forth)
      fRun->AddTask(kalman1);


      cout<<"Setting up Reclusterizer"<<endl;
      TpcReClusterizerTask* tpcReCl = new TpcReClusterizerTask();
      tpcReCl->SetTrackBranchName("TrackPrePostFit");
      tpcReCl->SetClusterBranchName("TpcPreCluster");
      tpcReCl->SetRecoHitBranchName("TpcPreSPHit");
      tpcReCl->SetReClusterBranchName("TpcCluster");
      //tpcReCl->SetAlign(kFALSE);
      tpcReCl->SetStepSize(1);
      tpcReCl->SetMinClusterSize(2);
      //tpcReCl->SetVerbose();
      //tpcReCl->SetTimeInfo();
      //tpcReCl->SetPersistence();
      //tpcReCl->SetAllRecoHitPersistence();
      tpcReCl->SetUseFirstDigiPos();
      //tpcReCl->SetUseChamberEdge();
      tpcReCl->SetUseCosmics();
      //tpcReCl->SetUseAllDigis();
      tpcReCl->SetTrackInit();
      fRun->AddTask(tpcReCl);

    }
  else
    {
      TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask();
      //tpcCF->SetDigiPersistence(); // keep Digis copys in clusters
      //tpcCF->SetPersistence(); // keep Clusters
      tpcCF->timeslice(6); //in samples
      // tpcCF->SetSingleDigiClusterAmpCut(15);
      // tpcCF->SetClusterAmpCut(9.1); // cut on mean digi amplitude
      tpcCF->SetErrorPars(600.,400.);
      tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple
      //      fRun->AddTask(tpcCF);
    

      //find TpcRiemannTracks in the TPC alone
      TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask();
      tpcSPR->SetClusterBranchName("TpcCluster");
      tpcSPR->SetRiemannTrackBranchName("RiemannTrack");
      //tpcSPR->SetPersistence();
      //tpcSPR->SetVerbose();
      tpcSPR->SetSortingParameters(true,1,0); //check if really needed
      //fRun->AddTask(tpcSPR);

      //build GFTracks from TpcRiemannTracks
      TpcTrackInitTask* trackInit=new TpcTrackInitTask();
      trackInit->SetRiemannBranchName("RiemannTrack");
      trackInit->SetClusterBranchName("TpcPreCluster");
      trackInit->SetTrackOutBranchName("TpcTrackPreFit");
      trackInit->SetRecoHitOutBranchName("TpcSPHit");
      //trackInit->SetPersistence();
      trackInit->SetCosmicSorting();
      trackInit->SetVerbose(0);
      //trackInit->SetMCPid(); // use ideal particle identification
      trackInit->SetPDG(211);
      //trackInit->useGeane(); // uses RKTrackrep and GeaneTrackrep
      if (smoothing==1)
	trackInit->SetSmoothing(true);
      fRun->AddTask(trackInit);


      // cout<<"Setting up AlignmentTask"<<endl;
      // TpcAlignmentTask* align = new TpcAlignmentTask();
      // align->SetRecoHitBranchName("TpcSPHit");
      // fRun->AddTask(align);
    }

  KalmanTask* kalman =new KalmanTask();
  kalman->SetTpcClusterBranchName("TpcSPHit");
  kalman->SetTrackBranchName("TpcTrackPreFit");
  kalman->SetOutBranchName("TrackPostFit");
  //kalman->SetPersistence();
  //kalman->SetClusterBranchName("TpcCluster_cut");
  kalman->SetNumIterations(3); // number of fitting iterations (back and forth)
  fRun->AddTask(kalman);

  cout<<"Setting up CosmicsTask"<<endl;
  TpcCosmicsTask* Cos = new TpcCosmicsTask();
  if (reclust)
    Cos->SetClusterBranchName("TpcCluster");
  else
    Cos->SetClusterBranchName("TpcPreCluster");
  Cos->SetPersistence();
  Cos->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();)
  Cos->SetDatafile(filename);
  if (smoothing==1)
    Cos->SetUnbiased();
  //Cos->SetVerbose();
  fRun->AddTask(Cos);

  cout<<"Setting up Track Cleaner"<<endl;
  TpcTrackCleaner* tpcTrCl = new TpcTrackCleaner();
  tpcTrCl->SetCosmics();
  tpcTrCl->SetThetaCut(40,140);
  tpcTrCl->SetClMeanCut(400);
  //tpcTrCl->SetPhiCut(70,110);
  //tpcTrCl->SetVerbose();
  fRun->AddTask(tpcTrCl);


  // TpcResidualTask* Res = new TpcResidualTask();
  // Res->SetPersistence();
  // Res->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();)
  // if (smoothing==1)
  //   Res->SetUnbiased();
  // fRun->AddTask(Res);
  

  // -----   Intialise and run   --------------------------------------------
 
  cout<<"************start ini***************"<<endl;
  fRun->Init();
  //Initialize te Digimapper:
  TpcDigiPar* tpcpar = FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar");
  TpcDigiMapper::getInstance()->init(tpcpar);
  //  TpcDigiMapper::getInstance()->forceManualDriftVel(kTRUE);
  TpcAlignmentManager::init(tpcpar->getAlignmentFile());
  cout<<"********************************************"<<endl;
  std::cout<<"pads:"<<TpcDigiMapper::getInstance()->getPadPlane()->GetNPads()<<std::endl;

  //  if (gaincorr==1)
    //    tpcCalib->SetExternalInput("./tpc/FOPI/par/gainmaps/Krypton_GEM85p_Nov_2.csv");
  rtdb->print();

  std::cout<<"start:"<<nEvStart<<" end:"<<nEvents<<std::endl;

  //  startClustVis(filename);
  fRun->Run(nEvStart,nEvents);

  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;


  //return 0;
}