#include "CbmLitCheckEnergyLossMuons.h"
#include "CbmLitMaterialInfo.h"
#include "CbmLitMaterialEffectsImp.h"
#include "CbmLitTrackParam.h"
#include "CbmLitUtils.h"
#include "CbmLitDrawHist.h"

#include "TCanvas.h"

#include <cstdlib>

CbmLitCheckEnergyLossMuons::CbmLitCheckEnergyLossMuons():
	fMat("iron")
{
	fMom[0] = 47.04;
	fMom[1] = 56.16;
	fMom[2] = 68.02;
	fMom[3] = 85.09;
	fMom[4] = 100.3;
	fMom[5] = 152.7;
	fMom[6] = 176.4;
	fMom[7] = 221.8;
	fMom[8] = 286.8;
	fMom[9] = 364.2;
	fMom[10] = 391.7;
	fMom[11] = 494.5;
	fMom[12] = 899.5;
	fMom[13] = 1101.;
	fMom[14] = 1502.;
	fMom[15] = 2103.;
	fMom[16] = 3104.;
	fMom[17] = 4104.;
	fMom[18] = 8105.;
	fMom[19] = 10110.;
	fMom[20] = 14110.;
	fMom[21] = 20110.;
	fMom[22] = 30110.;
	fMom[23] = 40110.;
	fMom[24] = 80110.;
	fMom[25] = 100100.;
}

CbmLitCheckEnergyLossMuons::~CbmLitCheckEnergyLossMuons()
{
}

void CbmLitCheckEnergyLossMuons::Check()
{
	SetStyles();

	CreateGraphs();

	if (fMat == "iron") FillTableIron(); else
	if (fMat == "tungsten") FillTableTungsten(); else
	if (fMat == "carbon") FillTableCarbon();

	CalcEloss();

	DrawGraphs();
}

void CbmLitCheckEnergyLossMuons::DrawGraphs()
{
	TCanvas *c1 = new TCanvas("mean_energy_loss_muons","mean_energy_loss_muons",700, 500);

	fTable[0]->GetXaxis()->SetLimits(45, 100102);
	fTable[0]->SetMinimum(0.1);
	fTable[0]->SetMaximum(10.);

	DrawGraph(fTable[0], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR1, LIT_LINE_WIDTH, LIT_LINE_STYLE1, LIT_MARKER_SIZE,
			kOpenCircle, true, true, "ALP");

	DrawGraph(fTable[1], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR1, LIT_LINE_WIDTH, LIT_LINE_STYLE1, LIT_MARKER_SIZE,
			kOpenSquare, true, true, "LP");

	DrawGraph(fTable[2], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR1, LIT_LINE_WIDTH, LIT_LINE_STYLE1, LIT_MARKER_SIZE,
			kOpenTriangleUp, true, true, "LP");

	DrawGraph(fTable[3], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR1, LIT_LINE_WIDTH, LIT_LINE_STYLE1, LIT_MARKER_SIZE,
			kOpenDiamond, true, true, "LP");


	DrawGraph(fCalc[0], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR2, LIT_LINE_WIDTH, LIT_LINE_STYLE2, LIT_MARKER_SIZE,
			kOpenCircle, true, true, "LP");

	DrawGraph(fCalc[1], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR2, LIT_LINE_WIDTH, LIT_LINE_STYLE2, LIT_MARKER_SIZE,
			kOpenSquare, true, true, "LP");

	DrawGraph(fCalc[2], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR2, LIT_LINE_WIDTH, LIT_LINE_STYLE2, LIT_MARKER_SIZE,
			kOpenTriangleUp, true, true, "LP");

	DrawGraph(fCalc[3], "momentum [MeV/c]", "energy loss [Mev*cm^2/g]",
			LIT_COLOR2, LIT_LINE_WIDTH, LIT_LINE_STYLE2, LIT_MARKER_SIZE,
			kOpenDiamond, true, true, "LP");

//	TLegend* l1 = new TLegend(0.20, 0.97, 0.9, 0.7);
//	l1->SetFillColor(kWhite);
//	l1->SetHeader("Energy losses for muons in iron vs. momentum");
//	l1->AddEntry(fTable[0],"total (table)","lp");
//	l1->AddEntry(fCalc[0],"total (calculation)","lp");
//	l1->AddEntry(fTable[1],"ionization (table)","lp");
//	l1->AddEntry(fCalc[1],"ionization (calculation)","lp");
//	l1->AddEntry(fTable[2],"bremsstrahlung (table)","lp");
//	l1->AddEntry(fCalc[2],"bremsstrahlung (calculation)","lp");
//	l1->AddEntry(fTable[3],"pair production (table)","lp");
//	l1->AddEntry(fCalc[3],"pair production (calculation)","lp");
//	l1->Draw();

	SaveCanvasAsImage(c1, fOutputDir);
}

void CbmLitCheckEnergyLossMuons::CalcEloss()
{
	CbmLitMaterialInfo material;
	if (fMat == "tungsten") {
		material.SetA(183.84);
		material.SetZ(74.0);
		material.SetRho(19.3);
		material.SetRL(0.35);
	} else
	if (fMat == "iron") {
		material.SetA(55.85);
		material.SetZ(26.0);
		material.SetRho(7.87);
		material.SetRL(1.757622);
	} else
	if (fMat == "carbon") {
		material.SetA(12.0107);
		material.SetZ(6);
		material.SetRho(2.265);
		material.SetRL(18.8);
	} else exit(0);

	CbmLitMaterialEffectsImp me;
	CbmLitTrackParam par;
	for (Int_t i = 0; i < 26; i++) {
		par.SetQp(1./(fMom[i]*1e-3));
		fCalc[1]->SetPoint(i, fMom[i], me.BetheBloch(&par, &material) * 1e3);
		fCalc[2]->SetPoint(i, fMom[i], me.BetheHeitler(&par, &material) * 1e3);
		fCalc[3]->SetPoint(i, fMom[i], me.PairProduction(&par, &material) * 1e3);
		fCalc[0]->SetPoint(i, fMom[i], me.dEdx(&par, &material) * 1e3);
	}
}

void CbmLitCheckEnergyLossMuons::CreateGraphs()
{
	for (int i = 0; i < 4; i++) {
	    fTable[i] = new TGraph();
	    fCalc[i] = new TGraph();
	}
}

void CbmLitCheckEnergyLossMuons::FillTableIron()
{
	Double_t table_iron[26][4] = {
		//ion    //brems   //pair   //total
		{5.494,   0.000,    0.000,   5.494},
		{4.321,   0.000,    0.000,   4.321},
		{3.399,   0.000,    0.000,   3.399},
		{2.654,   0.000,    0.000,   2.654},
		{2.274,   0.000,    0.000,   2.274},
		{1.717,   0.000,    0.000,   1.717},

		{1.616,   0.000,    0.000,   1.616},
		{1.516,   0.000,    0.000,   1.516},
		{1.463,   0.000,    0.000,   1.463},
		{1.451,   0.000,    0.000,   1.451},
		{1.452,   0.000,    0.000,   1.453},
		{1.467,   0.000,    0.000,   1.467},
		{1.548,   0.000,    0.000,   1.548},

		{1.581,   0.001,    0.000,   1.582},
		{1.635,   0.001,    0.000,   1.637},
		{1.694,   0.002,    0.001,   1.697},
		{1.760,   0.003,    0.002,   1.767},
		{1.806,   0.004,    0.004,   1.816},
		{1.911,   0.010,    0.011,   1.936},

		{1.942,   0.014,    0.015,   1.975},
		{1.987,   0.021,    0.024,   2.039},
		{2.032,   0.033,    0.039,   2.113},
		{2.080,   0.054,    0.068,   2.214},
		{2.112,   0.076,    0.099,   2.303},
		{2.184,   0.171,    0.236,   2.623},
		{2.207,   0.221,    0.310,   2.777}
	};

	for (Int_t i = 0; i < 26; i++) {
		fTable[1]->SetPoint(i, fMom[i], table_iron[i][0]);
		fTable[2]->SetPoint(i, fMom[i], table_iron[i][1]);
		fTable[3]->SetPoint(i, fMom[i], table_iron[i][2]);
		fTable[0]->SetPoint(i, fMom[i], table_iron[i][3]);
	}
}

void CbmLitCheckEnergyLossMuons::FillTableTungsten()
{
	Double_t table_tungsten[26][4] = {
		//ion    //brems   //pair   //total
		{4.000,   0.000,    0.000,   4.00},
		{3.185,   0.000,    0.000,   3.185},
		{2.534,   0.000,    0.000,   2.534},
		{2.000,   0.000,    0.000,   2.000},
		{1.726,   0.000,    0.000,   1.726},
		{1.323,   0.000,    0.000,   1.323},

		{1.251,   0.000,    0.000,   1.251},
		{1.182,   0.000,    0.000,   1.182},
		{1.149,   0.000,    0.000,   1.149},
		{1.145,   0.000,    0.000,   1.145},
		{1.150,   0.000,    0.000,   1.150},
		{1.168,   0.000,    0.000,   1.168},
		{1.247,   0.001,    0.000,   1.249},

		{1.279,   0.001,    0.000,   1.281},
		{1.329,   0.002,    0.001,   1.333},
		{1.384,   0.004,    0.002,   1.391},
		{1.446,   0.006,    0.005,   1.459},
		{1.489,   0.009,    0.009,   1.509},
		{1.587,   0.023,    0.026,   1.640},

		{1.616,   0.031,    0.036,   1.687},
		{1.658,   0.047,    0.057,   1.768},
		{1.700,   0.073,    0.092,   1.874},
		{1.745,   0.120,    0.159,   2.035},
		{1.774,   0.169,    0.231,   2.190},
		{1.840,   0.383,    0.548,   2.800},
		{1.860,   0.496,    0.717,   3.110}
	};

	for (Int_t i = 0; i < 26; i++) {
		fTable[1]->SetPoint(i, fMom[i], table_tungsten[i][0]);
		fTable[2]->SetPoint(i, fMom[i], table_tungsten[i][1]);
		fTable[3]->SetPoint(i, fMom[i], table_tungsten[i][2]);
		fTable[0]->SetPoint(i, fMom[i], table_tungsten[i][3]);
	}
}

void CbmLitCheckEnergyLossMuons::FillTableCarbon()
{
	Double_t table_carbon[26][4] = {
		//ion    //brems   //pair   //total
		{7.119,   0.000,    0.000,   7.119},
		{5.551,   0.000,    0.000,   5.551},
		{4.332,   0.000,    0.000,   4.332},
		{3.356,   0.000,    0.000,   3.356},
		{2.861,   0.000,    0.000,   2.861},
		{2.127,   0.000,    0.000,   2.127},

		{1.992,   0.000,    0.000,   1.992},
		{1.854,   0.000,    0.000,   1.854},
		{1.775,   0.000,    0.000,   1.775},
		{1.745,   0.000,    0.000,   1.745},
		{1.745,   0.000,    0.000,   1.745},
		{1.751,   0.000,    0.000,   1.751},
		{1.819,   0.000,    0.000,   1.820},

		{1.850,   0.000,    0.000,   1.851},
		{1.900,   0.000,    0.000,   1.901},
		{1.955,   0.000,    0.000,   1.957},
		{2.018,   0.001,    0.001,   2.021},
		{2.062,   0.001,    0.001,   2.066},
		{2.161,   0.003,    0.003,   2.171},

		{2.191,   0.004,    0.004,   2.204},
		{2.234,   0.006,    0.007,   2.254},
		{2.278,   0.010,    0.011,   2.308},
		{2.326,   0.016,    0.019,   2.374},
		{2.359,   0.022,    0.028,   2.427},
		{2.434,   0.050,    0.068,   2.587},
		{2.458,   0.065,    0.089,   2.655}
	};

	for (Int_t i = 0; i < 26; i++) {
		fTable[1]->SetPoint(i, fMom[i], table_carbon[i][0]);
		fTable[2]->SetPoint(i, fMom[i], table_carbon[i][1]);
		fTable[3]->SetPoint(i, fMom[i], table_carbon[i][2]);
		fTable[0]->SetPoint(i, fMom[i], table_carbon[i][3]);
	}

}

ClassImp(CbmLitCheckEnergyLossMuons);