/**
 * @author Elena Lebedeva <e.lebedeva@gsi.de>
 * @since 2015
 * @version 1.0
 **/


#ifndef CBM_LMVM_UTILS_H
#define CBM_LMVM_UTILS_H

#include "CbmMCTrack.h"
#include "CbmLmvmCandidate.h"
#include "CbmStsTrack.h"
#include "CbmKFVertex.h"
#include "CbmL1PFFitter.h"
#include "TMCProcess.h"

class CbmLmvmUtils{
public:

	/*
	 * Calculates and set track parameters to CbmLmvmCandidate.
	 * The following parameters are set: fChi2sts, fChi2Prim, fPosition, fMomentum, fMass, fCharge, fEnergy, fRapidity
	 */
	static void CalculateAndSetTrackParamsToCandidate(
		CbmLmvmCandidate* cand,
		CbmStsTrack* stsTrack,
		CbmKFVertex& kfVertex)
	{
		CbmL1PFFitter fPFFitter;
		vector<CbmStsTrack> stsTracks;
		stsTracks.resize(1);
		stsTracks[0] = *stsTrack;
		vector<L1FieldRegion> vField;
		vector<float> chiPrim;
		fPFFitter.GetChiToVertex(stsTracks, vField, chiPrim, kfVertex, 3e6);
		cand->fChi2sts = stsTracks[0].GetChiSq() / stsTracks[0].GetNDF();
		cand->fChi2Prim = chiPrim[0];
		const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();

		vtxTrack->Position(cand->fPosition);
		vtxTrack->Momentum(cand->fMomentum);

		cand->fMass = TDatabasePDG::Instance()->GetParticle(11)->Mass();
		cand->fCharge = (vtxTrack->GetQp() > 0) ?1 :-1;
		cand->fEnergy = sqrt(cand->fMomentum.Mag2() + cand->fMass * cand->fMass);
		cand->fRapidity = 0.5*TMath::Log((cand->fEnergy + cand->fMomentum.Z()) / (cand->fEnergy - cand->fMomentum.Z()));
	}

	/*
	 * Armenteros - Podolansky plot
	 */
    static void CalculateArmPodParams(
    	CbmLmvmCandidate* cand1,
		CbmLmvmCandidate* cand2,
        Double_t &alpha,
        Double_t &ptt)
    {
		alpha = ptt = 0.;
		Double_t spx = cand1->fMomentum.X() + cand2->fMomentum.X();
		Double_t spy = cand1->fMomentum.Y() + cand2->fMomentum.Y();
		Double_t spz = cand1->fMomentum.Z() + cand2->fMomentum.Z();
		Double_t sp  = sqrt(spx*spx + spy*spy + spz*spz);

		if( sp == 0.0) return;
		Double_t pn, /*pp,*/ pln, plp;
		if(cand1->fCharge < 0.) {
			pn = cand1->fMomentum.Mag();
			//pp = cand2->fMomentum.Mag();
			pln  = (cand1->fMomentum.X()*spx + cand1->fMomentum.Y()*spy + cand1->fMomentum.Z()*spz)/sp;
			plp  = (cand2->fMomentum.X()*spx + cand2->fMomentum.Y()*spy + cand2->fMomentum.Z()*spz)/sp;
		}else{
			pn = cand2->fMomentum.Mag();
			//pp = cand1->fMomentum.Mag();
			pln  = (cand2->fMomentum.X()*spx + cand2->fMomentum.Y()*spy + cand2->fMomentum.Z()*spz)/sp;
			plp  = (cand1->fMomentum.X()*spx + cand1->fMomentum.Y()*spy + cand1->fMomentum.Z()*spz)/sp;
		}
		if( pn == 0.0) return;
		Double_t ptm  = (1.-((pln/pn)*(pln/pn)));
		ptt= (ptm>=0.)?  pn*sqrt(ptm) :0;
		alpha = (plp-pln)/(plp+pln);
    }

    /*
	 * \brief Return true if MC track is signal primary electron.
	 */
	static Bool_t IsMcSignalElectron(
			CbmMCTrack* mctrack)
	{
		if (mctrack == NULL) return false;
		Int_t pdg = TMath::Abs(mctrack->GetPdgCode());
		if (mctrack->GetGeantProcessId() == kPPrimary && pdg == 11) return true;
		return false;
	}

	/*
	 * \brief Return true if MC track is electron from gamma conversion.
	 */
	static Bool_t IsMcGammaElectron(
			CbmMCTrack* mctrack,
			TClonesArray* mcTracks)
	{
		if (mctrack == NULL) return false;
		Int_t pdg = TMath::Abs(mctrack->GetPdgCode());
		if (pdg != 11) return false;
		Int_t motherId = mctrack->GetMotherId();
		if (motherId < 0){
			return false;
		} else {
			CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(mcTracks->At(motherId));
			Int_t motherPdg = mct1->GetPdgCode();
			if (mct1 != NULL && motherPdg == 22 && pdg == 11){
				return true;
			}
		}
		return false;
	}

	/*
	 * \brief Return true if MC track is electron from Pi0 dalitz decay.
	 */
	static Bool_t IsMcPi0Electron(
			CbmMCTrack* mctrack,
			TClonesArray* mcTracks)
	{
		if (mctrack == NULL) return false;
		Int_t pdg = TMath::Abs(mctrack->GetPdgCode());
		if (pdg != 11) return false;
		Int_t motherId = mctrack->GetMotherId();
		if (motherId < 0){
			return false;
		} else {
			CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(mcTracks->At(motherId));
			Int_t motherPdg = mct1->GetPdgCode();
			if (mct1 != NULL && motherPdg == 111 && pdg == 11 && mctrack->GetGeantProcessId() != kPPrimary){
				return true;
			}
		}
		return false;
	}

	/*
	 * \brief Return true if MC track is electron from Eta decay.
	 */
	static Bool_t IsMcEtaElectron(
			CbmMCTrack* mctrack,
			TClonesArray* mcTracks)
	{
		if (mctrack == NULL) return false;
		Int_t pdg = TMath::Abs(mctrack->GetPdgCode());
		if (pdg != 11) return false;
		Int_t motherId = mctrack->GetMotherId();
		if (motherId < 0){
			return false;
		} else {
			CbmMCTrack* mct1 = static_cast<CbmMCTrack*>(mcTracks->At(motherId));
			Int_t motherPdg = mct1->GetPdgCode();
			if (mct1 != NULL && motherPdg == 221 && pdg == 11){
				return true;
			}
		}
		return false;
	}

};

#endif