// *** routine to only keep PID matched candidates in list
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
{
	int removed = 0;
	
	for (int ii=l.GetLength()-1;ii>=0;--ii)
	{
		if ( !(ana->McTruthMatch(l[ii])) )
		{
			l.Remove(l[ii]);
			removed++;
		}
	}
	
	return removed;
}

void ana_ste(TString InFile="test_fast.root", int nevts=0, double pbarmom = 6.231552)
{
	// *** 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 muplus, muminus, piplus, piminus, jpsi, psi2s, all, mclist;
	
	// *** Mass selector for the jpsi cands
	double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass();   // Get nominal PDG mass of the J/psi
	RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,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);
	
	// ***
	// 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);
		nmc->DumpData();
		
		
		// *** Setup event shape object
		theAnalysis->FillList(all,  "All", pidalg);
		PndEventShape evsh(all, ini, 0.05, 0.1);
			
		
		// *** Select with no PID info ('All'); type and mass are set 		
		theAnalysis->FillList(muplus,  "MuonAllPlus", pidalg);
		theAnalysis->FillList(muminus, "MuonAllMinus", pidalg);
		theAnalysis->FillList(piplus,  "PionAllPlus", pidalg);
		theAnalysis->FillList(piminus, "PionAllMinus", pidalg);
		
		// ***
		// *** TRUE PID combinatorics
		// ***
		
		// *** do MC truth match for PID type
		SelectTruePid(theAnalysis, muplus);
		SelectTruePid(theAnalysis, muminus);
		SelectTruePid(theAnalysis, piplus);
		SelectTruePid(theAnalysis, piminus);
				
		// *****
		// *** combinatorics for J/psi -> mu+ mu-
		// *****
		
		jpsi.Combine(muplus, muminus);		
		jpsi.SetType(443);
		int njmct = theAnalysis->McTruthMatch(jpsi);
				
		for (j=0;j<jpsi.GetLength();++j) 
		{
		    // 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);
			
			PndKinFitter mfitter(jpsi[j]);		// instantiate the PndKinFitter in psi(2S)
			mfitter.AddMassConstraint(m0_jpsi);	// add the mass constraint
			mfitter.Fit();						// do fit
			
			Float_t  chi2_m = mfitter.GetChi2();	// get chi2 of fit
			Float_t  prob_m = mfitter.GetProb();	// access probability of fit
			Int_t  ndf_m = mfitter.GetNdf();
			RhoCandidate *jfit = jpsi[j]->GetFit();	// access the fitted cand
			Float_t  mass_m = jfit->M();
			
			ntp1->Column("kfchi2",    (Float_t) chi2_m);
			ntp1->Column("kfprob",    (Float_t) prob_m);
			ntp1->Column("kfm",    (Float_t) mass_m);
			ntp1->Column("kfndf",    (Int_t) ndf_m);
			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);
			
			Float_t  chi2_m = kinfit.GetChi2();	// get chi2 of fit
			Float_t  prob_m = kinfit.GetProb();	// access probability of fit
			Int_t  ndf_m = kinfit.GetNdf();
			Float_t  mass_m = psifit->M();
			
			ntp2->Column("kfchi2",    (Float_t) chi2_m);
			ntp2->Column("kfprob",    (Float_t) prob_m);
			ntp2->Column("kfm",    (Float_t) mass_m);
			ntp2->Column("kfndf",    (Int_t) ndf_m);
			ntp2->DumpData();
		}			
		
		muplus.Cleanup();
		muminus.Cleanup();
		piplus.Cleanup();
		piminus.Cleanup();
		jpsi.Cleanup();
		psi2s.Cleanup();
		all.Cleanup();
		mclist.Cleanup();
		
	}
	
	// *** write out all the histos
	out->cd();
	
	ntp1->GetInternalTree()->Write();
	ntp2->GetInternalTree()->Write();
	nmc->GetInternalTree()->Write();
		
	out->Save();
	
}