void analyse(TString filename="data/test_sim.root", int nevts=0)
{
  TStopwatch timer;
  timer.Start();
  int verbose=0;
  //int canvasplot=1; //whether to plot into a ps, too.
	gSystem->Load("libEGPythia6");  // needed for TDatabasePDG
	gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
	gSystem->Load("libRho");
	gSystem->Load("libPid");
	gSystem->Load("libAnalysisTools");
  
  gROOT->LoadMacro("$VMCWORKDIR/macro/mvd/Tools.C");
  LoadPandaStyle();
    
  int bins=100;
  double min=0,max=2.5; // GeV/c
  double min2=0.13,max2=0.15;
  double min3=1.,max3=1.5;
  double max4=0.3;//cm
  double maxchi = 50;
  
  TList* histolist = new TList();
  
  
  TH1F *hMompip=new TH1F("momppip","P(#pi+);p / GeV/c",bins,min,max);
  TH1F *hMompim=new TH1F("momppim","P(#pi-);p / GeV/c",bins,min,max);
  TH1F *hMomKp=new TH1F("mompkp","P(K+);p / GeV/c",bins,min,max);
  TH1F *hMomKm=new TH1F("mompkm","P(K-);p / GeV/c",bins,min,max);
  TH1I *hTFStat=new TH1I("tfstat","Track Fit Status Flag",10,-5,5);
  
  
  // Add histogram pointers to a list for iterative saving & plotting
  //---newpage
  histolist->Add(hMompip);
  histolist->Add(hMompim);
  histolist->Add(hMomKp);
  histolist->Add(hMomKm);
  histolist->Add(hTFStat);
  
  // Set up the Analysis
  TString histoFile = InitRun(filename);
  PndAnalysis* theAnalysis = new PndAnalysis();
  
  // the candidates lists we need
  TCandList pionp, pionm, kaonp, kaonm ,mctracks;
  TCandidate *tmpcand=0, *recocand=0, *mccand=0;
  
  int evts = theAnalysis->GetEntries();
  if (nevts>0 && nevts<evts) evts=nevts; 
  
  double tfcut=-1; // trackfit status cut
  
  TLorentzVector tmpMom;
  double mass=0,momentum=0,theta=0,chi2,energy=0;
  int ievt=0,j=0,mcidx=-1,fitstat=0;
  bool tester=false;
  
  
  // *************
  // this is the loop through the events ... as simple as this...
  // ****************
  std::cout<<"   *****   Start the eventloop with "<<evts<<" events.   *****"<<std::endl;
  while (theAnalysis->GetEvent() && ievt++<evts)
  {    
    theAnalysis->FillList(mctracks,"McTruth");
    theAnalysis->FillList(pionp,"PionLoosePlus");
    theAnalysis->FillList(pionm,"PionLooseMinus");
    theAnalysis->FillList(kaonp,"KaonLoosePlus");
    theAnalysis->FillList(kaonm,"KaonLooseMinus");
    if(verbose>0) printf("Event %i: TCandList sizes: pionp:%i, pionm:%i, kaonp:%i, kaonm:%i, mctracks:%i\n",ievt,pionp.GetLength(),pionm.GetLength(),kaonp.GetLength(),kaonm.GetLength(),mctracks.GetLength());
    
    for (j=0;j<pionp.GetLength();++j) {
      recocand=pionp.Get(j);
      if(!recocand) {if(verbose>0) cout<<"############## pionp\n#############"<<endl;continue;}
      if(verbose>1) std::cout<<" Pion+: "<<*recocand<<std::endl;
      fitstat=recocand->GetMicroCandidate().GetFitStatus();
      hTFStat->Fill(fitstat);
      if(fitstat<tfcut) {
        pionm.Remove(*recocand);
        continue;
      }
      tmpMom=recocand->P4();
      momentum=tmpMom.P();
      hMompip->Fill(momentum);
    }
    for (j=0;j<pionm.GetLength();++j) {
      recocand=pionm.Get(j);
      if(!recocand) {if(verbose>0) cout<<"############## pionm\n#############"<<endl;continue;}
      if(verbose>1) std::cout<<" Pion-: "<<*recocand<<std::endl;
      fitstat=recocand->GetMicroCandidate().GetFitStatus();
      hTFStat->Fill(fitstat);
      if(fitstat<tfcut) {
        pionm.Remove(*recocand);
        continue;
      }
      tmpMom=recocand->P4();
      momentum=tmpMom.P();
      hMompim->Fill(momentum);
    }
    for (j=0;j<kaonp.GetLength();++j) {
      recocand=kaonp.Get(j);
      if(!recocand) {if(verbose>0) cout<<"############## kaonp\n#############"<<endl;continue;}
      if(verbose>1) std::cout<<" K+: "<<*recocand<<std::endl;
      fitstat=recocand->GetMicroCandidate().GetFitStatus();
      hTFStat->Fill(fitstat);
      if(fitstat<tfcut) {
        kaonp.Remove(*recocand);
        continue;
      }
      tmpMom=recocand->P4();
      momentum=tmpMom.P();
      hMomKp->Fill(momentum);
    }
    for (j=0;j<kaonm.GetLength();++j) {
      recocand=kaonm.Get(j);
      if(!recocand) {if(verbose>0) cout<<"############## kaonm\n#############"<<endl;continue;}
      if(verbose>1) std::cout<<" K-: "<<*recocand<<std::endl;
      fitstat=recocand->GetMicroCandidate().GetFitStatus();
      hTFStat->Fill(fitstat);
      if(fitstat<tfcut) {
        kaonm.Remove(*recocand);
        continue;
      }
      tmpMom=recocand->P4();
      momentum=tmpMom.P();
      hMomKm->Fill(momentum);
    }
    
    if(verbose>0) cout << " ----- event done -----"<<endl;
  } // end eventloop
  
  // Save all histograms to a file
  TFile* histofile = new TFile(histoFile.Data(),"RECREATE");
  histofile->cd();
  
  TListIter iter(histolist);
  TH1* tmph=NULL;
  while( tmph=(TH1*)iter() ) tmph->Write();
  histofile->Close();
  
  //if(canvasplot){
    //plothistos(histoFile);
  //}
    
  // now take the time...
  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();
  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
}


TString InitRun(TString filename)
{
  FairRunAna* fRun = new FairRunAna();
  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
  
  filename.ReplaceAll("_sim.root",".root");
  PndFileNameCreator namecreator(filename.Data());
  TString simFile = namecreator.GetSimFileName();
  TString parFile = namecreator.GetParFileName();
  TString recoFile = namecreator.GetRecoFileName();
  TString tracksFile = namecreator.GetCustomFileName("pidtracks");
  TString histoFile = namecreator.GetCustomFileName("histos");
  TString evrdummy = namecreator.GetCustomFileName("evrdummy");
  cout<<"simFile="<<simFile.Data()<<endl;
  cout<<"parFile="<<parFile.Data()<<endl;
  cout<<"recoFile="<<recoFile.Data()<<endl;
  cout<<"trkFile="<<tracksFile.Data()<<endl;
  fRun->SetInputFile(simFile);
  fRun->AddFriend(recoFile);//,"rec");
  fRun->AddFriend(tracksFile);//,"trk");
  FairParRootFileIo* parIO = new FairParRootFileIo();
  parIO->open(parFile.Data());
  rtdb->setFirstInput(parIO);
  rtdb->setOutput(parIO);  
  
  TString outFile = evrdummy;
  TString sysFile = gSystem->Getenv("VMCWORKDIR");  
  
  fRun->SetOutputFile(outFile.Data());
  FairGeane* geane = new FairGeane();
  fRun->AddTask(geane);
  fRun->Init();  
  return histoFile;
}