#include <algorithm>
void ana(TString InFile = "test_fast.root", int nevts = 0, double pbarmom = 6.5) {
	// *** 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", "Dp Analysis");
	RhoTuple *ntp2 = new RhoTuple("ntp2", "Dm 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 kminus, piplus, piplus2, kplus, piminus, piminus2, dpluslist, dminuslist, all, mclist;
	
	// *** Mass selector for the D+
	double m0_dplus = TDatabasePDG::Instance()->GetParticle("D+")->Mass();   // Get nominal PDG mass of D+
	std::cout << "mass d0 = " << m0_dplus << std::endl;
	RhoMassParticleSelector * dPlusMassSelector = new RhoMassParticleSelector("dp", m0_dplus, 0.3);
	
	// *** the lorentz vector of the initial p
	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(kplus,  "KaonBestPlus", pidalg);
		theAnalysis->FillList(kminus, "KaonBestMinus", pidalg);
		theAnalysis->FillList(piplus,  "PionBestPlus", pidalg);
		theAnalysis->FillList(piminus, "PionBestMinus", pidalg);
		
		
		// *****
		// *** combinatorics for D+ -> K- Pi+ Pi+ and the antiparticle's process
		// *****
		
		// dpluslist.Cleanup();
		dpluslist.Combine(kminus, piplus, piplus);
		dpluslist.SetType(411);
		int ndplusmct = theAnalysis->McTruthMatch(dpluslist);

		// dminuslist.Cleanup();
		dminuslist.Combine(kplus, piminus, piminus);
		dminuslist.SetType(-411);
		int ndminusmct = theAnalysis->McTruthMatch(dminuslist);
				
		// D+
		for (j = 0; j < dpluslist.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("evI",		 (Int_t) i);
			ntp1->Column("candI",	 (Int_t) j);
			ntp1->Column("ncandDplus",  (Float_t) dpluslist.GetLength());
			ntp1->Column("ncandDplusI",  (Int_t) dpluslist.GetLength());
			ntp1->Column("nmctDplus",   (Float_t) ndplusmct);
			ntp1->Column("nmctDplusI",   (Int_t) ndplusmct);
			
			// store info about initial 4-vector
			qa.qaP4("beam", ini, ntp1);
			
			// store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
			qa.qaComp("Dplus", dpluslist[j], ntp1);
			
			// store info about event shapes
			qa.qaEventShapeShort("es", &evsh, ntp1);

			PndKinFitter kinfit(dpluslist[j]);  // should actually be a vertex fit
			kinfit.AddMassConstraint(m0_dplus);
			kinfit.Fit();
			// std::cout << "Fit Chi2/NDF: " << kinfit.GetChi2() << "/" << kinfit.GetNdf() << ", Fit Prob: " << kinfit.GetProb() << std::endl;
			RhoCandidate * dplusfit = dpluslist[j]->GetFit();
			qa.qaComp("fDplus", dplusfit, ntp1);
			ntp1->Column("fDpluschi2", (Float_t) kinfit.GetChi2());
			ntp1->Column("fDplusndf", (Float_t) kinfit.GetNdf());
			ntp1->Column("fDplusprob", (Float_t) kinfit.GetProb());

			qa.qaMcDiff("Dplus", dpluslist[j], ntp1);
			
			// store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
			RhoCandidate *truth = dpluslist[j]->GetMcTruth();		
			TLorentzVector lv;
			if (truth) lv = truth->P4();
			qa.qaP4("trDplus", lv, ntp1);
			
			ntp1->DumpData();
		}
		// D-
		for (j = 0; j < dminuslist.GetLength(); ++j) {
			ntp2->Column("ev",		 (Float_t) i);
			ntp2->Column("cand",	 (Float_t) j);
			ntp2->Column("evI",		 (Int_t) i);
			ntp2->Column("candI",	 (Int_t) j);
			ntp2->Column("ncandDminus", (Float_t) dminuslist.GetLength());
			ntp2->Column("ncandDminusI", (Int_t) dminuslist.GetLength());
			ntp2->Column("nmctDminus",  (Float_t) ndminusmct);
			ntp2->Column("nmctDminusI",  (Int_t) ndminusmct);

			qa.qaP4("beam", ini, ntp2);

			qa.qaComp("Dminus", dminuslist[j], ntp2);
			qa.qaEventShapeShort("es", &evsh, ntp2);

			RhoCandidate *truth = dminuslist[j]->GetMcTruth();		
			TLorentzVector lv;
			if (truth) lv = truth->P4();
			qa.qaP4("trDminus", lv, ntp2);

			ntp2->DumpData();
		}				
	}
	
	// *** write out all the histos
	out->cd();
	
	ntp1->GetInternalTree()->Write();
	ntp2->GetInternalTree()->Write();
	nmc->GetInternalTree()->Write();
		
	out->Save();
}