void runCosmicAna_mb(TString recofile,int run, unsigned int nEvents = 0, unsigned int nEvStart = 0)
{
// ========================================================================
// Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
  Int_t iVerbose = 1;

  TStopwatch timer;
  timer.Start();
  
  TString basedir = gSystem->Getenv("VMCWORKDIR");
  FairRunAna* fRun = new FairRunAna();
  //  fRun->SetInputFile(recofile);
  fRun->SetInputFile("outfiles/parts/runC_3888_smoothed_000.reco.root");
  fRun->AddFile("outfiles/parts/runC_3888_smoothed_001.reco.root");

  //create output file -----------------------------------------------  
  TString outName(recofile);
  outName.ReplaceAll(".reco","");
  outName.ReplaceAll(".root","_ana.root");
  cout<<outName<<"\n";
  fRun->SetOutputFile(outName);
  //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);
  fRun->SetGeomFile(geoFile);

  FOPIField *fMagField = new FOPIField(0.616);
  fMagField->SetTargetOffset(taroff);
  fRun->SetField(fMagField);
  
  // cout<<"filename: "<<recofile<<endl;
  // TFile testFile(recofile);
  // TTree* testtree = testFile.Get("tpcEvent");
  // unsigned int numEvents = testtree->GetEntries();
  // std::cout<<"Found "<<numEvents<<" events in input data file"<<std::endl;
  // if(nEvents==0 || nEvents>numEvents) nEvents=numEvents;
  // testFile.Close();

  //--------------------SET UP TASKS ------------------------------
  TpcEventCounter* evCount =  new TpcEventCounter();
  //evCount->SetnEvts(nEvents);
  //evCount->SetStep(1);
  fRun->AddTask(evCount); 

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

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

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

  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: "<<outName<<std::endl;
}