void anatut_psi2s(int nevts=0)
{
	// *** some variables
	int i=0,j=0, k=0, l=0;
	
	TString OutFile="output.root";  
					
	// *** the files coming from the simulation
	TString inPidFile  = "pid_complete.root";    // this file contains the PndPidCandidates
	TString inRecoFile  = "reco_complete.root";
	TString inSimFile = "sim_complete.root";  // this file contains the MC truth
	TString inParFile = "simparams.root";
	
	gStyle->SetOptFit(1011);
	gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
	
	FairLogger::GetLogger()->SetLogToFile(kFALSE);
	
	// *** initialization
	FairRunAna* fRun = new FairRunAna();
	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
	fRun->SetInputFile(inSimFile);
	fRun->AddFriend(inPidFile);
	fRun->AddFriend(inRecoFile);
	
	FairParRootFileIo* parIO = new FairParRootFileIo();
	parIO->open(inParFile);
	rtdb->setFirstInput(parIO);
	rtdb->setOutput(parIO);  
	
	fRun->SetOutputFile(OutFile);
	fRun->Init(); 
	
	// *** create an output file for all histograms
	TFile *out = TFile::Open("output_psi2sana.root","RECREATE");
	
	// *** create some histograms
	TH1F *hdplus_nopid = new TH1F("hdplus_nopid","J/#psi mass (no pid)",200,0,2.5);	
	TH1F *hdplus_ftm = new TH1F("hdplus_ftm","J/#psi mass (full truth match)",200,0,2.5);
	TH1F *hdplus_nm = new TH1F("hdplus_nm","J/#psi mass (no truth match)",200,0,2.5);
	TH1F *hdplus_pid = new TH1F("hdplus_pid","J/#psi mass (no pid)",200,0,2.5);	
	TH1F *hmcdplus_ftm = new TH1F("hmcdplus_ftm","J/#psi mass (full truth match)",200,0,2.5);
	TH1F *hmcdplus_nm = new TH1F("hmcdplus_nm","J/#psi mass (full truth match)",200,0,2.5);
	TH1F *hdminus_nopid = new TH1F("hdminus_nopid","J/#psi mass (no pid)",200,0,2.5);	
	TH1F *hdminus_ftm = new TH1F("hdminus_ftm","J/#psi mass (full truth match)",200,0,2.5);
	TH1F *hdminus_nm = new TH1F("hdminus_nm","J/#psi mass (no truth match)",200,0,2.5);
	TH1F *hdminus_pid = new TH1F("hdminus_pid","J/#psi mass (no pid)",200,0,2.5);
	TH1F *hmcdminus_ftm = new TH1F("hmcdminus_ftm","J/#psi mass (full truth match)",200,0,2.5);
	TH1F *hmcdminus_nm = new TH1F("hmcdminus_nm","J/#psi mass (full truth match)",200,0,2.5);

	// *** the data reader object
	PndAnalysis* theAnalysis = new PndAnalysis("SttMvdGemGenTrack","FtsIdealGenTrack");
	if (nevts==0) nevts= theAnalysis->GetEntries();
	
	// *** TCandLists for the analysis
	TCandList mctrk, kplus, kminus, piplus, piminus, dplus, dminus;
	TCandList mckplus, mckminus, mcpiplus, mcpiminus, mcdplus, mcdminus;
	// *** the MC truth matcher
	PndMcTruthMatch mcm;
	
	// *** the lorentz vector of the initial psi(3770)
	TLorentzVector ini(0, 0, 6.5788, 7.58364);
	
	// ***
	// the event loop
	// ***
	while (theAnalysis->GetEvent() && i++<nevts)
	{
		if ((i%100)==0) cout<<"evt " << i << endl;
		
		mctrk.Cleanup(); kplus.Cleanup(); kminus.Cleanup(); piplus.Cleanup(); piminus.Cleanup(); dplus.Cleanup(); dminus.Cleanup();
		mckplus.Cleanup(); mckminus.Cleanup(); mcpiplus.Cleanup(); mcpiminus.Cleanup(); mcdplus.Cleanup(); mcdminus.Cleanup();
		
		// *** the MC Truth objects
		theAnalysis->FillList(mctrk,"McTruth");
			
		// *** Select with no PID info ('All'); type and mass are set 		
		theAnalysis->FillList(kplus, "KaonAllPlus");
		theAnalysis->FillList(kminus, "KaonAllMinus");
		theAnalysis->FillList(piplus, "PionAllPlus");
		theAnalysis->FillList(piminus, "PionAllMinus");
	       
		dplus.Combine(kminus, piplus, piplus);
		dminus.Combine(kplus, piminus, piminus);
		
		// *** do the truth match for jpsi
		mcm.SetType(dplus,"D+");
		mcm.SetType(dminus,"D-");
		
		if (dplus.GetLength()>0)
		  for (j=0;j<dplus.GetLength();++j) 
		    {
		      hdplus_nopid->Fill( dplus[j].M() );
		      if (mcm.MctMatch(dplus[j], mctrk)) hdplus_ftm->Fill( dplus[j].M() );
		      else hdplus_nm->Fill( dplus[j].M() );
		    }
		if (dminus.GetLength()>0)
		  for (j=0;j<dminus.GetLength();++j) 
		    {
		      hdminus_nopid->Fill( dminus[j].M() );
		      if (mcm.MctMatch(dminus[j], mctrk)) hdminus_ftm->Fill( dminus[j].M() );
		      else hdminus_nm->Fill( dminus[j].M() );
		    }
		
		mcm.SetType(kplus,"k+");
		mcm.SetType(kminus,"k-");
		mcm.SetType(piplus,"pi+");
		mcm.SetType(piminus,"pi-");
		if (kplus.GetLength()>0)
		  for (j=0;j<kplus.GetLength();++j) 
		    {
		      if (mcm.MctMatch(kplus[j], mctrk))  mckplus.Put( kplus[j] ); 
		    }
		if (kminus.GetLength()>0)
		  for (j=0;j<kminus.GetLength();++j) 
		    {
		      if (mcm.MctMatch(kminus[j], mctrk))  mckminus.Put( kminus[j] ); 
		    }
		if (piplus.GetLength()>0)
		  for (j=0;j<piplus.GetLength();++j) 
		    {
		      if (mcm.MctMatch(piplus[j], mctrk))  mcpiplus.Put( piplus[j] ); 
		    }
		if (piminus.GetLength()>0)
		  for (j=0;j<piminus.GetLength();++j) 
		    {
		      if (mcm.MctMatch(piminus[j], mctrk))  mcpiminus.Put( piminus[j] ); 
		    }
		
		mcdplus.Combine(mckminus, mcpiplus, mcpiplus);
		mcdminus.Combine(mckplus, mcpiminus, mcpiminus);
		
		mcm.SetType(mcdplus,"D+");
		mcm.SetType(mcdminus,"D-");
		if (mcdplus.GetLength()>0)
		  for (j=0;j<mcdplus.GetLength();++j) 
		    {
		      hdplus_pid->Fill( mcdplus[j].M() );
		      if (mcm.MctMatch(mcdplus[j], mctrk)) hmcdplus_ftm->Fill( mcdplus[j].M() );
		      else hmcdplus_nm->Fill( mcdplus[j].M() );
		    }
		if (mcdminus.GetLength()>0)
		  for (j=0;j<mcdminus.GetLength();++j) 
		    {
		      hdminus_pid->Fill( mcdminus[j].M() );
		      if (mcm.MctMatch(mcdminus[j], mctrk)) hmcdminus_ftm->Fill( mcdminus[j].M() );
		      else hdminus_nm->Fill( mcdminus[j].M() );
		    }
	}
	
	// *** write out all the histos
	out->cd();
	
	hdplus_nopid->Write();
	hdplus_ftm->Write();
	hdplus_nm->Write();
	hdplus_pid->Write();
	hmcdplus_ftm->Write();
	hdminus_nopid->Write();
	hdminus_ftm->Write();
	hdminus_nm->Write();
	hdminus_pid->Write();
	hmcdminus_ftm->Write();
	
	out->Save();
	
}