///////////////////////////////////////////////////////////////////////////////
//
// $Id: LogWeights.hpp,v 1.7 2002/04/16 16:55:26 buescher Exp $
// 
// File: LogWeight.cpp
//
// Purpose:  supply information for EM cluster
//
// Created:    04-May-2001 Volker Buescher based on em_evt packages.
//
// History:  06-OCT-2001 DooKee Cho added w_cc(pt in cc) 
//                                  and Z correction(in cc) using NN (calo+cps)
// History:  18-OCT-2001 DooKee Cho added w_cc(pt in cc) 
//                                  and Z correction(in cc) using NN (caloonly)
// History:  14-JAN-2002 DooKee Cho small correction(thanks to Scott)
//
// History:  1-APR-2002 DooKee Cho fixed a bug(pt>155.0, use same correction as pt=155.0)
//                                  
//             
///////////////////////////////////////////////////////////////////////////////
#include "caladdress/Floor.hpp"
#include <cmath> 

namespace emid{

float weight_cut=5.5;     // weight cut for log energy weights
                          // used in run I
struct CCZCORR{
  float w0[5];
  float zcoeff[6];
  float ptcoeff[4];
};

struct ECRCORR{
  float w0[5];
  float rcoeff[6];
  float ptcoeff[4];
};


  CCZCORR cczcorr[4]={
  {{4.8709,-0.82507E-03,0.0,0.0,4.6},
   {-0.1020,0.1080E-01,0.1416E-04,0.3094E-06,0.0, 0.0},
   {-0.4917E-02, 0.7795E-04, 0.0, 0.0}},
  {{5.5500,0.30263E-02,0.0,0.0,6.0},
   {-0.58428E-01, 0.13349E-01, 0.89920E-05, -0.27759E-06,0.0, 0.0},
   {-0.4056E-02, 0.6342E-04, 0.0, 0.0}},
  {{4.2190,0.41364E-01,-0.39400E-03,0.11923E-05,6.0},
   {-0.10680E-01, 0.22974E-02, 0.27091E-04, -0.84591E-06,-0.29821E-08, 0.0},
   {-0.3649E-02, 0.5445E-04, 0.0, 0.0}},
  {{4.2182,0.18061E-01,-0.77280E-04,0.0,4.0},
   {-0.26135E-01, -0.18860E-01, 0.78556E-04, -0.32333E-05,-0.66184E-08,  0.53690E-10},
   {-0.4583E-02, 0.4794E-04, 0.0, 0.0}}
  }; 


  ECRCORR ecrcorr[4]={
  {{3.7,0.0,0.0,0.0,0.0},
   {-0.88547E-01,0.11576E-01,0.35242E-03,-0.65128E-05,0.0,0.0},
   {0.0, 0.0, 0.0, 0.0}},
  {{5.1,0.0,0.0,0.0,0.0},
   {0.46620E-01,0.76784E-02,0.25547E-03,-0.41281E-05,0.0,0.0},
   {0.0, 0.0, 0.0, 0.0}},
  {{4.0,0.0,0.0,0.0,0.0},
   {0.34149,-0.16705E-01,0.53222E-03,-0.38603E-05,0.0,0.0},
   {0.0, 0.0, 0.0, 0.0}},
  {{3.7,0.0,0.0,0.0,0.0},
   {-0.15962,0.24181E-01,-0.41633E-03,0.13662E-05,0.0,0.0},
   {0.0, 0.0, 0.0, 0.0}}
   }; 


  inline float LogWeight(int icorr,float z,Floor flr,float ECell,float EFloor,float pt);
  inline float LogZCorrCC(int icorr,int i,float z,float emf1,float emf2,float emf3,float emf4,float pt);
  inline float LogRCorrEC(int icorr,int i,float z,float ecr,float emf1,float emf2,float emf3,float emf4,float pt);

  inline float sigmoid(float x);
  inline float corr_cc(int i,float z,float emf1,float emf2,float emf3,float emf4,float pt);
}

inline float emid::LogWeight(int icorr, float z, Floor flr, float ECell, float EFloor, float pt)
{
  float logE=log(ECell/EFloor);
  float wt=emid::weight_cut+logE;
  if(std::abs(z)<=130.&& flr>0 && flr<5){ // cc recalculate weight
    if(std::abs(pt)>155.0) pt=155.0;
      if(icorr==1) {         // Fixed W_0 based on 50 GeV
	  wt=cczcorr[flr-1].w0[4]+logE;
      }else if(icorr==2) {  // W_0(pt)
	  wt=cczcorr[flr-1].w0[0]+cczcorr[flr-1].w0[1]*pt
	      +cczcorr[flr-1].w0[2]*pt*pt+cczcorr[flr-1].w0[3]*pt*pt*pt+logE;
      }
  }else if(std::abs(z)>=160.0&& flr>0 && flr<5){    //ec recalculate weight
    if(std::abs(pt)>155.0) pt=155.0;
      wt=ecrcorr[flr-1].w0[0]+logE;
  }
  return wt;
}

inline float emid::LogZCorrCC(int icorr,int i,float z,float emf1,float emf2,float emf3,float emf4,float pt)
{
  float z_corrected=z;
  if(std::abs(z)<=130.){
    if(std::abs(pt)>155.0) pt=155.0;
    if(icorr==1) {
      float dz=0;
      float zpower=1.0;
    
      // z dependent correction
      for(int k=0; k<6; ++k){
	dz+=cczcorr[i-1].zcoeff[k]*zpower;
	zpower=zpower*z;
      }
      z_corrected=z-dz;

      // pt dependent correction 
      float ptpower=1.0;
      dz=0;
      for(int k=0; k<4; ++k){
	dz+=cczcorr[i-1].ptcoeff[k]*ptpower;
	ptpower=ptpower*pt;
      }
      z_corrected=z_corrected*(1.0-dz);
    }else if(icorr==2) {
      z_corrected=z+corr_cc(i,z,emf1,emf2,emf3,emf4,pt)*5.0;
    }
  }
  return z_corrected;
}


inline float emid::LogRCorrEC(int icorr,int i,float z,float ecr,float emf1,float emf2,float emf3,float emf4,float pt)
{
  float r_corrected=ecr;
  if(std::abs(z)>160.0){
    if(std::abs(pt)>155.0) pt=155.0;

    float dr=0;
    float rpower=1.0;

    // r dependent correction
    for(int k=0; k<6; ++k){
      dr+=ecrcorr[i-1].rcoeff[k]*rpower;
      rpower=rpower*ecr;      
    }
    r_corrected=ecr-dr;
  }
  return r_corrected;
}

inline float emid::sigmoid(float x)
{
  if(x>37.) return 1.;
  else if(x<-37.) return 0.;
  else return 1./(1.+exp(-x));
}

inline float emid::corr_cc(int i,float z,float emf1,float emf2,float emf3,float emf4,float pt)
{
  float out1,out2,out3,out4,out5,out6,out7,out8;
  float rin1,rin2,rin3,rin4,rin5,rin6,rin7,rin8;

  out1 = z/125.;
  out2 = emf1;
  out3 = emf2;
  out4 = pt/150.;

  if(i==1) {

//     calo+cps weight
//     rin1 = 1.097293e+00+(-7.590874e-02)*out1+(6.889727e+01)*out2+(-2.680817e-01)*out3+(1.063754e-01)*out4;
//     rin2 = -3.117768e+00+(2.662246e-01)*out1+(5.318740e+01)*out2+(2.470878e+01)*out3+(2.342771e+00)*out4;
//     rin3 = 2.013535e+01+(1.406912e+01)*out1+(-3.808360e-01)*out2+(-8.853319e+00)*out3+(1.669219e+00)*out4;
//     rin4 = 2.980056e+00+(-1.504779e-01)*out1+(-7.245010e+01)*out2+(-1.455479e+01)*out3+(-2.337015e+00)*out4;
//     rin5 = 1.944731e+01+(-1.605394e+01)*out1+(-3.864603e+00)*out2+(5.376479e+00)*out3+(4.673327e-01)*out4;
//     rin6 = 2.299450e+00+(-3.831722e-02)*out1+(-8.884430e+01)*out2+(-4.329238e+00)*out3+(-1.640589e+00)*out4;
//     rin7 = 5.673327e+00+(-4.684875e+00)*out1+(1.588438e+01)*out2+(-1.594008e-01)*out3+(3.950433e-01)*out4;
//     rin8 = 2.608943e+00+(5.446163e-01)*out1+(1.185540e-01)*out2+(-2.052808e-01)*out3+(-3.912658e-02)*out4;
    
//     out1 = sigmoid(rin1);
//     out2 = sigmoid(rin2);
//     out3 = sigmoid(rin3);
//     out4 = sigmoid(rin4);
//     out5 = sigmoid(rin5);
//     out6 = sigmoid(rin6);
//     out7 = sigmoid(rin7);
//     out8 = sigmoid(rin8);

//     return  -5.099165e+00+(4.095054e+01)*out1+(-3.124696e+01)*out2+(4.999425e+00)*out3+(-5.421374e+01)*out4
//       +(-4.399975e+00)*out5+(3.546380e+01)*out6+(3.456229e+00)*out7+(-9.328331e+00)*out8;

//   calo only opimization
    rin1 = -5.360177e-01+(-2.453901e+00)*out1+(2.854906e+00)*out2+(-6.467948e+00)*out3+(-7.058127e-01)*out4;
    rin2 = 7.384829e+00+(2.474857e+00)*out1+(5.685580e+01)*out2+(-2.697409e+01)*out3+(-1.310651e+00)*out4;
    rin3 = 1.260522e+01+(5.548471e+00)*out1+(3.421428e+01)*out2+(-2.710610e+01)*out3+(-1.596936e+00)*out4;
    rin4 = -8.975015e-01+(1.391096e-02)*out1+(-3.744968e+01)*out2+(-8.128343e-01)*out3+(-1.744908e-01)*out4;
    rin5 = 1.846038e+01+(-1.209350e+01)*out1+(9.755953e+00)*out2+(-7.330227e+00)*out3+(7.784411e-01)*out4;
    rin6 = -7.261182e-01+(-6.468563e-01)*out1+(-6.709912e+01)*out2+(1.112497e+01)*out3+(1.229036e-01)*out4;
    rin7 = 2.790842e+00+(9.816161e-01)*out1+(9.051273e+01)*out2+(-2.326896e+01)*out3+(-6.385700e-01)*out4;
    rin8 = -7.101524e-01+(-2.439599e+00)*out1+(3.097952e+00)*out2+(-6.685354e+00)*out3+(-8.163578e-01)*out4;

    out1 = sigmoid(rin1);
    out2 = sigmoid(rin2);
    out3 = sigmoid(rin3);
    out4 = sigmoid(rin4);
    out5 = sigmoid(rin5);
    out6 = sigmoid(rin6);
    out7 = sigmoid(rin7);
    out8 = sigmoid(rin8);

    return  9.330931e-01+(7.433985e+00)*out1+(-6.739972e+01)*out2+(4.588412e+01)*out3+(-1.733942e+01)*out4
      +(-2.162184e+01)*out5+(2.300186e+01)*out6+(4.184828e+01)*out7+(-6.916628e+00)*out8;

  }else if(i==2){

//     calo+cps weight
//     rin1 = 2.725636e+00+(-1.568131e+00)*out1+(-4.500310e-01)*out2+(3.231739e+00)*out3+(3.339067e-01)*out4;
//     rin2 = 4.560314e+00+(3.125498e+00)*out1+(1.838387e+01)*out2+(1.828984e-01)*out3+(1.049720e+00)*out4;
//     rin3 = 6.348769e+00+(3.839253e+00)*out1+(1.547823e+00)*out2+(5.586311e+00)*out3+(1.305110e+01)*out4;
//     rin4 = 1.412610e+00+(6.054317e-02)*out1+(-2.003106e+01)*out2+(-4.482923e+01)*out3+(-2.678894e+01)*out4;
//     rin5 = -2.469605e+01+(-1.757508e+01)*out1+(6.385236e+00)*out2+(7.487655e+00)*out3+(3.979059e-01)*out4;
//     rin6 = -1.385140e+00+(2.867984e+00)*out1+(1.323049e+00)*out2+(-7.356154e+00)*out3+(-1.006133e+00)*out4;
//     rin7 = -1.807696e+01+(1.267538e+01)*out1+(5.846549e+00)*out2+(-9.746763e-01)*out3+(-5.681171e-01)*out4;
//     rin8 = -2.630656e-01+(1.674820e-01)*out1+(2.203393e+01)*out2+(3.936907e+01)*out3+(4.342471e+01)*out4;
    
//     out1 = sigmoid(rin1);
//     out2 = sigmoid(rin2);
//     out3 = sigmoid(rin3);
//     out4 = sigmoid(rin4);
//     out5 = sigmoid(rin5);
//     out6 = sigmoid(rin6);
//     out7 = sigmoid(rin7);
//     out8 = sigmoid(rin8);

//     return 9.547780e+00+(1.780360e+01)*out1+(-8.031423e+00)*out2+(1.379051e+01)*out3+(-1.253791e+01)*out4
//       +(-3.783483e+00)*out5+(3.619195e+00)*out6+(6.567922e+00)*out7+(-3.273851e+01)*out8;

//   calo only opimization
    rin1 = 5.104730e+00+(-3.400891e+00)*out1+(1.359165e+01)*out2+(1.478298e+01)*out3+(6.799255e-01)*out4;
    rin2 = 2.877062e+01+(2.054329e+01)*out1+(-6.851477e+00)*out2+(-1.117992e+01)*out3+(-3.407728e-01)*out4;
    rin3 = 6.081908e+00+(3.362017e+00)*out1+(1.686553e+01)*out2+(-1.341661e+00)*out3+(1.224359e-01)*out4;
    rin4 = -3.044222e+00+(1.867185e+00)*out1+(-5.561286e+00)*out2+(1.181619e+00)*out3+(3.745009e+00)*out4;
    rin5 = -2.355149e+00+(-3.571762e-01)*out1+(-1.217368e+00)*out2+(4.735327e+00)*out3+(9.456027e-01)*out4;
    rin6 = -1.327034e+00+(-1.444217e+00)*out1+(4.861928e+00)*out2+(4.191449e+00)*out3+(-2.363513e+00)*out4;
    rin7 = 2.063106e+01+(-1.491199e+01)*out1+(-5.377828e+00)*out2+(2.946555e+00)*out3+(4.551947e-01)*out4;
    rin8 = 1.918961e+00+(-6.971245e-01)*out1+(2.438949e+00)*out2+(-2.175864e+00)*out3+(-1.589770e+00)*out4;

    out1 = sigmoid(rin1);
    out2 = sigmoid(rin2);
    out3 = sigmoid(rin3);
    out4 = sigmoid(rin4);
    out5 = sigmoid(rin5);
    out6 = sigmoid(rin6);
    out7 = sigmoid(rin7);
    out8 = sigmoid(rin8);

    return -1.343235e+00+(9.939797e+00)*out1+(4.240392e+00)*out2+(-1.292941e+01)*out3+(2.004552e+00)*out4
      +(3.644026e+00)*out5+(-7.628740e-01)*out6+(-5.744442e+00)*out7+(6.484482e+00)*out8;


  }else if(i==3){

//     calo+cps weight
//     rin1 = 1.309300e+01+(9.397373e+00)*out1+(1.878087e+01)*out2+(3.076998e+00)*out3+(4.468628e-01)*out4;
//     rin2 = 2.267605e+00+(1.811185e+00)*out1+(-3.570715e+00)*out2+(8.918338e-01)*out3+(1.697065e-01)*out4;
//     rin3 = 1.461070e+01+(-1.274931e+01)*out1+(-1.529288e+00)*out2+(1.027187e+01)*out3+(2.249881e+00)*out4;
//     rin4 = -1.075816e+01+(7.422731e+00)*out1+(-1.871771e+01)*out2+(-5.157375e+00)*out3+(-5.079110e-01)*out4;
//     rin5 = 1.438973e+01+(1.246405e+01)*out1+(5.893971e+00)*out2+(8.753778e+00)*out3+(1.301357e+00)*out4;
//     rin6 = 1.352880e+00+(1.180933e+00)*out1+(-3.285885e-01)*out2+(7.042210e-01)*out3+(7.823139e-02)*out4;
//     rin7 = -1.695699e+01+(-7.238909e+00)*out1+(9.000217e+00)*out2+(9.360111e+00)*out3+(1.483872e+00)*out4;
//     rin8 = -2.703729e+00+(6.909974e-01)*out1 +(2.623371e+00)*out2+(4.555530e-01)*out3+(2.601189e-02)*out4;

//     out1 = sigmoid(rin1);
//     out2 = sigmoid(rin2);
//     out3 = sigmoid(rin3);
//     out4 = sigmoid(rin4);
//     out5 = sigmoid(rin5);
//     out6 = sigmoid(rin6);
//     out7 = sigmoid(rin7);
//     out8 = sigmoid(rin8);

//     return 1.241866e+01+(-1.442610e+01)*out1+(2.708457e+00)*out2+(-5.149567e+00)*out3+(-1.257937e+01)*out4
//       +(7.936903e+00)*out5+(-4.372071e+00)*out6+(-4.341796e+00)*out7+(4.042630e+00)*out8;

//   calo only opimization
    rin1 = 6.687513e+00+(4.742553e+00)*out1+(-1.921852e-01)*out2+(8.735341e+00)*out3+(5.669459e+00)*out4;
    rin2 = 2.060688e+00+(-8.148183e-01)*out1+(-1.166612e+00)*out2+(-2.448084e-01)*out3+(-7.259993e-02)*out4;
    rin3 = 1.454462e+01+(-1.189735e+01)*out1+(-1.613441e+00)*out2+(8.752237e+00)*out3+(1.699585e+00)*out4;
    rin4 = -9.160851e+00+(6.286569e+00)*out1+(-1.462017e+01)*out2+(-5.548660e+00)*out3+(-5.991488e-01)*out4;
    rin5 = 3.102560e+01+(2.036634e+01)*out1+(-1.164477e+01)*out2+(-2.456155e+00)*out3+(-3.053949e-01)*out4;
    rin6 = 9.800562e-02+(1.507541e+00)*out1+(-4.630475e+00)*out2+(1.796990e+00)*out3+(-1.994557e-01)*out4;
    rin7 = 1.752046e-02+(-2.086730e-02)*out1+(2.491200e+00)*out2+(-3.771323e-01)*out3+(1.263022e-01)*out4;
    rin8 = -1.898897e+01+(-1.171588e+01)*out1+(5.774021e+00)*out2+(-2.175084e+01)*out3+(-2.542364e+00)*out4;

    out1 = sigmoid(rin1);
    out2 = sigmoid(rin2);
    out3 = sigmoid(rin3);
    out4 = sigmoid(rin4);
    out5 = sigmoid(rin5);
    out6 = sigmoid(rin6);
    out7 = sigmoid(rin7);
    out8 = sigmoid(rin8);

    return 7.352498e+00+(2.893817e+00)*out1+(-5.417335e+00)*out2+(-5.839607e+00)*out3+(-1.426468e+01)*out4
      +(3.781558e+00)*out5+(-1.993675e+00)*out6+(-4.615290e+00)*out7+(2.887536e+01)*out8;

  }else if(i==4){

//     calo+cps weight
//     rin1 = 1.739504e+00+(-1.359310e+00)*out1+(4.871366e-01)*out2+(-3.675268e+00)*out3+(9.390654e+00)*out4;
//     rin2 = 1.101727e+00+(-3.695814e+00)*out1+(-3.519250e+00)*out2+(1.012465e+01)*out3+(1.632931e+01)*out4;
//     rin3 = 5.529250e+00+(-1.447094e+00)*out1+(-2.710102e+00)*out2+(-6.364409e+00)*out3+(-4.101287e-01)*out4;
//     rin4 = 2.466821e+00+(5.632797e+00)*out1+(1.077799e+00)*out2+(-1.082192e+01)*out3+(7.127756e+00)*out4;
//     rin5 = 5.977946e+00+(1.193340e+00)*out1+(-3.412404e+00)*out2+(-6.202793e+00)*out3+(-7.271314e-01)*out4;
//     rin6 = 1.311999e+00+(-4.860522e-01)*out1+(-1.350748e+00)*out2+(-5.562873e+00)*out3+(1.142923e+01)*out4;
//     rin7 = -1.012606e+00+(-2.307605e-01)*out1+(9.952796e-01)*out2+(5.846318e+00)*out3+(-1.269713e+01)*out4;
//     rin8 = -5.015484e-01+(-9.930616e-01)*out1+(-2.611544e+00)*out2+(6.926847e+00)*out3+(-1.343878e+01)*out4;

//     out1 = sigmoid(rin1);
//     out2 = sigmoid(rin2);
//     out3 = sigmoid(rin3);
//     out4 = sigmoid(rin4);
//     out5 = sigmoid(rin5);
//     out6 = sigmoid(rin6);
//     out7 = sigmoid(rin7);
//     out8 = sigmoid(rin8);

//     return 5.408648e+00+(8.039118e+00)*out1+(-1.769201e+00)*out2+(-9.886955e+00)*out3+(3.607326e-01)*out4
//       +(1.094678e+01)*out5+(-1.309899e+01)*out6+(-1.241196e+01)*out7+(4.139146e+00)*out8;

//   calo only opimization
    rin1 = 5.512694e+00+(2.526613e+00)*out1+(7.916696e+00)*out2+(-2.725401e+01)*out3+(-1.080395e-01)*out4;
    rin2 = 2.469562e+00+(1.889429e+00)*out1+(3.185034e-01)*out2+(-7.544770e+00)*out3+(5.015396e+00)*out4;
    rin3 = 7.734093e+00+(-1.149609e+00)*out1+(-5.505675e+00)*out2+(-6.492834e+00)*out3+(-6.539473e-01)*out4;
    rin4 = -5.437393e+00+(2.204851e+00)*out1+(-5.148233e+00)*out2+(2.505908e+01)*out3+(-3.637297e-01)*out4;
    rin5 = 7.113125e+00+(1.226893e+00)*out1+(-4.843938e+00)*out2+(-5.999518e+00)*out3+(-6.492181e-01)*out4;
    rin6 = 1.896361e+00+(-6.758293e-01)*out1+(-3.629415e+00)*out2+(-3.283566e+00)*out3+(5.963047e+00)*out4;
    rin7 = -3.464513e+00+(2.527691e+00)*out1+(-6.837227e-01)*out2+(1.235217e+01)*out3+(-5.135345e+00)*out4;
    rin8 = 1.768892e+00+(5.876383e-01)*out1+(-2.999926e+00)*out2+(-3.505533e+00)*out3+(5.737184e+00)*out4;

    out1 = sigmoid(rin1);
    out2 = sigmoid(rin2);
    out3 = sigmoid(rin3);
    out4 = sigmoid(rin4);
    out5 = sigmoid(rin5);
    out6 = sigmoid(rin6);
    out7 = sigmoid(rin7);
    out8 = sigmoid(rin8);

    return 4.031643e+00+(2.153386e+00)*out1+(-5.510784e+00)*out2+(-1.852509e+01)*out3+(2.804891e+00)*out4
      +(1.635435e+01)*out5+(-8.254666e+00)*out6+(-3.871061e+00)*out7+(9.779259e+00)*out8;

  }else return 0.0;
}