#include<iostream>
#include<cmath>
#include<iomanip>
#include<string>
#include "functions_hadcalib.h"
#include "TMath.h"

using namespace std;

//////////////////////////////////////////////
// Function file for hadronic calibration   //
//        using the numeric method          //
//////////////////////////////////////////////


void printvector(double N_vector[], TString ttext, const int N)
{
 cout<<ttext<<" = {";
 int i;
 for(i=0;i<(N-1);i++)
   { cout<<N_vector[i]<<", ";   }
cout<<N_vector[N-1]<<"};"<<endl;
}

void printvector(int N_vector[], TString ttext, const int N)
{
 cout<<ttext<<" = {";
 int i;
 for(i=0;i<(N-1);i++)
   { cout<<N_vector[i]<<", ";   }
cout<<N_vector[N-1]<<"};"<<endl;
}


void logaxisbins(const Int_t nbins, Double_t lowbin, Double_t upbin, Double_t* logbins)
{
  // Gives visually equidistant bins in a logarithmic histogram

  Double_t ln10=log(10.0);
  Double_t Lmin =log(lowbin)/ln10;
  Double_t Lmax =log(upbin)/ln10;
  Double_t dL = (Lmax-Lmin)/(Double_t(nbins));
  Double_t L;

  logbins[0] = lowbin;
  logbins[nbins] = upbin;
  for(int i=1;i<nbins;i++)
    { L=Lmin+i*dL;
      logbins[i]=TMath::Power(10.0,L);
    }
}

Int_t find_seedcluster(Int_t Nclusters, Float_t* cluster_energy)
{
  // Find maximum-energy cluster (seed cluster for the cone)
  // Syntax:
  // seed=find_seedcluster(<number of clusters>, <cluster energy array>)

  Double_t Emax = -100000., Ecluster;
  Int_t seedcluster, s;
  for(int i=0;i<Nclusters;i++) 
    { Ecluster = cluster_energy[i];
    if(Ecluster > Emax)
      { Emax = Ecluster;
        seedcluster = i;
      }
    }
  s=seedcluster;
  return s;
}

Double_t make_cone(Int_t Nclusters, Float_t* cluster_energy, Float_t* eta_cl, Float_t* phi_cl, Double_t coneR, Int_t* coneclusters, Int_t* cellcluster_topo, Int_t* Nclusters_in_cone)
{
  //Make the cone and return the cone energy (energy of clusters within the cone)

  Int_t seedcluster = find_seedcluster(Nclusters, cluster_energy);
  Double_t eta0=eta_cl[seedcluster]; //eta and phi for the seed cluster
  Double_t phi0=phi_cl[seedcluster];
  Double_t eta, phi, e, p, testR;
  Double_t Econe=0., Ecluster; 
  Int_t n=0;
  for(int i=0;i<Nclusters;i++) //Make a CONE
    { Ecluster = cluster_energy[i];
    //Eallclusters=Eallclusters+Ecluster;
      eta=eta_cl[i];
      phi=phi_cl[i];
      e=eta-eta0;
      p=phi-phi0;
      testR=sqrt(e*e+p*p);
      if((testR<coneR) || (testR==coneR))
        // cone condition
        { Econe=Econe+Ecluster;
	  coneclusters[n]=i;
	  //*(coneclusters+n)=cellcluster_topo[i];
          n++;
	}
    } //end for(int i=0;i<Nclusters;i++)
  *Nclusters_in_cone=n;
  return Econe;
}


Double_t make_cones(Int_t Nclusters, Float_t* cluster_energy, Float_t* eta_cl, Float_t* phi_cl, Double_t coneR, Int_t* coneclusters, Int_t* cellcluster_topo, Int_t* Nclusters_in_cone, Int_t* excluded_clusters)
{
  //Make the cone and return the cone energy 
  //(energy of clusters within the cone)
  //Doesn't include clusters in other cones


  /////Find the seedcluster (cluster with maximum energy)
  // of the allowed clusters
  Double_t Emax = -100000., Ecluster;
  Int_t seedcluster, ncl, im;
  Int_t coneveto=1, ec, ecstart=0;
  for(int ncl=0;ncl<Nclusters;ncl++) 
    { 
      coneveto=1;
      for(im=0;im<Nclusters;im++)
	{ ec = excluded_clusters[im];
	  if(ec==ncl) //true if this cluster belongs to excluded_clusters
	    { coneveto=0; 
	    }
	}//end for(im=0;im<Nclusters;im++)

      if(coneveto!=0)
	{ //Do not search for seedcluster if cluster
	  // belongs to excluded_clusters
	  Ecluster = cluster_energy[ncl];
	  if(Ecluster > Emax)
	    { Emax = Ecluster;
	      seedcluster = ncl;
	    }
	}//end if(coneveto!=0)
    }//end for(int ncl=0;ncl<Nclusters;ncl++)
  //cout<<"  seedcluster = "<<seedcluster<<endl;

  //Find how many clusters there are in excluded_clusters
  for(im=0;im<Nclusters;im++)
    { ec = excluded_clusters[im];
      if(ec!=-1)
	{ ecstart=im+1;
	}
    }

  //////Making the cone
  Double_t eta0=eta_cl[seedcluster]; //eta and phi for the seed cluster
  Double_t phi0=phi_cl[seedcluster];
  Double_t eta, phi, e, p, testR;
  Double_t Econe=0.; 
  Int_t n=0;
  for(int ncl=0;ncl<Nclusters;ncl++)
    { 
      coneveto=1;
      for(im=0;im<Nclusters;im++)
	{ ec = excluded_clusters[im];
	  if(ncl==ec)
	    { coneveto=0; }
	}
 
      if(coneveto!=0)
	{ //Do not test for cone if cluster
	  // belongs to excluded_clusters
	  Ecluster = cluster_energy[ncl];
	  eta=eta_cl[ncl];
	  phi=phi_cl[ncl];
	  e=eta-eta0;
	  p=phi-phi0;
	  testR=sqrt(e*e+p*p);
	  if((testR<coneR) || (testR==coneR))   // cone condition
	    { Econe=Econe+Ecluster;
	      coneclusters[n]=ncl;
	      n++;
	    }
        }//end if(coneveto!=0)
    }//end for(int ncl=0;ncl<Nclusters;ncl++)


  //Add the found cone clusters to excluded_clusters
  for(im=0;im<n;im++)
    { excluded_clusters[im+ecstart]=coneclusters[im];
    }


  *Nclusters_in_cone=n;
  return Econe;
}


Double_t get_weight(TProfile2D *weights_table, Double_t Econe, Double_t E_dens)
{
  
  //Find histo table limits
  Double_t xlow, ylow, xup, yup;
  Int_t nbinsx, nbinsy;
  Double_t xbins[50], ybins[50];
  nbinsx = weights_table->GetXaxis()->GetNbins();
  nbinsy = weights_table->GetYaxis()->GetNbins();
  TAxis *xaxis = weights_table->GetXaxis();
  TAxis *yaxis = weights_table->GetYaxis();

  for(int i=1;i<nbinsx+2;i++)//Bin 0 contains underflow, start at 1
    { xlow = xaxis->GetBinLowEdge(i);  
      xbins[i-1] = xlow;
    }
  for(int i=1;i<nbinsy+2;i++)
    { ylow = yaxis->GetBinLowEdge(i);
      ybins[i-1] = ylow;
    }
 

  //Find weight
  Double_t w=0., weight=0.; 
  Int_t binl;
  Double_t x_low, x_up, y_low, y_up;
  for(int nx=0; nx<nbinsx+1; nx++)
    { x_low=xbins[nx]; x_up=xbins[nx+1];
      if(Econe>=x_low && Econe<x_up)
	{ for(int ny=0; ny<nbinsy+1; ny++)
	    { y_low=ybins[ny]; y_up=ybins[ny+1];
	      if(E_dens>=y_low && E_dens<y_up)
		{ binl=weights_table->GetBin(nx+1, ny+1);  
                    // "global" bin no
		  w=weights_table->GetBinContent(binl);
		}//end if(E_dens>=y_low...)
	    }//end for(int i=0;i<nbinsy;i++)
	}//end if(Econe>=x_low && Econe<x_up)
    }//end for(int i=0;i<nbinsx+1;i++)

  if(w!=0)
    {weight = w;}
  else
    {weight = 1;}

return weight;
}


Double_t get_weight_old(TH2F *weights_table, Double_t Econe, Double_t E_dens)
{
  
  //Find histo table limits
  Double_t xlow, ylow, xup, yup;
  Int_t nbinsx, nbinsy;
  Double_t xbins[50], ybins[50];
  nbinsx = weights_table->GetXaxis()->GetNbins();
  nbinsy = weights_table->GetYaxis()->GetNbins();
  TAxis *xaxis = weights_table->GetXaxis();
  TAxis *yaxis = weights_table->GetYaxis();

  for(int i=1;i<nbinsx+2;i++)//Bin 0 contains underflow, start at 1
    { xlow = xaxis->GetBinLowEdge(i);  
      xbins[i-1] = xlow;
    }
  for(int i=1;i<nbinsy+2;i++)
    { ylow = yaxis->GetBinLowEdge(i);
      ybins[i-1] = ylow;
    }
 

  //Find weight
  Double_t w=0., weight=0.; 
  Int_t binl;
  Double_t x_low, x_up, y_low, y_up;
  for(int nx=0; nx<nbinsx+1; nx++)
    { x_low=xbins[nx]; x_up=xbins[nx+1];
      if(Econe>=x_low && Econe<x_up)
	{ for(int ny=0; ny<nbinsy+1; ny++)
	    { y_low=ybins[ny]; y_up=ybins[ny+1];
	      if(E_dens>=y_low && E_dens<y_up)
		{ binl=weights_table->GetBin(nx+1, ny+1);  
                    // "global" bin no
		  w=weights_table->GetBinContent(binl);
		}//end if(E_dens>=y_low...)
	    }//end for(int i=0;i<nbinsy;i++)
	}//end if(Econe>=x_low && Econe<x_up)
    }//end for(int i=0;i<nbinsx+1;i++)

  if(w!=0)
    {weight = w;}
  else
    {weight = 1;}

return weight;
}

void estimate_EMfraction(Float_t* cl_EM_frac, Float_t* cl_e_topo, Float_t* cl_et_topo, Float_t* cl_m2_lambda_topo, Float_t* cl_m2_r_topo, Float_t* cl_center_lambda_topo, Int_t Nclusters, Double_t Ecl_cut, Double_t* c)
{
  //The variable cl_EM_frac will be filled with the estimated
  // electromagnetic fractions of the cluster energies.
  //  If Et_cluster < Ecl_cut then this indicates that 
  //  cl_EM_frac[i] = -100

  //Double_t par_l2[2] = {100.96, -89.21};
  //Double_t par_r2[2] = {21.09, -15.83};
  //Double_t par_lc[2] = {203.57, -188.72};
  //Double_t par_l2[2] = {44.62, -30.87};
  //Double_t par_r2[2] = {22.08, -14.33};
  //Double_t par_lc[2] = {132.82, -106.83};
  //Double_t par_l2[2] = {88.32, -75.92};
  //Double_t par_r2[2] = {19.77, -13.54};
  //Double_t par_lc[2] = {169.97, -147.67}; 

  //Double_t par_l2[2] = {77.18, -66.79};//Et_cluster cut > 1000 MeV
  //Double_t par_r2[2] = {13.92, -7.15};
  //Double_t par_lc[2] = {160.97, -148.63};
  Double_t par_l2[2] = {82.29, -72.72}; //Et_cluster cut > 2000 MeV
  Double_t par_r2[2] = {15.02, -8.76};
  Double_t par_lc[2] = {163.00, -150.15};
  //Double_t par_l2[2] = {84.79, -75.38};//Et_cluster cut > 2500 MeV
  //Double_t par_r2[2] = {15.36, -9.32};
  //Double_t par_lc[2] = {159.07, -146.69};

  Double_t Ecluster,Etcluster, Ecllog, lambda2, r2, lambdac;
  Double_t EMest, EMestm, EMl2, EMr2, EMlc;
  //Double_t c[3]={0.186,0.151,0.663};//lc, l2, r2

  for(int ncl=0;ncl<Nclusters;ncl++)
    { cl_EM_frac[ncl]=-100.;
      Ecluster  = cl_e_topo[ncl];
      Etcluster = cl_et_topo[ncl];
      //Eemcluster= cl_Eem[ncl];
      //Eemfrac   = Eemcluster/Ecluster;
      Ecllog    = log(Ecluster);
      lambda2 = cl_m2_lambda_topo[ncl];
      r2      = cl_m2_r_topo[ncl];
      lambdac = cl_center_lambda_topo[ncl];

      if(Etcluster>Ecl_cut)
	{
	  Double_t l2E = sqrt(lambda2)/Ecllog;
	  Double_t r2E = sqrt(r2)/Ecllog;
	  Double_t lcE = lambdac/Ecllog;
	  EMl2 = (l2E - par_l2[0])/par_l2[1];
	  EMr2 = (r2E - par_r2[0])/par_r2[1];
	  EMlc = (lcE - par_lc[0])/par_lc[1];
	  EMest = c[0]*EMlc + c[1]*EMl2 + c[2]*EMr2;
	  cl_EM_frac[ncl] = EMest;
	  //EMest_vs_EMfrac->Fill(Eemfrac,EMest);
	  //EMest_vs_EMfrac_mean->Fill(Eemfrac,EMestm);

	}

    }//end for(int ncl=0;ncl<Nclusters;ncl++)

}

void CalibHit_cells_to_cluster(Float_t* cell_var, Float_t* cellweight_topo, Double_t* cluster_var, Int_t* cellcluster_topo,Int_t Ncells, Int_t Nclusters)
{ //Loops over the entries in the variable cell_var, and sums them 
  //according to cluster.

  Int_t cluster_id;
  Double_t cellv;
  for(int i=0;i<Nclusters;i++)
    { cluster_var[i]=0.;
    }

  for(int i=0;i<Ncells;i++)
    { cluster_id = cellcluster_topo[i]-1;
      cellv=cell_var[i]*cellweight_topo[i];
      cluster_var[cluster_id]+=cellv;
    }
}

Double_t EMfrac_from_histo(TProfile2D *EM_histo, Double_t log_densmom, Double_t log_lambdac, char opt, Double_t EMnoise)
{
  
  //Find EM_histo axes
  Double_t xlow, ylow, xup, yup;
  Int_t nbinsx, nbinsy;
  Double_t xbins[50], ybins[50];
  nbinsx = EM_histo->GetXaxis()->GetNbins();
  nbinsy = EM_histo->GetYaxis()->GetNbins();
  TAxis *xaxis = EM_histo->GetXaxis();
  TAxis *yaxis = EM_histo->GetYaxis();

  for(int i=1;i<nbinsx+2;i++)//Bin 0 contains underflow, start at 1
    { xlow = xaxis->GetBinLowEdge(i);  
      xbins[i-1] = xlow;
    }
  for(int i=1;i<nbinsy+2;i++)
    { ylow = yaxis->GetBinLowEdge(i);
      ybins[i-1] = ylow;
    }
 

  //Find EM fraction
  //Double_t EMf=0., EMfrac=0., sEMf=0., mEMf; 
  Double_t EMf=-0.15, EMfrac=0., sEMf=0., mEMf;
  Int_t binl;
  Double_t x_low, x_up, y_low, y_up;
  for(int nx=0; nx<nbinsx+1; nx++)
    { x_low=xbins[nx]; x_up=xbins[nx+1];
      if(log_densmom>=x_low && log_densmom<x_up)
	{ for(int ny=0; ny<nbinsy+1; ny++)
	    { y_low=ybins[ny]; y_up=ybins[ny+1];
	      if(log_lambdac>=y_low && log_lambdac<y_up)
		{ binl=EM_histo->GetBin(nx+1, ny+1); // "global" bin no
		  EMf=EM_histo->GetBinContent(binl);
		  EM_histo->SetErrorOption("s");
		  sEMf=EM_histo->GetBinError(binl);
		}//end if(log_lambdac>=y_low...)
	    }//end for(int i=0;i<nbinsy;i++)
	}//end if(log_densmom>=x_low && log_densmom<x_up)
    }//end for(int i=0;i<nbinsx+1;i++)

  /*switch(opt)
    { case 's': EMfrac = EMf+sEMf; break;
	//case 'm': EMfrac = EMf+mEMf; break;
      default: EMfrac = EMf; break;
    }
  */
  if(opt=='s')
    {
      EMfrac = EMf+sEMf;
      //cout<<"opt 's'!"<<endl;
    }
  else
    {
      EMfrac = EMf;
    }
  
return EMfrac;
}
/*
Double_t EMfrac_from_histospread(TProfile2D *EM_histo, Double_t log_densmom, Double_t log_lambdac)
{
  
  //Find histo table limits
  Double_t xlow, ylow, xup, yup;
  Int_t nbinsx, nbinsy;
  Double_t xbins[50], ybins[50];
  nbinsx = EM_histo->GetXaxis()->GetNbins();
  nbinsy = EM_histo->GetYaxis()->GetNbins();
  TAxis *xaxis = EM_histo->GetXaxis();
  TAxis *yaxis = EM_histo->GetYaxis();

  for(int i=1;i<nbinsx+2;i++)//Bin 0 contains underflow, start at 1
    { xlow = xaxis->GetBinLowEdge(i);  
      xbins[i-1] = xlow;
    }
  for(int i=1;i<nbinsy+2;i++)
    { ylow = yaxis->GetBinLowEdge(i);
      ybins[i-1] = ylow;
    }
 

  //Find EM fraction
  Double_t EMf=0., EMfrac=0., sEMf=0.; 
  Int_t binl;
  Double_t x_low, x_up, y_low, y_up;
  for(int nx=0; nx<nbinsx+1; nx++)
    { x_low=xbins[nx]; x_up=xbins[nx+1];
      if(log_densmom>=x_low && log_densmom<x_up)
	{ for(int ny=0; ny<nbinsy+1; ny++)
	    { y_low=ybins[ny]; y_up=ybins[ny+1];
	      if(log_lambdac>=y_low && log_lambdac<y_up)
		{ binl=EM_histo->GetBin(nx+1, ny+1);// "global" bin no
		  EMf=EM_histo->GetBinContent(binl);
		  EM_histo->SetErrorOption("s");
		  sEMf=EM_histo->GetBinError(binl);
		}//end if(log_lambdac>=y_low...)
	    }//end for(int i=0;i<nbinsy;i++)
	}//end if(log_densmom>=x_low && log_densmom<x_up)
    }//end for(int i=0;i<nbinsx+1;i++)

  EMfrac = EMf+sEMf;

return EMfrac;
}
*/

void cellcluster_array(Int_t* cellcluster, Int_t Nclusters, Int_t* ncells_topoC)
{
  int Ncells_icluster;
  int Nn=0;

  for(int i=0;i<Nclusters;i++)
    { Ncells_icluster = ncells_topoC[i]; //Number of cells in this cluster
      for(int n=0; n<Ncells_icluster; n++)
	{ cellcluster[Nn+n] = i+1;
	}
      Nn+=Ncells_icluster;

    }
}