class RhoCandList;
class RhoCandidate;
class PndAnaPidSelector;
class PndAnaPidCombiner;
class PndAnalysis;
class RhoTuple;
class PndPidCandidate;


void psi4160at8GeV10000eventsAnalysis(TString InFile="psi4160_fast.root", int nevts=0, double pbarmom = 8.00000)
{
	// *** some variables
	int i=0,j=0, k=0, l=0;
	//gStyle->SetOptFit(1011);
	
	// *** the output file for FairRunAna
	TString OutFile="psi4160_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 psi4160.root
	fRun->Init(); 
	
	// *** take constant field; needed for PocaVtx
	RhoCalculationTools::ForceConstantBz(20.0);
	
	// *** create an output file for all histograms
	TFile *out = TFile::Open("psi4160.root","RECREATE");

	// *** create some ntuples
	RhoTuple *npsi4160 = new RhoTuple("npsi4160", "psi4160 analysis");
	RhoTuple *nd0 = new RhoTuple("nd0", "d0 analysis");
	RhoTuple *nantid0 = new RhoTuple("nantid0", "antid0 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
	PndRhoTupleQA qa(theAnalysis, pbarmom); 

	// *** RhoCandLists for the analysis
	RhoCandList psi4160, d0, antid0, kplus, kminus, piplus, piminus, all, mclist;
	
	// *** Mass selector for the d0 cands
	double m0_d0 = TDatabasePDG::Instance()->GetParticle("D0")->Mass();   // Get nominal PDG mass of the D0/anti-D0
	RhoMassParticleSelector *d0MassSel=new RhoMassParticleSelector("d0",m0_d0,1.0);
	
	// *** the lorentz vector of the initial psi(4160)
	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,  "KaonLoosePlus", pidalg);
		theAnalysis->FillList(kminus, "KaonLooseMinus", pidalg);
		theAnalysis->FillList(piplus,  "PionLoosePlus", pidalg);
		theAnalysis->FillList(piminus, "PionLooseMinus", pidalg);
		
		
/*		// some mass pre selection
		d0.Select(d0MassSel);
*/			
		// *** combinatorics for d0 -> k- pi+ *** //
		d0.Combine(kminus, piplus);
		d0.SetType(421);
		int nd0mct = theAnalysis->McTruthMatch(d0);
				
		for (j=0;j<d0.GetLength();++j) 
		{
		    // some general info about event (actually written for each candidate!)
			nd0->Column("ev",	(Float_t) i);
			nd0->Column("cand",	(Float_t) j);
			nd0->Column("ncand",    (Float_t) d0.GetLength());
			nd0->Column("nmct",     (Float_t) nd0mct);
			
			// store info about initial 4-vector
			qa.qaP4("beam", ini, nd0);
			
			// store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
			qa.qaComp("d0", d0[j], nd0);
			
			// store info about event shapes
			qa.qaEventShapeShort("es",&evsh, nd0);
		
			
			// store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
			RhoCandidate *truth = d0[j]->GetMcTruth();		
			TLorentzVector lv;
			if (truth) lv = truth->P4();
			qa.qaP4("trd0", lv, nd0);
			
			nd0->DumpData();
		}
		
/*		
		// some mass pre selection
		antid0.Select(d0MassSel);		
*/
		// *** combinatorics for antid0 -> k+ pi- *** //
		antid0.Combine(kplus, piminus);
		antid0.SetType(-421);
		int nantid0mct = theAnalysis->McTruthMatch(antid0);

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

		// *** combinatorics for psi4160 -> d0 anti-d0 *** //
		
		psi4160.Combine(d0, antid0);
		psi4160.SetType(88880);
		int npsi4160mct = theAnalysis->McTruthMatch(psi4160);
		
		cout << psi4160.GetLength() << endl;

		for (j=0;j<psi4160.GetLength();++j) 
		{
		    // some general info about event (actually written for each candidate!)
			npsi4160->Column("ev",	(Float_t) i);
			npsi4160->Column("cand",	(Float_t) j);
			npsi4160->Column("ncand",   (Float_t) psi4160.GetLength());
			npsi4160->Column("nmct",    (Float_t) npsi4160mct);

			PndKinFitter *fitter= new PndKinFitter(psi4160[j]);	// instantiate the kin fitter in psi(4160)
			fitter->Add4MomConstraint(ini);	// set 4 constraint
			fitter->Fit();				
			RhoCandidate *psif = psi4160[j]->GetFit();// get fitted psi
	
			npsi4160->Column("psi4160fit",	(Float_t) psif->M());
			
			// store info about initial 4-vector
			qa.qaP4("beam", ini, npsi4160);
			
			// store information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
			qa.qaComp("psi4160", psi4160[j], npsi4160);

			// store info about event shapes
			qa.qaEventShapeShort("es",&evsh, npsi4160);
			
			// *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
			RhoCandidate *truth = psi4160[j]->GetMcTruth();
			TLorentzVector lv;
			if (truth) lv = truth->P4();
			qa.qaP4("trpsi", lv, npsi4160);
			
			npsi4160->DumpData();
			delete fitter;	
		}
			
	}
	
	// *** write out all the histos
	out->cd();
	
	nd0->GetInternalTree()->Write();
	nantid0->GetInternalTree()->Write();
	npsi4160->GetInternalTree()->Write();
	nmc->GetInternalTree()->Write();
		
	out->Save();
	
}