// // Simone Bianco 14/07/2010

// This Class' Header ------------------
#include "TtTotAliTask.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 <fstream>

using namespace ROOT::Math;
using namespace std;

TtTotAliTask::TtTotAliTask()
  : FairTask("Alignment")
{
  fTCandBranchName = "MVDHitsStrip";

  /*for (Int_t gg = 0 ; gg < 4 ; gg++)
    {
      sX[gg] = 0.;
      //if (gg==2) sX[gg] = -0.0453456;
      //if (gg==2) sX[gg] = -0.0254522;
      sY[gg] = 0.;
      //if (gg==2) sY[gg] = +0.130958;
      sigX[gg] = 0.;
      sigY[gg] = 0.;
      m_X[gg] = 0.;
      m_Y[gg] = 0.;
    }*/

  fPrint = 0;
  fAxis = 1;
  fEvent = 0;
  fSmallRange = kFALSE;
  fMean = kFALSE;
  fDoubleGaus = kFALSE;
  fUseGeo = kFALSE;

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

  rrX = new TH2F("rrX","rrX",5000,-5.,+5.,100,-2.,+2.);
  rrY = new TH2F("rrY","rrY",5000,-5.,+5.,100,-2.,+2.);

}

TtTotAliTask::TtTotAliTask(Int_t ExcludedBox)
  : FairTask("3D-Straight-Line-Fit")
{
  fTCandBranchName = "MVDHitsStrip";
 
  if (ExcludedBox < 1 || ExcludedBox > 4) 
    {
      std::cout << "Excluded box: Wrong value, setting as default 2!" << std::endl; 
      ExcludedBox = 2;
    }

  fExclBox = ExcludedBox;

  /*for (Int_t gg = 0 ; gg < 4 ; gg++)
    {
      sX[gg] = 0.;
      //if (gg==2) sX[gg] = -0.0453456;
      //if (gg==0) sX[gg] = +0.184629;
      //if (gg==2) sX[gg] = -0.0254522;
      sY[gg] = 0.;
      //if (gg==2) sY[gg] = +0.130958;
      //if (gg==2) sY[gg] = -0.0320571;
      sigX[gg] = 0.;
      sigY[gg] = 0.;
      m_X[gg] = 0.;
      m_Y[gg] = 0.;
    }*/

  fPrint = 0;

  fAxis = 1;
  fEvent = 0;
  fSmallRange = kFALSE;
  fMean = kFALSE;
  fDoubleGaus = kFALSE;
  fUseGeo = kFALSE;

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

  rrX = new TH2F("rrX","rrX",5000,-5.,+5.,100,-2.,+2.);
  rrY = new TH2F("rrY","rrY",5000,-5.,+5.,100,-2.,+2.);

}

TtTotAliTask::~TtTotAliTask()
{

}

InitStatus TtTotAliTask::Init()
{
  
  //Get ROOT Manager
  FairRootManager* ioman= FairRootManager::Instance();
  
  if(ioman==0)
    {
      Error("TtTotAliTask::Init","RootManager not instantiated!");
      return kERROR;
    }
  
  //Get input collection
  
  fTCandArray=(TClonesArray*) ioman->GetObject(fTCandBranchName);

  if(fTCandArray==0)
    {
      Error("TtTotAliTask::Init","trackcand-array not found!");
      return kERROR;
    }
  
  std::cout << "-I- TtTotAliTask: Initialisation successfull" << std::endl;

  fEvent = 0;

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


  /*
  if (!hx) 
  {
 	hx = new TH1F("hx","hx",50000,-5.,+5.);
	std::cout << "New hx" << std::endl;
  }  
  else 
  {
        hx->Reset();
	std::cout << "hx reset" << std::endl;
  }

  if (!hy)  
  {
	hy = new TH1F("hy","hy",50000,-5.,+5.);
	std::cout << "New hy" << std::endl;
  }
  else 
  {
	hy->Reset();
	std::cout << "hx reset" << std::endl;
  }
  */

  for (Int_t gg = 0 ; gg < fNbox ; 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.);
  }


  return kSUCCESS;
}

Double_t TtTotAliTask::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 TtTotAliTask::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 TtTotAliTask::Exec(Option_t* opt)
{

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


  if (fEvent == 0)
  {

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

	delete rrX;
	delete rrY;
	rrX = new TH2F("rrX","rrX",5000,-5.,+5.,100,-2.,+2.);
  	rrY = new TH2F("rrY","rrY",5000,-5.,+5.,100,-2.,+2.);

  }

  Int_t ntcand=fTCandArray->GetEntriesFast();

//  std::cout << "Event: " << fEvent << ", size: " << ntcand << std::endl; 

  if (ntcand >= fNbox)
    {
      
      //  std::cout << "Entries: " << ntcand << std::endl;

 
      std::map<Double_t,Int_t> SensorsPos;

       
      //    std::cout<<"TtTotAliTask::Exec"<<std::endl;
      // Reset output Array
      if(fTCandArray==0) Fatal("TtTotAliTask::Exec","No Points");
   

      Int_t name;

    
      // Detailed output
      //   std::cout<<" -I- TtTotAliTask: contains "<<ntcand<<" TCandidates"<<std::endl;
	
      //std::cout<< " Detailed Debug info on the candidates:"<<std::endl;
    
      //for(Int_t itr=0;itr<ntcand;++itr){

      Int_t itrr = 0;

      while (SensorsPos.size() < fNbox && itrr < ntcand)
	{
		PndSdsHit *theHitP = (PndSdsHit*)   fTCandArray->At(itrr);


		if(theHitP)
		{

			//theHitP->Print();

			name = theHitP->GetSensorID();
	
			SensorsPos[theHitP->GetZ()] = name;
		}

		itrr++;

	}
    
	Bool_t test[fNbox];

	for (Int_t vv = 0 ; vv < fNbox ; vv++)
	{
		test[vv] = kFALSE;
	}

      for(Int_t itr=0;itr<ntcand;++itr){

	//PndSdsHit *theHit = (PndSdsHit*)   fTCandArray->At(itr);
	PndSdsHit *theHit = (PndSdsHit*)   fTCandArray->At(itr);
	
	Int_t buffDet = theHit->GetSensorID();
	test[buffDet] = kTRUE;
      }


      //  const Int_t sizeMap = Sensors.size();
      
      const Int_t sizeMap = SensorsPos.size();
  
      Int_t DetNames[sizeMap];
      Double_t Pos[sizeMap];
      Int_t jj = 0;

      if (SensorsPos.size()>=fNbox)
	{
	  //std::cout << "Number of Sensors:" << SensorsPos.size() << std::endl;
	  for (std::map<Double_t,Int_t>::iterator it=SensorsPos.begin();it!=SensorsPos.end();++it)
	    {
	      if (fEvent < 10) std::cout << "position: " << it->first << " name: " << (it->second) << std::endl;
	  
	      DetNames[jj] = it->second;
	      Pos[jj] = it -> first;
	      jj++;
	    }

	}

	//std::cout << "Event: " << fEvent << ", size: " << SensorsPos.size() << ", t1: " << test1 << ", t2: " << test2 << ", t3: " << test3 << std::endl;  
 
     
    //   if (fEvent < 5)
// 	{
      
// 	  for (Int_t l = 0 ; l < sizeMap ; l++)
// 	    {
// 	      std::cout << DetNames[l] << std::endl;
// 	    }
// 	}
    

      ///////////////
      
	  
      Bool_t sat=kTRUE;

      for (Int_t kk = 0 ; kk < fNbox ; kk++)
      {
	sat = sat * test[kk];
      }
	
      
      if ((SensorsPos.size()>=fNbox) && sat )
	{
	  Int_t DetRem[sizeMap-1];
	  for (Int_t ss = 0 ; ss < fNbox ; ss++) // just for this test, restore the previous value (sizeMap-1) at the end
	    {

	      DetRem[ss] = 0;

	    }
  
  
	  Int_t counter = 0;

	  Double_t BuffZ = -9999.;

	  for (Int_t kk = 0 ; kk < (sizeMap) ; kk++)
	    {

	      if ((fExclBox-1) == kk) 
		{
		  BuffZ = Pos[kk];
		  continue;
		}
	      else
		{
		  DetRem[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;

	
	  Int_t ihit = 0;

	  //std::cout << "Event: " << fEvent << std::endl;
	
	  for (Int_t it1 = 0 ; it1 < ntcand ; it1++)
	    
	    {
		
	      //PndSdsHit *theHit = (PndSdsHit*)   fTCandArray->At(it1);
	      PndSdsHit *theHit = (PndSdsHit*)   fTCandArray->At(it1);
	      
              if (theHit)
	      {	


	      for(Int_t ww = 0 ; ww < fNbox ; ww++)
	      {
		if ((theHit->GetSensorID()) == DetNames[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;
		}
	      }
	      
	      ihit++;
	      }
		
	    } // end loop on points
	  
	  
	  // setting the points
	  
	  Int_t uu = 0;

	  // to check alignment...

	  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];
		}
	    }

	
	// until here


	Double_t y1_low = -0.3;
	Double_t y1_hi = +0.2;
	Double_t x1_low = -0.5;
	Double_t x1_hi = +0.;

	//Double_t y2_low = -0.7;
	//Double_t y2_hi = +0.3;
	//Double_t x2_low = -0.7;
	//Double_t x2_hi = +0.7;

	Bool_t spot1 = kFALSE;
	//Bool_t spot2 = kFALSE;
	 
	if ( (x[0] > x1_low) && (x[0] < x1_hi) &&  (y[0] > y1_low) && (y[0] < y1_hi))  // here the excluded box has not yet been removed, so this check is always on box1
		{	
			spot1 = kTRUE;
		}
	
	//	if ( (x2 > x2_low) && (x2 < x2_hi) &&  (y2 > y2_low) && (y2 < y2_hi) )
	//	{
	//		spot2 = kTRUE;
	//	}

	//if (spot1)// && spot2)
	{

	
 	  Double_t DX,DY;

	  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);

	  rrX->Fill(DX,x[fExclBox-1]);
	  rrY->Fill(DY,y[fExclBox-1]);
      
	} // if spot1 and 2

	}// if 3 hits
      
      
      
    }
  

  fEvent++;
  return;
   
}


void TtTotAliTask::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++)
    {
      //std::cout << j << "@(" << Erx[j] << "," <<  Ery[j] << "," <<  Erz[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);


  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 TtTotAliTask::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 TtTotAliTask::FinishTask()
{
 
  Double_t mmX,mmY,siX,siY;

  if (fMean)
  {

	// too sensible to binning

	hx->GetXaxis()->SetRangeUser(hx->GetMean() - 0.01,hx->GetMean() + 0.01);
	hy->GetXaxis()->SetRangeUser(hy->GetMean() - 0.01,hy->GetMean() + 0.01);
	mmX = hx->GetMean();
	mmY = hy->GetMean();

	//mX = hx->GetBinCenter(hx->GetMaximumBin());

  	siX = hx->GetRMS();

	//mY = hy->GetBinCenter(hy->GetMaximumBin());

  	siY = hy->GetRMS();

	sX[fExclBox-1] += mmX;
	m_X[fExclBox-1] = mmX;
	sigX[fExclBox-1] = siX;
	sY[fExclBox-1] += mmY;
	m_Y[fExclBox-1] = mmY;
	sigY[fExclBox-1] = siY;

	hx->GetXaxis()->SetRangeUser(-5.,+5.);
	hy->GetXaxis()->SetRangeUser(-5.,+5.);


  }
  else
  {

		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.)
  //if (fAxis==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 TtTotAliTask::PrintVal()
{

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




void TtTotAliTask::PrintMeanResiduals()
{

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

void TtTotAliTask::PrintSigmaResiduals()
{

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

TVector2  TtTotAliTask::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 TtTotAliTask::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(-5,+5);

  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()));
  hx->GetXaxis()->SetRangeUser(-5,+5);
 
  std::string nameY = Form("HistResidualsY_%d.png",fExclBox);

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

  std::string name2DX = Form("HistResiduals2DX_%d.png",fExclBox);
  
  can->cd();

  rrX->Draw("COLZ");

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

  std::string name2DY = Form("HistResiduals2DY_%d.png",fExclBox);
  
  can->cd();

  rrY->Draw("COLZ");

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

}



ClassImp(TtTotAliTask);