class TCandList;
class TCandidate;
class TFitParams;
class PndAnaPidSelector;
class PndAnaPidCombiner;


void anatut_psi2s(int nevts=0)
{
	TString OutFile="output.root";

	// *** some variables
	int ievts=0, pionskicked=0, mcpions=0, detectpions=0;
	const double val=15.0; // will be pbar momentum

	//val is the momentum of the pbar beam
	const double mp=0.93827;
	double P = val;
	double E = mp+sqrt(P*P+mp*mp);
	TLorentzVector ini(0, 0, P, E);  // *** the lorentz vector of the initial pbar

// *** 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";

gROOT->SetStyle("Plain");

gStyle->SetOptFit(1011);
gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
// make StatBox in histograms show more useful information
gStyle->SetOptStat(1111110); // with integral
//gStyle->SetOptStat(111110); // without integral

FairLogger::GetLogger()->SetLogToFile(kFALSE);

// *** initialization
FairRunAna* fRun = new FairRunAna();
FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
// *** add signal files
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");

	// *** creating tuple
	TNtuple *truetuple = new TNtuple("truetuple","truetuple","rm:pocadistance:pocaVtxX:pocaVtxY:pocaVtxZ:radial:massdifftrue:PRGvtxX:PRGvtxY:PRGvtxZ:PRGradial");
	TNtuple *combtuple = new TNtuple("combtuple","combtuple","rm:pocadistance:pocaVtxX:pocaVtxY:pocaVtxZ:radial:PRGvtxX:PRGvtxY:PRGvtxZ:PRGradial");
	TNtuple *pidtuple = new TNtuple("pidtuple","pidtuple","rm:pocadistance:pocaVtxX:pocaVtxY:pocaVtxZ:radial:PRGvtxX:PRGvtxY:PRGvtxZ:PRGradial");
	TNtuple *truemctuple = new TNtuple("truemctuple","truemctuple","rm:pocadistance:pocaVtxX:pocaVtxY:pocaVtxZ:radial:massdifftrue:PRGvtxX:PRGvtxY:PRGvtxZ:PRGradial");

	// *** the data reader object
	PndAnalysis* theAnalysis = new PndAnalysis("SttMvdGemGenTrack","FtsIdealGenTrack");
	if (nevts==0) nevts= theAnalysis->GetEntries();

	// *** PID algorithm
	TString pidalgos("PidAlgoStt;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
//	TString pidalgos("PidAlgoMvd;PidAlgoStt;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");

	// *** TCandLists for the analysis
	TCandList mctrk, piplus, piminus, comball, eplus, eminus, muplus, muminus, kplus,kminus;
	TCandList pidpiplus, pidpiminus, combpid;
	TCandList truepiplus, truepiminus, combtrue;
	TCandList mcpiplus, mcpiminus, truemcpiplus, truemcpiminus, combmctrue;

	// *** the MC truth matcher variable
	PndMcTruthMatch mcm;

	// ***
	// the event loop
	// ***
	while (theAnalysis->GetEvent() && ievts++<nevts)
	{
		if ((ievts%100)==0)
			cout << "Evt: " << ievts << endl;

		// *** the MC Truth objects
		theAnalysis->FillList(mctrk,"McTruth");

		// *** Select with PID info; type and mass are set
		theAnalysis->FillList(eplus, "ElectronVeryTightPlus",pidalgos);
		theAnalysis->FillList(eminus, "ElectronVeryTightMinus",pidalgos);
		theAnalysis->FillList(muplus, "MuonVeryTightPlus",pidalgos);
		theAnalysis->FillList(muminus, "MuonVeryTightMinus",pidalgos);
		theAnalysis->FillList(kminus, "KaonVeryTightMinus",pidalgos);
		theAnalysis->FillList(kplus, "KaonVeryTightPlus",pidalgos);
		theAnalysis->FillList(piplus, "PionVeryLoosePlus",pidalgos);
		theAnalysis->FillList(piminus, "PionVeryLooseMinus",pidalgos);
		theAnalysis->FillList(mcpiplus, "PionVeryLoosePlus");
		theAnalysis->FillList(mcpiminus, "PionVeryLoosePlus");

			detectpions+=piplus.GetLength(); detectpions+=piminus.GetLength();
			mcpions+=mcpiplus.GetLength(); mcpions+=mcpiminus.GetLength();

		// Remove electrons and muons in pi lists
		for(Int_t l=0;l<eplus.GetLength();++l) {int k=piplus.Remove(eplus[l]); if(k!=0) pionskicked+=k;}
		for(Int_t l=0;l<muplus.GetLength();++l) {int k=piplus.Remove(muplus[l]); if(k!=0) pionskicked+=k;}
		for(Int_t l=0;l<kplus.GetLength();++l) {int k=piplus.Remove(kplus[l]); if(k!=0) pionskicked+=k;}
		for(Int_t l=0;l<eminus.GetLength();++l) {int k=piminus.Remove(eminus[l]); if(k!=0) pionskicked+=k;}
		for(Int_t l=0;l<muminus.GetLength();++l) {int k=piminus.Remove(muminus[l]); if(k!=0) pionskicked+=k;}
		for(Int_t l=0;l<kminus.GetLength();++l) {int k=piminus.Remove(kminus[l]); if(k!=0) pionskicked+=k;}

		pidpiplus.Cleanup();
		pidpiminus.Cleanup();
		truepiplus.Cleanup();
		truepiminus.Cleanup();

// truth match for pi +
for (Int_t ipiplus=0;ipiplus<piplus.GetLength();++ipiplus){
	mcm.SetType(piplus,"pi+");
	if (mcm.MctMatch(piplus[ipiplus], mctrk)){
		pidpiplus.Add(piplus[ipiplus]);

		// TRUE PRIMARY PION SELECTION
		if(mctrk[mctrk[piplus[ipiplus].GetMcIdx()].GetMcMotherIdx()].PdgCode()==88888)
			truepiplus.Add(piplus[ipiplus]);  // PdgCode 88888 is pbarpSystem
	}
}
// truth match for pi -
for (Int_t ipiminus=0;ipiminus<piminus.GetLength();++ipiminus){
	mcm.SetType(piminus,"pi-");
	if (mcm.MctMatch(piminus[ipiminus], mctrk)){
		pidpiminus.Add(piminus[ipiminus]);

		// TRUE PRIMARY PION SELECTION
		if(mctrk[mctrk[piminus[ipiminus].GetMcIdx()].GetMcMotherIdx()].PdgCode()==88888)
			truepiminus.Add(piminus[ipiminus]);
	}
}

		// COMBINATION
		comball.Combine(piplus,piminus);
		combpid.Combine(pidpiplus,pidpiminus);
		combtrue.Combine(truepiplus,truepiminus);


		// ANALYSE of PRIMARY combined PIONS
		Float_t massdifftrue=10000; // mass difference between best missing mass and true mass for h_c^'
		Int_t minmassdiffIdx=0;// stores index of best combtrue found
		for (int l=0;l<combtrue.GetLength();++l) {
		// recoil mass calculation
			TLorentzVector mm = ini-combtrue[l].P4();
			Float_t rm=mm.M();
			if (fabs(rm-mctrk[1].Mass())<massdifftrue){
				massdifftrue= fabs(rm-mctrk[1].Mass());
				minmassdiffIdx=l;
			}
		} // analyse for just the one combined particle
		if(combtrue.GetLength()>0){
		for(int l=0;l<1;l++){
		// recoil mass calculation
			TLorentzVector mm = ini-combtrue[minmassdiffIdx].P4();
			Float_t rm=mm.M();

		// PndVtxPRG Vertex Fitter
			PndVtxPRG prgfitter(combtrue[minmassdiffIdx]);
			TVector3 PRGvtx;
			TMatrixD vtxcov(3,3);
			prgfitter.FitVertexFast(PRGvtx,vtxcov);
			const TCandidate *combtruefit = combtrue[minmassdiffIdx].GetFit();
			Float_t PRGvtxX=PRGvtx.X();
			Float_t PRGvtxY=PRGvtx.Y();
			Float_t PRGvtxZ=PRGvtx.Z();
			Float_t PRGradial;
			PRGradial=PRGvtxX*PRGvtxX+PRGvtxY*PRGvtxY;

		// determine pseudo vertex *** POCA
			PndVtxPoca poca(combtrue[minmassdiffIdx]);
			TVector3 pocavtx;
			Float_t pocadistance = poca.GetPocaVtx(pocavtx);
			Float_t pocaVtxX=pocavtx.X();
			Float_t pocaVtxY=pocavtx.Y();
			Float_t pocaVtxZ=pocavtx.Z();
			Float_t radial;
			radial=pocaVtxX*pocaVtxX+pocaVtxY*pocaVtxY;

		// TUPLE FILL
			Float_t combtuplearray[] = {rm, pocadistance, pocaVtxX, pocaVtxY, pocaVtxZ, radial, massdifftrue, PRGvtxX, PRGvtxY, PRGvtxZ, PRGradial};
			truetuple->Fill(combtuplearray);
		}}

		// ANALYSE of ALL combined "PIONS"
		for (int l=0;l<comball.GetLength();++l) {
		// recoil mass calculation
			TLorentzVector mm = ini-comball[l].P4();
			Float_t rm=mm.M();

		// PndVtxPRG Vertex Fitter
			PndVtxPRG prgfitter(comball[l]);
			TVector3 PRGvtx;
			TMatrixD vtxcov(3,3);
			prgfitter.FitVertexFast(PRGvtx,vtxcov);
			Float_t PRGvtxX=PRGvtx.X();
			Float_t PRGvtxY=PRGvtx.Y();
			Float_t PRGvtxZ=PRGvtx.Z();
			Float_t PRGradial;
			PRGradial=PRGvtxX*PRGvtxX+PRGvtxY*PRGvtxY;

		// determine pseudo vertex *** POCA
			PndVtxPoca poca(comball[l]);
			TVector3 pocavtx;
			Float_t pocadistance = poca.GetPocaVtx(pocavtx);
			Float_t pocaVtxX=pocavtx.X();
			Float_t pocaVtxY=pocavtx.Y();
			Float_t pocaVtxZ=pocavtx.Z();
			Float_t radial;
			radial=pocaVtxX*pocaVtxX+pocaVtxY*pocaVtxY;

			Float_t combtuplearray2[]={rm,pocadistance, pocaVtxX, pocaVtxY, pocaVtxZ, radial, PRGvtxX, PRGvtxY, PRGvtxZ, PRGradial};
			combtuple->Fill(combtuplearray2);
		}

		// combined recoil mass of mc truth (pid) pions
		for (int l=0;l<combpid.GetLength();++l) {
		// recoil mass calculation
			TLorentzVector mm = ini-combpid[l].P4();
			Float_t rm=mm.M();

		// PndVtxPRG Vertex Fitter
			PndVtxPRG prgfitter(combpid[l]);
			TVector3 PRGvtx;
			TMatrixD vtxcov(3,3);
			prgfitter.FitVertexFast(PRGvtx,vtxcov);
			Float_t PRGvtxX=PRGvtx.X();
			Float_t PRGvtxY=PRGvtx.Y();
			Float_t PRGvtxZ=PRGvtx.Z();
			Float_t PRGradial;
			PRGradial=PRGvtxX*PRGvtxX+PRGvtxY*PRGvtxY;

		// determine pseudo vertex *** POCA
			PndVtxPoca poca(combpid[l]);
			TVector3 pocavtx;
			Float_t pocadistance = poca.GetPocaVtx(pocavtx);
			Float_t pocaVtxX=pocavtx.X();
			Float_t pocaVtxY=pocavtx.Y();
			Float_t pocaVtxZ=pocavtx.Z();
			Float_t radial;
			radial=pocaVtxX*pocaVtxX+pocaVtxY*pocaVtxY;

			Float_t combtuplearray3[]={rm,pocadistance, pocaVtxX, pocaVtxY, pocaVtxZ, radial, PRGvtxX, PRGvtxY, PRGvtxZ, PRGradial};
			pidtuple->Fill(combtuplearray3);
		}

		// truth match for mc-pi +
		for (Int_t ipiplus=0;ipiplus<mcpiplus.GetLength();++ipiplus){
			mcm.SetType(mcpiplus,"pi+");
			if (mcm.MctMatch(mcpiplus[ipiplus], mctrk)){ if(mctrk[mctrk[mcpiplus[ipiplus].GetMcIdx()].GetMcMotherIdx()].PdgCode()==88888) truemcpiplus.Add(mcpiplus[ipiplus]);}
		}
		// truth match for mc-pi -
		for (Int_t ipiminus=0;ipiminus<mcpiminus.GetLength();++ipiminus){
			mcm.SetType(mcpiminus,"pi-");
			if (mcm.MctMatch(mcpiminus[ipiminus], mctrk)){ if(mctrk[mctrk[mcpiminus[ipiminus].GetMcIdx()].GetMcMotherIdx()].PdgCode()==88888) truemcpiminus.Add(mcpiminus[ipiminus]);}
		}
		combmctrue.Combine(mcpiplus,mcpiminus);
		// ANALYSE of PRIMARY combined PIONS
		massdifftrue=10000; // mass difference between best missing mass and true mass for h_c^'
		minmassdiffIdx=0;// stores index of best combtrue found
		for (int l=0;l<combmctrue.GetLength();++l) {
		// recoil mass calculation
			TLorentzVector mm = ini-combmctrue[l].P4();
			Float_t rm=mm.M();
			if (fabs(rm-mctrk[1].Mass())<massdifftrue){
				massdifftrue= fabs(rm-mctrk[1].Mass());
				minmassdiffIdx=l;
			}
		} // analyse for just the one combined particle
		if(combmctrue.GetLength()>0){
		for(int l=0;l<1;l++){
		// recoil mass calculation
			TLorentzVector mm = ini-combmctrue[minmassdiffIdx].P4();
			Float_t rm=mm.M();

		// PndVtxPRG Vertex Fitter
			PndVtxPRG prgfitter(combmctrue[minmassdiffIdx]);
			TVector3 PRGvtx;
			TMatrixD vtxcov(3,3);
			prgfitter.FitVertexFast(PRGvtx,vtxcov);
			const TCandidate *combtruefit = combmctrue[minmassdiffIdx].GetFit();
			Float_t PRGvtxX=PRGvtx.X();
			Float_t PRGvtxY=PRGvtx.Y();
			Float_t PRGvtxZ=PRGvtx.Z();
			Float_t PRGradial;
			PRGradial=PRGvtxX*PRGvtxX+PRGvtxY*PRGvtxY;

		// determine pseudo vertex *** POCA
			PndVtxPoca poca(combmctrue[minmassdiffIdx]);
			TVector3 pocavtx;
			Float_t pocadistance = poca.GetPocaVtx(pocavtx);
			Float_t pocaVtxX=pocavtx.X();
			Float_t pocaVtxY=pocavtx.Y();
			Float_t pocaVtxZ=pocavtx.Z();
			Float_t radial;
			radial=pocaVtxX*pocaVtxX+pocaVtxY*pocaVtxY;

		// TUPLE FILL
			Float_t combtuplearray4[] = {rm, pocadistance, pocaVtxX, pocaVtxY, pocaVtxZ, radial, massdifftrue, PRGvtxX, PRGvtxY, PRGvtxZ, PRGradial};
			truemctuple->Fill(combtuplearray4);
		}}
	}

	cout << "Number of MC pions:       " << mcpions << endl;
	cout << "Number of kicked pions:  - " << pionskicked << endl;
	cout << "                         ------" << endl;
	cout << "Number of detected pions: " << detectpions << endl;
	out->cd();
	truetuple->Write();
	combtuple->Write();
	pidtuple->Write();
	out->Save();
}