void ana_lambdac(TString InFile="test_fast.root", int nevts=0, double pbarmom = 15.)
{
	// *** some variables
	int i=0,j=0, k=0, l=0;
	gStyle->SetOptFit(1011);

	if (!InFile.EndsWith(".root")) InFile+="_fast.root";
	
	// *** the output file for FairRunAna
	TString OutFile="dummy_out.root";  
					
	// *** initialization
	FairLogger::GetLogger()->SetLogToFile(kFALSE);
	
	FairRunAna* fRun = new FairRunAna();
	fRun->SetWriteRunInfoFile(kFALSE);
	fRun->SetInputFile(InFile);
	fRun->SetOutputFile(OutFile); // only dummy; the real output is 
	
	fRun->Init(); 
	
	// *** take constant field; needed for PocaVtx
	RhoCalculationTools::ForceConstantBz(20.0);
	
	// *** create an output file for all histograms
	TFile *out = TFile::Open(InFile+"_ana.root","RECREATE");
	
	// *** create some ntuples
	//	RhoTuple *ntp1 = new RhoTuple("ntp1", "jpsi analysis");
	//	RhoTuple *ntp2 = new RhoTuple("ntp2", "psi(2S) analysis");
	RhoTuple *nmc  = new RhoTuple("nmc",  "mctruth info");
	
	//
	// Now the analysis stuff comes...
	//
	
	// *** the data reader object
	PndAnalysis* theAnalysis = new PndAnalysis();
	if (nevts==0) nevts= theAnalysis->GetEntries();
	
	// *** name of the only PidAlgo TClonesArray in fsim
	TString pidalg = "PidChargedProbability";

	// *** QA tool for simple dumping of analysis results in RhoRuple
	// *** WIKI: https://panda-wiki.gsi.de/foswiki/bin/view/Computing/PandaRootAnalysisJuly13#PndRhoTupleQA
	PndRhoTupleQA qa(theAnalysis, pbarmom); 

	// *** RhoCandLists for the analysis
	RhoCandList pplus, pminus, piplus, piminus, kaon, kplus, kminus, jpsi, psi2s, all, mclist, blah, lambdac, lambdacb, gamma, beam;

	RhoCandidate *cand;

	// *** Mass selector for the jpsi cands
	double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass();   // Get nominal PDG mass of the J/psi
	double m0_Kplus = TDatabasePDG::Instance()->GetParticle("K+")->Mass();   
	double m0_Kminus = TDatabasePDG::Instance()->GetParticle("K-")->Mass();   

	double m0_piplus = TDatabasePDG::Instance()->GetParticle("pi+")->Mass();   
	double m0_piminus = TDatabasePDG::Instance()->GetParticle("pi-")->Mass();   

	double m0_pplus = TDatabasePDG::Instance()->GetParticle("proton")->Mass();   
	double m0_pminus = TDatabasePDG::Instance()->GetParticle("antiproton")->Mass();   

	double m0_lplus = TDatabasePDG::Instance()->GetParticle("Lambda_c+")->Mass();   
	double m0_lminus = TDatabasePDG::Instance()->GetParticle("Lambda_c-")->Mass();   
	

	RhoMassParticleSelector *lpMassSel=new RhoMassParticleSelector("lambdac",m0_lplus,0.2);
	RhoMassParticleSelector *lmMassSel=new RhoMassParticleSelector("lambdacb",m0_lminus,1.0);
	
	// *** the lorentz vector of the initial psi(2S)
	double m0_p    = TDatabasePDG::Instance()->GetParticle("proton")->Mass();   // Get nominal PDG mass of the proton
	TLorentzVector ini(0, 0, pbarmom, sqrt(m0_p*m0_p + pbarmom*pbarmom) + m0_p);
	

	TH2F *hlm = new TH2F("hlm","",200,1.9,2.7,200,1.9,2.7);
	TH1F *hb = new TH1F("hb","",200,4,8);

	TH2F *hlmf = new TH2F("hlmf","",200,1.9,2.7,200,1.9,2.7);
	TH1F *hbf = new TH1F("hbf","",200,4,8);

	TH1F *hchisq = new TH1F("hchisq","",200,0,100.);
    

	// ***
	// the event loop
	// ***
		
	while (theAnalysis->GetEvent() && i++<nevts)
	{
		if ((i%100)==0) cout<<"evt " << i << endl;
		
		// *** get MC list and store info in ntuple

		
	 	theAnalysis->FillList(mclist, "McTruth");
		
		nmc->Column("ev", (Int_t) i);
		qa.qaMcList("",   mclist, nmc);
		//cout << mclist.GetNumberOfTracks() << endl;
		nmc->DumpData();
		
		// *** Setup event shape object
		theAnalysis->FillList(all,  "All", pidalg);
		
		PndEventShape evsh(all, ini, 0.05, 0.1);
			
		
		theAnalysis->FillList(kplus, "KaonLoosePlus", "DrcBarrelProbability;DrcDiscProbability");
		theAnalysis->FillList(kminus, "KaonLooseMinus", "DrcBarrelProbability;DrcDiscProbability");

		theAnalysis->FillList(piminus, "PionLooseMinus", "DrcBarrelProbability;DrcDiscProbability");
		theAnalysis->FillList(piplus, "PionLoosePlus", "DrcBarrelProbability;DrcDiscProbability");

		theAnalysis->FillList(pminus, "ProtonLooseMinus", "DrcBarrelProbability;DrcDiscProbability");
		theAnalysis->FillList(pplus, "ProtonLoosePlus", "DrcBarrelProbability;DrcDiscProbability");

		theAnalysis->FillList(gamma, "Neutral");

		//	cout << "******************" << endl;
		//cout << "All:" << all.GetNumberOfTracks() << endl;
		//cout << "K+:" << kplus.GetNumberOfTracks() << endl;		
		//cout << "K-:" << kminus.GetNumberOfTracks() << endl;		
		//cout << "pi+:" << piplus.GetNumberOfTracks() << endl;		
		//cout << "pi-:" << piminus.GetNumberOfTracks() << endl;		
		//cout << "p+:" << pplus.GetNumberOfTracks() << endl;		
		//cout << "p-:" << pminus.GetNumberOfTracks() << endl;		
		//cout << "g:" << gamma.GetNumberOfTracks() << endl;		


		lambdac.Combine(pplus,kminus,piplus);
		lambdacb.Combine(pminus,gamma);

		lambdac.Select(lpMassSel);
		lambdacb.Select(lmMassSel);

		beam.Combine(lambdac,lambdacb);
		
		for (j=0;j<beam.GetLength();j++)
		  {

		    Pnd4CFitter *fitter=new Pnd4CFitter(beam[j],ini);
		    fitter->FitConserveMasses();
		    hchisq->Fill(fitter->GetChi2());

		    hb->Fill(beam[j]->M());
		    hlm->Fill(beam[j]->Daughter(0)->M(),
			      beam[j]->Daughter(1)->M());

		    if (fitter->GetChi2()<20.)
		      {
			hbf->Fill(beam[j]->GetFit()->M());
			hlmf->Fill(beam[j]->Daughter(0)->GetFit()->M(),
				   beam[j]->Daughter(1)->GetFit()->M());
		      }
		    delete fitter;
		  }
		
		
		/*		
		jpsi.Combine(kplus, kminus);		
		jpsi.SetType(443);
		int njmct = theAnalysis->McTruthMatch(jpsi);
				
		for (j=0;j<jpsi.GetLength();++j) 
		{

		  cand = jpsi.Get(j);
		  cout << j << "/" << cand->M() << endl;

		    // some general info about event (actually written for each candidate!)
			ntp1->Column("ev",		(Float_t) i);
			ntp1->Column("cand",	(Float_t) j);
			ntp1->Column("ncand",   (Float_t) jpsi.GetLength());
			ntp1->Column("nmct",    (Float_t) njmct);
			
			// store info about initial 4-vector
			qa.qaP4("beam", ini, ntp1);
			
			// store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
			qa.qaComp("j", jpsi[j], ntp1);
			
			// store info about event shapes
			qa.qaEventShapeShort("es",&evsh, ntp1);
		
			
			// store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
			RhoCandidate *truth = jpsi[j]->GetMcTruth();		
			TLorentzVector lv;
			if (truth) lv = truth->P4();
			qa.qaP4("trj", lv, ntp1);
			
			ntp1->DumpData();
		}
		
		
		// *****
		// *** combinatorics for psi(2S) -> J/psi pi+ pi-
		// *****

		// *** some rough mass selection on J/psi before
		jpsi.Select(jpsiMassSel);
		
		psi2s.Combine(jpsi, piplus, piminus);
		psi2s.SetType(100443);
		int npsimct = theAnalysis->McTruthMatch(psi2s);

		for (j=0;j<psi2s.GetLength();++j) 
		{
		    // some general info about event (actually written for each candidate!)
			ntp2->Column("ev",		(Float_t) i);
			ntp2->Column("cand",	(Float_t) j);
			ntp2->Column("ncand",   (Float_t) psi2s.GetLength());
			ntp2->Column("nmct",    (Float_t) npsimct);
			
			PndKinFitter kinfit(psi2s[j]);
			kinfit.Add4MomConstraint(ini);
			kinfit.Fit();
			
			RhoCandidate *psifit=psi2s[j]->GetFit();
			
			// store info about initial 4-vector
			qa.qaP4("beam", ini, ntp2);
			
			// store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
			qa.qaComp("psi", psi2s[j], ntp2);
			qa.qaComp("fpsi",psifit, ntp2);
			
			// store info about event shapes
			qa.qaEventShapeShort("es",&evsh, ntp2);
			
			// *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
			RhoCandidate *truth = psi2s[j]->GetMcTruth();
			TLorentzVector lv;
			if (truth) lv = truth->P4();
			qa.qaP4("trpsi", lv, ntp2);
			
			ntp2->DumpData();
		}			
		*/	
	}
	
		
	// *** write out all the histos
	out->cd();
	
	//	ntp1->GetInternalTree()->Write();
	//	ntp2->GetInternalTree()->Write();
	nmc->GetInternalTree()->Write();
		
	out->Save();
	
}