/* 
 * File:   CheckAlign.cc
 * Author: daniel
 * 
 * Created on 17. Februar 2014, 19:42
 */

#include "CheckAlign.h"

void CheckAlign::init(TBCore* core)
{ 
	for(auto dut: core->usedDUT)
	{
		int iden = dut->iden;
      int nCols = dut->getNcols();
      int nRows = dut->getNrows();
		
		CheckAlign::matchX[iden] = std::fabs(dut->getMatchX());
		CheckAlign::matchY[iden] = std::fabs(dut->getMatchY());
      
      CheckAlign::h_hitxVresx[iden] = new TProfile("", ";Hit Position [col];X Residual [#mum]", nCols, -.5, nCols-.5);
      CheckAlign::h_hityVresx[iden] = new TProfile("", ";Hit Position [row];X Residual [#mum]", nRows, -.5, nRows-.5);
      CheckAlign::h_hitxVresy[iden] = new TProfile("", ";Hit Position [col];Y Residual [#mum]", nCols, -.5, nCols-.5);
      CheckAlign::h_hityVresy[iden] = new TProfile("", ";Hit Position [row];Y Residual [#mum]", nRows, -.5, nRows-.5);
		
		std::string param = core->getParam(-1, CheckAlign::name, "scaleMatchX");
		if(param.compare("") != 0)
		{
			CheckAlign::scaleMatchX[iden] = std::fabs(std::stod(param));
		}
		else // default
		{
			CheckAlign::scaleMatchX[iden] = 1.0;
		}
		
		param = core->getParam(-1, CheckAlign::name, "scaleMatchY");
		if(param.compare("") != 0)
		{
			CheckAlign::scaleMatchY[iden] = std::fabs(std::stod(param));
		}
		else // default
		{
			CheckAlign::scaleMatchY[iden] = 1.0;
		}
		
		CheckAlign::matchXscaled[iden] = CheckAlign::matchX[iden] * CheckAlign::scaleMatchX[iden];
		CheckAlign::matchYscaled[iden] = CheckAlign::matchY[iden] * CheckAlign::scaleMatchY[iden];
   }
	
	CheckAlign::doCuts = true;
}

void CheckAlign::initRun(TBCore* core)
{  
	for(auto dut: core->usedDUT)
	{
		int iden = dut->iden;
		CheckAlign::h_resX[iden] = new TH1D("", ";X Residual [#mum]", 2000, -1000, 1000);
		CheckAlign::h_resY[iden] = new TH1D("", ";Y Residual [#mum]", 2000, -1000, 1000);
	}
}

// check it is correct?
void CheckAlign::buildEvent(TBCore* core, TBEvent* event)
{
	DUT* dut = event->dut;
	int iden = event->iden;
	
	// do cuts
	if(CheckAlign::doCuts == true)
	{
		if(event->fTracks != kGood)
		{
			return;
		}
	}
	
	for(auto tbtrack: event->tracks)
	{	
		// match cluster with higher mathing value
		TBCluster* matchedCluster = TBCluster::getMatchedCluster(&event->clusters, tbtrack, event, CheckAlign::matchXscaled[iden], CheckAlign::matchYscaled[iden]);
		
		if(matchedCluster == NULL)
		{
			continue;
		}
		
		int sumToT = TBCluster::getSumToT(matchedCluster);
		double hitX = TBCluster::getChargeWeightedX(matchedCluster, event);
		double hitY = TBCluster::getChargeWeightedY(matchedCluster, event);
		
		if(sumToT == 0)
		{
			continue;
		}
		
		if(matchedCluster->hits.size() == 1 and sumToT >= 3)
		{
			if(false) // does not use Etacorrection
			{
				/* Eta corrected */
			}
			else
			{
				CheckAlign::h_resX[iden]->Fill(tbtrack->trackX - hitX); 
				CheckAlign::h_resY[iden]->Fill(tbtrack->trackY - hitY);
			}
		}
		
		double resX = tbtrack->trackX - hitX;
		double resY = tbtrack->trackY - hitY;
		
		if(CheckAlign::doCuts == false or (CheckAlign::doCuts == true and tbtrack->ndof >= 2))
		{
			int hitCol;
			int hitRow;
			
			dut->getColRow(hitX, hitY, &hitCol, &hitRow);
			
			h_hitxVresx[iden]->Fill(hitCol, resX);
			h_hityVresx[iden]->Fill(hitRow, resX);
			h_hitxVresy[iden]->Fill(hitCol, resY);
			h_hityVresy[iden]->Fill(hitRow, resY);
		}
	}
}

void CheckAlign::finalizeRun(TBCore* core)
{
	core->output->processName = CheckAlign::name;
	
	bool badRun = false;
	char* histoTitle = new char[500];
	for(auto dut: core->usedDUT)
	{
		int iden = dut->iden;
		
		double resXmean = CheckAlign::h_resX[iden]->GetMean();
		double resYmean = CheckAlign::h_resY[iden]->GetMean();
		
		if(std::fabs(resXmean) > 10 or std::fabs(resYmean) > 10)
		{
			badRun = true;
			
			// set up cuts
			core->output->currentPreprocessCuts = "only matched cluster with size 1 and sum ToT >= 3";
			// -----------
			
			std::sprintf(histoTitle, "Residuals X DUT %i Run %i", iden, core->currentRun);
			CheckAlign::h_resX[iden]->SetTitle(histoTitle);
			std::sprintf(histoTitle, "res_X_dut_%i_run_%i", iden, core->currentRun);
			CheckAlign::h_resX[iden]->SetName(histoTitle);
			core->output->drawAndSave(CheckAlign::h_resX[iden], "", "emr");
			
			std::sprintf(histoTitle, "Residuals Y DUT %i Run %i", iden, core->currentRun);
			CheckAlign::h_resY[iden]->SetTitle(histoTitle);
			std::sprintf(histoTitle, "res_Y_dut_%i_run_%i", iden, core->currentRun);
			CheckAlign::h_resY[iden]->SetName(histoTitle);
			core->output->drawAndSave(CheckAlign::h_resY[iden], "", "emr");
		}
		
		CheckAlign::deltaX[iden].push_back(resXmean);
		CheckAlign::deltaY[iden].push_back(resYmean);
	 
		// set values in TBCore
		core->translations[iden][core->currentRun] = new TBTranslation(resXmean, resYmean);

		delete CheckAlign::h_resX[iden];
		delete CheckAlign::h_resY[iden];
	 }
	
	 CheckAlign::run.push_back(core->currentRun);
	 if(badRun == true)
	 {
		CheckAlign::badRuns.push_back(core->currentRun);
	 }
	 
	 CheckAlign::h_resX.clear();
	 CheckAlign::h_resY.clear();
	 delete[] histoTitle;
}

void CheckAlign::finalize(TBCore* core)
{
	core->output->processName = CheckAlign::name;
	
	char* histoTitle = new char[500];
	for(auto dut: core->usedDUT)
	{
		int iden = dut->iden;
		
		// set up cuts
			core->output->currentPreprocessCuts = "only matched cluster, sum ToT > 0 and ndof >= 2";
		// -----------
		
		std::sprintf(histoTitle, "Hit X Vs. Residuals X DUT %i", iden);
		CheckAlign::h_hitxVresx[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitXresX_dut_%i", iden);
		CheckAlign::h_hitxVresx[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hitxVresx[iden], "", "emr");
		
		std::sprintf(histoTitle, "Hit Y Vs. Residuals X DUT %i", iden);
		CheckAlign::h_hityVresx[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitYresX_dut_%i", iden);
		CheckAlign::h_hityVresx[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hityVresx[iden], "" , "emr");
		
		std::sprintf(histoTitle, "Hit X Vs. Residuals Y DUT %i", iden);
		CheckAlign::h_hitxVresy[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitXresY_dut_%i", iden);
		CheckAlign::h_hitxVresy[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hitxVresy[iden], "", "emr");

		std::sprintf(histoTitle, "Hit Y Vs. Residuals Y DUT %i", iden);
		CheckAlign::h_hityVresy[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitYresY_dut_%i", iden);
		CheckAlign::h_hityVresy[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hityVresy[iden], "", "emr");
   
		//Fit a pol 1 to get a measure for the misalignment and the tilt of the sample
		//Prepare the functions and fit
		TF1* pol1_xx = new TF1("pol1_xx", "pol1", 3, 15);
		CheckAlign::h_hitxVresx[iden]->Fit("pol1_xx", "+ QR");
		TF1* pol1_yx = new TF1("pol1_yx", "pol1", 15, 135);
		CheckAlign::h_hityVresx[iden]->Fit("pol1_yx", "+ QR");
		TF1* pol1_xy = new TF1("pol1_xy", "pol1", 3, 15);
		CheckAlign::h_hitxVresy[iden]->Fit("pol1_xy", "+ QR");
		TF1* pol1_yy = new TF1("pol1_yy", "pol1", 15, 135);
		CheckAlign::h_hityVresy[iden]->Fit("pol1_yy", "+ QR");
		
    //Pipe out the plots
		std::sprintf(histoTitle, "Hit X Vs. Residuals X Fitted DUT %i", iden);
		CheckAlign::h_hitxVresx[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitXresX_fitted_dut_%i", iden);
		CheckAlign::h_hitxVresx[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hitxVresx[iden], "", "emr");
		
		std::sprintf(histoTitle, "Hit Y Vs. Residuals X Fitted DUT %i", iden);
		CheckAlign::h_hityVresx[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitYresX_fitted_dut_%i", iden);
		CheckAlign::h_hityVresx[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hityVresx[iden], "", "emr");
		
		std::sprintf(histoTitle, "Hit X Vs. Residuals Y Fitted DUT %i", iden);
		CheckAlign::h_hitxVresy[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitXresY_fitted_dut_%i", iden);
		CheckAlign::h_hitxVresy[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hitxVresy[iden], "", "emr");
  
		std::sprintf(histoTitle, "Hit Y Vs. Residuals Y Fitted DUT %i", iden);
		CheckAlign::h_hityVresy[iden]->SetTitle(histoTitle);
		std::sprintf(histoTitle, "hitYresY_fitted_dut_%i", iden);
		CheckAlign::h_hityVresy[iden]->SetName(histoTitle);
		core->output->drawAndSave(CheckAlign::h_hityVresy[iden], "", "emr");

		TBLOG(kINFO, "DUT " << iden);
		TBLOG(kINFO, "hitXresX Offset: "
						<< pol1_xx->GetParameter(0) << " +- " << pol1_xx->GetParError(0) << " µm, tilt: "
						<< pol1_xx->GetParameter(1) << " +- " << pol1_xx->GetParError(1));
		TBLOG(kINFO, "hitYresX Offset: "
						<< pol1_yx->GetParameter(0) << " +- " << pol1_yx->GetParError(0) << " µm, tilt: "
						<< pol1_yx->GetParameter(1) << " +- " << pol1_yx->GetParError(1));
		TBLOG(kINFO, "hitXresY Offset: "
						<< pol1_xy->GetParameter(0) << " +- "<< pol1_xy->GetParError(0) << " µm, tilt: "
						<< pol1_xy->GetParameter(1) << " +- " << pol1_xy->GetParError(1));
		TBLOG(kINFO, "hitYresY Offset: "
						<< pol1_yy->GetParameter(0) << " +- " << pol1_yy->GetParError(0) << " µm, tilt: "
						<< pol1_yy->GetParameter(1) << " +- " << pol1_yy->GetParError(1));
  
		if(CheckAlign::badRuns.size() > 0)
		{
			TBLOG(kINFO, "Found " << CheckAlign::badRuns.size() << " suspicious runs:");
			for(int i = 0; i < badRuns.size(); i++)
			{
				TBLOG(kINFO, "Bad Run: " << CheckAlign::badRuns[i]);
				TBLOG(kINFO, "Delta X [um] = " << CheckAlign::deltaX[iden][i]);
				TBLOG(kINFO, "Delta Y [um] = " << CheckAlign::deltaY[iden][i]);
			}
			
		}
		else
		{
			TBLOG(kINFO, "Found 0 suspicious runs!");
		}
		
		
		// set up cuts
			core->output->currentPreprocessCuts = "only matched cluster with size 1 and sum ToT >= 3";
		// -----------
		
		// save deltaX and deltaY in histograms
		TH1D* h_deltaX = new TH1D("", ";Run Number;delta X",
										   CheckAlign::run.size(),
										   CheckAlign::run.front(),
										   CheckAlign::run.back());
		
		TH1D* h_deltaY = new TH1D("", ";Run Number;delta Y",
										   CheckAlign::run.size(),
										   CheckAlign::run.front(),
										   CheckAlign::run.back());
		
		for(int i = 0; i < CheckAlign::run.size(); i++)
		{
			h_deltaX->Fill(CheckAlign::run[i], CheckAlign::deltaX[iden][i]);
			h_deltaY->Fill(CheckAlign::run[i], CheckAlign::deltaY[iden][i]);
		}
		std::sprintf(histoTitle, "Translation Delta X DUT %i", iden);
		h_deltaX->SetTitle(histoTitle);
		std::sprintf(histoTitle, "translation_deltaX_dut_%i", iden);
		h_deltaX->SetName(histoTitle);
		core->output->drawAndSave(h_deltaX);
		
		std::sprintf(histoTitle, "Translation Delta Y DUT %i", iden);
		h_deltaY->SetTitle(histoTitle);
		std::sprintf(histoTitle, "translation_deltaY_dut_%i", iden);
		h_deltaY->SetName(histoTitle);
		core->output->drawAndSave(h_deltaY);
		
		
		delete pol1_xx;
		delete pol1_yx;
		delete pol1_xy;
		delete pol1_yy;
		delete h_deltaX;
		delete h_deltaY;
		
		delete CheckAlign::h_hitxVresx[iden];
		delete CheckAlign::h_hityVresx[iden];
		delete CheckAlign::h_hitxVresy[iden];
		delete CheckAlign::h_hityVresy[iden];
		delete CheckAlign::h_resX[iden];
		delete CheckAlign::h_resY[iden];
		
		CheckAlign::deltaX[iden].clear();
		CheckAlign::deltaY[iden].clear();
	}
	
	CheckAlign::h_hitxVresx.clear();
	CheckAlign::h_hityVresx.clear();
	CheckAlign::h_hitxVresy.clear();
	CheckAlign::h_hityVresy.clear();
	CheckAlign::h_resX.clear();
	CheckAlign::h_resY.clear();
	CheckAlign::deltaX.clear();
	CheckAlign::deltaY.clear();
	CheckAlign::badRuns.clear();
	CheckAlign::run.clear();
	delete[] histoTitle;
}