//Elisabetta: testing genfit2....

class RhoCandList;
class RhoCandidate;
class PndAnaPidSelector;
class PndAnaPidCombiner;
class PndAnalysis;

// *** 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 printCand(RhoCandidate *c)
{
	TLorentzVector lv=c->P4();
	
	cout <<c->PdgCode()<<" ("<<lv.X()<<"/"<<lv.Y()<<"/"<<lv.Z()<<"/"<<lv.E()<<")"<<endl;
}

void countDoubles(RhoCandList &l, int &n1, int &n2, int &n3)
{
	int n_smc  = 0;
	int n_strk = 0;
	int n_both = 0;
	double d = 0.00001;
	
	for (int i=0;i<l.GetLength()-1;++i)
	{
		for (int j=i+1;j<l.GetLength();++j)
		{
			TLorentzVector dl = l[i]->P4() - l[j]->P4();
		
			bool chkmc = (l[i]->GetMcTruth()==l[j]->GetMcTruth());
			bool chktrk = (fabs(dl.X())<d) && (fabs(dl.Y())<d) && (fabs(dl.Z())<d) && (fabs(dl.E())<d);
			if (chkmc) n_smc++;
			if (chktrk) n_strk++;
			if (chktrk && chkmc) n_both++;
		}	
	}
	n1 = n_strk;
	n2 = n_smc;
	n3 = n_both;
}

void ana_complete_test(int nevts=0)
{
  
  double pbarmom = 6.231552;

      TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",1.9,kFALSE,0.1,0,"",88888);
        TStopwatch timer;
	// *** some variables
	int i=0,j=0, k=0, l=0;
	gStyle->SetOptFit(1011);
	
	// *** the output file for FairRunAna
	TString OutFile="output.root";  
					
	// *** the files coming from the simulation
	TString inPidFile  = "pid_complete.root";    // this file contains the PndPidCandidates and McTruth
	TString inParFile  = "simparams.root";
	
	// *** PID table with selection thresholds; can be modified by the user
	TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par";	
	
	// *** initialization
	FairLogger::GetLogger()->SetLogToFile(kFALSE);
	FairRunAna* fRun = new FairRunAna();
	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
	fRun->SetInputFile(inPidFile);
	
	// *** setup parameter database 	
	FairParRootFileIo* parIO = new FairParRootFileIo();
	parIO->open(inParFile);
	FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo();
	parIOPid->open(pidParFile.Data(),"in");
	
	rtdb->setFirstInput(parIO);
	rtdb->setSecondInput(parIOPid);
	rtdb->setOutput(parIO);  
	
	fRun->SetOutputFile(OutFile);
	fRun->Init(); 
	
        // *** create an output file for all histograms
	TFile *out = TFile::Open("output_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();
	
	PndRhoTupleQA qa(theAnalysis, pbarmom); 

	// *** RhoCandLists for the analysis
	RhoCandList mclist, chrg, muplus, muminus, piplus, piminus, jpsi, psi2s;
	
	// *** 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)
	TLorentzVector ini(0, 0, 6.231552, 7.240065);
	
	// ***
	// the event loop
	// ***
	
	int cntdbltrk=0, cntdblmc=0, cntdblboth=0, cnttrk=0, cnt_dbl_jpsi=0, cnt_dbl_psip=0;
	
	while (theAnalysis->GetEvent() && i++<nevts)
	{
		if ((i%100)==0) cout<<"evt " << i << endl;
				
		// *** Select with no PID info ('All'); type and mass are set 		

		
		theAnalysis->FillList(mclist,  "McTruth");


		
		
		nmc->Column("ev", (Int_t) i);
		qa.qaMcList("",   mclist, nmc);
		nmc->DumpData();
		


		theAnalysis->FillList(chrg,    "Charged");
		theAnalysis->FillList(muplus,  "MuonAllPlus");
		theAnalysis->FillList(muminus, "MuonAllMinus");
		theAnalysis->FillList(piplus,  "PionAllPlus");
		theAnalysis->FillList(piminus, "PionAllMinus");

		

		// *** combinatorics for J/psi -> mu+ mu-
		jpsi.Combine(muplus, muminus);
		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("jpsi_", 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("true_jpsi_", lv, ntp1);
			
			ntp1->DumpData();
		}
		
			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("fit_to_psi_",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("true_psi_", lv, ntp2);
			
			ntp2->DumpData();
		}			
	}

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