/* 
 * File:   qEfficiency.cc
 * Author: daniel
 * 
 * Created on 8. März 2014, 15:28
 */

#include "qEfficiency.h"

void qEfficiency::init(const TBCore* core)
{
	for(auto dut: core->usedDUT)
	{
		int iden = dut->iden;
		
		std::string param = core->getParam(iden, qEfficiency::name, "electrodeWidth");
		if(param.compare("") != 0)
		{
			qEfficiency::electrodeWidth[iden] = std::stod(param);
		}
		else
		{
			qEfficiency::electrodeWidth[iden] = 5.0;
		}
		
		param = core->getParam(iden, qEfficiency::name, "xy_totCutoff");
		if(param.compare("") != 0)
		{
			qEfficiency::xy_totCutoff[iden] = std::stod(param);
		}
		else
		{
			qEfficiency::xy_totCutoff[iden] = 25.0;
		}
		
		//Make histograms
		qEfficiency::xRes = 4;
		qEfficiency::yRes = 4;
		qEfficiency::totRes = 1;
		qEfficiency::qRes = 750;
		for(auto geometry: dut->pixelGeometries)
		{
			double maxPixelPitchX = geometry->getMaxPitchX();
			double maxPixelPitchY = geometry->getMaxPitchY();
			
			TH3D* totScatter = new TH3D("h_totScatter", ";Xmod position [#mum];Ymod position [#mum];TOT",
						(int) (maxPixelPitchX / qEfficiency::xRes), 0, maxPixelPitchX,
						(int) (maxPixelPitchY / qEfficiency::yRes), 0, maxPixelPitchY,
						(int) (150 / qEfficiency::totRes), 0, 150);
			
			TH3D* chargeScatter = new TH3D("h_chargeScatter", ";Xmod position [#mum];Ymod position [#mum];Charge [electrons]",
						(int) (maxPixelPitchX / qEfficiency::xRes), 0, maxPixelPitchX,
						(int) (maxPixelPitchY / qEfficiency::yRes), 0, maxPixelPitchY,
						(int) (60e3 / qEfficiency::qRes), 0, 60e3);
			
			TH1D* chargeLandau = new TH1D("h_chargeLandau", ";Charge [ke];N", (int) (60e3 / qEfficiency::qRes), 0, 60);
			
			qEfficiency::h_totScatter[iden].push_back(totScatter); 
			qEfficiency::h_chargeScatter[iden].push_back(chargeScatter);
			qEfficiency::h_chargeLandau[iden].push_back(chargeLandau);
		}
		
	}
}

void qEfficiency::event(const TBCore* core, const TBEvent* event)
{
	DUT* dut = event->dut;
	int iden = event->iden;
	
	if(event->fTracks != kGood)
	{
		return;
	}
	
	if(event->fClusters != kGood)
	{
		return;
	}
	
	if(event->fHits != kGood)
	{
		return;
	}
	
	for(auto track: event->tracks)
	{
		if(track->fTrackMaskedRegion == kGood)
		{
			continue;
		}
		
		if(track->fMatchedCluster != kGood)
		{
			continue;
		}
		
		PixelGeometry* pixelGeometry = dut->pixelArray[track->trackCol][track->trackRow].getGeometry();
		int geometryNumber = dut->getGeometryNumber(pixelGeometry);
		double foldedX;
		double foldedY;
		dut->getFoldedXY(track->trackCol, track->trackRow, track->trackX, track->trackY, &foldedX, &foldedY);
		
		int sumToT = TBCluster::getSumToT(track->matchedCluster);
		int sumCharge = TBCluster::getSumCharge(track->matchedCluster, event);
		
		qEfficiency::h_totScatter[iden][geometryNumber]->Fill(foldedX, foldedY, sumToT);
		
		if(dut->totcalib->hasPerPixelCalib())
		{
			qEfficiency::h_chargeScatter[iden][geometryNumber]->Fill(foldedX, foldedY, sumCharge);
			qEfficiency::h_chargeLandau[iden][geometryNumber]->Fill(sumCharge / 1000.);
		}
		else if(dut->totcalib->hasCalib())
		{
			qEfficiency::h_chargeScatter[iden][geometryNumber]->Fill(foldedX, foldedY, sumCharge);
			qEfficiency::h_chargeLandau[iden][geometryNumber]->Fill(sumCharge / 1000.);
		}
	}
}

void qEfficiency::finalize(const TBCore* core)
{
	core->output->processName = qEfficiency::name;
	
	char* histoTitle = new char[500];
	for(auto dut: core->usedDUT)
	{
		int iden = dut->iden;
		
		for(auto geometry: dut->pixelGeometries)
		{
			int geometryNumber = dut->getGeometryNumber(geometry);
			double maxPixelPitchX = geometry->getMaxPitchX();
			double maxPixelPitchY = geometry->getMaxPitchY();
			
			std::string edge = "";
			if(dut->pixelGeometries[geometryNumber]->isEdgePixel() == true)
			{
				edge = "(Edge Pixel)";
			}
			
			// set up cuts
			if(dut->totcalib->hasPerPixelCalib())
			{
				core->output->cuts = "only matched tracks with per pixel calib";
			}
			else if(dut->totcalib->hasCalib())
			{
				core->output->cuts = "only matched tracks with calib";
			}
			else
			{
				core->output->cuts = "only matched tracks";
			}
			core->output->currentPreprocessCuts = "";
			// -----------
			
			std::sprintf(histoTitle, "Charge Scatter DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			qEfficiency::h_chargeScatter[iden][geometryNumber]->SetTitle(histoTitle);
			std::sprintf(histoTitle, "chargeScatter_dut_%i_geometry_%i", iden, geometryNumber);
			qEfficiency::h_chargeScatter[iden][geometryNumber]->SetName(histoTitle);
			core->output->drawAndSave(qEfficiency::h_chargeScatter[iden][geometryNumber], "", "e");
			
			std::sprintf(histoTitle, "Charge Landau DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			qEfficiency::h_chargeLandau[iden][geometryNumber]->SetTitle(histoTitle);
			std::sprintf(histoTitle, "chargeLandau_dut_%i_geometry_%i", iden, geometryNumber);
			qEfficiency::h_chargeLandau[iden][geometryNumber]->SetName(histoTitle);
			core->output->drawAndSave(qEfficiency::h_chargeLandau[iden][geometryNumber], "", "e");
			
			
			
			// set up cuts
			core->output->cuts = "only matched tracks";
			core->output->currentPreprocessCuts = "";
			// -----------
			
			
			//Save raw data to file
			std::sprintf(histoTitle, "ToT Scatter DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			qEfficiency::h_totScatter[iden][geometryNumber]->SetTitle(histoTitle);
			std::sprintf(histoTitle, "totScatter_dut_%i_geometry_%i", iden, geometryNumber);
			qEfficiency::h_totScatter[iden][geometryNumber]->SetName(histoTitle);
			core->output->drawAndSave(qEfficiency::h_totScatter[iden][geometryNumber], "", "e");
			
			TF1* ffit = new TF1("Langaus", langaufunction, (qEfficiency::h_chargeLandau[iden][geometryNumber]->GetMean()) / 2., 40, 4);
			double startvalues[4] = { 0.65, qEfficiency::h_chargeLandau[iden][geometryNumber]->GetMean(), qEfficiency::h_chargeLandau[iden][geometryNumber]->Integral((qEfficiency::h_chargeLandau[iden][geometryNumber]->GetMean()) / 2., 40), 1.6 };
			ffit->SetParameters(startvalues);
			ffit->SetParNames("Width", "MPV", "Area", "GSigma");
			TBLOG(kINFO, "DUT " << iden << " Geometry: " << geometryNumber);
			qEfficiency::h_chargeLandau[iden][geometryNumber]->Fit("Langaus", "LLR");
			
			std::sprintf(histoTitle, "Charge Landau Fitted DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			qEfficiency::h_chargeLandau[iden][geometryNumber]->SetTitle(histoTitle);
			std::sprintf(histoTitle, "chargeLandau_fitted_dut_%i_geometry_%i", iden, geometryNumber);
			qEfficiency::h_chargeLandau[iden][geometryNumber]->SetName(histoTitle);
			core->output->drawAndSave(qEfficiency::h_chargeLandau[iden][geometryNumber], "", "e");
			
			TBLOG(kINFO, "DUT " << iden << " Geometry: " << geometryNumber);
			TBLOG(kINFO, "Langaus MPV:    " << ffit->GetParameter(1) << " +- " << ffit->GetParError(1) << " ke");
			TBLOG(kINFO, "Langaus Width:  " << ffit->GetParameter(0) << " +- " << ffit->GetParError(0));
			TBLOG(kINFO, "Langaus Area:   " << ffit->GetParameter(2) << " +- " << ffit->GetParError(2));
			TBLOG(kINFO, "Langaus GSigna: " << ffit->GetParameter(3) << " +- " << ffit->GetParError(3));
			
			delete ffit;
			
			char* histoT = new char[500];
			std::sprintf(histoT, "ToT Scatter DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			std::sprintf(histoTitle, "totScater_dut_%i_geometry_%i", iden, geometryNumber);
			qEfficiency::zxProject(core, histoT, histoTitle, qEfficiency::h_totScatter[iden][geometryNumber]);
			
			std::sprintf(histoT, "Charge Scatter DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			std::sprintf(histoTitle, "chargeScater_dut_%i_geometry_%i", iden, geometryNumber);
			qEfficiency::zxProject(core, histoT, histoTitle, qEfficiency::h_chargeScatter[iden][geometryNumber]);
		

			//Project out y-axis to make x-z scatterplot
			TH2D* h_totProject = (TH2D*) qEfficiency::h_totScatter[dut->iden][geometryNumber]->Project3D("zx");
			TH2D* h_chargeProject = (TH2D*) qEfficiency::h_chargeScatter[dut->iden][geometryNumber]->Project3D("zx");
			h_totProject->SetTitle("");
			h_chargeProject->SetTitle("");
			
			std::sprintf(histoTitle, "ToT Project DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			h_totProject->SetTitle(histoTitle);
			std::sprintf(histoTitle, "totProject_dut_%i_geometry_%i", iden, geometryNumber);
			h_totProject->SetName(histoTitle);
			core->output->drawAndSave(h_totProject, "", "e");
			
			std::sprintf(histoTitle, "Charge Project DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			h_chargeProject->SetTitle(histoTitle);
			std::sprintf(histoTitle, "totProject_dut_%i_geometry_%i", iden, geometryNumber);
			h_chargeProject->SetName(histoTitle);
			core->output->drawAndSave(h_chargeProject, "", "e");
			
			//Make normal (mean/spread) 2D and 1D profiles
			TProfile2D* p_totProfile = h_totScatter[dut->iden][geometryNumber]->Project3DProfile("yx");
			
			std::sprintf(histoTitle, "ToT Profile DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			p_totProfile->SetTitle(histoTitle);
			std::sprintf(histoTitle, "totProfile_dut_%i_geometry_%i", iden, geometryNumber);
			p_totProfile->SetName(histoTitle);
			core->output->drawAndSave(p_totProfile, "colz", "e");
		
			// BJD: this is a dirty hack for IBL sensor review request, 4 July 2011
			// Matthias George made me do it!!!
			/*TH2D* h2_totProfile = (TH2D*) p_totProfile->ProjectionXY();
			h2_totProfile->SetTitle(";Long Pixel Pitch [#mum];Short Pixel Pitch [#mum];TOT");
			double k = maxPixelPitchY / maxPixelPitchX;
			TCanvas* c_aspect = new TCanvas("", "", 1200, (int) floor(1200 * k));
			c_aspect->cd();
			gStyle->SetPalette(1, 0);
			h2_totProfile->SetStats(false);
			h2_totProfile->GetYaxis()->SetTickLength(0.01);
			// Tinker with presentation, since ROOT tries so hard to make these plots unreadable
			gPad->SetRightMargin(0.7 * k);
			gPad->SetLeftMargin(0.5 * k);
			gPad->SetBottomMargin(0.22);
			h2_totProfile->GetXaxis()->SetNdivisions((int) maxPixelPitchX / 50.0);
			h2_totProfile->GetYaxis()->SetNdivisions(5);
			h2_totProfile->GetXaxis()->SetLabelOffset(0.005);
			h2_totProfile->GetYaxis()->SetLabelOffset(0.005);
			h2_totProfile->GetZaxis()->SetLabelOffset(0.005);
			h2_totProfile->GetXaxis()->SetLabelSize(0.1);
			h2_totProfile->GetYaxis()->SetLabelSize(0.1);
			h2_totProfile->GetZaxis()->SetLabelSize(0.08);
			h2_totProfile->GetXaxis()->SetTitleOffset(1.0);
			h2_totProfile->GetYaxis()->SetTitleOffset(2.0 * k);
			h2_totProfile->GetZaxis()->SetTitleOffset(2.0 * k);
			h2_totProfile->GetXaxis()->SetTitleSize(0.1);
			h2_totProfile->GetYaxis()->SetTitleSize(0.07);
			h2_totProfile->GetZaxis()->SetTitleSize(0.07);
			h2_totProfile->Draw("colz");
			gPad->Update();
			
			
			char* fileName = core.tboutput->buildHistName(this->name, ("totProfileAspect_dut"+strIden+"_geometry_"+std::to_string(geometryNumber)).c_str());
			c_aspect->SaveAs(fileName);
			delete c_aspect;
			 */
			
			TProfile2D* p_chargeProfile = h_chargeScatter[dut->iden][geometryNumber]->Project3DProfile("yx");
			
			std::sprintf(histoTitle, "Charge Profile DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			p_chargeProfile->SetTitle(histoTitle);
			std::sprintf(histoTitle, "chargeProfile_dut_%i_geometry_%i", iden, geometryNumber);
			p_chargeProfile->SetName(histoTitle);
			core->output->drawAndSave(p_chargeProfile, "colz", "e");
			
			
			TProfile* p_totProfile_all = h_totProject->ProfileX();
			// TODO
			TProfile* p_chargeProfile_all = h_chargeProject->ProfileX();
			
			std::sprintf(histoTitle, "ToT Profile All DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			p_totProfile_all->SetTitle(histoTitle);
			std::sprintf(histoTitle, "totProfileAll_dut_%i_geometry_%i", iden, geometryNumber);	
			p_totProfile_all->SetName(histoTitle);
			core->output->drawAndSave(p_totProfile_all, "", "e");
			
			std::sprintf(histoTitle, "Charge Profile All DUT %i Geometry %i %s", iden, geometryNumber, edge.c_str());
			p_chargeProfile_all->SetTitle(histoTitle);
			std::sprintf(histoTitle, "chargeProfileAll_dut_%i_geometry_%i", iden, geometryNumber);
			p_chargeProfile_all->SetName(histoTitle);
			core->output->drawAndSave(p_chargeProfile_all, "", "e");
			
				
			delete h_totProject;
			delete h_chargeProject;
			delete p_totProfile;
			delete p_chargeProfile;
			delete p_totProfile_all;
			//delete p_chargeProfile_all;
			delete qEfficiency::h_totScatter[iden][geometryNumber];
			delete qEfficiency::h_chargeScatter[iden][geometryNumber];
			delete qEfficiency::h_chargeLandau[iden][geometryNumber];
			
			delete [] histoT;
		}
		
		qEfficiency::h_totScatter[iden].clear();
		qEfficiency::h_chargeScatter[iden].clear();
		qEfficiency::h_chargeLandau[iden].clear();
	}
	
	qEfficiency::h_totScatter.clear();
	qEfficiency::h_chargeScatter.clear();
	qEfficiency::h_chargeLandau.clear();
}


void qEfficiency::zxProject(const TBCore* core, const char* histoTitle, const char* histoName, TH3D* h1, TH3D* h2)
{
	char* buffer = new char[500];

	double range_ylow = h1->GetZaxis()->GetBinLowEdge(h1->GetZaxis()->GetFirst());
	double range_yhi = h1->GetZaxis()->GetBinUpEdge(h1->GetZaxis()->GetLast());

	TH2D* h1_project = (TH2D*) h1->Project3D("zx");
	TH2D* h2_project = NULL;
	h1_project->SetTitle(histoTitle);
	if(h2 != NULL)
	{
		TH2D* h2_project = (TH2D*) h2->Project3D("zx");
		h1_project->Add(h2_project);
	}
	std::sprintf(buffer, "%s %s", histoTitle, "Project");
	h1_project->SetTitle(buffer);
	std::sprintf(buffer, "%s_%s", histoName, "project");
	h1_project->GetYaxis()->SetRangeUser(range_ylow, range_yhi);
	h1_project->SetName(buffer);
	core->output->drawAndSave(h1_project, "", "e");

	TProfile* p_profile = h1_project->ProfileX();
	std::sprintf(buffer, "%s %s", histoTitle, "Profile");
	p_profile->SetTitle(buffer);
	std::sprintf(buffer, "%s_%s", histoName, "profile");
	p_profile->GetYaxis()->SetRangeUser(range_ylow, range_yhi);
	p_profile->SetName(buffer);
	core->output->drawAndSave(p_profile, "", "e");
	
	delete h1_project;
	if(h2 != NULL)
	{
		delete h2_project;
	}
	delete p_profile;
	delete [] buffer;
}

Double_t langaufunction(Double_t *x, Double_t *par)
{
	//Fit parameters:
	//par[0]=Width (scale) parameter of Landau density
	//par[1]=Most Probable (MP, location) parameter of Landau density
	//par[2]=Total area (integral -inf to inf, normalization constant)
	//par[3]=Width (sigma) of convoluted Gaussian function
	//
	//In the Landau distribution (represented by the CERNLIB approximation), 
	//the maximum is located at x=-0.22278298 with the location parameter=0.
	//This shift is corrected within this function, so that the actual
	//maximum is identical to the MP parameter.

	// Numeric constants
	Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
	Double_t mpshift = -0.22278298; // Landau maximum location

	// Control constants
	Double_t np = 100.0; // number of convolution steps
	Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas

	// Variables
	Double_t xx;
	Double_t mpc;
	Double_t fland;
	Double_t sum = 0.0;
	Double_t xlow, xupp;
	Double_t step;
	Double_t i;

	// MP shift correction
	mpc = par[1] - mpshift * par[0];

	// Range of convolution integral
	xlow = x[0] - sc * par[3];
	xupp = x[0] + sc * par[3];

	step = (xupp - xlow) / np;

	// Convolution integral of Landau and Gaussian by sum
	for(i = 1.0; i <= np / 2; i++)
	{
		xx = xlow + (i - .5) * step;
		fland = TMath::Landau(xx, mpc, par[0]) / par[0];
		sum += fland * TMath::Gaus(x[0], xx, par[3]);

		xx = xupp - (i - .5) * step;
		fland = TMath::Landau(xx, mpc, par[0]) / par[0];
		sum += fland * TMath::Gaus(x[0], xx, par[3]);
	}

	return (par[2] * step * sum * invsq2pi / par[3]);
}