#include "MessageService/MsgService.h"


#include "ZdetWeight.h"
#include "Parsers.h"
#include "FitBeam.h"
#include "FitData.h"
#include "FitPar.h"
#include "GenericFactory.h"

#include <cmath>

#include <boost/algorithm/string/split.hpp>
#include <boost/algorithm/string/classification.hpp>
#include <boost/algorithm/string/predicate.hpp>
#include <boost/algorithm/string.hpp>

REGISTER_OBJECT(SkzpWeight,ZdetWeight)

using namespace boost;
using namespace Parsers;
CVSID("");


ZdetWeight::ZdetWeight()
  :fName("ZdetWeight"),
   fAntiNuXsError(0.3)
{
  
  
}

ZdetWeight::~ZdetWeight()
{
  
}

Bool_t ZdetWeight::SetProperty(std::string prop) 
{
  
  string property=prop.substr(0,prop.find("="));
  string val=prop.substr(prop.find("=")+1);

  if(property=="Parameters") {
    vector<std::string> params;

    parse_val(val,params);
    std::vector<std::string>::iterator iter=params.begin();
    std::vector<std::string>::iterator iter_end=params.end();
    for( ; iter!=iter_end ; ++iter) {
      //cout<<*iter<<endl;
      std::vector<string> vals;
      split(vals, *iter, is_any_of("() "), token_compress_on);
      
      FitPar par(vals);
      fParameters.push_back(par);
    }
  }
  else if(property=="AntiNuXsError") parse_val(val,fAntiNuXsError);
  else if(property=="UseShwScale") parse_val(val,fUseShwScale);
  else {
    MSG("ZdetWeight",Msg::kWarning)<<"Cannot set property "<<property
				    <<" for ZdetWeight function "<<fName<<endl;
    return false;
  } 
   
  return true;
  
}


double ZdetWeight::GetWeight(FitData &data)
{
  
  double weight=1.;

  // Shower offset
  if(data.RecoShwEnergy()>0.) {
    data.RecoShwEnergy(data.RecoShwEnergy()+fPars[2]);
  }
  if( data.RecoShwEnergy() < 0.) data.RecoShwEnergy(0.);

  // Shower energy scale
  if(!fUseShwScale) {
    data.RecoShwEnergy(data.RecoShwEnergy()*(1-fPars[8]));
  }
  else {
    data.RecoShwEnergy(data.RecoShwEnergy()*(1-fPars[8]*data.ShwError()));
  }
  
  // Muon energy scale
  data.RecoTrkEnergy(data.RecoTrkEnergy()*(1.-fPars[9]));
  
   // NC scaling
  if(data.InteractionType()==0) {
    //  if(data.NuType()==14)       weight*=(1-fPars[3]);
    //     else if(data.NuType()==-14) weight*=(1-fPars[4]);
    //     else if(data.NuType()==13)  weight*=(1-fPars[5]);
    
    if      (data.Beam()>=5600&&data.Beam()<5700) weight*=(1- fPars[3]);
    else if (data.Beam()>=5500&&data.Beam()<5600) weight*=(1- fPars[4]);
    else if (data.Beam()>=5300&&data.Beam()<5400) weight*=(1- fPars[5]);
  }
  // Xs scaling
  else {
    if(data.InteractionProcess()==1001)  {
      weight*=(1+fPars[10]*data.GeneratorError());
    }
    else if(data.InteractionProcess()==1002) {
      weight*=(1+fPars[11]*data.GeneratorError());
    }
    else if(data.InteractionProcess()==1003) {
      weight*=(1+fPars[12]*data.GeneratorError());
    }
  }

  
  // Energy Scaling
  if( fabs(data.NuType())==14 || data.NuType()==55 || data.NuType()==56) {
    data.RecoEnergy( (data.RecoShwEnergy()+data.RecoTrkEnergy())*(1-fPars[0]) );
  }
  else if( fabs(data.NuType())==12 || data.NuType()==52 || data.NuType()==53) {
    data.RecoEnergy( (data.RecoShwEnergy()+data.RecoTrkEnergy())*(1-fPars[1]) );
  }
  else data.RecoEnergy( data.RecoShwEnergy()+data.RecoTrkEnergy());



  // Other 
  if( (data.NuType()==-14||data.NuType()==55) && data.TrueEnergy() < fPars[7] ) {
    double a=fPars[6];
    double xp=fPars[7];
    weight*=a+2.*(1.-a)*data.TrueEnergy()/xp+(a-1.)*
      data.TrueEnergy()*data.TrueEnergy()/(xp*xp);
  }

  MSG("ZdetWeight",Msg::kVerbose)
    <<"Det Weight="<<weight<<endl;
  return weight;
}

double ZdetWeight::GetPT() const
{
  double penalty=0;
  for (int j=0;j<4;j++) {
    double a=fPars[6];
    double xp=fPars[7];
    double true_enu=(2*j+1)*xp/8.;
    penalty+=pow((a+2.*(1.-a)*true_enu/xp+(a-1.)*
		 true_enu*true_enu/(xp*xp)-1.),2)/(fAntiNuXsError*fAntiNuXsError)/4.;
  }

  return penalty;

}