void ana_hits(Int_t nEvents=10, Int_t iSel=1, Int_t iGenCor=1, char *cFileId="", Int_t iSet=0, Int_t iSel2=0) 
{
   Int_t iVerbose = 1;
   // Specify log level (INFO, DEBUG, DEBUG1, ...)
   TString logLevel = "FATAL";
   //TString logLevel = "ERROR";
   TString logLevel = "INFO";
   //TString logLevel = "DEBUG"; 
   //TString logLevel = "DEBUG1";
   //TString logLevel = "DEBUG2";
   //TString logLevel = "DEBUG3";
   FairLogger* log;  

   TString workDir       = gSystem->Getenv("VMCWORKDIR");
   TString paramDir      = workDir  + "/macro/tof/beamtime/thuqa";
   TString ParFile       = paramDir + "/unpack_" + cFileId + ".params.root";
   TString InputFile     = paramDir + "/unpack_" + cFileId + ".out.root";
   TString InputDigiFile = paramDir + "/digi_" + cFileId + Form("_%03d",iSet) + ".out.root";
   TString OutputFile    = paramDir + "/hits_" + cFileId + Form("_%03d%1d",iSet,iSel2) + ".out.root";
   TString cAnaFile=Form("%s_%03d%1d_%1d_tofAnaTestBeam.hst.root",cFileId,iSet,iSel2,iSel);

   cout << " InputDigiFile = "
	<< InputDigiFile
	<< endl;

   TList *parFileList = new TList();

   TObjString mapParFile = paramDir + "/parMapCosmicsThuMay2015.txt";
   parFileList->Add(&mapParFile);

   TString TofGeo="v14d";
   TString FId=cFileId;
   TString FPar="";
   if(FId.Contains("tsu")) {
     cout << FId << ": Analyse TSU counter "<< endl;
     //TofGeo="v14c_tsu";
     //FPar="tsu.";
   }

   TObjString tofDigiFile = workDir + "/parameters/tof/tof_" + TofGeo + ".digi.par"; // TOF digi file
   parFileList->Add(&tofDigiFile);   

   TObjString tofDigiBdfFile =  paramDir + "/tof.digibdf.par";
   //   TObjString tofDigiBdfFile =  paramDir + "/tof." + FPar + "digibdf.par";
   parFileList->Add(&tofDigiBdfFile);

   TString geoDir  = gSystem->Getenv("VMCWORKDIR");
   TString geoFile = geoDir + "/geometry/tof/geofile_tof_" + TofGeo + ".root";
   TFile* fgeo = new TFile(geoFile);
   TGeoManager *geoMan = (TGeoManager*) fgeo->Get("FAIRGeom");
   if (NULL == geoMan){
     cout << "<E> FAIRGeom not found in geoFile"<<endl;
     return;
   }
   if(0){
   TGeoVolume* master=geoMan->GetTopVolume();
   master->SetVisContainers(1); 
   master->Draw("ogl"); 
   }

  // -----   Reconstruction run   -------------------------------------------
  FairRunAna *run= new FairRunAna();
  cout << "InputFile:     "<<InputFile.Data()<<endl;
  cout << "InputDigiFile: "<<InputDigiFile.Data()<<endl;

  run->SetInputFile(InputFile.Data());
  run->AddFriend(InputDigiFile.Data());
  //run->SetInputFile(InputDigiFile.Data());
  //run->AddFriend(InputFile.Data());
  run->SetOutputFile(OutputFile);
  gLogger->SetLogScreenLevel(logLevel.Data());

   // =========================================================================
   // ===                         Mapping                                   ===
   // =========================================================================
   TMbsMappingTof* tofMapping = new TMbsMappingTof("Tof Mapping", iVerbose);
   //   run->AddTask(tofMapping);

   // =========================================================================
   // ===                       Tracking                                    ===
   // =========================================================================
   CbmTofTrackFinder* tofTrackFinder= new CbmTofTrackFinderNN();
   CbmTofFindTracks* tofFindTracks = new CbmTofFindTracks("TOF Track Finder");
   tofFindTracks->UseFinder(tofTrackFinder);
   run->AddTask(tofFindTracks);

   // =========================================================================
   // ===                       Analysis                                    ===
   // =========================================================================
   CbmTofAnaTestbeam* tofAnaTestbeam = new CbmTofAnaTestbeam("TOF TestBeam Analysis",iVerbose);
   tofAnaTestbeam->SetCorMode(iGenCor); // 1 - DTD4, 2 - X4, 3 - Y4, 4 - Texp 
   tofAnaTestbeam->SetHitDistMin(30.);  // initialization

   //CbmTofAnaTestbeam defaults  
   tofAnaTestbeam->SetDXMean(0.);
   tofAnaTestbeam->SetDYMean(0.);
   tofAnaTestbeam->SetDTMean(0.);      // in ps
   tofAnaTestbeam->SetDXWidth(0.5);
   tofAnaTestbeam->SetDYWidth(0.5);
   tofAnaTestbeam->SetDTWidth(90.);    // in ps
   tofAnaTestbeam->SetCalParFileName(cAnaFile);
   tofAnaTestbeam->SetPosY4Sel(0.5);   // Y Position selection in fraction of strip length
   tofAnaTestbeam->SetDTDia(0.);       // Time difference to additional diamond
   tofAnaTestbeam->SetMul0Max(20);     // Max Multiplicity in dut 
   tofAnaTestbeam->SetMul4Max(10);      // Max Multiplicity in Ref - RPC 
   tofAnaTestbeam->SetMulDMax(20);      // Max Multiplicity in Diamond    
   tofAnaTestbeam->SetTOffD4(10000.);  // initialization
   tofAnaTestbeam->SetDTD4MAX(6000.);     // initialization of Max time difference Ref - BRef

   //tofAnaTestbeam->SetTShift(-28000.);  // initialization
   tofAnaTestbeam->SetPosYS2Sel(0.5);   // Y Position selection in fraction of strip length
   tofAnaTestbeam->SetChS2Sel(0.);      // Center of channel selection window
   tofAnaTestbeam->SetDChS2Sel(100.);   // Width  of channel selection window
   tofAnaTestbeam->SetSel2TOff(450.);   // Shift Sel2 time peak to 0 

   Int_t iRSel=0;
   if(iSel2>=0){
     iRSel=2;   // use diamond											//change PW
     tofAnaTestbeam->SetTShift(-2000.);  // initialization
   }else{
     iSel2=-iSel2;
     iRSel=iSel2;
     tofAnaTestbeam->SetTShift(50.);  // initialization
   }
   tofAnaTestbeam->SetMrpcSel2(iSel2); // initialization of second selector Mrpc 
   tofAnaTestbeam->SetChi2Lim(10000.);     // initialization of Chi2 selection limit  
   tofAnaTestbeam->SetBeamRefSmType(iRSel); // common reaction reference 
   tofAnaTestbeam->SetBeamRefSmId(0);

  switch(iSet){
  case 920:
  case 921:
    gInterpreter->ProcessLine("pl_over_clu(2,0)");
	gInterpreter->ProcessLine("pl_over_clu(2,1)");
    gInterpreter->ProcessLine("pl_over_clu(9)");
    gInterpreter->ProcessLine("pl_over_cluSel(0,9)"); //Just plastic SM 0 Selector?!?
    gInterpreter->ProcessLine("pl_all_dTSel()");
    break;

  case 2021:
  case 2120:
    gInterpreter->ProcessLine("pl_over_clu(2,0)");
	gInterpreter->ProcessLine("pl_over_clu(2,1)");
    gInterpreter->ProcessLine("pl_over_clu(2)");
    gInterpreter->ProcessLine("pl_over_cluSel(0,2)"); //Just plastic SM 0 Selector?!?
    gInterpreter->ProcessLine("pl_all_dTSel()");
    break;

  default:
    ;
  }


   }
   cout << " Initialize TSHIFT to "<<tofAnaTestbeam->GetTShift()<<endl;
   run->AddTask(tofAnaTestbeam);

   // =========================================================================
   /*
   CbmTofOnlineDisplay* display = new CbmTofOnlineDisplay();
   display->SetUpdateInterval(1000);
   run->AddTask(display);   
   */
  // -----  Parameter database   --------------------------------------------
   FairRuntimeDb* rtdb = run->GetRuntimeDb();
   Bool_t kParameterMerged = kTRUE;
   FairParRootFileIo* parIo2 = new FairParRootFileIo(kParameterMerged);
   parIo2->open(ParFile.Data(), "UPDATE");
   parIo2->print();
   rtdb->setFirstInput(parIo2);
   
   FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
   parIo1->open(parFileList, "in");
   parIo1->print();
   rtdb->setSecondInput(parIo1);
   rtdb->print();
   rtdb->printParamContexts();

   //  FairParRootFileIo* parInput1 = new FairParRootFileIo();
   //  parInput1->open(ParFile.Data());
   //  rtdb->setFirstInput(parInput1);

  // -----   Intialise and run   --------------------------------------------
  run->Init();
  cout << "Starting run" << endl;
  run->Run(0, nEvents);
  // ------------------------------------------------------------------------
  // default display 
  TString Display_Status = "pl_over_Mat04D4best.C";
  TString Display_Funct = "pl_over_Mat04D4best()";  
  gROOT->LoadMacro(Display_Status);
  gInterpreter->ProcessLine(Display_Funct);

  gROOT->LoadMacro("pl_over_MatD4sel.C");
  gInterpreter->ProcessLine("pl_over_MatD4sel()");
}