// // Simone Bianco 14/07/2010

// This Class' Header ------------------
#include "TtTotAli.h"

// C/C++ Headers ----------------------
#include <iostream>

// Collaborating Class Headers --------
#include "FairRootManager.h"
#include "TClonesArray.h"

#include "PndTrackCand.h"
#include "PndSdsHit.h"
#include "PndTrackCandHit.h"

#include "TROOT.h"
#include "TFile.h"
#include "TGeoTrack.h"
#include "TGeoManager.h"
#include "TLorentzVector.h"
#include "TVector2.h"

// Fit Classes -----------
#include <TMath.h>
#include <TVector3.h>
#include <TRandom2.h>
//#include <TStyle.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TVirtualFitter.h>
#include <TPolyLine3D.h>
#include <Math/Vector3D.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TChain.h>

#include <fstream>

using namespace ROOT::Math;
using namespace std;

TtTotAli::TtTotAli()
{

  fnamePixIn = "MVDHitsPixel";
  fnameStrIn = "MVDHitsStrip";
  fTreeName = "cbmsim";

  fPrint = 0;
  fEvent = 0;
  fSmallRange = kFALSE;
  fDoubleGaus = kFALSE;
  fUseGeo = kFALSE;

  hx = new TH1F("hx","hx",1000,-5.,+5.);
  hy = new TH1F("hy","hy",1000,-5.,+5.);

  t = new TTree();

  fExcludeStrips = kFALSE;
  fExcludePixels = kFALSE;
  fOnlyRes = kFALSE;
  fIsAChain = kFALSE;
  fPixZOff = 0.;

}

TtTotAli::~TtTotAli()
{
 
}  

void TtTotAli::Initialize()
{
    std::cout << "hi!" << std::endl;	

    std::cout << "Input file:" << std::endl;	
    std::cout << fInFileName.Data() << std::endl;	

    std::cout << "Branches:" << std::endl;	
    std::cout << fnameStrIn.Data() << std::endl;	
    std::cout << fnamePixIn.Data() << std::endl;
    
    

    if (!fIsAChain)
    {
    	TFile *f = new TFile(fInFileName.Data());
    	t = (TTree*) f->Get(fTreeName.Data());
    }
    else
    {
	TChain *ch = new TChain(fTreeName.Data());
	std::cout << "Number of files to include in the chain: ";
	Int_t nF;
	std::cin >> nF;
	std::cout << std::endl;
	for (Int_t ww = 0 ; ww < nF ; ww++)
	{
		std::cout << ww <<"° file: ";
		TString nameF;
		std::cin >> nameF;
		ch->Add(nameF.Data());
	}
	t = ch;
    }

    


}

void TtTotAli::Run()
{
 
    if (!fExcludeStrips)
    {
    	strInput = new TClonesArray("PndSdsHit",10);
    	t->SetBranchAddress(fnameStrIn.Data(),&strInput);
    }
 
    if (!fExcludePixels)
    {
	pixInput = new TClonesArray("PndSdsHit",10);
    	t->SetBranchAddress(fnamePixIn.Data(),&pixInput);
    }

    fEvent = 0;
    fNbox = fNStrip+fNPixel;

    if (fExclBoxStr>0) 
    {
	fExclBox = fExclBoxStr;
    }
    else
    {
	fExclBox = fExclBoxPix + fNStrip;
    }

    if (sX.size() == 0) 
    {  // this means it's the first execution of Run(), so we "define" the sizes
    	for (Int_t gg = 0 ; gg < (fNStrip+fNPixel) ; gg++)
    	{
		sX.push_back(0.);
		sY.push_back(0.);
		sigX.push_back(0.);
		sigY.push_back(0.);
		m_X.push_back(0.);
		m_Y.push_back(0.);
    	}
    }


    for (Int_t ie = 0 ; ie < t->GetEntries() ; ie++)
    {
	
	t->GetEntry(ie);
	if (ie % 10000 == 0) std::cout << "Event: " << ie << std::endl;

	AnaEvent();
	
    }

    FinishRun();

}


Double_t TtTotAli::distance2(double x,double y,double z,double *par) { 
   // distance line point is D= | (xp-x0) cross  ux | 
   // where ux is direction of line and x0 is a point in the line (like t = 0) 
   TVector3 xp(x,y,z); 
   TVector3 x0(par[0],par[2],0); 
   TVector3 x1(par[0]+par[1],par[2]+par[3],1); 
   TVector3 u = (x1-x0).Unit(); 
   double d2 = ((xp-x0).Cross(u)) .Mag2(); 

   //std::cout << "Root of the distance: " << TMath::Sqrt(d2) << std::endl;

   return d2; 
}

void TtTotAli::SumDistance2(int &, double *, double & sum, double * par, int ) {

   TGraph2D * gr = dynamic_cast<TGraph2D*>( (TVirtualFitter::GetFitter())->GetObjectFit() );
   assert(gr != 0);
   double * x = gr->GetX();
   double * y = gr->GetY();
   double * z = gr->GetZ();
   double * Ex = gr->GetEX();
   double * Ey = gr->GetEY();
   double * Ez = gr->GetEZ();

   int npoints = gr->GetN();
   sum = 0;
   double d = 0.;

   //for (int i  = 0; i < 3 ; ++i) {
   for (int i  = 0; i < npoints ; ++i) {
       //std::cout << "Sensor: " << i;
      d = distance2(x[i],y[i],z[i],par);
      sum += d;
   }
}


void TtTotAli::AnaEvent()
{

  //  std::cout << "Event: " << fEvent << std::endl;


  if (fEvent == 0)
  {

  	delete hx;
	delete hy;
	hx = new TH1F("hx","hx",1000,-5.,+5.);
	hy = new TH1F("hy","hy",1000,-5.,+5.);

  }


  Int_t ntcand = 0;

  Int_t sizS = 0;
  Int_t sizP = 0;

  if (!fExcludeStrips) 
  {
	ntcand +=strInput->GetEntries();
	sizS = strInput->GetEntries();
  }
  if (!fExcludePixels) 
  {
	ntcand +=pixInput->GetEntries();
	sizP = pixInput->GetEntries();
  }

  if ((sizS >= fNStrip) && (sizP >= fNPixel))
    {
      
      	std::map<Double_t,Int_t> SensorsPosStr;
      	std::map<Double_t,Int_t> SensorsPosPix;

	std::map<Int_t,Int_t> SensorsStrCount;
      	std::map<Int_t,Int_t> SensorsPixCount;

      	Int_t name;


      	Int_t itrr = 0;

      // strips

      	if (!fExcludeStrips)
      	{
      		while (SensorsPosStr.size() < strInput->GetEntries() && itrr < strInput->GetEntries())
		{
			PndSdsHit *theHitP = (PndSdsHit*)   strInput->At(itrr);
			if(theHitP)
			{
				name = theHitP->GetSensorID();
				SensorsPosStr[theHitP->GetZ()] = name;
			}
			itrr++;
		}

      	}
      	// pixels

      	itrr = 0;

      	if (!fExcludePixels)
      	{
      	while (SensorsPosPix.size() < pixInput->GetEntries() && itrr < pixInput->GetEntries())
		{
			PndSdsHit *theHitP = (PndSdsHit*)   pixInput->At(itrr);
			if(theHitP)
			{
				name = theHitP->GetSensorID();
				SensorsPosPix[theHitP->GetZ()] = name;
			}
			itrr++;
		}
      	}
    
	Bool_t testStrip[sizS];
	Bool_t testPixel[sizP];

	if (!fExcludeStrips)
      	{
		for (Int_t vv = 0 ; vv < strInput->GetEntries() ; vv++)
		{
			testStrip[vv] = kFALSE;
		}
	}

	if (!fExcludePixels)
      	{
		for (Int_t vv = 0 ; vv < pixInput->GetEntries() ; vv++)
		{
			testPixel[vv] = kFALSE;
		}
	}

	if (!fExcludeStrips)
      	{
		for(Int_t itr1=0;itr1<strInput->GetEntries();++itr1){
			PndSdsHit *theHit = (PndSdsHit*)   strInput->At(itr1);
			Int_t buffDet = theHit->GetSensorID();
			if (!SensorsStrCount[buffDet]) SensorsStrCount[buffDet] = 1;
			else SensorsStrCount[buffDet]++;
			if (SensorsStrCount[buffDet] == 1) testStrip[buffDet] = kTRUE;
			else testStrip[buffDet] = kFALSE;
      		}
	}

	if (!fExcludePixels)
      	{
		for(Int_t itr2=0;itr2<pixInput->GetEntries();++itr2){
			PndSdsHit *theHit2 = (PndSdsHit*)   pixInput->At(itr2);
			Int_t buffDet = theHit2->GetSensorID();
			if (!SensorsPixCount[buffDet]) SensorsPixCount[buffDet] = 1;
			else SensorsPixCount[buffDet]++;
			if (SensorsPixCount[buffDet] == 1) testPixel[buffDet] = kTRUE;
			else testPixel[buffDet] = kFALSE;
			//std::cout << "Pix hit " << itr2 << ", sensID " << buffDet << ", z " << theHit2->GetZ() << std::endl;
      		}
	}

	const Int_t sizeMapStr = SensorsPosStr.size();
	Int_t DetNamesStr[sizeMapStr];
	Double_t PosStr[sizeMapStr];

	const Int_t sizeMapPix = SensorsPosPix.size();
	Int_t DetNamesPix[sizeMapPix];
	Double_t PosPix[sizeMapPix];

	Int_t jj = 0;

	if ((sizeMapStr >= fNStrip) && (sizeMapPix >= fNPixel))
      	{
		//std::cout << "Strip Planes" << std::endl;
  	  	for (std::map<Double_t,Int_t>::iterator it=SensorsPosStr.begin();it!=SensorsPosStr.end();++it)
	    	{
	      		if (fEvent < 10)std::cout << "position: " << it->first << " name: " << (it->second) << std::endl;
			DetNamesStr[jj] = it->second;
			PosStr[jj] = it -> first;
			jj++;
	    	}
		jj = 0;
		//std::cout << "Pixel Planes" << std::endl;
  	  	for (std::map<Double_t,Int_t>::iterator it=SensorsPosPix.begin();it!=SensorsPosPix.end();++it)
	    	{
	      		if (fEvent < 10)std::cout << "position: " << it->first << " name: " << (it->second) << std::endl;
			DetNamesPix[jj] = it->second;
			PosPix[jj] = it -> first;
			jj++;
	    	}
	}

	Bool_t sat=kTRUE;

	for (Int_t kk = 0 ; kk < fNStrip ; kk++)
	{
		sat = sat * testStrip[kk];
		//std::cout << "Strip test - " << kk << " : " << testStrip[kk] << " --> " << sat << std::endl;
	}
	for (Int_t kk = 0 ; kk < fNPixel ; kk++)
	{
		sat = sat * testPixel[kk];
		//std::cout << "Pixel test - " << kk << " : " << testPixel[kk] << " --> " << sat << std::endl;
	}
	
	//std::cout << "sizS " << sizS << ", fNStrip " << fNStrip << ", sizP " << sizP << ", fNPixel " << fNPixel << ", sat " << sat << std::endl;
     
	if ((sizS>=fNStrip) && (sizP>=fNPixel) && sat )
	{

		Int_t maxStrip,maxPixel;

		if (fExclBoxStr <= 0)	
		{	
			maxStrip = sizeMapStr;
		}
		if (fExclBoxStr >  0)
		{	
			maxStrip = sizeMapStr-1;
		}

		if (fExclBoxPix <= 0)
		{	
			maxPixel = sizeMapPix;
		}
		if (fExclBoxPix >  0)
		{	
			maxPixel = sizeMapPix-1;
		}

		
		Int_t DetRemStr[maxStrip];		
		Int_t DetRemPix[maxPixel];


		if (!fExcludeStrips)
      		{
			for (Int_t ss = 0 ; ss < maxStrip ; ss++) 
	    		{
				DetRemStr[ss] = 0;
			}
		}

		if (!fExcludePixels)
      		{
			for (Int_t ss = 0 ; ss < maxPixel ; ss++) 
		    	{
				DetRemPix[ss] = 0;
			}
		}

		Int_t counter = 0;

		Double_t BuffZ = -9999.;

		if (!fExcludeStrips)
      		{
			for (Int_t kk = 0 ; kk < (maxStrip) ; kk++)
			{	
				if ((fExclBoxStr-1) == kk) 
				{
		  			BuffZ = PosStr[kk];
		  			continue;
				}
				else
				{
		  			DetRemStr[counter] = kk;
		  			counter++;
				}
      			}
	    	}

		counter = 0;

		if (!fExcludePixels)
      		{
			for (Int_t kk = 0 ; kk < (maxPixel) ; kk++)
			{	
				if ((fExclBoxPix-1) == kk) 
				{
		  			BuffZ = PosPix[kk];
		  			continue;
				}
				else
				{
		  			DetRemPix[counter] = kk;
		  			counter++;
				}
      		    	}
		}


	  	Double_t x[fNbox];
	  	Double_t y[fNbox];
	  	Double_t z[fNbox];

	  	Double_t Erx[fNbox];
	  	Double_t Ery[fNbox];
	  	Double_t Erz[fNbox];

	  	Int_t track = 0;

      
	  	Double_t RealX = -1999999., RealY = -1999999;
      

	  	Double_t buffErrX = -0.099999;
	  	Double_t buffErrY = -0.099999;

		if (!fExcludeStrips)
      		{
			for (Int_t it1 = 0 ; it1 < strInput->GetEntries(); it1++)
			{
		
				PndSdsHit *theHit = (PndSdsHit*)   strInput->At(it1);

				for(Int_t ww = 0 ; ww < fNStrip ; ww++)
				{
					if ((theHit->GetSensorID()) == DetNamesStr[ww])
					{
		  				x[ww] = theHit->GetX();
		  				y[ww] = theHit->GetY();
		  				z[ww] = theHit->GetZ();
		  				Erx[ww] = theHit->GetDx();
		  				Ery[ww] = theHit->GetDy();
		  				Erz[ww] = theHit->GetDz();
		  				//std::cout << ww+1 << "° det: " << x[ww] << "," << y[ww] << "," << z[ww] << std::endl;
					}
				}
			} // end loop on strip hits
		}

		if (!fExcludePixels)
      		{
			for (Int_t it1 = 0 ; it1 < pixInput->GetEntries(); it1++)
			{

				PndSdsHit *theHit = (PndSdsHit*)   pixInput->At(it1);

				for(Int_t ww = 0 ; ww < fNPixel ; ww++)
				{
					if ((theHit->GetSensorID()) == DetNamesPix[ww])
					{
		  				x[fNStrip+ww] = -1.*theHit->GetX();
		  				y[fNStrip+ww] = theHit->GetY();
		  				z[fNStrip+ww] = theHit->GetZ()+fPixZOff;
		  				Erx[fNStrip+ww] = theHit->GetDx();
		  				Ery[fNStrip+ww] = theHit->GetDy();
		  				Erz[fNStrip+ww] = theHit->GetDz();
		  				//std::cout << fNStrip+ww+1 << "° det: " << x[fNStrip+ww] << "," << y[fNStrip+ww] << "," << z[fNStrip+ww] << std::endl;
					}
				}
			} // end loop on pixel hits

		}

		// now the array contain the strip hits in z-order and then the pixel hits in z-order
	  
		Int_t uu = 0;

		// to check alignment...

		if (!fOnlyRes)
		{
			for (Int_t ww = 0 ; ww < fNbox ; ww++)
		    	{
	      			if (TMath::Abs(Erx[ww]) < 0.5) 
				{
		  			x[ww]+=sX[ww];
				}
	      			if (TMath::Abs(Ery[ww]) < 0.5) 
				{
					y[ww]+=sY[ww];
				}
	    		}
		}

		Double_t DX,DY;

		// to check

		//for (Int_t ww = 0 ; ww < fNbox ; ww++)
	    	//{
		//	std::cout << "# (" << x[ww] << "," << y[ww] << "," << z[ww] << ")"<< std::endl;
	    	//}



		if (!fUseGeo)
	  	{
	 		MyFit(x,y,z,Erx,Ery,Erz,x[fExclBox-1],y[fExclBox-1],z[fExclBox-1],DX,DY);	
	  	}
	  	else
	  	{
	 		MyFit3D(x,y,z,Erx,Ery,Erz,x[fExclBox-1],y[fExclBox-1],z[fExclBox-1],DX,DY);	
		}

		Double_t pointX,pointY,pointZ;
				
		hx->Fill(DX);
		hy->Fill(DY);

	}// if one hit per box

    }
  

  fEvent++;
  return;
   
}


void TtTotAli::MyFit(Double_t *x,Double_t *y,Double_t *z,Double_t *Erx,Double_t *Ery,Double_t *Erz,Double_t realX, Double_t realY, Double_t realZ, Double_t &DELTAX, Double_t &DELTAY)
{

  TGraphErrors grX;
  TGraphErrors grY;

  Int_t ix=0,iy=0;

  for (Int_t j = 0 ; j < fNbox ; j++)
    {
      	//if (j == 0 ) std::cout << "Excluded plane: " << fExclBox << std::endl;
	//std::cout << j << "@(" << x[j] << "," <<  y[j] << "," <<  z[j] << ")" << std::endl;
      if (j == fExclBox-1) continue;

	//if (j==0 || j==1 || j==3 || j==5)
      if (TMath::Abs(Erx[j])<0.5) 
      {
	grX.SetPoint(ix,z[j],x[j]);
	ix++;
      }
      //if(j==0 || j==2 || j==4 || j==5)
      if (TMath::Abs(Ery[j])<0.5) 
      {
	grY.SetPoint(iy,z[j],y[j]);
	iy++;
      }

    }

  //  std::cout << "POINTS  " << ix << "   " << iy << std::endl;


  Double_t mmx,nx,mmy,ny;


  grX.Fit("pol1","Q");
  
  nx = (grX.GetFunction("pol1"))->GetParameter(0);
  mmx = (grX.GetFunction("pol1"))->GetParameter(1);

  grX.SetMarkerStyle(20);
  //grX.Draw("APL");

  grY.Fit("pol1","Q");

  ny = (grY.GetFunction("pol1"))->GetParameter(0);
  mmy = (grY.GetFunction("pol1"))->GetParameter(1);

  Double_t buffX, buffY;

  buffX = mmx * realZ + nx;
  buffY = mmy * realZ + ny;

  DELTAX = buffX - realX;
  DELTAY = buffY - realY;



}

void TtTotAli::MyFit3D(Double_t *x,Double_t *y,Double_t *z,Double_t *Erx,Double_t *Ery,Double_t *Erz,Double_t realX, Double_t realY, Double_t realZ, Double_t &DELTAX, Double_t &DELTAY)
{

  TGraph2D *gr2 = new TGraph2D();

  Int_t ixy=0;

  for (Int_t j = 0 ; j < fNbox ; j++)
    {
      //std::cout << j << "@(" << Erx[j] << "," <<  Ery[j] << "," <<  Erz[j] << ")" << std::endl;
      if (j == fExclBox-1) continue;

      gr2->SetPoint(ixy,x[j],y[j],z[j]);
      ixy++;
    }

    

  ///////////////////////////////////////

  Int_t Npoint = gr2->GetN();
  Double_t *xx = gr2->GetX();
  Double_t *yy = gr2->GetY();
  Double_t *zz = gr2->GetZ();

  TVector3 dir;
  dir.SetXYZ(xx[Npoint-1]-xx[0],yy[Npoint-1]-yy[0],zz[Npoint-1]-zz[0]);

  TVirtualFitter::SetDefaultFitter("Minuit");
  TVirtualFitter *min = TVirtualFitter::Fitter(0,4);
  min->SetObjectFit(gr2);
  min->SetFCN(*SumDistance2);//using local coordinate in FCN

  Double_t arglist[100];

  arglist[0] = 0;
  min->ExecuteCommand("SET NOWarnings",arglist,1);
  arglist[0] = -1;
  min->ExecuteCommand("SET PRINT",arglist,1);

  double pStart[4] = {xx[0],(dir.X())/(dir.Z()),yy[0],(dir.Y())/(dir.Z())};

  min->SetParameter(0,"X0",pStart[0],1,0,0);
  min->SetParameter(1,"A",pStart[1],1,0,0);
  min->SetParameter(2,"Y0",pStart[2],1,0,0);
  min->SetParameter(3,"B",pStart[3],1,0,0);

  
/*	  min->SetParameter(0,"y0",0,1,0,0);
	  min->SetParameter(1,"A",0,1,0,0);
	  min->SetParameter(2,"z0",0,1,0,0);
	  min->SetParameter(3,"B",0,1,0,0);*/


  arglist[0] = 1000; //number of functiona calls
  arglist[1] = 0.001; //tolerance
  min->ExecuteCommand("MIGRAD", arglist ,2);

   int nvpar,nparx;
   double amin,edm, errdef;
   min->GetStats(amin,edm,errdef,nvpar,nparx);
   //if(fVerbose>1) min->PrintResults(1,amin);

  Double_t fitPar[4];
  Double_t fitParErr[4];
  Double_t chi2 = amin/(2.*Npoint-4);

  fitPar[0] = min->GetParameter(0);
  fitPar[1] = min->GetParameter(1);
  fitPar[2] = min->GetParameter(2);
  fitPar[3] = min->GetParameter(3);

  fitParErr[0] = min->GetParError(0);
  fitParErr[1] = min->GetParError(1);
  fitParErr[2] = min->GetParError(2);
  fitParErr[3] = min->GetParError(3);

  Double_t mmx,nx,mmy,ny;
  Double_t buffX, buffY;

  buffX = fitPar[1] * realZ + fitPar[0];
  buffY = fitPar[3] * realZ + fitPar[2];

  DELTAX = buffX - realX;
  DELTAY = buffY - realY;



}

void TtTotAli::FinishRun()
{
 
  Double_t mmX,mmY,siX,siY;


		TF1 *fun = (TF1*) gROOT->GetFunction("gaus");
		TF1 *funSum = new TF1("funSum","gaus+gaus(3)");
		
		if (fExclBox-1 != 2)
		{
		funSum->SetParameter(0,100.);
 		funSum->SetParameter(1,hx->GetMean());
 		funSum->SetParameter(2,0.01);
 		funSum->SetParLimits(2,0.,10.);
	 	funSum->SetParameter(3,10.);
 		funSum->SetParameter(4,hx->GetMean());
 		funSum->SetParameter(5,0.1);
 		funSum->SetParLimits(5,0.,10.);
		}
		else
		{
		funSum->SetParameter(0,10.);
 		funSum->SetParameter(1,0.);
 		funSum->SetParameter(2,0.05);
 		funSum->SetParLimits(2,0.,10.);
	 	funSum->SetParameter(3,1.);
 		funSum->SetParameter(4,0.);
 		funSum->SetParameter(5,1.);
 		funSum->SetParLimits(5,0.,10.);
		}


 	 if (fSmallRange)
 	 {
		if (fDoubleGaus)
		{
			hx->Fit(funSum,"Q","",(hx->GetMean() - 0.1),(hx->GetMean() + 0.1));
		}
		else
		{
  			hx->Fit(fun,"Q","",(hx->GetMean() - 0.1),(hx->GetMean() + 0.1));
		}
 	}
  	else
  	{
		if (fDoubleGaus)
		{
			hx->Fit(funSum,"Q","",(hx->GetMean() - 2*(hx->GetRMS())),(hx->GetMean() + + 2*(hx->GetRMS())));
		}
		else
		{
			hx->Fit(fun,"Q","",(hx->GetMean() - 2*(hx->GetRMS())),(hx->GetMean() + 2*(hx->GetRMS())));
		}
  	}
  	//hx->Fit(fun,"Q");

  	if (fDoubleGaus)
  	{
  		mmX = funSum->GetParameter(1);
  		siX = funSum->GetParameter(2);
  	}
  	else
  	{
  		mmX = fun->GetParameter(1);
  		siX = fun->GetParameter(2);
	}

  	cout << "RMS: " << hx->GetRMS() << " sigma " << siX << endl;

		if (fExclBox-1 != 2)
		{
		funSum->SetParameter(0,100.);
 		funSum->SetParameter(1,hy->GetMean());
 		funSum->SetParameter(2,0.01);
 		funSum->SetParLimits(2,0.,10.);
	 	funSum->SetParameter(3,10.);
 		funSum->SetParameter(4,hy->GetMean());
 		funSum->SetParameter(5,0.1);
 		funSum->SetParLimits(5,0.,10.);
		}
		else
		{
		funSum->SetParameter(0,10.);
 		funSum->SetParameter(1,0.);
 		funSum->SetParameter(2,0.05);
 		funSum->SetParLimits(2,0.,10.);
	 	funSum->SetParameter(3,1.);
 		funSum->SetParameter(4,0.);
 		funSum->SetParameter(5,1.);
 		funSum->SetParLimits(5,0.,10.);
		}



  	if (fSmallRange)
  	{
		if (fDoubleGaus)
		{
        		hy->Fit(funSum,"Q","",(hy->GetMean() - 0.1),(hy->GetMean() + 0.1));
		}
		else
		{
			hy->Fit(fun,"Q","",(hy->GetMean() - 0.1),(hy->GetMean() + 0.1));
		}
  	}
  	else
  	{
		if (fDoubleGaus)
		{
        		hy->Fit(funSum,"Q","",(hy->GetMean() - 2*(hy->GetRMS())),(hy->GetMean() + 2*(hy->GetRMS())));
		}
		else
		{
			hy->Fit(fun,"Q","",(hy->GetMean() - 2*(hy->GetRMS())),(hy->GetMean() + 2*(hy->GetRMS())));
		}
  	}
  	//hy->Fit(fun,"Q");

  	if (fDoubleGaus)
  	{
  		mmY = funSum->GetParameter(1);
  		siY = funSum->GetParameter(2);
  	}
  	else
  	{
  		mmY = fun->GetParameter(1);
  		siY = fun->GetParameter(2);
  	}

  	cout << "X, mean: " << mmX << ", sig: " << siX << endl;
  	cout << "Y, mean: " << mmY << ", sig: " << siY << endl;

  	
    	if (fDoubleGaus)
	{
		delete funSum;
	}
	else
	{
		delete fun;
	}

  //if ((TMath::Abs(siX)) < 1.)
      	
     	 sX[fExclBox-1] += mmX;
     	 m_X[fExclBox-1] = mmX;
     	 sigX[fExclBox-1] = siX;
    		
  //else
  //  {
  //    sX[fExclBox-1] = 0.;
  //    m_X[fExclBox-1] = 0.;
  //    sigX[fExclBox-1] = 0.;
  //  }
  
  //if ((TMath::Abs(siY)) < 1.)
  //if (fAxis==2)
   
   	sY[fExclBox-1] += mmY;
      	m_Y[fExclBox-1] = mmY;
      	sigY[fExclBox-1] = siY;
   
  //else
  //  {
  //    sY[fExclBox-1] = 0.;
  //    m_Y[fExclBox-1] = 0.;
  //    sigY[fExclBox-1] = 0.;
  //  }
      
 
  
  std::cout << "Printing histos, excluded box: " << fExclBox << std::endl;  

  if (fPrint == 1)
    {
      std::cout << "fPrint == 1" << std::endl;  
      PrintHistos();
    }
  


//  hx->Reset();
//  hy->Reset();

   

  fEvent = 0;
}


void TtTotAli::PrintVal()
{

  cout << "SHIFTS" << endl;
  for (Int_t k = 0 ; k < fNbox ; k++)
    {
      cout << "X: " << sX[k] << "   Y: " << sY[k] << endl;
    }
}




void TtTotAli::PrintMeanResiduals()
{

  cout << "RESIDUALS" << endl;
  for (Int_t k = 0 ; k < fNbox ; k++)
    {
      cout << "X: " << m_X[k] << "   Y: " << m_Y[k] << endl;
    }
}

void TtTotAli::PrintSigmaResiduals()
{

  cout << "SIGMA-RESIDUALS" << endl;
  for (Int_t k = 0 ; k < fNbox ; k++)
    {
      cout << "sigX: " << sigX[k] << "   sigY: " << sigY[k] << endl;
    }
}

TVector2  TtTotAli::GetRes()
{
  
  Double_t resx,resy;
  Double_t bufff = 1.;
  for (Int_t aa = 0 ; aa < fNbox ; aa++) bufff=bufff*sigX[aa];
  resx = TMath::Power(bufff,1./(fNbox));
  bufff = 1.;
  for (Int_t bb = 0 ; bb < fNbox ; bb++) bufff=bufff*sigY[bb];
  resy = TMath::Power(bufff,1./(fNbox));
  
  return TVector2(resx,resy);

}

void TtTotAli::PrintHistos()
{

  gStyle->SetOptFit(1);

  TCanvas *can = new TCanvas();

  can->cd();

  hx->Draw();
  
  //hx->GetXaxis()->SetRangeUser((hx->GetMean())-8*(hx->GetRMS()),(hx->GetMean())+8*(hx->GetRMS()));
  hx->GetXaxis()->SetRangeUser(-1,+1);

  std::string nameX = Form("HistResidualsX_%d.png",fExclBox);

  can->Print(nameX.c_str(),"png");


  hy->Draw();
  
  //hx->GetXaxis()->SetRangeUser((hx->GetMean())-8*(hx->GetRMS()),(hx->GetMean())+8*(hx->GetRMS()));
  hy->GetXaxis()->SetRangeUser(-1,+1);
 
  std::string nameY = Form("HistResidualsY_%d.png",fExclBox);

  can->Print(nameY.c_str(),"png");

}



ClassImp(TtTotAli);