#include "BlastWaveLongitudinal.h"
#include <iostream>

namespace BlastWaveLongitudinalNamespace {

#include "../NumericalIntegration.h"

}

BlastWaveLongitudinal::BlastWaveLongitudinal(double mass, int PDGID, bool UseAcc, double ymin, double ymax, double ycm, double T):
		xlag32(), wlag32(), xleg32(), wleg32(), xlegeta(), wlegeta(),
		fMass(mass), fPDGID(PDGID), fUseAcceptance(UseAcc), fYmin(ymin), fYmax(ymax), fYcm(ycm), fT(T),
		fY2Min(), fY2Max(), fEtaMin(), fEtaMax(), fAcceptance(), fReconstructionEfficiency(), fetaay2(), fNormeta(), fNormeta4pi()
//		fNorm(1.), fY2Min(), fY2Max(), fEtaMin(), fEtaMax(), fAcceptance(), fReconstructionEfficiency(), fetaay2(), fNormeta(), fNormeta4pi()
{
	BlastWaveLongitudinalNamespace::GetCoefsIntegrateLegendre32(0., 3., xlag32, wlag32);
	BlastWaveLongitudinalNamespace::GetCoefsIntegrateLegendre32(fYmin, fYmax, xleg32, wleg32);
	{
	    TString work = getenv("VMCWORKDIR");
		TString dir  = work + "/KF/KFModelParameters/common/";
		char spdg[300];
		sprintf(spdg, "%d", fPDGID);
		if (!Acceptance::ReadAcceptanceFunction(fAcceptance, dir + "pty_acc_" + spdg + ".txt")) {
			fUseAcceptance = false;
			return;
		}
		fReconstructionEfficiency = Acceptance::ReconstructionEfficiencyFunction(0.99, 0.98, 0.135);
		std::vector<double> etavec(0), y2acc(0), norm(0), norm4pi(0);
		double deta = 0.01;
		double etamax = ycm;
		for(int i=0;;++i) {
		  double tmpeta = 1.e-5 + deta*(i);
		  if (i==0) { fY2Min = y2Av(tmpeta); fEtaMin = tmpeta; }
		  if (tmpeta>etamax) { fY2Max = y2Av(tmpeta - deta); fEtaMax = tmpeta - deta; break; }
		  etavec.push_back(tmpeta);
		  y2acc.push_back(y2Av(tmpeta));
		  norm.push_back(Normalization(tmpeta));
		  norm4pi.push_back(Normalization4pi(tmpeta));
		  // std::cout << "Calc: " << etavec[i] << " " << y2acc[i] << " " << norm[i] << "\n";
		}
		fetaay2 = TSpline3("fetaay2acc", &y2acc[0], &etavec[0], etavec.size());
		fNormeta = TSpline3("fNormacc", &etavec[0], &norm[0], etavec.size());
		fNormeta4pi = TSpline3("fNorm4pi", &etavec[0], &norm4pi[0], etavec.size());
		for(unsigned int i=0;i<etavec.size();++i) {
		  // std::cout << "ACalc: " << fetaay2.Eval(y2acc[i]) << " " << y2acc[i] << " " << fNormeta.Eval(etavec[i]) << "\n";
		//  std::cout << "Avs: " << mtAv(Tvec[i]) << " " << mtAv2(Tvec[i]) << "\n";
		}
	}
}

double BlastWaveLongitudinal::y2Av(double eta) {
		BlastWaveLongitudinalNamespace::GetCoefsIntegrateLegendre32(-eta, eta, xlegeta, wlegeta);
		double ret1 = 0., ret2 = 0.;
		for(Int_t ie = 0 ; ie < 32 ; ie++){
			for(Int_t i = 0 ; i < 32 ; i++){
				for(Int_t j = 0 ; j < 32 ; j++){
					double tmpf = 0.;
					double ty = xleg32[j] - xlegeta[ie];
					double tp = sqrt(xlag32[i]*xlag32[i] + fMass*fMass) * TMath::CosH(xleg32[j] + fYcm);
					tp = sqrt(tp*tp - fMass*fMass);
					double tmt = sqrt(xlag32[i]*xlag32[i] + fMass*fMass);
					tmpf = wlegeta[ie] * wlag32[i] * wleg32[j] * xlag32[i] * tmt * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / fT);
					if (fUseAcceptance)	tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i]) * fReconstructionEfficiency.f(tp);
					ret1 += tmpf;
					ret2 += tmpf * xleg32[j] * xleg32[j];
				}
			}
		}
	    return ret2 / ret1;
}

double BlastWaveLongitudinal::y2Av2(double eta) {
		double ret1 = 0., ret2 = 0.;
		for(Int_t j = 0 ; j < 32 ; j++) {
			ret1 += wleg32[j] * xleg32[j] * xleg32[j] * fy(xleg32[j], eta);
			ret2 += wleg32[j] * fy(xleg32[j], eta);
			}
		return ret1 / ret2;
}

double BlastWaveLongitudinal::fy(double y, double eta) {
		double ret = 0.;
		BlastWaveLongitudinalNamespace::GetCoefsIntegrateLegendre32(-eta, eta, xlegeta, wlegeta);
		for(Int_t ie = 0 ; ie < 32 ; ie++){
			for(Int_t i = 0 ; i < 32 ; i++){
				double tmpf = 0.;
				double ty = y - xlegeta[ie];
				double tmt = sqrt(xlag32[i]*xlag32[i] + fMass*fMass);
				double tp = tmt * TMath::CosH(y + fYcm);
				tp = sqrt(tp*tp - fMass*fMass);
				tmpf = wlegeta[ie] * wlag32[i] * xlag32[i] * tmt * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / fT);
				if (fUseAcceptance)	tmpf *= fAcceptance.getAcceptance(y + fYcm, xlag32[i]) * fReconstructionEfficiency.f(tp);
				ret += tmpf;
			}
		}
		return ret / fNormeta.Eval(eta);
}

double BlastWaveLongitudinal::Normalization(double eta) {
	double ret = 0.;
	//if (fUseAcceptance) 
	{
		BlastWaveLongitudinalNamespace::GetCoefsIntegrateLegendre32(-eta, eta, xlegeta, wlegeta);
		for(Int_t ie = 0 ; ie < 32 ; ie++){
			for(Int_t i = 0 ; i < 32 ; i++){
				for(Int_t j = 0 ; j < 32 ; j++){
					double tmpf = 0.;
					double ty = xleg32[j] - xlegeta[ie];
					double tp = sqrt(xlag32[i]*xlag32[i] + fMass*fMass) * TMath::CosH(xleg32[j] + fYcm);
					tp = sqrt(tp*tp - fMass*fMass);
					double tmt = sqrt(xlag32[i]*xlag32[i] + fMass*fMass);
					tmpf = wlegeta[ie] * wlag32[i] * wleg32[j] * xlag32[i] * tmt * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / fT);
					if (fUseAcceptance)	tmpf *= fAcceptance.getAcceptance(xleg32[j] + fYcm, xlag32[i]) * fReconstructionEfficiency.f(tp);
					ret += tmpf;
				}
			}
		}
	}
	return ret;
}

double BlastWaveLongitudinal::Normalization4pi(double eta) {
	double ret = 0.;
	//if (fUseAcceptance) 
	{
		BlastWaveLongitudinalNamespace::GetCoefsIntegrateLegendre32(-eta, eta, xlegeta, wlegeta);
		for(Int_t ie = 0 ; ie < 32 ; ie++){
			for(Int_t i = 0 ; i < 32 ; i++){
				for(Int_t j = 0 ; j < 32 ; j++){
					double tmpf = 0.;
					double ty = xleg32[j] - xlegeta[ie];
					double tp = sqrt(xlag32[i]*xlag32[i] + fMass*fMass) * TMath::CosH(xleg32[j] + fYcm);
					tp = sqrt(tp*tp - fMass*fMass);
					double tmt = sqrt(xlag32[i]*xlag32[i] + fMass*fMass);
					tmpf = wlegeta[ie] * wlag32[i] * wleg32[j] * xlag32[i] * tmt * TMath::CosH(ty) * TMath::Exp(-tmt * TMath::CosH(ty) / fT);
					ret += tmpf;
				}
			}
		}
	}
	return ret;
}