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


void ana_pionangle(int nevts=0)
{
	TString OutFile="output_macrotest.root";

	// *** some variables
	int ievts=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);
	// load Root classes
	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_test.root","RECREATE");

	// ** createing some histograms
	TH1F *pionangle = new TH1F("pionangle","angle of pions", 360, 0, 20);
	TH1F *pionanglemc = new TH1F("pionanglemc","angle of mc pions", 360, 0, 20);
	pionanglemc->SetLineColor(4);
	TH2F *pionanglemom = new TH2F("pionanglemom", "angle to momentum (true)", 180,0,180,50,0,10);
	TH2F *pionanglemom2 = new TH2F("pionanglemom2", "angle to momentum (false)", 180,0,180,50,0,10);

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

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

	TCandList mctrk, pions, comb;
	PndMcTruthMatch mcm;

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

		Float_t sroot, tan;
		for(Int_t i=0;i<pions.GetLength();i++) {
			if (!mcm.MctMatch(pions[i], mctrk)){
//			if (pions[i].PdgCode()!=211 && pions[i].PdgCode()!=-211) {
//				if(mctrk[pions[i].GetMcIdx()].GetMcMotherIdx()<0 || mctrk[mctrk[pions[i].GetMcIdx()].GetMcMotherIdx()].PdgCode()==88888){
//					sroot=TMath::Sqrt(pions[i].Px()*pions[i].Px()+pions[i].Py()*pions[i].Py());
//					tan=TMath::ATan2(sroot,pions[i].Pz());
//	//				pionangle->Fill(tan*180/(TMath::Pi()));
////					if(tan*180/(TMath::Pi())<20)pionangle->Fill(pions[i].GetMomentum().Mag());
//					pionanglemom->Fill(tan*180/(TMath::Pi()),pions[i].GetMomentum().Mag());
//				} else {
//					sroot=TMath::Sqrt(pions[i].Px()*pions[i].Px()+pions[i].Py()*pions[i].Py());
//					tan=TMath::ATan2(sroot,pions[i].Pz());
//					pionanglemom2->Fill(tan*180/(TMath::Pi()),pions[i].GetMomentum().Mag());
//				}
//			} else {
//			sroot=TMath::Sqrt(pions[i].Px()*pions[i].Px()+pions[i].Py()*pions[i].Py());
//			tan=TMath::ATan2(sroot,pions[i].Pz());
//			pionanglemom2->Fill(tan*180/(TMath::Pi()),pions[i].GetMomentum().Mag());
//			}
//			if(mcm.MctMatch(pions[i],mctrk)) {
				cout << "McIdx: " << pions[i].GetMcIdx() << " pdg: " << pions[i].PdgCode() << endl;
				cout << "McPdg: " << mctrk[pions[i].GetMcIdx()].PdgCode() << endl << endl;
			}
		}
		// McTruth particles
		for(Int_t i=0;i<mctrk.GetLength();i++) {
			if(mctrk[i].PdgCode()==211 || mctrk[i].PdgCode()==-211){
				if(mctrk[i].GetMcMotherIdx()<0 || mctrk[mctrk[i].GetMcMotherIdx()].PdgCode()==88888) {
					sroot=TMath::Sqrt(mctrk[i].Px()*mctrk[i].Px()+mctrk[i].Py()*mctrk[i].Py());
					tan=TMath::ATan2(sroot,mctrk[i].Pz());
	//				pionanglemc->Fill(tan*180/(TMath::Pi()));
					if(tan*180/(TMath::Pi())<20)pionanglemc->Fill(mctrk[i].GetMomentum().Mag());
				}
			}
		}

	} // end of event loop
	out->cd();
//	pionanglemc->Draw();
//	pionangle->Draw("same");
	TCanvas *ch1 = new TCanvas("ch1","sg1");
	pionanglemom->Draw("COLZ");
	TCanvas *ch2 = new TCanvas("ch2","sg2");
	ch2->cd();
	pionanglemom2->Draw("COLZ");
}
