SKZPWeightCalculator.cxx

Go to the documentation of this file.

#include <cassert>
#include <map>
#include "DatabaseInterface/DbiConfigStream.h"
#include "MCReweight/Zbeam.h"
#include "MCReweight/Zfluk.h"
#include "MCNtuple/NtpMCFluxInfo.h"
#include "MCNtuple/NtpMCTruth.h"
#include "MCReweight/SKZPWeightCalculator.h"
#include "MessageService/MsgService.h"

#include "Conventions/BeamType.h"
#include "MCReweight/BeamSys.h"

#include "TMath.h"

using namespace std;

CVSID("$Id: SKZPWeightCalculator.cxx,v 1.28 2010/02/06 11:12:03 med Exp $");

const double SKZPWeightCalculator::XP=25.;

SKZPWeightCalculator::SKZPWeightCalculator(Registry *stdconfig):
  WeightCalculator(stdconfig),
  MAXBSPARS(20),
  MAXHPPARS(20),
  MAXDETPARS(20),
  read_pars_from_db(false),
  cfgname(""),
  cfgnum(kUnknown),
  sampsel(kNotSpecified),
  runfrac()
{
  vtsRunIEnd = new VldTimeStamp(2006,3,1,0,0,0);
  vtsRunIIEnd = new VldTimeStamp(2007,8,1,0,0,0);
  vtsRunIIIEnd = new VldTimeStamp(2009,8,1,0,0,0);
  zbeam = new Zbeam();
  zfluk = new Zfluk();
  nwc = 0; // hold off until needed // new NeugenWeightCalculator();
  mcr = &MCReweight::Instance();
  Config();
}

SKZPWeightCalculator::SKZPWeightCalculator(string cn, bool readfromdb, bool beamonly):
  WeightCalculator(0),
  MAXBSPARS(20),
  MAXHPPARS(20),
  MAXDETPARS(20),
  read_pars_from_db(readfromdb),
  cfgname(cn),
  cfgnum(kUnknown),
  sampsel(kNotSpecified),
  runfrac()
{
  fBeamOnly = beamonly;

  vtsRunIEnd = new VldTimeStamp(2006,4,1,0,0,0);
  vtsRunIIEnd = new VldTimeStamp(2007,8,1,0,0,0);
  vtsRunIIIEnd = new VldTimeStamp(2009,8,1,0,0,0);
  zbeam = new Zbeam();
  zfluk = new Zfluk();
  nwc = 0; // hold off until needed // new NeugenWeightCalculator();
  mcr = &MCReweight::Instance();
  Config();
}

SKZPWeightCalculator::~SKZPWeightCalculator()
{

  delete vtsRunIEnd;
  vtsRunIEnd=0;
  delete vtsRunIIEnd;
  vtsRunIIEnd=0;
  delete vtsRunIIIEnd;
  vtsRunIIIEnd=0;
  delete zbeam;
  zbeam=0;
  delete zfluk;
  zfluk=0;
  delete nwc;
  nwc=0;

}

void SKZPWeightCalculator::Config()
{
  if(cfgname==""){
    cfgname="PRL";
    MSG("SKZPWeightCalculator",Msg::kInfo)<<"You did not specify a "
                                          <<"SKZPWeight configuration."
                                          <<" Setting ConfigName to :"
                                          <<cfgname<<endl;
    cfgnum=kPRL;
  }
  SetConfigNum(cfgname);


  if(!read_pars_from_db){
    MSG("SKZPWeightCalculator",Msg::kInfo)<<"You did not select to read in "
                                          <<"SKZP Fit parameters from the DB."
                                          <<" You must specify your "
                                          <<"own parameters"<<endl;
    return;
  }

  //retrieve appropriate constants from DB
  Registry reg;
  
  //  cout<<"Trying to get parameters from DB"<<endl;
  DbiConfigStream stream("SKZPWeights",cfgname.c_str());
  if(stream.IsEmpty()){
    MSG("SKZPWeightCalculator",Msg::kWarning)
      << "Could not find SKZP configuration: "
      << cfgname
      << " in database.  " << endl
      << "If you want to run without database access"
      << " construct the weight calculator"
      << " with false as second arg." << endl;
    MSG("SKZPWeightCalculator",Msg::kFatal) << "Now assert()" << endl;
    assert(false);
  }
  //  cout<<"made dbiconfigstream"<<endl;
  stream >> &reg;

  //  cout<<"Read registry from db, and it contains:"<<endl;
  //  reg.PrettyPrint(cout);



  //load information from registry to private data members
  double tmpd;
  for(int i=0;i<MAXBSPARS;i++){
    if(reg.Get(Form("bs_par_%d",i),tmpd)){
      bs_pars.push_back(tmpd);
    }
  }

  for(int i=0;i<MAXHPPARS;i++){
    if(reg.Get(Form("hp_par_%d",i),tmpd)){
      hp_pars.push_back(tmpd);
    }
  }

  for(int i=0;i<MAXDETPARS;i++){
    if(reg.Get(Form("det_par_%d",i),tmpd)){
      det_pars.push_back(tmpd);
    }
  }

  if(cfgnum==kPiMinus){
    zfluk->SetParameters(hp_pars);
  }
  if(cfgnum==kDetXs || cfgnum==kDogwood1_Daikon07 || cfgnum==kDogwood1_Daikon07_v2){
    zfluk->SetParameters(hp_pars);

    if(!fBeamOnly) {
      nwc = new NeugenWeightCalculator();
      nwc->SetOptimization(true);
      Registry *nwc_reg=new Registry;
      nwc_reg->UnLockValues();
      nwc_reg->UnLockKeys();
      nwc_reg->Clear();
      
      nwc_reg->Set("neugen:config_name","MODBYRS");
      nwc_reg->Set("neugen:config_no",4);
      nwc_reg->Set("neugen:use_scale_factors",1);
      nwc_reg->Set("neugen:ma_qe",1.1);
      nwc_reg->Set("neugen:ma_res",1.1); 
      
      // This is specific to MODBYRS v4 where we weight them 
      // by their 1 sigma error as given in docdb-2989
      
      double kno1_val=0.1; //for ri12 and ri42
      double kno2_val=0.3; //for ri22 and ri32
      double kno3_val=1.0; //for rij3 
    
      double kno1_err=0.1; //for rij2 
      double kno2_err=0.2; //for rij3
      
      double kno_shift=1;
      
      nwc_reg->Set("neugen:kno_r112",
                 (kno1_val+kno_shift*kno1_err)/(kno1_val));
      nwc_reg->Set("neugen:kno_r142",
                   (kno1_val+kno_shift*kno1_err)/(kno1_val));
      nwc_reg->Set("neugen:kno_r122",
                   (kno2_val+kno_shift*kno1_err)/(kno2_val));
      nwc_reg->Set("neugen:kno_r132",
                   (kno2_val+kno_shift*kno1_err)/(kno2_val));
      nwc_reg->Set("neugen:kno_r113",
                   (kno3_val+kno_shift*kno2_err)/(kno3_val));
      nwc_reg->Set("neugen:kno_r123",
                   (kno3_val+kno_shift*kno2_err)/(kno3_val));
      nwc_reg->Set("neugen:kno_r133",
                   (kno3_val+kno_shift*kno2_err)/(kno3_val));
      nwc_reg->Set("neugen:kno_r143",
                   (kno3_val+kno_shift*kno2_err)/(kno3_val));
      
      nwc_reg->LockValues();
      nwc_reg->LockKeys();
      nwc->SetReweightConfig(nwc_reg);
      mcr->AddWeightCalculator(nwc);
    }
  }
  string tmpcfgname=cfgname;
  if(tmpcfgname.find("PRL_Syst")!=string::npos){
    tmpcfgname="PRL";
  }
  zbeam->SetReweightConfig(tmpcfgname);
}

void SKZPWeightCalculator::SetSampleSelection(std::string ss)
{
  if(ss=="numu"||ss=="NuMu"||ss=="CC"){
    sampsel=kNuMuCC;
  }
  else if(ss=="anumu"||ss=="antinumu"||ss=="numubar"||ss=="AntiNuMu"||ss=="NuMuBar"){
    sampsel=kAntiNuMu;
  }
  else if(ss=="nue"||ss=="Nue"){
    sampsel=kNue;
  }

}

double SKZPWeightCalculator::GetWeight(Registry *eventinfo)
{
  //don't use me I'm slow.
  //the GetWeight function this version calls has been tested,
  //but reading the event info from the registry has not 
  //been tested (it should work, I just haven't tried it yet)
  bool gotallinfo=true;
  int det,ibeam,tptype,ccnc,inu;
  double pt, pz, true_enu, reco_emu, reco_eshw;
  double new_reco_emu, new_reco_eshw, beamweight, detweight;
  
  if(!eventinfo->Get("evt:det",det)) gotallinfo=false;
  if(!eventinfo->Get("evt:ibeam",ibeam)) gotallinfo=false;
  if(!eventinfo->Get("evt:tptype",tptype)) gotallinfo=false;
  if(!eventinfo->Get("evt:ccnc",ccnc)) gotallinfo=false;
  if(!eventinfo->Get("evt:inu",inu)) gotallinfo=false;
  if(!eventinfo->Get("evt:pt",pt)) gotallinfo=false;
  if(!eventinfo->Get("evt:pz",pz)) gotallinfo=false; 
  if(!eventinfo->Get("evt:true_enu",true_enu)) gotallinfo=false;  
  if(!eventinfo->Get("evt:reco_emu",reco_emu)) gotallinfo=false;  
  if(!eventinfo->Get("evt:reco_eshw",reco_eshw)) gotallinfo=false;

  if(!gotallinfo){
    MSG("SKZPWeightCalculator",Msg::kError)<<"Could not retrieve all necessary"
                                           <<" event information from "
                                           <<"eventinfo registry, returning 1."
                                           <<endl;
    return 1.;
  }

  GetWeight(det,ibeam,tptype,pt,pz,ccnc,
            true_enu,inu,reco_emu,reco_eshw,
            new_reco_emu,new_reco_eshw,beamweight,detweight);

  //put new reco energies and weights into registry;
  eventinfo->UnLockKeys();
  eventinfo->Set("evt:new_reco_emu",new_reco_emu);
  eventinfo->Set("evt:new_reco_eshw",new_reco_eshw);
  eventinfo->Set("evt:beamweight",beamweight);
  eventinfo->Set("evt:detweight",detweight);
  eventinfo->LockKeys();

  // this function returns only the beam weight, as that is 
  //what most people are encouraged to use, the detweight is
  //returned through the registry
  return beamweight;
}

double SKZPWeightCalculator::GetWeight(const NtpMCTruth *mc, int det, int Ibeam,
                                       double reco_emu, double reco_eshw, 
                                       double &new_reco_emu, double &new_reco_eshw,
                                       double &beamweight, double &detweight,
                                       RunPeriod_t rp,
                                       double true_eshw,MCEventInfo *ei)
{


  //retrieve truth information on the neutrino interaction from
  //the NtpMCTruth object
  int ccnc = mc->iaction;
  double true_enu = 1.*mc->p4neu[3];
  int inu = mc->inu;


  //pass the truth info and the fluxinfo object to the next
  //form of GetWeight
  return GetWeight(&(mc->flux),det,Ibeam,
                   ccnc,true_enu,inu,
                   reco_emu,reco_eshw,
                   new_reco_emu,new_reco_eshw,
                   beamweight, detweight,rp,true_eshw,ei);

}

double SKZPWeightCalculator::GetWeight(const NtpMCFluxInfo *fi, int det, int Ibeam,
                                       int ccnc, double true_enu,int inu,
                                       double reco_emu, double reco_eshw, 
                                       double &new_reco_emu, double &new_reco_eshw,
                                       double &beamweight, double &detweight,
                                       RunPeriod_t rp,
                                       double true_eshw,MCEventInfo *ei)
{
  int parenttype=-1;
  double parentpz=0.;
  double parentpx=0.;
  double parentpy=0.;

  //extract the flux weighting variables from the fluxinfo object
  //due to bug in v18 flux files, the tp* variables
  //(i.e. info pertaining to the particle created in the target)
  //are the same as the pp* variables
  //(i.e. the info pertaining to the particle leaving the target)
  //in v19 flux and beyond, this has been fixed.
  //In the future, we will probably want to weight based on the tp variables,
  //but to compare with what was done before, I want the option to 
  //use the pp* variables, which is why this if statement is here.

  if(cfgname=="v19fit"){  //this isn't a real config yet
    parenttype = fi->ptype;
    parentpz = 1.*fi->pppz;
    parentpy = 1.*fi->ppdydz*parentpz;
    parentpx = 1.*fi->ppdxdz*parentpz;
  }
  else{
    parenttype = fi->tptype;
    parentpz = 1.*fi->tpz;
    parentpy = 1.*fi->tpy;
    parentpx = 1.*fi->tpx;
  }

  double pt = sqrt(parentpy*parentpy+parentpx*parentpx);

  //pass the variables on to the next version of GetWeight
  return GetWeight(det,Ibeam,parenttype,pt,parentpz,ccnc,true_enu,inu,
                   reco_emu,reco_eshw,new_reco_emu,new_reco_eshw,beamweight,detweight,
                   rp,true_eshw,ei);

}
  
double SKZPWeightCalculator::GetWeight(int det, int Ibeam,
                                       int tptype, double pt, double pz,    
                                       int ccnc, double true_enu,int inu,
                                       double reco_emu, double reco_eshw,
                                       double &new_reco_emu, double &new_reco_eshw,
                                       double &beamweight, double &detweight,
                                       RunPeriod_t rp,
                                       double true_eshw,MCEventInfo *ei)
{
  
  //do the beam weighting
  beamweight=GetBeamWeight(det,Ibeam,tptype,pt,pz,true_enu,inu,rp);
  //do the detector weighting
  detweight=GetDetWeight(ccnc,true_enu,inu,
                         reco_emu,reco_eshw,new_reco_emu,new_reco_eshw,true_eshw,ei);

  //return just the beam weight, as this is what 
  //most people will want
  //the value of the detweight is passed by reference
  //do not use the new values of new_reco_emu and new_reco_eshw
  //in the case of beam weighting only!!!!!
  return beamweight;

}

double SKZPWeightCalculator::GetBeamWeight(int det, int Ibeam,
                                           int tptype, double pt,
                                           double pz, double true_enu, int inu, 
                                           RunPeriod_t rp)
{
  //do the flux reweighting


  double beamsysw=1.;
  double hadpw=1.;
  if(cfgnum==kPRL){
    //to let zbeam know I want the PRL version of 
    //of it's weights, I'm passing the beam systematic
    //parameters parameters one by one

    //reweight for spot size in me and he
    //now done in RunPeriodWeight
    //    beamsysw*=zbeam->GetWeight(det,Ibeam,6,0,true_enu);   
    //horn 1 offset
    beamsysw*=zbeam->GetWeight(det,Ibeam,1,bs_pars[0],true_enu);
    //baffle scraping
    beamsysw*=zbeam->GetWeight(det,Ibeam,2,bs_pars[1],true_enu);
    //pot
    beamsysw*=zbeam->GetWeight(det,Ibeam,3,bs_pars[2],true_enu);
    //horn current miscalibration
    beamsysw*=zbeam->GetWeight(det,Ibeam,4,bs_pars[3],true_enu);
    //horn current distribution
    beamsysw*=zbeam->GetWeight(det,Ibeam,5,bs_pars[4],true_enu);

    //now change pt and pz
    double hp[7]={hp_pars[0],hp_pars[1],hp_pars[2],hp_pars[3],
                  hp_pars[4],hp_pars[5],hp_pars[6]};
    //to let zfluk know I want the PRL version of it's weights
    //I'm passing the hadron production parameters as an array
    hadpw = zfluk->GetWeightPRL(tptype,pt,pz,hp);
  }  
  else if(cfgnum==kBoston){ 
    //Boston version of zbeam takes a vector of beam systematic parameters
     //first 5 are for numu's
     //second 5 are for antinumu's
     vector<double> tmp_bs_pars;
     if(inu>0){
        for(int i=0;i<5;i++){
           tmp_bs_pars.push_back(bs_pars[i]);
        }
     }
     else{
        for(int i=0;i<5;i++){
           tmp_bs_pars.push_back(bs_pars[i+5]);
        }
     }
     beamsysw=zbeam->GetWeight(det,Ibeam,tmp_bs_pars,inu,true_enu);     

     //reweight for spot size in me and he
     //now done in RunPeriodWeight
     //     beamsysw*=zbeam->GetWeight(det,Ibeam,6,0,true_enu);   
     
     //updated version of zbeam takes a vector of hadron production parameters
     hadpw=zfluk->GetWeightBoston(tptype,pt,pz,hp_pars);
  }
  else if(cfgnum==kPiMinus){
    //updated version of zbeam takes a vector of beam systematic parameters
    beamsysw=zbeam->GetWeight(det,Ibeam,bs_pars,inu,true_enu);
    //reweight for spot size in me and he
    //now done in RunPeriodWeight
    //    beamsysw*=zbeam->GetWeight(det,Ibeam,6,0,true_enu);   

    //pi minus version of fit has parameters set one time at the beginning,
    //the parameters are not reset each time in getweight
    hadpw = zfluk->GetWeight(tptype,pt,pz);
  }
  else if(cfgnum==kDetXs || cfgnum==kDogwood1_Daikon07){

    // Hadron Production weight
    hadpw = zfluk->GetWeight(tptype,pt,pz);
 
  }
  else if(cfgnum==kDogwood1_Daikon07_v2){

    Detector::Detector_t detector;
    if(det==1) detector=Detector::kNear;
    else if(det==2) detector=Detector::kFar;
    else if(det==3) detector=Detector::kUnknown;
    
    Zbeam::ZbeamData_t zdata;
    zdata.ntype    = inu;
    zdata.true_enu = true_enu;
    zdata.detector = detector;
    zdata.beam     =  BeamType::FromZarko(Ibeam);
    
    // Focusing alignment effects
    beamsysw *= zbeam->GetWeight(zdata,BeamSys::kHorn1Offset   ,bs_pars[0]);
    beamsysw *= zbeam->GetWeight(zdata,BeamSys::kBaffleScraping,bs_pars[1]);
    beamsysw *= zbeam->GetWeight(zdata,BeamSys::kPOT           ,bs_pars[2]);
    beamsysw *= zbeam->GetWeight(zdata,BeamSys::kHornIMiscal   ,bs_pars[3]);
    beamsysw *= zbeam->GetWeight(zdata,BeamSys::kHornIDist     ,bs_pars[4]);
    
    // Hadron production weight
    hadpw = zfluk->GetWeight(tptype,pt,pz);

  }
  else{ 
    //updated version of zbeam takes a vector of beam systematic parameters
    beamsysw=zbeam->GetWeight(det,Ibeam,bs_pars,inu,true_enu);
    //updated version of zbeam takes a vector of hadron production parameters
    hadpw=zfluk->GetWeight(tptype,pt,pz);
  }

  double rpw = GetRunPeriodWeight(rp,det,Ibeam,inu,true_enu);
  //  cout<<"run period "<<rp<<" weight "<<rpw<<endl;
  //return the product of beam systematic weights and had. prod. weights
  //  cout<<"inu "<<inu<<"hadpw "<<hadpw<<" beamsysw "<<beamsysw<<endl;
  return hadpw*beamsysw*rpw;
}

double SKZPWeightCalculator::GetDetWeight(int ccnc, double true_enu, int inu,
                                          double reco_emu, double reco_eshw,
                                          double &new_reco_emu, double &new_reco_eshw,
                                           double true_eshw,MCEventInfo *ei)
{

  //computes the effect of the changing the detector parameters in the skzp fit
  //this function changes the reconstructed values of muon and shower energy,
  //so the new values are passed back to the calling function by reference as
  //part of the function parameter list
  double dw = 1.;
  new_reco_emu=reco_emu;
  new_reco_eshw=reco_eshw;


  if(cfgnum==kPRL){
    //right now this is the only version of detector weighting,
    //but this may evolve in the future

    //prl version of weights weights the NC component up or down by 1-fitparam
    if(ccnc==0){
      dw = 1.-det_pars[2];
    }
    
    //prl version of skzp changes the shower energy by an offset,
    //then changes the total neutrino energy (eshw+emu) by 
    //a multiplicative constant that is 1-fit parameter
    if(reco_eshw>0){
      new_reco_eshw = reco_eshw+det_pars[1];
    }
    if(new_reco_eshw<0){
      new_reco_eshw=0.;
    }
    
    new_reco_eshw*=(1-det_pars[0]);
    new_reco_emu=reco_emu*(1-det_pars[0]);
  }
  else if(cfgnum==kBoston){
    //boston version of weights weights the numu's and anumus separately
    if(ccnc==0){
       if(inu>0){
          dw = 1.-det_pars[2];
       }
       else{
          dw = 1.-det_pars[5];
       }
    }
    
    //boston version of skzp changes the shower energy by an offset,
    //then changes the total neutrino energy (eshw+emu) by 
    //a multiplicative constant that is 1-fit parameter
    if(reco_eshw>0){
       if(inu>0){
          new_reco_eshw = reco_eshw+det_pars[1];
       }
       else{
          new_reco_eshw = reco_eshw+det_pars[4];
       }
    }
    if(new_reco_eshw<0){
       new_reco_eshw=0.;
    }
    
    if(inu>0){
       new_reco_eshw*=(1-det_pars[0]);
       new_reco_emu=reco_emu*(1-det_pars[0]);
    }
    else{
       new_reco_eshw*=(1-det_pars[3]);
       new_reco_emu=reco_emu*(1-det_pars[3]);
    }
  }
  else if(cfgnum==kPiMinus){
    //Piminus version of weights weights the numu's and anumus separately
    if(ccnc==0){
      if(sampsel==kNuMuCC){
        //      cout<<"Found an NC, and I'm using the CC parameter"<<endl;
        dw = 1-det_pars[2];
      }
      else if(sampsel==kAntiNuMu){
        //      cout<<"Found an NC, and I'm using the Anti parameter"<<endl;
        dw = 1-det_pars[3];
      }
      else{//if a sample isn't specified, 
           //do it like I think you should.
       if(inu>0){
         dw = 1.-det_pars[2];
       }
       else{
         dw = 1.-det_pars[3];
       }
      }
    }
    
    //Piminus version of skzp changes the shower energy by an offset,
    //then changes the total neutrino energy (eshw+emu) by 
    //a multiplicative constant that is 1-fit parameter
    //it also treats nu's and anu's in the same way
    new_reco_eshw=reco_eshw;
    if(reco_eshw>0){
      new_reco_eshw = reco_eshw+det_pars[1];
    }
    if(new_reco_eshw<0){
      new_reco_eshw=0.;
    }
    
    if(inu==14||inu==-14){
      new_reco_eshw*=(1-det_pars[0]);
      new_reco_emu=reco_emu*(1-det_pars[0]);
    }
    else if(inu==12||inu==-12){
      new_reco_eshw*=(1-0.);
      new_reco_emu=reco_emu*(1-0.);
    }
    else{
      new_reco_eshw*=(1-0.);
      new_reco_emu=reco_emu*(1-0.);
    }

    //Piminus version of the fit has a scale for numubar xsec
    if(inu==-14){
      double xsw=1.;
      if(true_enu<XP){
        xsw = det_pars[4]+2.*(1.-det_pars[4])*true_enu/XP
          +(det_pars[4]-1.)*true_enu*true_enu/(XP*XP);
      }
      dw*=xsw;
    }
  }
  else if(cfgnum==kDetXs || cfgnum==kDogwood1_Daikon07){

    // You must give an MCEventInfo object to get the
    // correct detector weight
    if(fBeamOnly||ei==0) return dw;

    // Shower offset
    if(reco_eshw>0.) {
      new_reco_eshw = reco_eshw+det_pars[0];
    }
    if(new_reco_eshw<0){
      new_reco_eshw=0.;
    }

    // Shower energy scale from functional form
    double x0=10;
    double x = TMath::Min((double)true_eshw,x0);
    double err_a = 1.04/1.0 - 1.0;
    double err_b = (1.1+2.*(1-1.1)*(x/x0)+(1.1-1)*(x/x0)*(x/x0))/(1.0) - 1.0;
    double err_c = 0.057;
    double shw_scale=sqrt(pow(err_a,2)+pow(err_b,2)+pow(err_c,2));
    
    new_reco_eshw*=(1-det_pars[4]*shw_scale);

    // Scale muon energy 
    new_reco_emu=reco_emu*(1.-det_pars[5]);
    
    // NC parameters
    if(ei->iaction==0){
      if(sampsel==kNuMuCC)         dw*=(1-det_pars[1]);
      else if(sampsel==kAntiNuMu)  dw*=(1-det_pars[2]);
      else{
        if(inu==14)  dw = 1.-det_pars[1];
        if(inu==-14) dw = 1.-det_pars[2];
      }
    }
    else {
      // Cross Section Parameters
      if(ei!=0) {
        double gw=1;
        if(ei->iresonance!=1005) gw=mcr->ComputeWeight(ei)-1.0;
        
        if(ei->iresonance==1001) dw*=(1+det_pars[6]*gw);
        if(ei->iresonance==1002) dw*=(1+det_pars[7]*gw);
        //if(ei->iresonance==1003) dw*=(1+det_pars[8]*gw);
      }
    }
        
    //numubar xsec
    if(inu==-14){
      double xsw=1.;
      if(true_enu<XP){
        xsw = det_pars[3]+2.*(1.-det_pars[3])*true_enu/XP
          +(det_pars[3]-1.)*true_enu*true_enu/(XP*XP);
      }
      dw*=xsw;
    }
   
  }
  else if(cfgnum==kDogwood1_Daikon07_v2){

    // Similar to kDetXs but with parameters incremented
    // by 1 as new parameter 0 is the overall neutrino 
    // energy scale.

    // You must give an MCEventInfo object to get the
    // correct detector weight
    if(fBeamOnly||ei==0) return dw;

    // Shower offset
    if(reco_eshw>0.) {
      new_reco_eshw = reco_eshw+det_pars[1];
    }
    if(new_reco_eshw<0){
      new_reco_eshw=0.;
    }

    // Shower energy scale from functional form
    double x0=10;
    double x = TMath::Min((double)true_eshw,x0);
    double err_a = 1.04/1.0 - 1.0;
    double err_b = (1.1+2.*(1-1.1)*(x/x0)+(1.1-1)*(x/x0)*(x/x0))/(1.0) - 1.0;
    double err_c = 0.057;
    double shw_scale=sqrt(pow(err_a,2)+pow(err_b,2)+pow(err_c,2));
    
    new_reco_eshw*=(1-det_pars[5]*shw_scale);

    // Scale muon energy 
    new_reco_emu=reco_emu*(1.-det_pars[6]);
    
    // Now need to scale the overall neutrino energy BUT the current
    // weighting functions don't allow me to pass back a new value 
    // for the neutrino energy, just the muon and shower separately.
    // So I'm going to assume everyone will just add these two to
    // get the new neutrino energy and hence implement the neutrino
    // energy scale parameter by applying it to both the muon and
    // shower energy seperately:

    new_reco_eshw*=(1.-det_pars[0]);
    new_reco_emu*=(1.-det_pars[0]);

    // NC parameters
    if(ei->iaction==0){
      if(sampsel==kNuMuCC)         dw*=(1-det_pars[2]);
      else if(sampsel==kAntiNuMu)  dw*=(1-det_pars[3]);
      else{
        if(inu==14)  dw = 1.-det_pars[2];
        if(inu==-14) dw = 1.-det_pars[3];
      }
    }
    else {
      // Cross Section Parameters
      if(ei!=0) {
        double gw=1;
        if(ei->iresonance!=1005) gw=mcr->ComputeWeight(ei)-1.0;
        
        if(ei->iresonance==1001) dw*=(1+det_pars[7]*gw);
        if(ei->iresonance==1002) dw*=(1+det_pars[8]*gw);
      }
    }
        
    //numubar xsec
    if(inu==-14){
      double xsw=1.;
      if(true_enu<XP){
        xsw = det_pars[4]+2.*(1.-det_pars[4])*true_enu/XP
          +(det_pars[4]-1.)*true_enu*true_enu/(XP*XP);
      }
      dw*=xsw;
    }

  }
  else{
    //this does nothing for now, and 
    //should never be called,
    //but I want the structure
    MSG("SKZPWeightCalculator",Msg::kError)<<"You are in the else of DetWeighting!"
                                           <<"this should never happen!"
                                           <<"How did you get here?"<<endl;
  }

  /*
  if(new_reco_eshw+new_reco_emu<1.5){
    cout<<"Detetor weight "<<inu<<" ccnc "<<ccnc<<"true enu "<<true_enu<<" reco emu "<<reco_emu<<" new_reco_emu "<<new_reco_emu<<" reco eshw "<<reco_eshw<<" new_reco_eshw "<<new_reco_eshw<<" weight "<<dw<<endl;
  }
  */
  return dw;

}

    

void SKZPWeightCalculator::PrintReweightConfig(ostream &o)
{
  o<<"Reweighting to : "<<cfgname<<endl;
  for(unsigned int i=0;i<bs_pars.size();i++){
    o<<"Beam Sys par "<<i<<" "<<bs_pars[i]<<endl;
  }
  for(unsigned int i=0;i<hp_pars.size();i++){
    o<<"Had. Prod par "<<i<<" "<<hp_pars[i]<<endl;
  }
  for(unsigned int i=0;i<det_pars.size();i++){
    o<<"Det par "<<i<<" "<<det_pars[i]<<endl;
  }
}


void SKZPWeightCalculator::ChangeParameters(vector<double> *bs, vector<double> *hp, vector<double> *det)
{
  //this function is provided so that a user can compute weights for
  //many different szkp parameters.
  //it is provided for special studies, and I don't expect
  //most people will need this functionality.

  //clear anything that's already stored in the vectors
  bs_pars.clear();
  hp_pars.clear();
  det_pars.clear();

  //copy the function parameters into the private data members
  for(unsigned int i=0;i<bs->size();i++){
    double b=(*bs)[i];
    bs_pars.push_back(b);
  }
  for(unsigned int i=0;i<hp->size();i++){
    double h=(*hp)[i];
    hp_pars.push_back(h);
  }
  for(unsigned int i=0;i<det->size();i++){
    double d=(*det)[i];
    det_pars.push_back(d);
  }

  if(cfgnum==kPiMinus || cfgnum==kDetXs || cfgnum==kDogwood1_Daikon07){
    zfluk->SetParameters(hp_pars);
  }
  
}

void SKZPWeightCalculator::GetParameters(vector<double> &bs, vector<double> &hp, vector<double> &det)
{
  //this function is provided so that a user can compute weights for
  //many different szkp parameters.
  //it is provided for special studies, and I don't expect
  //most people will need this functionality.

  //clear anything that's already stored in the vectors
  bs.clear();
  hp.clear();
  det.clear();

  //copy the private data members into the function parameters
  for(unsigned int i=0;i<bs_pars.size();i++){
    bs.push_back(bs_pars[i]);
  }
  for(unsigned int i=0;i<hp_pars.size();i++){
    hp.push_back(hp_pars[i]);
  }
  for(unsigned int i=0;i<det_pars.size();i++){
    det.push_back(det_pars[i]);
  }
}

void SKZPWeightCalculator::SetConfigNum(string c)
{
  if(c=="PRL") cfgnum = kPRL;
  else if(c.find("PRL_Syst")!=string::npos) cfgnum = kPRL;
  else if(c=="Boston") cfgnum = kBoston;
  else if(c=="PiMinus") cfgnum = kPiMinus;
  else if(c=="PiMinus_CedarDaikon") cfgnum = kPiMinus;
  else if(c=="DetXs") cfgnum = kDetXs;
  else if(c=="Dogwood1_Daikon07") cfgnum = kDogwood1_Daikon07;
  else if(c=="Dogwood1_Daikon07_v2") cfgnum = kDogwood1_Daikon07_v2;
  else cfgnum = kUnknown;

}

double SKZPWeightCalculator::GetFluxError(int det, int Ibeam, 
                                          int inu, double true_enu, 
                                          FluxErrEffect_t neffect,double nsigma)
{
 zbeam->SetReweightConfig(cfgname);
  
  // Add target z with different than nsigma weight
  if(cfgname=="DetXs" && neffect==kTotalErrorAfterTune) {
  
    double e=0;
 
    double weight=zbeam->GetWeight(det,Ibeam,BeamSys::kHorn1Offset,nsigma,inu,true_enu)-1;
    e+=weight*weight;
    weight=zbeam->GetWeight(det,Ibeam,BeamSys::kBaffleScraping,nsigma,inu,true_enu)-1;
    e+=weight*weight;
    weight=zbeam->GetWeight(det,Ibeam,BeamSys::kPOT,nsigma,inu,true_enu)-1;
    e+=weight*weight;
    weight=zbeam->GetWeight(det,Ibeam,BeamSys::kHornIMiscal,nsigma,inu,true_enu)-1;
    e+=weight*weight;
    weight=zbeam->GetWeight(det,Ibeam,BeamSys::kHornIDist,nsigma,inu,true_enu)-1;
    e+=weight*weight;
    weight=zbeam->GetWeight(det,Ibeam,BeamSys::kBeamWidth,nsigma,inu,true_enu)-1;
    e+=weight*weight;
    weight=zbeam->GetWeight(det,Ibeam,BeamSys::kHadProdAfter,nsigma,inu,true_enu)-1;
    e+=weight*weight;
        
    // add error from best fit for target z

    // le010 configurations
    if(Ibeam==2 || Ibeam==6 || Ibeam==7) {
      double targetz=TMath::Max(TMath::Abs(bs_pars[0]),TMath::Abs(bs_pars[4]));
      weight=zbeam->GetWeight(det,Ibeam,BeamSys::kTargetZ,targetz*nsigma,inu,true_enu)-1;
      e+=weight*weight;
    }
    // le100
    else if(Ibeam==3) {
      double targetz=TMath::Abs(bs_pars[1]);
      weight=zbeam->GetWeight(det,Ibeam,BeamSys::kTargetZ,targetz*nsigma,inu,true_enu)-1;
      e+=weight*weight;
    }
    // le250
    else if(Ibeam==4) {
      double targetz=TMath::Max(TMath::Abs(bs_pars[3]),TMath::Abs(bs_pars[5]));
      weight=zbeam->GetWeight(det,Ibeam,BeamSys::kTargetZ,targetz*nsigma,inu,true_enu)-1;
      e+=weight*weight;
    }
    // le150
    else if(Ibeam==9) {
      double targetz=TMath::Abs(bs_pars[2]);
      weight=zbeam->GetWeight(det,Ibeam,BeamSys::kTargetZ,targetz*nsigma,inu,true_enu)-1;
      e+=weight*weight;
    }

    return 1+sqrt(e);
  }
  else return zbeam->GetWeight(det,Ibeam,neffect,nsigma,inu,true_enu);

}

int SKZPWeightCalculator::RunPeriodIntFromValidity(VldContext vc)
{
  int rp = 1; // default Run I
  VldTimeStamp ts = vc.GetTimeStamp();
  if(ts.GetDate()<(*vtsRunIEnd).GetDate()) rp = 1;
  else if(ts.GetDate()<(*vtsRunIIEnd).GetDate()) rp = 2;
  else if(ts.GetDate()<(*vtsRunIIIEnd).GetDate()) rp = 3;
  else rp = 4;
  return rp;
}

double SKZPWeightCalculator::GetRunPeriodWeight(RunPeriod_t rp,
                                                int det, 
                                                int Ibeam, 
                                                int ntype, 
                                                double Enu)
{
  double weight=1.;

  if(cfgnum==kDetXs || cfgnum==kDogwood1_Daikon07){

    // Target z weighting comes from the fit parameters
    Detector::Detector_t detector = Detector::kUnknown;
    if      (det==1) detector=Detector::kNear;
    else if (det==2) detector=Detector::kFar;
    
    Zbeam::ZbeamData_t zdata;
    zdata.ntype    = ntype;
    zdata.true_enu = Enu;
    zdata.detector = detector;
    zdata.beam     = BeamType::FromZarko(Ibeam);
    

    if(rp==kRunI){ // both weighting schemes have same parameters
      
      // target z effects
      if ( zdata.beam == BeamType::kL010z170i || 
           zdata.beam == BeamType::kL010z185i || 
           zdata.beam == BeamType::kL010z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[0]);
      } else if (zdata.beam == BeamType::kL100z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[1]);
      } else if (zdata.beam == BeamType::kL150z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[2]);
      } else if (zdata.beam == BeamType::kL250z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[3]);
      }
      
      // spot size
      if (zdata.beam == BeamType::kL250z200i ||
          zdata.beam == BeamType::kL100z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kBeamWidth,0.);
      }
    }
    else if(rp==kRunII) {
      
      if (zdata.beam == BeamType::kL010z185i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[4]);
      }
      else if (zdata.beam == BeamType::kL250z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[5]);
      }

      if(cfgnum==kDogwood1_Daikon07){ // apply target decay weight
        weight *= GetTargetDecayWeight(det,Ibeam,ntype,Enu,bs_pars[7]);
      }

    }

    else if(rp==kRunIII && cfgnum==kDetXs) { // use Run II beam position and Helium weight

      if (zdata.beam == BeamType::kL010z185i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[4]);
      }
      weight *= GetHeliumWeight(det,Ibeam,ntype,Enu);

      // DetXs target decay weight to be added when finalised

    }

    else if(rp==kRunIII && cfgnum==kDogwood1_Daikon07){

      if (zdata.beam == BeamType::kL010z185i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[6]);
      }
      weight *= GetTargetDecayWeight(det,Ibeam,ntype,Enu,bs_pars[8]);
      // Don't need Helium weight for Daikon07

    }

  } // end if(DetXs or Dogwood1_Daikon07)
  else if(cfgnum==kDogwood1_Daikon07_v2){

    // Target z weighting comes from the fit parameters
    Detector::Detector_t detector = Detector::kUnknown;
    if      (det==1) detector=Detector::kNear;
    else if (det==2) detector=Detector::kFar;
    
    Zbeam::ZbeamData_t zdata;
    zdata.ntype    = ntype;
    zdata.true_enu = Enu;
    zdata.detector = detector;
    zdata.beam     = BeamType::FromZarko(Ibeam);
    
    if(rp==kRunI){
      
      // target z effects
      if ( zdata.beam == BeamType::kL010z170i || 
           zdata.beam == BeamType::kL010z185i || 
           zdata.beam == BeamType::kL010z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[5]);
      } else if (zdata.beam == BeamType::kL100z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[6]);
      } else if (zdata.beam == BeamType::kL150z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[7]);
      } else if (zdata.beam == BeamType::kL250z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[8]);
      }
      
      // spot size
      if (zdata.beam == BeamType::kL250z200i ||
          zdata.beam == BeamType::kL100z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kBeamWidth,0.);
      }
    }
    else if(rp==kRunII) {
      
      if (zdata.beam == BeamType::kL010z185i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[9]);
      }
      else if (zdata.beam == BeamType::kL250z200i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[10]);
      }

      // apply target decay weight
      weight *= GetTargetDecayWeight(det,Ibeam,ntype,Enu,bs_pars[12]);

    }

    else if(rp==kRunIII){

      if (zdata.beam == BeamType::kL010z185i) {
        weight *= zbeam->GetWeight(zdata,BeamSys::kTargetZ ,bs_pars[11]);
      }

      // apply target decay weight
      weight *= GetTargetDecayWeight(det,Ibeam,ntype,Enu,bs_pars[13]);

      // Don't need Helium weight for Daikon07

    }

  }
  else {
    if(rp==kRunI){
      // cout<<"on run I"<<endl;
      //RUN I
      //ME and HE runs were lower intensity, smaller spot size,  
      if(BeamType::FromZarko(Ibeam)==BeamType::kL100z200i||
         BeamType::FromZarko(Ibeam)==BeamType::kL250z200i){
        //GetWeight with effect 6 already takes care of the beam type "if's"
        //but I thought I should just make them explicit here too
        weight=zbeam->GetWeight(det,Ibeam,6,0,ntype,Enu);
        //      cout<<"in run I weight"<<endl;
      }
    }
    else if(rp==kRunII || rp==kRunIII){
      //RUN II
      //reweight le10 beam to le09 beam
      //    cout<<"in run II weight "<<endl;
      if(BeamType::FromZarko(Ibeam)==BeamType::kL010z185i||
         BeamType::FromZarko(Ibeam)==BeamType::kL010z000i){
        weight=zbeam->GetWeight(det,Ibeam,15,-1.,ntype,Enu);
        //      cout<<"in run II weight "<<endl;
      }

      if(rp==kRunII){
        // target decay weight to be added when finalised, 0.6 should give ~average effect
        weight *= GetTargetDecayWeight(det,Ibeam,ntype,Enu,0.6);
      }
      if(rp==kRunIII){
        weight *= GetHeliumWeight(det,Ibeam,ntype,Enu);
        // target decay weight to be added when finalised, 1.25 should give ~average effect
        weight *= GetTargetDecayWeight(det,Ibeam,ntype,Enu,1.25);
      }

    }
    //  cout<<"run perioud weight "<<weight<<endl;
  }
  return weight;
}

double SKZPWeightCalculator::GetRFWWeight(const NtpMCTruth *mc, int det, int Ibeam,
                                          double reco_emu, double reco_eshw, 
                                          double &new_reco_emu, double &new_reco_eshw,
                                          double &beamweight, double &detweight)
{
  //retrieve truth information on the neutrino interaction from
  //the NtpMCTruth object
  int ccnc = mc->iaction;
  double true_enu = 1.*mc->p4neu[3];
  int inu = mc->inu;

  //pass the truth info and the fluxinfo object to the next
  //form of GetWeight
  return GetRFWWeight(&(mc->flux),det,Ibeam,
                      ccnc,true_enu,inu,
                      reco_emu,reco_eshw,
                      new_reco_emu,new_reco_eshw,
                      beamweight, detweight);
  
}

double SKZPWeightCalculator::GetRFWWeight(const NtpMCFluxInfo *fi, int det, int Ibeam,
                                          int ccnc, double true_enu, int inu,
                                          double reco_emu, double reco_eshw, 
                                          double &new_reco_emu, double &new_reco_eshw,
                                          double &beamweight, double &detweight)
{
  int parenttype=-1;
  double parentpz=0.;
  double parentpx=0.;
  double parentpy=0.;

  parenttype = fi->tptype;
  parentpz = 1.*fi->tpz;
  parentpy = 1.*fi->tpy;
  parentpx = 1.*fi->tpx;

  double pt = sqrt(parentpy*parentpy+parentpx*parentpx);

  //pass the variables on to the next version of GetWeight
  return GetRFWWeight(det,Ibeam,parenttype,pt,parentpz,ccnc,true_enu,inu,
                      reco_emu,reco_eshw,new_reco_emu,
                      new_reco_eshw,beamweight,detweight);

}

double SKZPWeightCalculator::GetRFWWeight(int det, int Ibeam,
                                          int tptype, double pt, double pz,    
                                          int ccnc, double true_enu,int inu,
                                          double reco_emu, double reco_eshw,
                                          double &new_reco_emu, double &new_reco_eshw,
                                          double &beamweight, double &detweight)
{

  //do the detector weighting
  detweight=GetDetWeight(ccnc,true_enu,inu,
                         reco_emu,reco_eshw,new_reco_emu,new_reco_eshw);
  beamweight=1.;
  if(runfrac.size()==0){
    MSG("SKZPWeightCalculator",Msg::kError)<<"You did not set the "
                                           <<"Run Period exposures."<<endl
                                           <<"I can not compute a "
                                           <<" fractional weight without "
                                           <<"knowing the relative exposures"<<endl;
    return 1.;
  }
    

  map<RunPeriod_t,double> beamc(runfrac.find(Ibeam)->second);
  map<RunPeriod_t,double>::iterator bcit(beamc.begin());
  double num=0.;
  double denom=0.;
  int i=0;
  while(bcit!=beamc.end()){
    RunPeriod_t rp = bcit->first;
    double pot = bcit->second;
    double w=GetBeamWeight(det,Ibeam,tptype,pt,pz,true_enu,inu,rp);
    num+=w*pot;
    denom+=pot;
    bcit++;
    //        cout<<"bcit->first "<<bcit->first<<endl;
    //    cout<<"bcit->second "<<bcit->second<<endl;

    //    cout<<"i "<<i<<" w "<<w<<" pot "<<pot<<" w*pot "<<w*pot<<" num "<<num<<" denom "<<denom<<endl;
    i++;
  }
    
    
  if(denom!=0){
    beamweight*=num/denom;
  }


  //return just the beam weight, as this is what 
  //most people will want
  //the value of the detweight is passed by reference
  //do not use the new values of new_reco_emu and new_reco_eshw
  //in the case of beam weighting only!!!!!
  return beamweight;
}

void SKZPWeightCalculator::SetRunFrac(int Ibeam, RunPeriod_t rp, double pot)
{

  map<int, map<RunPeriod_t,double> >::iterator it(runfrac.find(Ibeam));
  if(it!=runfrac.end()){
    //already have an entry for this beam type, just insert (runnum,pot)
    if((it->second).find(rp)==it->second.end()){
      (it->second)[rp]=pot;
    }
    else{
      MSG("SKZPWeightCalculator",Msg::kError)<<"You already specified pot for "
                                             <<BeamType::AsString(BeamType::FromZarko(Ibeam))
                                             <<" beam, with run period "<<rp<<endl
                                             <<"I will not overwrite!!!"<<endl;
    } 
  }
  else{
    //don't have this Ibeam yet, make a new map
    map<RunPeriod_t, double> m;
    m[rp]=pot;
    runfrac[Ibeam]=m;
  }


}

//-----------------------------------------------------------------------------------
Zbeam* SKZPWeightCalculator::GetZbeam()
{
   return zbeam;
}

//-----------------------------------------------------------------------------------
Zfluk* SKZPWeightCalculator::GetZfluk()
{
   return zfluk;
}



// Target decay weighting function

double SKZPWeightCalculator::GetTargetDecayWeight(int det, int Ibeam, int Ntype, double true_enu, double fin78fraction)
{

  double weight78 = 1.0;

  if(det==1){ // ND
    if(Ibeam==2){ // L010z185i
      if(Ntype==56 || Ntype==14) weight78 = GetMissingFin78Weight_L010z185i_ND_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) weight78 = GetMissingFin78Weight_L010z185i_ND_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) weight78 = GetMissingFin78Weight_L010z185i_ND_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) weight78 = GetMissingFin78Weight_L010z185i_ND_Nuebar(true_enu);
      else weight78 = 1.0;
    }
    else if(Ibeam==8){ // L010z000i
      if(Ntype==56 || Ntype==14) weight78 = GetMissingFin78Weight_L010z000i_ND_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) weight78 = GetMissingFin78Weight_L010z000i_ND_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) weight78 = GetMissingFin78Weight_L010z000i_ND_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) weight78 = GetMissingFin78Weight_L010z000i_ND_Nuebar(true_enu);
      else weight78 = 1.0;
    }
    else weight78 = 1.0;
  }
  else if(det==2){ // FD
    if(Ibeam==2){ // L010z185i
      if(Ntype==56 || Ntype==14) weight78 = GetMissingFin78Weight_L010z185i_FD_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) weight78 = GetMissingFin78Weight_L010z185i_FD_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) weight78 = GetMissingFin78Weight_L010z185i_FD_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) weight78 = GetMissingFin78Weight_L010z185i_FD_Nuebar(true_enu);
      else weight78 = 1.0;
    }
    else if(Ibeam==8){ // L010z000i
      if(Ntype==56 || Ntype==14) weight78 = GetMissingFin78Weight_L010z000i_FD_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) weight78 = GetMissingFin78Weight_L010z000i_FD_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) weight78 = GetMissingFin78Weight_L010z000i_FD_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) weight78 = GetMissingFin78Weight_L010z000i_FD_Nuebar(true_enu);
      else weight78 = 1.0;
    }
    else weight78 = 1.0;
  }
  else weight78 = 1.0;

  // Now convert to fractional weight
  // eg: if want to remove 4 fins then double the distance from unity for the weight

  double weight = 1.;

  if(weight78<0.){
    weight = 1.0-(fin78fraction*(1.0-weight78));
  }
  else if(weight78>0.){
    weight = 1.0+(fin78fraction*(weight78-1.0));
  }

  return weight;

}



// Helium weighting function

double SKZPWeightCalculator::GetHeliumWeight(int det, int Ibeam, int Ntype, double true_enu)
{
  if(det==1){ // ND
    if(Ibeam==2){ // L010z185i
      if(Ntype==56 || Ntype==14) return GetHeliumWeight_L010z185i_ND_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) return GetHeliumWeight_L010z185i_ND_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) return GetHeliumWeight_L010z185i_ND_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) return GetHeliumWeight_L010z185i_ND_Nuebar(true_enu);
      else return 1.0;
    }
    else if(Ibeam==8){ // L010z000i
      if(Ntype==56 || Ntype==14) return GetHeliumWeight_L010z000i_ND_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) return GetHeliumWeight_L010z000i_ND_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) return GetHeliumWeight_L010z000i_ND_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) return GetHeliumWeight_L010z000i_ND_Nuebar(true_enu);
      else return 1.0;
    }
    else return 1.0;
  }
  else if(det==2){ // FD
    if(Ibeam==2){ // L010z185i
      if(Ntype==56 || Ntype==14) return GetHeliumWeight_L010z185i_FD_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) return GetHeliumWeight_L010z185i_FD_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) return GetHeliumWeight_L010z185i_FD_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) return GetHeliumWeight_L010z185i_FD_Nuebar(true_enu);
      else return 1.0;
    }
    else if(Ibeam==8){ // L010z000i
      if(Ntype==56 || Ntype==14) return GetHeliumWeight_L010z000i_FD_Numu(true_enu);
      else if(Ntype==55 || Ntype==-14) return GetHeliumWeight_L010z000i_FD_Numubar(true_enu);
      else if(Ntype==53 || Ntype==12) return GetHeliumWeight_L010z000i_FD_Nue(true_enu);
      else if(Ntype==52 || Ntype==-12) return GetHeliumWeight_L010z000i_FD_Nuebar(true_enu);
      else return 1.0;
    }
    else return 1.0;
  }
  else return 1.0;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_ND_Numu(double enu)
{
  double weight = 1.;
  double par[19] = {
    1.10472,-0.0861419,0.0116673,0.0193203,0.230526,0.323048,-0.0434165,0.00189932,1.1483,
    -0.0322421,0.00224179,-4.51111e-05,0.818716,0.0177147,-0.000366594,2.1493e-06,1.20406,
    -0.00411679,2.18602e-05
  };
  if(enu<5.15){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*TMath::Log(enu);
  }
  else if(enu<7.25){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu;
  }
  else if(enu<18.0){
    weight = par[8]+par[9]*enu+par[10]*enu*enu+par[11]*enu*enu*enu;
  }
  else if(enu<53.5){
    weight = par[12]+par[13]*enu+par[14]*enu*enu+par[15]*enu*enu*enu;
  }
  else{
    weight = par[16]+par[17]*enu+par[18]*enu*enu;
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_ND_Numubar(double enu)
{
  double weight = 1.;
  double par[15] = {
    1.06239,0.0342807,-0.00879299,0.000479629,-0.346922,0.267095,-0.016762,0.00033883,
    0.8,0.01,-0.0005,1e-06,1.03,0.0042195,0.000108851
  };
  if(enu<11.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu;
  }
  else if(enu<21.4){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu;
  }
  else if(enu<37.){
    weight = par[8]+par[9]*enu+par[10]*(enu-20.)*(enu-20.)+par[11]*(enu-20.)*(enu-20.)*(enu-20.);
  }
  else{
    weight = par[12]-par[13]*(enu-37.)+par[14]*(enu-37.)*(enu-37.);
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_ND_Nue(double enu)
{
  double weight = 1.;
  double par[11] = {
    0.97709,0.0939192,-0.0536218,0.00805987,1.03,-0.0163451,0.00263497,
    -0.000118506,1.57655e-06,0.998,600
  };
  if(enu<2.43){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu;
  }
  else if(enu<23.0){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu+par[8]*enu*enu*enu*enu;
  }
  else{
    weight = par[9]+(par[10]/(enu*enu*enu));
  }
  return weight;

}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_ND_Nuebar(double enu)
{
  double weight = 1.;
  double par[7] = {
    1.02113,0.0140663,-0.000791038,1.37577e-05,-7.48181e-08,0.9995,1150
  };
  if(enu<30.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_FD_Numu(double enu)
{
  double weight = 1.;
  double par[14] = {
    1.04327,-0.0352661,0.0047907,0.00388763,0.88061,0.0321005,-0.00164144,
    -4.5869e-05,1.05352,-0.00869623,0.000390208,-4.45885e-06,0.997,5000
  };
  if(enu<4.4){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*TMath::Log(enu);
  }
  else if(enu<8.9){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu;
  }
  else if(enu<48.0){
    weight = par[8]+par[9]*enu+par[10]*enu*enu+par[11]*enu*enu*enu;
  }
  else{
    weight = par[12]+(par[13]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_FD_Numubar(double enu)
{
  double weight = 1.;
  double par[10] = {
    1.04689,0.00226136,-0.000790701,2.61669e-05,-0.0744887,
    0.0847917,-0.00190558,1.21492e-05,0.99,4000
  };
  if(enu<21.9){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu;
  }
  else if(enu<38.0){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu;
  }
  else{
    weight = par[8]+(par[9]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_FD_Nue(double enu)
{
  double weight = 1.;
  double par[8] = {
    1.00099,0.022962,-0.0106913,0.00153978,-6.74999e-05,1.021,-0.0018,1.5e-05
  };
  if(enu<11.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]*(enu-11.))+(par[7]*(enu-11.)*(enu-11.));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z185i_FD_Nuebar(double enu)
{
  double weight = 1.;
  double par[7] = {
    1.02439,-0.000679498,0.000218119,-6.49776e-06,4.58303e-08,0.997,3300
  };
  if(enu<43.1){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_ND_Numu(double enu)
{
  double weight = 1.;
  double par[10] = {
    1.08341,-0.00598166,-0.000644125,0.0153799,1.16251,
    -0.0259465,0.00130708,-1.8318e-05,1.03,1500
  };
  if(enu<7.4){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*TMath::Log(enu);
  }
  else if(enu<38.4){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu;
  }
  else{
    weight = par[8]+(par[9]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_ND_Numubar(double enu)
{
  double weight = 1.;
  double par[15] = {
    1.06489,0.0137669,-0.00415539,0.000218168,0.923382,0.0108634,-0.000376271,
    3.86015e-06,0.8,0.01,-0.0005,1e-06,1.03,0.00883329,0.000347229
  };
  if(enu<11.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu;
  }
  else if(enu<21.4){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu+0.0014;
  }
  else if(enu<37.){
    weight = par[8]+par[9]*enu+par[10]*(enu-20.)*(enu-20.)+par[11]*(enu-20.)*(enu-20.)*(enu-20.)+0.0097;
  }
  else{
    weight = par[12]-par[13]*(enu-37.)+par[14]*(enu-37.)*(enu-37.)+0.0097;
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_ND_Nue(double enu)
{
  double weight = 1.;
  double par[7] = {
    1.02614,0.0135056,-0.00151707,5.53927e-05,-6.55825e-07,0.978,1.
  };
  if(enu<20.45){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]/(enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_ND_Nuebar(double enu)
{
  double weight = 1.;
  double par[7] = {
    1.02609,0.0255862,-0.00296928,9.18715e-05,-7.22365e-07,1.101,-2500
  };
  if(enu<28.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_FD_Numu(double enu)
{
  double weight = 1.;
  double par[10] = {
    1.04313,-0.00624758,0.000104616,0.00872899,1.12567,
    -0.0171559,0.000639488,-6.72704e-06,1.0197,1500
  };
  if(enu<10.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*TMath::Log(enu);
  }
  else if(enu<46.0){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu;
  }
  else{
    weight = par[8]+(par[9]/(enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_FD_Numubar(double enu)
{
  double weight = 1.;
  double par[11] = {
    1.0378,-0.000996055,-0.000280837,0.00197533,1.25837,-0.0519768,
    0.00349608,-9.25331e-05,8.19759e-07,1,-1.5e+07
  };
  if(enu<9.5){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*TMath::Log(enu);
  }
  else if(enu<49.3){
    weight = par[4]+par[5]*enu+par[6]*enu*enu+par[7]*enu*enu*enu+par[8]*enu*enu*enu*enu;
  }
  else{
    weight = par[9]+(par[10]/(enu*enu*enu*enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_FD_Nue(double enu)
{
  double weight = 1.;
  double par[7] = {
    1.02653,-0.00287051,0.000666625,-4.19876e-05,7.69875e-07,0.99805,10
  };
  if(enu<20.0){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]/(enu*enu));
  }
  return weight;
}

double SKZPWeightCalculator::GetHeliumWeight_L010z000i_FD_Nuebar(double enu)
{
  double weight = 1.;
  double par[7] = {
    1.02206,0.00932646,-0.00100526,2.37525e-05,-4.71822e-08,1.2,-180
  };
  if(enu<33.45){
    weight = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else{
    weight = par[5]+(par[6]/(enu*enu));
  }
  return weight;
}



double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_ND_Numu(double enu)
{
  double y;

  double par[20];

  par[0] =  9.91823e-01;
  par[1] =  8.14956e-02;
  par[2] = -1.44191e-01;
  par[3] =  7.80651e-02;
  par[4] = -1.29277e-02;
  par[5] = -1.41734e-04;
  par[6] =  3.62036e+00;
  par[7] = -8.98865e-01;
  par[8] =  9.61367e-02;
  par[9] = -2.44885e+00;
  par[10] =  6.81085e-01;
  par[11] =  9.40087e-02;
  par[12] = -6.50227e-03;
  par[13] = -2.86720e+00;
  par[14] = -1.14010e-01;
  par[15] =  1.06866e-03;
  par[16] =  1.77433e+00;
  par[17] =  8.30457e+00;
  par[18] =  1.018;
  par[19] =  0.005;

  if(enu<2.875){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else if(enu<4.375){
    y = par[6]+par[7]*enu+par[8]*enu*enu+par[9]*(1./enu);
  }
  else if(enu<7.460){
    y = par[10]+par[11]*enu+par[12]*enu*enu;
  }
  else if(enu<15.15){
    y = par[13]+par[14]*enu+par[15]*enu*enu+par[16]*TMath::Log(enu)+par[17]*(1./enu);
  }
  else{
    y = par[18]+par[19]*(1./(enu-13.8));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_ND_Numubar(double enu)
{
  double y;

  double par[11];

  par[0] = 1.01079;
  par[1] = 0.0163014;
  par[2] = -0.00951575;
  par[3] = 0.00189497;
  par[4] = -0.000114420;
  par[5] = 0.745891;
  par[6] = 0.0198691;
  par[7] = -0.000440386;
  par[8] = 1.22518;
  par[9] = 1.04;
  par[10] = -0.2;

  if(enu<8.6){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else if(enu<15.5){
    y = par[5]+par[6]*enu+par[7]*enu*enu+par[8]*(1./enu);
  }
  else {
    y = par[9]+par[10]*(1./enu);
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_ND_Numu(double enu)
{
  double y;

  double par[7];

  par[0] = 1.01052;
  par[1] = 0.00135492;
  par[2] = 0.000352042;
  par[3] = -0.0000281668;
  par[4] = 0.000000465827;
  par[5] = 1.018;
  par[6] = 0.3;

  if(enu<17.2){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else {
    y = par[5]+par[6]*(1./enu);
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_ND_Numubar(double enu)
{
  double y;

  double par[7];

  par[0] = 1.00806;
  par[1] = 0.00509031;
  par[2] = -0.000422559;
  par[3] = 0.0000180352;
  par[4] = -0.000000272726;
  par[5] = 1.018;
  par[6] = 0.85;

  if(enu<30.5){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else {
    y = par[5]+par[6]*(1./enu);
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_FD_Numu(double enu)
{
  double y;

  double par[20];

  par[0] =   9.79176e-01;
  par[1] =   5.94877e-02;
  par[2] =  -6.33591e-02;
  par[3] =   1.24537e-02;
  par[4] =   6.85855e-03;
  par[5] =  -2.08938e-03;
  par[6] =   2.56824e+00;
  par[7] =  -5.16706e-01;
  par[8] =   5.21584e-02;
  par[9] =  -1.52990e+00;
  par[10] =   7.10491e-01;
  par[11] =   7.95342e-02;
  par[12] =  -5.11905e-03;
  par[13] =  -4.11699e+00;
  par[14] =  -1.27508e-01;
  par[15] =   1.08948e-03;
  par[16] =   2.22361e+00;
  par[17] =   1.17195e+01;
  par[18] = 1.018;
  par[19] = 0.005;

  if(enu<2.83){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else if(enu<4.925){
    y = par[6]+par[7]*enu+par[8]*enu*enu+par[9]*(1./enu);
  }
  else if(enu<9.0){
    y = par[10]+par[11]*enu+par[12]*enu*enu;
  }
  else if(enu<17.6){
    y = par[13]+par[14]*enu+par[15]*enu*enu+par[16]*TMath::Log(enu)+par[17]*(1./enu);
  }
  else{
    y = par[18]+par[19]*(1./(enu-13.8));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_FD_Numubar(double enu)
{
  double y;

  double par[11];

  par[0] = 0.992562;
  par[1] = 0.0301506;
  par[2] = -0.0129656;
  par[3] = 0.00215436;
  par[4] = -0.000115552;
  par[5] = 0.958894;
  par[6] = 0.00228454;
  par[7] = -0.0000164352;
  par[8] = 0.468619;
  par[9] = 1.04;
  par[10] = -0.2;

  if(enu<8.35){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else if(enu<38.75){
    y = par[5]+par[6]*enu+par[7]*enu*enu+par[8]*(1./enu);
  }
  else {
    y = par[9]+par[10]*(1./enu);
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_FD_Numu(double enu)
{
  double y;

  double par[7];

  par[0] = 0.991629;
  par[1] = 0.0150514;
  par[2] = -0.00358913;
  par[3] = 0.000384422;
  par[4] = -0.0000136394;
  par[5] = 1.008;
  par[6] = 0.35;

  if(enu<13.38){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else {
    y = par[5]+par[6]*(1./enu);
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_FD_Numubar(double enu)
{
  double y;

  double par[7];

  par[0] = 0.991063;
  par[1] = 0.0114890;
  par[2] = -0.00123288;
  par[3] = 0.0000558658;
  par[4] = -0.000000834802;
  par[5] = 1.056;
  par[6] = -0.35;

  if(enu<17.925){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu;
  }
  else {
    y = par[5]+par[6]*(1./enu);
  }

  return y;
}



double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_ND_Nue(double enu)
{
  double y;

  double par[12];

  par[0] = 0.964536;
  par[1] = 0.0341348;
  par[2] = -0.0167779;
  par[3] = 0.00275733;
  par[4] = -0.000173092;
  par[5] = 0.00000364374;
  par[6] = 65.6200;
  par[7] = 1.55324;
  par[8] = -0.0132149;
  par[9] = -27.7156;
  par[10] = -147.925;
  par[11] = 1.11675;

  if(enu<15.70){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else if(enu<29.8){
    y = -0.012+par[6]+par[7]*enu+par[8]*enu*enu+par[9]*TMath::Log(enu)+par[10]*(1./enu);
  }
  else{
    y = par[11];
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_ND_Nuebar(double enu)
{
  double y;

  double par[8];

  par[0] = 1.01535;
  par[1] = -0.0167946;
  par[2] = 0.00921187;
  par[3] = -0.00133391;
  par[4] = 0.0000684085;
  par[5] = -0.00000113447;
  par[6] = 0.9;
  par[7] = 1.0;

  if(enu<26.95){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else {
    y = par[6]+par[7]*(1./(enu-22.));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_ND_Nue(double enu)
{
  double y;

  double par[8];

  par[0] = 1.09927;
  par[1] = -0.0580560;
  par[2] = 0.00847919;
  par[3] = -0.000341542;
  par[4] = -0.00000545052;
  par[5] = 0.000000394652;
  par[6] = 1.15;
  par[7] = -1.0;

  if(enu<20.8){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else {
    y = par[6]+par[7]*(1./(enu-14.));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_ND_Nuebar(double enu)
{
  double y;

  double par[8];

  par[0] = 0.960938;
  par[1] = 0.169592;
  par[2] = -0.0611788;
  par[3] = 0.00782170;
  par[4] = -0.000416337;
  par[5] = 0.00000785160;
  par[6] = 1.12;
  par[7] = -1.0;

  if(enu<19.775){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else {
    y = par[6]+par[7]*(1./(enu-13.8));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_FD_Nue(double enu)
{
  double y;

  double par[8];

  par[0] = 0.976334;
  par[1] = 0.0146822;
  par[2] = -0.00805149;
  par[3] = 0.00132456;
  par[4] = -0.0000810696;
  par[5] = 0.00000166323;
  par[6] = 1.0851;
  par[7] = -1.0;

  if(enu<19.45){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else{
    y = par[6]+par[7]*(1./(enu-9.));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z185i_FD_Nuebar(double enu)
{
  double y;

  double par[9];

  par[0] = 0.912223;
  par[1] = 0.106252;
  par[2] = -0.0341165;
  par[3] = 0.00490576;
  par[4] = -0.000344942;
  par[5] = 0.0000114327;
  par[6] = -0.000000142208;
  par[7] = 0.9;
  par[8] = 1.0;

  if(enu<27.55){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu+par[6]*enu*enu*enu*enu*enu*enu;
  }
  else {
    y = par[7]+par[8]*(1./(enu-22.));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_FD_Nue(double enu)
{
  double y;

  double par[8];

  par[0] = 1.05046;
  par[1] = -0.0197051;
  par[2] = -0.000709986;
  par[3] = 0.000580179;
  par[4] = -0.0000443486;
  par[5] = 0.000000962613;
  par[6] = 0.92;
  par[7] = 0.6;

  if(enu<19.5){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else {
    y = par[6]+par[7]*(1./(enu-12.));
  }

  return y;
}

double SKZPWeightCalculator::GetMissingFin78Weight_L010z000i_FD_Nuebar(double enu)
{
  double y;

  double par[8];

  par[0] = 1.03867;
  par[1] = 0.0130506;
  par[2] = -0.00751797;
  par[3] = 0.00109883;
  par[4] = -0.0000639034;
  par[5] = 0.00000130481;
  par[6] = 1.04;
  par[7] = -0.2;

  if(enu<18.9){
    y = par[0]+par[1]*enu+par[2]*enu*enu+par[3]*enu*enu*enu+par[4]*enu*enu*enu*enu+par[5]*enu*enu*enu*enu*enu;
  }
  else {
    y = par[6]+par[7]*(1./(enu-12.));
  }

  return y;
}

---