NuZBeamReweight.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell July/2006       
//                                                                    
// Contact: j.j.hartnell@rl.ac.uk

#include "TClonesArray.h"
#include "TDirectory.h"

#include "Conventions/BeamType.h"
#include "MessageService/MsgService.h"
#include "MCNtuple/NtpMCTruth.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "TruthHelperNtuple/NtpTHEvent.h"
#include "StandardNtuple/NtpStRecord.h"
#include "Conventions/ReleaseType.h"

#include "MCReweight/MCReweight.h"
#include "MCReweight/NeugenWeightCalculator.h"
#include "MCReweight/ReweightHelpers.h"
#include "MCReweight/SKZPWeightCalculator.h"
#include "MCReweight/Zbeam.h"
#include "MCReweight/Zfluk.h"

#include "NtupleUtils/NuUtilities.h"
#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuZBeamReweight.h"

CVSID("$Id: NuZBeamReweight.cxx,v 1.32 2010/02/09 00:06:46 ahimmel Exp $");

//......................................................................

NuZBeamReweight::NuZBeamReweight()
{
  MSG("NuZBeamReweight",Msg::kDebug)
    <<"Running NuZBeamReweight Constructor..."<<endl;


  MSG("NuZBeamReweight",Msg::kDebug)
    <<"Finished NuZBeamReweight Constructor"<<endl;
}

//......................................................................

NuZBeamReweight::~NuZBeamReweight()
{
  MSG("NuZBeamReweight",Msg::kDebug)
    <<"Running NuZBeamReweight Destructor..."<<endl;
  

  MSG("NuZBeamReweight",Msg::kDebug)
    <<"Finished NuZBeamReweight Destructor"<<endl;
}

//......................................................................

void NuZBeamReweight::ExtractZBeamReweight(NuEvent& nu, bool useCustom) const
{

  MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"Running ExtractZBeamReweight..."<<endl;
  
  if (!this->EventPreCheck(nu)) return;
  
//  // Use NA49 NuMu-only weights for RHC
//  if (nu.hornIsReverse || 
//      nu.runPeriod == 4 || 
//      nu.beamType == BeamType::kL010z185i_rev) {
//    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
//    <<"Using numu-only weights for RHC"<<endl;
//    useCustom = true;
//  }
  
  //get a reference
  SKZPWeightCalculator* wc;
  if (useCustom) {
    wc = this->GetSKZPWeightCalculatorCustom();
  }
  else {
    wc = this->GetSKZPWeightCalculator(nu);
  }
    

  // If runPeriod == 0, then treat as Daikon04 and earlier
  if (nu.runPeriod == 0) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"Getting old-style weights for nu.runPeriod="<< nu.runPeriod <<endl;

    //calc the weights for runI
    GetWeights(nu, *wc, SKZPWeightCalculator::kRunI);
    nu.trkEnWeightRunI                = nu.trkEnWeight;
    nu.shwEnWeightRunI                = nu.shwEnWeight;
    nu.beamWeightRunI                 = nu.beamWeight;
    nu.fluxErrHadProdAfterTuneRunI    = nu.fluxErrHadProdAfterTune;
    nu.fluxErrTotalErrorPreTuneRunI   = nu.fluxErrTotalErrorPreTune;
    nu.fluxErrTotalErrorAfterTuneRunI = nu.fluxErrTotalErrorAfterTune;
    nu.detectorWeightNMBRunI          = nu.detectorWeightNMB;
    nu.detectorWeightNMRunI           = nu.detectorWeightNM;
    
    
    //calc the weights for runII
    GetWeights(nu, *wc, SKZPWeightCalculator::kRunII);
    nu.trkEnWeightRunII                = nu.trkEnWeight;
    nu.shwEnWeightRunII                = nu.shwEnWeight;
    nu.beamWeightRunII                 = nu.beamWeight;
    nu.fluxErrHadProdAfterTuneRunII    = nu.fluxErrHadProdAfterTune;
    nu.fluxErrTotalErrorPreTuneRunII   = nu.fluxErrTotalErrorPreTune;
    nu.fluxErrTotalErrorAfterTuneRunII = nu.fluxErrTotalErrorAfterTune;
    nu.detectorWeightNMBRunII          = nu.detectorWeightNMB;
    nu.detectorWeightNMRunII           = nu.detectorWeightNM;
    
    
    //calc the weights for the weighted average
    GetRFWeights(nu, *wc);
  }
  else {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"Getting new-style weights for nu.runPeriod="<< nu.runPeriod <<endl;

    GetWeights(nu, *wc);
  }
}

//......................................................................

void NuZBeamReweight::ExtractZBeamReweightCustom(NuEvent& nu) const
{ 
  this->ExtractZBeamReweight(nu, true);
}


//......................................................................

void NuZBeamReweight::ExtractZBeamReweight(const NtpStRecord& /*ntp*/,
                                           const NtpSREvent& /*evt*/,
                                           NuEvent& nu) const
{ 
  //get the weights
  this->ExtractZBeamReweight(nu);
}

//......................................................................

void NuZBeamReweight::ExtractZBeamReweight(const NtpStRecord& /*ntp*/,
                                           NuEvent& nu) const
{ 
  MAXMSG("NuZBeamReweight",Msg::kError,2000)
    <<"NuZBeamReweight::ExtractZBeamReweight(const NtpStRecord& ntp,NuEvent& nu) const"
    <<" is deprecated, passing off to ExtractZBeamReweight(NuEvent& nu)!"<<endl;
  this->ExtractZBeamReweight(nu);  
}


//......................................................................

void NuZBeamReweight::GetBeamWeightsOnly(NuEvent& nu) const
{
  
  MAXMSG("NuZBeamReweight",Msg::kInfo,1)
  <<"Running GetZBeamReweight..."<<endl;
  
  if (!this->EventPreCheck(nu)) return;
  
  //get a reference
  SKZPWeightCalculator* wc=this->GetSKZPWeightCalculator(nu);
  
  // If runPeriod == 0, then treat as Daikon04 and earlier
  if (nu.runPeriod == 0) {
    //calc the weights for runI
    GetWeights(nu, *wc, SKZPWeightCalculator::kRunI, true);
    nu.beamWeightRunI                 = nu.beamWeight;
    nu.fluxErrHadProdAfterTuneRunI    = nu.fluxErrHadProdAfterTune;
    nu.fluxErrTotalErrorPreTuneRunI   = nu.fluxErrTotalErrorPreTune;
    nu.fluxErrTotalErrorAfterTuneRunI = nu.fluxErrTotalErrorAfterTune;
    
    
    //calc the weights for runII
    GetWeights(nu, *wc, SKZPWeightCalculator::kRunII, true);
    nu.beamWeightRunII                 = nu.beamWeight;
    nu.fluxErrHadProdAfterTuneRunII    = nu.fluxErrHadProdAfterTune;
    nu.fluxErrTotalErrorPreTuneRunII   = nu.fluxErrTotalErrorPreTune;
    nu.fluxErrTotalErrorAfterTuneRunII = nu.fluxErrTotalErrorAfterTune;
    
    
    //calc the weights for the weighted average
    GetRFWeights(nu, *wc);
  }
  else {
    GetWeights(nu, *wc, true);
  }
}

//......................................................................

SKZPWeightCalculator* NuZBeamReweight::
GetSKZPWeightCalculator(const NuEvent& nu) const
{
  
  MAXMSG("NuZBeamReweight",Msg::kInfo,1)
  <<"Loading normal weights"<<endl;

  //make a weight calculator, defining the skzp version to use
  static SKZPWeightCalculator* pwc=0;
  
  //create the weight calculator on the first call
  if (pwc==0) {
    //get a string of the reweightVersion
    string sReweightVersion=this->AsString
      (static_cast<SKZPWeightCalculator::SKZPConfig_t>
       (nu.reweightVersion));
    MSG("NuZBeamReweight",Msg::kInfo)
      <<"To access database, using the string="<<sReweightVersion<<endl;

    TDirectory* tmpd = gDirectory;
    MSG("NuZBeamReweight",Msg::kDebug)
      <<"TDirectory before Reweight is:"<<endl;
    //tmpd->Print();
    
    //create the weight calculator
    //read constants from DB
    pwc=new SKZPWeightCalculator(sReweightVersion,true);

    // Can always do this since it will only be accessed when getting RFWeights
    //define and set the numbers of POTs in each run period
    //numbers from Andy Blakes minos-doc-3727
    Float_t potRunI=1.269e20;
    //runIIa=1.229
    //runIIb=0.721
    Float_t potRunII=1.950e20;
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"ExtractZBeamReweight: SetRunFrac w/ potRunI="<<potRunI
    <<", potRunII="<<potRunII<<endl;
    
    Int_t Ibeam;
    if (nu.beamType == BeamType::kL010z185i_rev) {
      MAXMSG("NuZBeamReweight",Msg::kInfo,1) 
      << "Changing beamtype from kL010z185i_rev to kL010z185i for SKZP weighting" << endl;
      Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(BeamType::kL010z185i));
    }
    else {
      Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(nu.beamType));
    }
    
    
    pwc->SetRunFrac(Ibeam,SKZPWeightCalculator::kRunI,potRunI);
    pwc->SetRunFrac(Ibeam,SKZPWeightCalculator::kRunII,potRunII);       
    
    //print the config
    pwc->PrintReweightConfig(cout);

    MSG("NuZBeamReweight",Msg::kDebug)
      <<"TDirectory after building Zbeam and Zflux is:"<<endl;
    //gDirectory->Print();
    MSG("NuZBeamReweight",Msg::kDebug)
      <<"After reseting gDirectory its location is:"<<endl;
    gDirectory=tmpd;
    //gDirectory->Print();
  }
  
  return pwc;
}

//......................................................................

SKZPWeightCalculator* NuZBeamReweight::
GetSKZPWeightCalculatorCustom() const
{
  MAXMSG("NuZBeamReweight",Msg::kInfo,1)
  <<"Loading numu-only weights"<<endl;

  //make a weight calculator, defining the skzp version to use
  static SKZPWeightCalculator* pwc=0;
  
  //create the weight calculator on the first call
  if (pwc==0) {
    
    TDirectory* tmpd = gDirectory;
    MSG("NuZBeamReweight",Msg::kDebug) <<"TDirectory before Reweight is:"<<endl;
    
    //create the weight calculator
    //read constants from DB
    pwc=new SKZPWeightCalculator("DetXs",true);
    
    pwc->PrintReweightConfig(cout);
    
    std::vector<double> beam;
    std::vector<double> fluk;
    std::vector<double> det;
    
    MSG("NuZBeamReweight",Msg::kInfo) << std::endl;
    MSG("NuZBeamReweight",Msg::kInfo) << "Setting new target positions." << std::endl;
    
    // target positions
    beam.push_back(1.52410);
    beam.push_back(-3.34533);
    beam.push_back(-3.68550);
    beam.push_back(0.177135);
    beam.push_back(-1.73917);
    beam.push_back(-0.0842292);
    
    MSG("NuZBeamReweight",Msg::kInfo) << std::endl;
    MSG("NuZBeamReweight",Msg::kInfo) << "Setting new hadron production parameters." << std::endl;
    
    // pi+
    fluk.push_back(0.53363);
    fluk.push_back(-4.09962);
    fluk.push_back(0.959969);
    fluk.push_back(0.464735);
    fluk.push_back(0.948545);
    fluk.push_back(1.22209);
    
    // K+
    fluk.push_back(5.69124);
    fluk.push_back(17.4408);
    fluk.push_back(1.49750);
    fluk.push_back(0.972495);
    fluk.push_back(0.968975);
    fluk.push_back(2.19318);
    
    // pi-
    fluk.push_back(1.05039);
    fluk.push_back(-0.867787);
    
    // K-
    fluk.push_back(0.842913);
    fluk.push_back(0.0595505);
    
    MSG("NuZBeamReweight",Msg::kInfo) << std::endl;
    MSG("NuZBeamReweight",Msg::kInfo) << "Setting new detector parameters." << std::endl;
    
    // detector parameters
    det.push_back(-0.0183836);
    det.push_back(-0.244669);
    det.push_back(0.0); // numubar NC at nominal
    det.push_back(1.0); // numubar CC xsec at nominal
    det.push_back(0.342342);
    det.push_back(0.00517944);
    det.push_back(0.389318);
    det.push_back(-0.136436);
    
    MSG("NuZBeamReweight",Msg::kInfo) << std::endl;
    MSG("NuZBeamReweight",Msg::kInfo) << "Calling change parameters." << std::endl;
    
    pwc->ChangeParameters(&beam,&fluk,&det);
    
    MSG("NuZBeamReweight",Msg::kInfo) << std::endl;
    MSG("NuZBeamReweight",Msg::kInfo) << "Changed parameters, new configuration:" << std::endl;
    MSG("NuZBeamReweight",Msg::kInfo) << std::endl;
      
    
    // Can always do this since iy will only be accessed when getting RFWeights
    //define and set the numbers of POTs in each run period
    //numbers from Andy Blakes minos-doc-3727
    Float_t potRunI=1.269e20;
    //runIIa=1.229
    //runIIb=0.721
    Float_t potRunII=1.950e20;
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"ExtractZBeamReweight: SetRunFrac w/ potRunI="<<potRunI
    <<", potRunII="<<potRunII<<endl;
    const Int_t Ibeam=BeamType::ToZarko(BeamType::kL010z185i);
    pwc->SetRunFrac(Ibeam,SKZPWeightCalculator::kRunI,potRunI);
    pwc->SetRunFrac(Ibeam,SKZPWeightCalculator::kRunII,potRunII);   
    
    
    //print the config
    pwc->PrintReweightConfig(cout);
    
    MSG("NuZBeamReweight",Msg::kDebug)
    <<"TDirectory after building Zbeam and Zflux is:"<<endl;
    //gDirectory->Print();
    MSG("NuZBeamReweight",Msg::kDebug)
    <<"After reseting gDirectory its location is:"<<endl;
    gDirectory=tmpd;
    //gDirectory->Print();
  }
  
  return pwc;
}

//......................................................................

void NuZBeamReweight::CalcGeneratorReweight(const NtpStRecord& ntp,
                                            const NtpSREvent& evt,
                                            NuEvent& nu) const
{

  //return if not MC
  if (nu.simFlag!=SimFlag::kMC) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
      <<"Not doing CalcGeneratorReweight"
      <<"for simFlag="<<nu.simFlag<<endl;
    
    //set the weights to 1 that would otherwise be set below
    nu.generatorWeight=1;//no reweighting
    return;
  }
 
  if (!nu.useGeneratorReweight) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
      <<"CalcGeneratorReweight: not calculating Generator reweight"
      <<" since nu.useGeneratorReweight="<<nu.useGeneratorReweight
      <<endl;
    return;
  }

  if (ReleaseType::IsDaikon(nu.releaseType)) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
      <<"Not calculating Generator weight for Daikon, not required"
      <<endl;
    return;
  }
  
  MAXMSG("NuZBeamReweight",Msg::kDebug,100)
    <<"CalcGeneratorReweight:"<<endl
    <<"  generator:config_name="<<nu.sGeneratorConfigName
    <<"  generator:config_no="<<nu.generatorConfigNo
    <<endl;
  
  //declare the registry and create it only once
  static Registry* rwtconfig=this->CreateNeugenRegistry(nu);

  //make a WeightCalculator and add it to the MCReweight singleton
  this->AddNeugenWeightCalculator();

  //get the MC info from the ntp
  TClonesArray& thevtTca=*ntp.thevt;
  const NtpTHEvent* the=
    dynamic_cast<NtpTHEvent*>(thevtTca[evt.index]);
  Int_t thidx=the->neumc;

  //fill MCEventInfo to do MODBYRS3;
  MCEventInfo ei;
  ei.UseStoredXSec(true);
  NtpStRecord* st=const_cast<NtpStRecord*>(&ntp);//nasty!
  ReweightHelpers::MCEventInfoFilla(&ei,st,thidx);

  //a second one for a sanity check on the values in the NuEvent
  MCEventInfo ei2;
  ei2.UseStoredXSec(true);
  this->MCEventInfoFilla(&ei2,nu);
  this->CompareMCEventInfo(ei,ei2);
  
  //Compute the MODBYRS3 weight
  Double_t generatorweight=1.;
  NuParent* np=0;
  if (nu.useGeneratorReweight){
    MCReweight* mcr=&MCReweight::Instance();
    generatorweight=mcr->ComputeWeight(&ei,np,rwtconfig);
  }
  
  //set the output of the function
  nu.generatorWeight=generatorweight;
}

//......................................................................

void NuZBeamReweight::CalcGeneratorReweight(NuEvent& nu) const
{
  if (!nu.useGeneratorReweight) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
      <<"CalcGeneratorReweight(nu): not calculating Generator reweight"
      <<" since nu.useGeneratorReweight="<<nu.useGeneratorReweight
      <<endl;
    return;
  }

  if (ReleaseType::IsDaikon(nu.releaseType)) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
      <<"Not doing CalcGeneratorReweight(nu) for Daikon, not required"
      <<endl;
    return;
  }

  MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"CalcGeneratorReweight(nu):"<<endl
    <<"  generator:config_name="<<nu.sGeneratorConfigName
    <<"  generator:config_no="<<nu.generatorConfigNo
    <<endl;
  
  //declare the registry and create it only once
  static Registry* rwtconfig=this->CreateNeugenRegistry(nu);

  //make a WeightCalculator and add it to the MCReweight singleton
  this->AddNeugenWeightCalculator();
  
  //fill MCEventInfo
  MCEventInfo ei;
  ei.UseStoredXSec(true);
  //ReweightHelpers::MCEventInfoFilla(&ei,st,thidx);
  this->MCEventInfoFilla(&ei,nu);
  
  //Compute the MODBYRS3 weight
  Double_t generatorweight=1.;
  NuParent* np=0;
  if (nu.useGeneratorReweight){
    MCReweight* mcr=&MCReweight::Instance();
    generatorweight=mcr->ComputeWeight(&ei,np,rwtconfig);
  }
  
  //set the output of the function
  nu.generatorWeight=generatorweight;
}


//......................................................................

double NuZBeamReweight::GetWeightHelium(NuEvent &nu) const
{
  if (ReleaseType::IsDaikon(nu.releaseType) && 
      ReleaseType::GetMCSubVersion(nu.releaseType) > 4) {
    MAXMSG("NuZBeamReweight",Msg::kWarning, 10) 
    << "Asking for He weights for releaseType = " 
    << NuUtilities::PrintRelease(nu.releaseType) 
    << " -- He should already be accounted for.  Returning 1." << endl;
    return 1;
  }
  
  // It doesn't matter which weight calculator we get
  SKZPWeightCalculator *wc = this->GetSKZPWeightCalculatorCustom();

  int inu = nu.inu;
  double E = nu.neuEnMC;
  int det = nu.detector;
  int Ibeam;
  if (nu.beamType == BeamType::kL010z185i_rev) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1) 
    << "Changing beamtype from kL010z185i_rev to kL010z185i for SKZP weighting" << endl;
    Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(BeamType::kL010z185i));
    inu *= -1;
  }
  else {
    Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(nu.beamType));
  } 
  
  
  return wc->GetHeliumWeight(det, Ibeam, inu, E);
}



//......................................................................

bool NuZBeamReweight::EventPreCheck(NuEvent &nu) const
{
  if (nu.simFlag!=SimFlag::kMC || 
      nu.reweightVersion==SKZPWeightCalculator::kUnknown ||
      nu.runPeriod == -1) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
    <<"Not doing ExtractZBeamReweight, reweightVersion="
    << nu.reweightVersion << ", simFlag=" << nu.simFlag
    << ", runPeriod=" << nu.runPeriod << endl;
    
    //set the weights to 1 that would otherwise be set below
    nu.trkEnWeight=1;//no reweighting
    nu.shwEnWeight=1;//no reweighting
    nu.beamWeight=1;//no reweighting
    nu.fluxErrHadProdAfterTune=-1;
    nu.fluxErrTotalErrorPreTune=-1;
    nu.fluxErrTotalErrorAfterTune=-1;
    nu.detectorWeightNMB=1;//no reweighting
    nu.detectorWeightNM=1;//no reweighting
    
    nu.trkEnWeightRunI=1;//no reweighting
    nu.shwEnWeightRunI=1;//no reweighting
    nu.beamWeightRunI=1;//no reweighting
    nu.fluxErrHadProdAfterTuneRunI=-1;
    nu.fluxErrTotalErrorPreTuneRunI=-1;
    nu.fluxErrTotalErrorAfterTuneRunI=-1;
    nu.detectorWeightNMBRunI=1;//no reweighting
    nu.detectorWeightNMRunI=1;//no reweighting
  
    nu.trkEnWeightRunII=1;//no reweighting
    nu.shwEnWeightRunII=1;//no reweighting
    nu.beamWeightRunII=1;//no reweighting
    nu.fluxErrHadProdAfterTuneRunII=-1;
    nu.fluxErrTotalErrorPreTuneRunII=-1;
    nu.fluxErrTotalErrorAfterTuneRunII=-1;
    nu.detectorWeightNMBRunII=1;//no reweighting
    nu.detectorWeightNMRunII=1;//no reweighting
    
    return false;
  }
  
  return true;
}


//......................................................................

void NuZBeamReweight::GetWeights(NuEvent &nu, SKZPWeightCalculator& wc, bool beamOnly) const
{  
  if (nu.runPeriod <= 0) {
    MSG("NuZBeamReweight",Msg::kError) << "Run period not specified in NuEvent -- GetWeights failing." << endl;
    return;
  }
  
  SKZPWeightCalculator::RunPeriod_t rp = SKZPWeightCalculator::kNone;
  if (nu.runPeriod == 1) rp = SKZPWeightCalculator::kRunI;
  if (nu.runPeriod == 2) rp = SKZPWeightCalculator::kRunII;
  if (nu.runPeriod == 3) rp = SKZPWeightCalculator::kRunIII;
  if (nu.runPeriod == 4) {
    rp = SKZPWeightCalculator::kRunI;
    MAXMSG("NuZBeamReweight",Msg::kInfo,1) << "Reweighting RunIV as if it were Run I." << endl;
  }
  
  this->GetWeights(nu, wc, rp, beamOnly); 
}

//......................................................................

void NuZBeamReweight::GetWeights(NuEvent &nu, SKZPWeightCalculator& wc, 
                                 SKZPWeightCalculator::RunPeriod_t rp,
                                 bool beamOnly) const
{
  //fill MCEventInfo
  MCEventInfo ei;
  if (!beamOnly) {
    ei.UseStoredXSec(true);
    this->MCEventInfoFilla(&ei,nu);
  }
  
  //convert the std BeamType into Zarko's convention
  //was 2=le010z185i (Zbeam.h)
  Int_t Ibeam;
  if (nu.beamType == BeamType::kL010z185i_rev) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1) 
    << "Changing beamtype from kL010z185i_rev to kL010z185i for SKZP weighting" << endl;
    Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(BeamType::kL010z185i));
  }
  else {
    Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(nu.beamType));
  }
  
  //copied from SKZPWeightCalculator
  Int_t ccnc=nu.iaction;
  Double_t true_enu=nu.neuEnMC;//p4neu[3]
  Double_t true_eshw=nu.shwEnMC;//p4shw[3]
  Int_t inu=nu.inu; 
  Float_t nsigma=1;
  //FUTURE VERSIONS may use this:
  //parenttype = fi->ptype;
  //parentpz = 1.*fi->pppz;
  //parentpy = 1.*fi->ppdydz*parentpz;
  //parentpx = 1.*fi->ppdxdz*parentpz;
  Int_t parenttype = nu.tptype;
  Double_t parentpz = nu.tpz;
  Double_t parentpy = nu.tpy;
  Double_t parentpx = nu.tpx;
  Double_t pt = sqrt(parentpy*parentpy+parentpx*parentpx);
  
  //create variables that are passed by reference to function to fill
  Double_t new_reco_emu=0;
  Double_t new_reco_eshw=0;
  Double_t beamweight=0;
  Double_t detweightNMB=0;
  Double_t detweightNM=0;
  Float_t trkEnWeight=1;
  Float_t shwEnWeight=1;
  
  if (beamOnly) {
    // Get only the beam weights
    beamweight = wc.GetBeamWeight(nu.detector, Ibeam, parenttype, pt,
                                  parentpz, true_enu, inu, rp);
  }
  else {
    // Get the skzp weights (beam and detector)
    wc.SetSampleSelection("NuMuBar");
    wc.GetWeight(nu.detector,Ibeam,
                 parenttype,pt,parentpz,ccnc,true_enu,inu,
                 nu.trkEn,nu.shwEn,new_reco_emu,new_reco_eshw,
                 beamweight,detweightNMB,rp,true_eshw,&ei);
    
    
    
    MAXMSG("NuZBeamReweight",Msg::kDebug,200)
    <<endl<<"iaction="<<nu.iaction
    <<", nu.trkEn="<<nu.trkEn<<", nu.shwEn="<<nu.shwEn<<endl;
    
    MAXMSG("NuZBeamReweight",Msg::kDebug,200)
    <<"NMB: newEmu="<<new_reco_emu
    <<", newEShw="<<new_reco_eshw
    <<", detNMB="<<detweightNMB
    <<", beamW="<<beamweight<<endl;
    
    //get the skzp weights (beam and detector)
    wc.SetSampleSelection("NuMu");
    wc.GetWeight(nu.detector,Ibeam,
                 parenttype,pt,parentpz,ccnc,true_enu,inu,
                 nu.trkEn,nu.shwEn,new_reco_emu,new_reco_eshw,
                 beamweight,detweightNM,rp,true_eshw,&ei);
    
    MAXMSG("NuZBeamReweight",Msg::kDebug,200)
    <<"NM:  newEmu="<<new_reco_emu
    <<", newEshw="<<new_reco_eshw
    <<",  detNM="<<detweightNM
    <<", beamW="<<beamweight<<endl;
    
    //new for kDetXs
    //double GetWeight(int det, int Ibeam,
    //     int tptype, double pt, double pz,    
    //     double reco_emu, double reco_eshw,
    //     double &new_reco_emu, double &new_reco_eshw,
    //     double &beamweight, double &detweight,
    //     RunPeriod_t rp=kRunI,double true_eshw=-1,MCEventInfo *ei=0);
    
    if (nu.trkEn) trkEnWeight=new_reco_emu/nu.trkEn;
    if (nu.shwEn) shwEnWeight=new_reco_eshw/nu.shwEn;
  }
  
  //get the flux errors
  Float_t fluxErrHadProdAfterTune = 
  wc.GetFluxError(nu.detector,Ibeam,inu,true_enu,
                  SKZPWeightCalculator::kHadProdAfterTune,
                  nsigma);
  Float_t fluxErrTotalErrorPreTune = 
  wc.GetFluxError(nu.detector,Ibeam,inu,true_enu,
                  SKZPWeightCalculator::kTotalErrorPreTune,
                  nsigma);
  Float_t fluxErrTotalErrorAfterTune = 
  wc.GetFluxError(nu.detector,Ibeam,inu,true_enu,
                  SKZPWeightCalculator::kTotalErrorAfterTune,
                  nsigma);
  
  
  
  MAXMSG("NuZBeamReweight",Msg::kDebug,300)
  << "Wgts for run " << (int)rp
  << ", b="   << beamweight 
  << ", dNMB="<< detweightNMB 
  << ", dNM=" << detweightNM
  << ", trk=" << trkEnWeight 
  << ", shw=" << shwEnWeight 
  << endl
  << "energyMC="  << nu.neuEnMC
  << ", iaction=" << nu.iaction
  << ", inu="     << nu.inu
  << endl
  << "fluxErrs: " << " HadProdAfterTune   =" << fluxErrHadProdAfterTune << endl
  << "          " << " TotalErrorPreTune  =" << fluxErrTotalErrorPreTune << endl
  << "          " << " TotalErrorAfterTune=" << fluxErrTotalErrorAfterTune << endl 
  << endl;  
  
  nu.beamWeight                 = beamweight;
  nu.fluxErrHadProdAfterTune    = fluxErrHadProdAfterTune;
  nu.fluxErrTotalErrorPreTune   = fluxErrTotalErrorPreTune;
  nu.fluxErrTotalErrorAfterTune = fluxErrTotalErrorAfterTune;
  
  if (!beamOnly) {
    nu.trkEnWeight                = trkEnWeight;
    nu.shwEnWeight                = shwEnWeight;
    nu.detectorWeightNMB          = detweightNMB;
    nu.detectorWeightNM           = detweightNM;
  }
  
}


//......................................................................

void NuZBeamReweight::GetRFWeights(NuEvent &nu, SKZPWeightCalculator& wc) const
{

  //fill MCEventInfo
  MCEventInfo ei;
  ei.UseStoredXSec(true);
  this->MCEventInfoFilla(&ei,nu);
  //convert the std BeamType into Zarko's convention
  //was 2=le010z185i (Zbeam.h)
  
  Int_t Ibeam;
  if (nu.beamType == BeamType::kL010z185i_rev) {
    MAXMSG("NuZBeamReweight",Msg::kInfo,1) 
    << "Changing beamtype from kL010z185i_rev to kL010z185i for SKZP weighting" << endl;
    Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(BeamType::kL010z185i));
  }
  else {
    Ibeam = BeamType::ToZarko(static_cast<BeamType::BeamType_t>(nu.beamType));
  }

  
  
  //copied from SKZPWeightCalculator
  Int_t ccnc=nu.iaction;
  Double_t true_enu=nu.neuEnMC;//p4neu[3]
  //Double_t true_eshw=nu.shwEnMC;//p4shw[3]
  Int_t inu=nu.inu; 
  Float_t nsigma=1;
  //FUTURE VERSIONS may use this:
  //parenttype = fi->ptype;
  //parentpz = 1.*fi->pppz;
  //parentpy = 1.*fi->ppdydz*parentpz;
  //parentpx = 1.*fi->ppdxdz*parentpz;
  Int_t parenttype = nu.tptype;
  Double_t parentpz = nu.tpz;
  Double_t parentpy = nu.tpy;
  Double_t parentpx = nu.tpx;
  Double_t pt = sqrt(parentpy*parentpy+parentpx*parentpx);

  //create variables that are passed by reference to function to fill
  Double_t new_reco_emu=0;
  Double_t new_reco_eshw=0;
  Double_t beamweight=0;
  Double_t detweightNMB=0;
  Double_t detweightNM=0;
  
  //get the skzp weights (beam and detector)
  wc.SetSampleSelection("NuMuBar");
  wc.GetRFWWeight(nu.detector,Ibeam,
                  parenttype,pt,parentpz,ccnc,true_enu,inu,
                  nu.trkEn,nu.shwEn,new_reco_emu,new_reco_eshw,
                  beamweight,detweightNMB);
  
  MAXMSG("NuZBeamReweight",Msg::kDebug,200)
  <<endl<<"iaction="<<nu.iaction
  <<", nu.trkEn="<<nu.trkEn<<", nu.shwEn="<<nu.shwEn<<endl;
  
  MAXMSG("NuZBeamReweight",Msg::kDebug,200)
  <<"NMB: newEmu="<<new_reco_emu
  <<", newEShw="<<new_reco_eshw
  <<", detNMB="<<detweightNMB
  <<", beamW="<<beamweight<<endl;
  
  //get the skzp weights (beam and detector)
  wc.SetSampleSelection("NuMu");
  wc.GetRFWWeight(nu.detector,Ibeam,
                  parenttype,pt,parentpz,ccnc,true_enu,inu,
                  nu.trkEn,nu.shwEn,new_reco_emu,new_reco_eshw,
                  beamweight,detweightNM);
  
  MAXMSG("NuZBeamReweight",Msg::kDebug,200)
  <<"NM:  newEmu="<<new_reco_emu
  <<", newEshw="<<new_reco_eshw
  <<",  detNM="<<detweightNM
  <<", beamW="<<beamweight<<endl;
  
  //new for kDetXs
  //double GetWeight(int det, int Ibeam,
  //       int tptype, double pt, double pz,    
  //       double reco_emu, double reco_eshw,
  //       double &new_reco_emu, double &new_reco_eshw,
  //       double &beamweight, double &detweight,
  //       RunPeriod_t rp=kRunI,double true_eshw=-1,MCEventInfo *ei=0);
  
  Float_t trkEnWeight=1;
  Float_t shwEnWeight=1;
  if (nu.trkEn) trkEnWeight=new_reco_emu/nu.trkEn;
  if (nu.shwEn) shwEnWeight=new_reco_eshw/nu.shwEn;
  
  //get the flux errors
  Float_t fluxErrHadProdAfterTune = 
  wc.GetFluxError(nu.detector,Ibeam,inu,true_enu,
                  SKZPWeightCalculator::kHadProdAfterTune,
                  nsigma);
  Float_t fluxErrTotalErrorPreTune = 
  wc.GetFluxError(nu.detector,Ibeam,inu,true_enu,
                  SKZPWeightCalculator::kTotalErrorPreTune,
                  nsigma);
  Float_t fluxErrTotalErrorAfterTune = 
  wc.GetFluxError(nu.detector,Ibeam,inu,true_enu,
                  SKZPWeightCalculator::kTotalErrorAfterTune,
                  nsigma);
  
  
  
  MAXMSG("NuZBeamReweight",Msg::kDebug,300)
  << "Wgts for run " << "I+II Avg"
  << ", b="   << beamweight 
  << ", dNMB="<< detweightNMB 
  << ", dNM=" << detweightNM
  << ", trk=" << trkEnWeight 
  << ", shw=" << shwEnWeight 
  << endl
  << "energyMC="  << nu.neuEnMC
  << ", iaction=" << nu.iaction
  << ", inu="     << nu.inu
  << endl
  << "fluxErrs: " << " HadProdAfterTune   =" << fluxErrHadProdAfterTune << endl
  << "          " << " TotalErrorPreTune  =" << fluxErrTotalErrorPreTune << endl
  << "          " << " TotalErrorAfterTune=" << fluxErrTotalErrorAfterTune << endl 
  << endl;  
  
  
  nu.trkEnWeight                = trkEnWeight;
  nu.shwEnWeight                = shwEnWeight;
  nu.beamWeight                 = beamweight;
  nu.fluxErrHadProdAfterTune    = fluxErrHadProdAfterTune;
  nu.fluxErrTotalErrorPreTune   = fluxErrTotalErrorPreTune;
  nu.fluxErrTotalErrorAfterTune = fluxErrTotalErrorAfterTune;
  nu.detectorWeightNMB          = detweightNMB;
  nu.detectorWeightNM           = detweightNM;
  
}




//......................................................................

void NuZBeamReweight::AddNeugenWeightCalculator() const
{
  Bool_t firstTime=true;
  if (firstTime) {
    firstTime=false;
    MAXMSG("NuZBeamReweight",Msg::kInfo,1)
      <<"Adding NeugenWeightCalculator to MCReweight singleton"<<endl;
    static NeugenWeightCalculator *n=new NeugenWeightCalculator();
    MCReweight* mcr=&MCReweight::Instance();    
    mcr->AddWeightCalculator(n);
  }
}

//......................................................................

Registry* NuZBeamReweight::CreateNeugenRegistry(const NuEvent& nu) const
{
  static Registry* rwtconfig=0;
  
  //create it the first time
  if (rwtconfig==0){
    
    //fill the registry to tell MCReweight to do MODBYRS3
    rwtconfig=new Registry();
    rwtconfig->UnLockValues();
    rwtconfig->UnLockKeys();
    rwtconfig->Clear();
    rwtconfig->Set("neugen:config_name",nu.sGeneratorConfigName.c_str());
    rwtconfig->Set("neugen:config_no",nu.generatorConfigNo);
    rwtconfig->LockValues();
    rwtconfig->LockKeys();
  }
  return rwtconfig;
}

//......................................................................

void NuZBeamReweight::CompareMCEventInfo(const MCEventInfo& ei,
                                         const MCEventInfo& ei2) const
{
  if (ei.initial_state!=ei2.initial_state ||
      ei.nucleus!=ei2.nucleus ||
      ei.had_fs!=ei2.had_fs){
    MAXMSG("NuZBeamReweight",Msg::kError,300)
      <<"Ahhh, CompareMCEventInfo difference found:"<<endl
      <<"1st: initial_state="<<ei.initial_state
      <<", nucleus="<<ei.nucleus
      <<", had_fs="<<ei.had_fs
      <<endl
      <<"2nd: initial_state="<<ei2.initial_state
      <<", nucleus="<<ei2.nucleus
      <<", had_fs="<<ei2.had_fs
      <<endl
      <<"NOTE: This discrepancy has to be understood if we"
      <<" are to trust the generator reweighting"
      <<endl;
  }
}

//......................................................................

void NuZBeamReweight::MCEventInfoFilla(MCEventInfo* ei,
                                       const NuEvent& nu) const
{

  ei->nuE=nu.neuEnMC;
  ei->nuPx=nu.neuPxMC;
  ei->nuPy=nu.neuPyMC;
  ei->nuPz=nu.neuPzMC;
  
  ei->tarE=nu.tgtEnMC;
  ei->tarPx=nu.tgtPxMC;
  ei->tarPy=nu.tgtPyMC;
  ei->tarPz=nu.tgtPzMC;
  
  ei->y=nu.yMC;
  ei->x=nu.xMC;
  ei->q2=nu.q2MC;
  ei->w2=nu.w2MC;
  
  ei->iaction=nu.iaction;
  ei->inu=nu.inu;
  ei->iresonance=nu.iresonance;
  ei->initial_state=nu.initialStateMC;
  
  ei->nucleus=nu.nucleusMC;
  
  ei->had_fs=nu.hadronicFinalStateMC;


  MAXMSG("NuZBeamReweight",Msg::kDebug,100)
    <<endl<<endl
    <<", ei->nuE="<<ei->nuE
    <<", ei->nuPx="<<ei->nuPx
    <<", ei->nuPy="<<ei->nuPy
    <<", ei->nuPz="<<ei->nuPz
    <<endl
    <<", ei->tarE="<<ei->tarE
    <<", ei->tarPx="<<ei->tarPx
    <<", ei->tarPy="<<ei->tarPy
    <<", ei->tarPz="<<ei->tarPz
    <<endl
    <<", ei->y="<<ei->y<<endl
    <<", ei->x="<<ei->x<<endl
    <<", ei->q2="<<ei->q2<<endl
    <<", ei->w2="<<ei->w2<<endl
    <<endl
    <<", ei->iaction="<<ei->iaction<<endl
    <<", ei->inu="<<ei->inu<<endl
    <<", ei->iresonance="<<ei->iresonance<<endl
    <<", ei->initial_state="<<ei->initial_state<<endl
    <<", ei->nucleus="<<ei->nucleus<<endl
    <<", ei->had_fs="<<ei->had_fs<<endl
    <<endl;
}

//......................................................................

const Char_t* NuZBeamReweight::
AsString(SKZPWeightCalculator::SKZPConfig_t c) const
{
  // return char string in cannonical offline form
  switch (c) {
    case SKZPWeightCalculator::kPRL:     return "PRL";
    case SKZPWeightCalculator::kBoston:  return "Boston";
    case SKZPWeightCalculator::kPiMinus: return "PiMinus_CedarDaikon";
    case SKZPWeightCalculator::kDetXs:   return "DetXs";
    case SKZPWeightCalculator::kDogwood1_Daikon07: return "Dogwood1_Daikon07";
    case SKZPWeightCalculator::kDogwood1_Daikon07_v2: return "Dogwood1_Daikon07_v2";
      case SKZPWeightCalculator::kUnknown: return "Unknown"; break;
    default:         return "!?Bad Value Passed?!"; break;
  }
}

//......................................................................

---