// ************************************************************************
//
// Template for an Analysis Task, 
// 
// for the December 2013 CM Tutorial 
//
// For further info see also
//
// http://panda-wiki.gsi.de/cgi-bin/viewauth/Computing/PandaRootRhoTutorial
// http://panda-wiki.gsi.de/cgi-bin/view/Computing/PandaRootAnalysisJuly13
//
// K.Goetzen 11/2013
//
// ************************************************************************


// The header file
#include "PndTutAnaTaskD0.h"

// C++ headers
#include <string>
#include <iostream>

// FAIR headers
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"

// ROOT headers
#include "TClonesArray.h"
#include "TVector3.h"
#include "TH1F.h"
#include "TH2F.h"

// RHO headers
#include "RhoCandidate.h"
#include "RhoHistogram/RhoTuple.h"
#include "RhoFactory.h"
#include "RhoMassParticleSelector.h"
#include "RhoMomentumParticleSelector.h"
#include "RhoEnergyParticleSelector.h"
#include "RhoTuple.h"

// Analysis headers
#include "PndAnalysis.h"
#include "Pnd4CFitter.h"
#include "PndKinVtxFitter.h"
#include "PndChiVtxFitter.h"
#include "PndVtxPRG.h"
#include "PndKinFitter.h"
#include "PndVtxPoca.h"
#include "PndPidCandidate.h"
#include "PndEventShape.h"
		
		
using std::cout;
using std::endl;


// -----   Default constructor   -------------------------------------------
PndTutAnaTaskD0::PndTutAnaTaskD0(double pbarmom) :
  FairTask("Panda Tutorial Analysis Task") 
{ 
	double mp=0.938272;
	fIni.SetXYZT(0,0,pbarmom, sqrt(pbarmom*pbarmom+mp*mp)+mp);
}
// -------------------------------------------------------------------------


// -----   Destructor   ----------------------------------------------------
PndTutAnaTaskD0::~PndTutAnaTaskD0() { }
// -------------------------------------------------------------------------


// -----   Public method Init   --------------------------------------------
InitStatus PndTutAnaTaskD0::Init() 
{		
	// *** initialize PndAnalysis object
	fAnalysis = new PndAnalysis();
	
	// *** reset the event counter
	fEvtCount = 0;

	// *******
	// ******* PREPARE/CREATE THE STUFF YOU NEED
	// *******
	fD0Sel 	 = new RhoMassParticleSelector("D0Sel",   1.864, 1.0);
	nd0 	 = new RhoTuple("nd0","my D0 tuple");
	fVtxPoca = new PndVtxPoca();
	
	return kSUCCESS;
}

// -------------------------------------------------------------------------
	
void PndTutAnaTaskD0::SetParContainers() 
{
  // *** Get run and runtime database
  FairRun* run = FairRun::Instance();
  if ( ! run ) Fatal("SetParContainers", "No analysis run");
}

// -------------------------------------------------------------------------


// -----   Public method Exec   --------------------------------------------
void PndTutAnaTaskD0::Exec(Option_t*)
{
	// *** some variables
	int i=0,j=0, k=0, l=0;
	
	// *** necessary to read the next event
	fAnalysis->GetEvent();
	
	// *** print event counter
	if (!(++fEvtCount%100)) cout << "evt "<<fEvtCount<<endl;
	
	// *******
	// ******* PUT ANALYSIS CODE HERE
	// *******
	RhoCandList kp, km, pip, pim, d0, d0b;
	TString pidsel = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts";
	
	fAnalysis->FillList( kp,  "KaonAllPlus",  pidsel);
	fAnalysis->FillList( km,  "KaonAllMinus", pidsel);
	fAnalysis->FillList( pip, "PionAllPlus",  pidsel);
	fAnalysis->FillList( pim, "PionAllMinus", pidsel);

	d0.Combine(km, pip);
	d0.SetType(421);
	d0b.Combine(kp, pim);
	d0b.SetType(-421);
	d0.Append(d0b);
	
	d0.Select(fD0Sel);
	
	for (j=0;j<d0.GetLength();++j)
	{
		RhoCandidate *dka = d0[j]->Daughter(0);
		RhoCandidate *dpi = d0[j]->Daughter(1);
		
		PndPidCandidate *dkar = (PndPidCandidate*)dka->GetRecoCandidate();
		PndPidCandidate *dpir = (PndPidCandidate*)dpi->GetRecoCandidate();

		TLorentzVector d0cm = d0[j]->P4();
		TLorentzVector kacm = dka->P4();
		TLorentzVector picm = dpi->P4();
		d0cm.Boost(-fIni.BoostVector());
		kacm.Boost(-fIni.BoostVector());
		picm.Boost(-fIni.BoostVector());
				
		bool mct = fAnalysis->McTruthMatch(d0[j]);
		
		TVector3 vtx;
		double qavtx = fVtxPoca->GetPocaVtx(vtx, d0[j]);
	
		nd0->Column("ev", 		(Float_t) i);
		nd0->Column("cand",		(Float_t) j);
		nd0->Column("pdg",		(Float_t) d0[j]->PdgCode());
		nd0->Column("mct",		(Float_t) mct);
		
		nd0->Column("d0m",		(Float_t) d0[j]->M());
		nd0->Column("d0p",		(Float_t) d0[j]->P());
		nd0->Column("d0tht",	(Float_t) d0[j]->P3().Theta());
		nd0->Column("d0phi",	(Float_t) d0[j]->P3().Phi());
		nd0->Column("d0pt",		(Float_t) d0[j]->P3().Pt());
		nd0->Column("d0pcm",	(Float_t) d0cm.P());
		nd0->Column("d0thtcm",	(Float_t) d0cm.Theta());
		nd0->Column("d0miss",	(Float_t) (fIni-(d0[j]->P4())).M());
		
		nd0->Column("d0vx",		(Float_t) vtx.X());
		nd0->Column("d0vy",		(Float_t) vtx.Y());
		nd0->Column("d0vz",		(Float_t) vtx.Z());
		nd0->Column("d0vqa",	(Float_t) qavtx);
		
		nd0->Column("kch",		(Float_t) dka->Charge());
		nd0->Column("kp",		(Float_t) dka->P());
		nd0->Column("ktht",		(Float_t) dka->P3().Theta());
		nd0->Column("kphi",		(Float_t) dka->P3().Phi());
		nd0->Column("kpt",		(Float_t) dka->P3().Pt());
		nd0->Column("kpid",		(Float_t) dka->GetPidInfo(3));
		nd0->Column("kapcm",	(Float_t) kacm.P());
		nd0->Column("kathtcm",	(Float_t) kacm.Theta());
		
		nd0->Column("pich",		(Float_t) dpi->Charge());
		nd0->Column("pip",		(Float_t) dpi->P());
		nd0->Column("pitht",	(Float_t) dpi->P3().Theta());
		nd0->Column("piphi",	(Float_t) dpi->P3().Phi());
		nd0->Column("pit",		(Float_t) dpi->P3().Pt());
		nd0->Column("pipid",	(Float_t) dpi->GetPidInfo(2));
		nd0->Column("pipcm",	(Float_t) picm.P());
		nd0->Column("pithtcm",	(Float_t) picm.Theta());
		
		nd0->DumpData();
	}
		
}


void PndTutAnaTaskD0::Finish()
{
	
	// *******
	// ******* STORE YOUR HISTOS AND TUPLES
	// *******
	nd0->GetInternalTree()->Write();
}

ClassImp(PndTutAnaTaskD0)