NCEnergyBin.cxx

Go to the documentation of this file.

//$Id: NCEnergyBin.cxx,v 1.41 2009/08/10 17:56:47 bckhouse Exp $
//
//NCEnergyBin.cxx
//
//B. Rebel 3/2005

#include "NCUtils/Extrapolation/NCEnergyBin.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"
#include "MessageService/MsgService.h"
#include "TMath.h"

#include <cassert>

CVSID("$Id: NCEnergyBin.cxx,v 1.41 2009/08/10 17:56:47 bckhouse Exp $");

//-------------------------------------------------------------------
NCEnergyBin::NCEnergyBin(double binCenter, double binWidth,
                         NCType::EEventType nccc) :
  fSignal(0.),
  fBinCenter(binCenter),
  fBinWidth(binWidth),
  fBinType(nccc)
{
  fNCBin = (nccc == NCType::kNC);

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

//----------------------------------------------------------------------
void NCEnergyBin::AddEventToBin(const ANtpRecoInfo* recoInfo,
                                const ANtpTruthInfoBeam* truthInfo,
                                NCType::EFileType fileType)
{
  const double weight = recoInfo->weight;
  // (PAR) We set track energy to zero for NC selected events, because
  // the total energy returned elsewhere is the track+shower.
  double trackEnergy=fNCBin ? 0 : recoInfo->muEnergy;

  if(!truthInfo){
    // data is easy...
    fSignal += weight;
    fDataEnergies.push_back(recoInfo->nuEnergy);
    fDataWeights.push_back(weight);

    return;
  }

  const bool electron = TMath::Abs(truthInfo->nuFlavor) == 12;

  //see if this is background event - either NC for CC bins or
  //nu_mu CC for NC bins
  if(fileType == NCType::kBeamFile && !electron &&
     ((fNCBin && truthInfo->interactionType == NCType::kCC) ||
      (!fNCBin && truthInfo->interactionType == NCType::kNC))){
    fMCBeamWeightBG.push_back(weight);
    fMCRecoShowerEBG.push_back(recoInfo->showerEnergy);
    fMCRecoMuonEBG.push_back(trackEnergy);
    fMCFlavorBG.push_back(truthInfo->nuFlavor);
    fMCEnergyBG.push_back(truthInfo->nuEnergy);
    return;
  }

  if(fileType == NCType::kTauFile &&
     truthInfo->interactionType == NCType::kCC){
    fMCBeamWeightNuTau.push_back(weight);
    fMCRecoShowerENuTau.push_back(recoInfo->showerEnergy);
    fMCRecoMuonENuTau.push_back(trackEnergy);
    fMCFlavorNuTau.push_back(truthInfo->nuFlavor);
    fMCEnergyNuTau.push_back(truthInfo->nuEnergy);
    return;
  }

  if(electron &&
     fileType == NCType::kBeamFile){
    fMCBeamWeightBeamNuE.push_back(weight);
    fMCRecoShowerEBeamNuE.push_back(recoInfo->showerEnergy);
    fMCRecoMuonEBeamNuE.push_back(trackEnergy);
    fMCFlavorBeamNuE.push_back(truthInfo->nuFlavor);
    fMCEnergyBeamNuE.push_back(truthInfo->nuEnergy);
    return;
  }

  if(electron &&
     fileType == NCType::kElectronFile &&
     truthInfo->interactionType == NCType::kCC){
    fMCBeamWeightOscNuE.push_back(weight);
    fMCRecoShowerEOscNuE.push_back(recoInfo->showerEnergy);
    fMCRecoMuonEOscNuE.push_back(trackEnergy);
    fMCFlavorOscNuE.push_back(truthInfo->nuFlavor);
    fMCEnergyOscNuE.push_back(truthInfo->nuEnergy);
    return;
  }

  fMCBeamWeight.push_back(weight);
  fMCRecoShowerE.push_back(recoInfo->showerEnergy);
  fMCRecoMuonE.push_back(trackEnergy);
  fMCFlavor.push_back(truthInfo->nuFlavor);
  fMCEnergy.push_back(truthInfo->nuEnergy);
}

//----------------------------------------------------------------------
int NCEnergyBin::GetDataVectorSize() const
{
  return int(fDataEnergies.size());
}

//----------------------------------------------------------------------
void NCEnergyBin::GetDataInformation(double &energy, double &weight, int eventNum) const
{
  assert(eventNum < int(fDataEnergies.size()));

  energy = fDataEnergies[eventNum];
  weight = fDataWeights[eventNum];
}

//----------------------------------------------------------------------
double NCEnergyBin::GetSignal() const
{
  return fSignal;
}

//----------------------------------------------------------------------
double NCEnergyBin::GetBinCentralValue() const
{
  return fBinCenter;
}

//----------------------------------------------------------------------
double NCEnergyBin::GetBinWidth() const
{
  return fBinWidth;
}

//----------------------------------------------------------------------
int NCEnergyBin::GetBinType() const
{
  return fBinType;
}

//----------------------------------------------------------------------
void NCEnergyBin::SetSignal(double signal)
{
  fSignal = signal;
}

//----------------------------------------------------------------------
int NCEnergyBin::GetMCSignalVectorSize() const
{
  return int(fMCBeamWeight.size());
}

//----------------------------------------------------------------------
void NCEnergyBin::GetMCInformation(double &trueEnergy,
                                   double &recoShowerE,
                                   double &recoMuonE,
                                   double &weight,
                                   int    &flavor,
                                   int    eventNum) const
{
  assert(eventNum < int(fMCEnergy.size()));

  trueEnergy = fMCEnergy[eventNum];
  recoShowerE = fMCRecoShowerE[eventNum];
  recoMuonE = fMCRecoMuonE[eventNum];
  weight = fMCBeamWeight[eventNum];
  flavor = fMCFlavor[eventNum];
}

//----------------------------------------------------------------------
int NCEnergyBin::GetMCBackgroundVectorSize() const
{
  return int(fMCBeamWeightBG.size());
}

//----------------------------------------------------------------------
void NCEnergyBin::GetMCBackgroundInformation(double &trueEnergy,
                                             double &recoShowerE,
                                             double &recoMuonE,
                                             double &weight,
                                             int    &flavor,
                                             int    eventNum) const
{

  assert(eventNum < int(fMCEnergyBG.size()));

  trueEnergy = fMCEnergyBG[eventNum];
  recoShowerE = fMCRecoShowerEBG[eventNum];
  recoMuonE = fMCRecoMuonEBG[eventNum];
  weight = fMCBeamWeightBG[eventNum];
  flavor = fMCFlavorBG[eventNum];
}

//----------------------------------------------------------------------
int NCEnergyBin::GetMCNuTauVectorSize() const
{
  return int(fMCBeamWeightNuTau.size());
}

//----------------------------------------------------------------------
void NCEnergyBin::GetMCNuTauInformation(double &trueEnergy,
                                      double &recoShowerE,
                                      double &recoMuonE,
                                      double &weight,
                                      int    &flavor,
                                      int    eventNum) const
{

  assert(eventNum < int(fMCEnergyNuTau.size()));

  trueEnergy = fMCEnergyNuTau[eventNum];
  recoShowerE = fMCRecoShowerENuTau[eventNum];
  recoMuonE = fMCRecoMuonENuTau[eventNum];
  weight = fMCBeamWeightNuTau[eventNum];
  flavor = fMCFlavorNuTau[eventNum];
}

//----------------------------------------------------------------------
int NCEnergyBin::GetMCOscNuEVectorSize() const
{
  return int(fMCBeamWeightOscNuE.size());
}

//----------------------------------------------------------------------
void NCEnergyBin::GetMCOscNuEInformation(double &trueEnergy,
                                         double &recoShowerE,
                                         double &recoMuonE,
                                         double &weight,
                                         int    &flavor,
                                         int    eventNum) const
{
  assert(eventNum < int(fMCEnergyOscNuE.size()));

  trueEnergy = fMCEnergyOscNuE[eventNum];
  recoShowerE = fMCRecoShowerEOscNuE[eventNum];
  recoMuonE = fMCRecoMuonEOscNuE[eventNum];
  weight = fMCBeamWeightOscNuE[eventNum];
  flavor = fMCFlavorOscNuE[eventNum];
}

//----------------------------------------------------------------------
int NCEnergyBin::GetMCBeamNuEVectorSize() const
{
  return int(fMCBeamWeightBeamNuE.size());
}

//----------------------------------------------------------------------
void NCEnergyBin::GetMCBeamNuEInformation(double &trueEnergy,
                                          double &recoShowerE,
                                          double &recoMuonE,
                                          double &weight,
                                          int    &flavor,
                                          int    eventNum) const
{
  assert(eventNum < int(fMCEnergyBeamNuE.size()));

  trueEnergy = fMCEnergyBeamNuE[eventNum];
  recoShowerE = fMCRecoShowerEBeamNuE[eventNum];
  recoMuonE = fMCRecoMuonEBeamNuE[eventNum];
  weight = fMCBeamWeightBeamNuE[eventNum];
  flavor = fMCFlavorBeamNuE[eventNum];
}

//----------------------------------------------------------------------
void NCEnergyBin::Add(NCEnergyBin *bin)
{

  //make sure you are adding the correct bin to this one
  if(bin->GetBinType() == fBinType
     && bin->GetBinCentralValue() == fBinCenter){

    double energy = 0.;
    double weight = 0.;
    double recoShowerE = 0.;
    double recoMuonE = 0.;
    int flavor = 0;

    int size = bin->GetDataVectorSize();
    for(int i = 0; i < size; ++i){
      bin->GetDataInformation(energy, weight, i);
      fDataEnergies.push_back(energy);
      fDataWeights.push_back(weight);
      fSignal += weight;
    }//end data

    size = bin->GetMCSignalVectorSize();
    for(int j = 0; j < size; ++j){
      bin->GetMCInformation(energy, recoShowerE, recoMuonE, weight, flavor, j);
      fMCBeamWeight.push_back(weight);
      fMCRecoShowerE.push_back(recoShowerE);
      fMCRecoMuonE.push_back(recoMuonE);
      fMCFlavor.push_back(flavor);
      fMCEnergy.push_back(energy);

    }//end signal

    size = bin->GetMCBackgroundVectorSize();
    for(int j = 0; j < size; ++j){
      bin->GetMCBackgroundInformation(energy, recoShowerE, recoMuonE, weight, flavor, j);
      fMCBeamWeightBG.push_back(weight);
      fMCRecoShowerEBG.push_back(recoShowerE);
      fMCRecoMuonEBG.push_back(recoMuonE);
      fMCFlavorBG.push_back(flavor);
      fMCEnergyBG.push_back(energy);

    }//end background

    size = bin->GetMCNuTauVectorSize();
    for(int j = 0; j < size; ++j){
      bin->GetMCNuTauInformation(energy, recoShowerE, recoMuonE, weight, flavor, j);
      fMCBeamWeightNuTau.push_back(weight);
      fMCRecoShowerENuTau.push_back(recoShowerE);
      fMCRecoMuonENuTau.push_back(recoMuonE);
      fMCFlavorNuTau.push_back(flavor);
      fMCEnergyNuTau.push_back(energy);
    }//end taus

    size = bin->GetMCOscNuEVectorSize();
    for(int j = 0; j < size; ++j){
      bin->GetMCOscNuEInformation(energy, recoShowerE, recoMuonE, weight, flavor, j);
      fMCBeamWeightOscNuE.push_back(weight);
      fMCRecoShowerEOscNuE.push_back(recoShowerE);
      fMCRecoMuonEOscNuE.push_back(recoMuonE);
      fMCFlavorOscNuE.push_back(flavor);
      fMCEnergyOscNuE.push_back(energy);
    }//end electrons

    size = bin->GetMCBeamNuEVectorSize();
    for(int j = 0; j < size; ++j){
      bin->GetMCBeamNuEInformation(energy, recoShowerE, recoMuonE, weight, flavor, j);
      fMCBeamWeightBeamNuE.push_back(weight);
      fMCRecoShowerEBeamNuE.push_back(recoShowerE);
      fMCRecoMuonEBeamNuE.push_back(recoMuonE);
      fMCFlavorBeamNuE.push_back(flavor);
      fMCEnergyBeamNuE.push_back(energy);
    }//end electrons

  }
  else
    MSG("NCEnergyBin", Msg::kWarning) << "trying to add incompatible bins "
                                      << " " << fBinType << "/" << bin->GetBinType()
                                      << " " << fBinCenter << "/" << bin->GetBinCentralValue()
                                      << endl;
}


template<class T> void ZeroVector(vector<T>& vec)
{
  for(unsigned int n = 0; n < vec.size(); ++n) vec[n] = 0;
}


void NCEnergyBin::Reset(bool data, bool mc)
{
  // TODO - I'm by no means sure this is the right set of things to reset

  if(data){
    fSignal = 0;

    ZeroVector(fDataEnergies);
    ZeroVector(fDataWeights);
  }

  if(mc){
    ZeroVector(fMCEnergy);
    ZeroVector(fMCEnergyBG);
    ZeroVector(fMCEnergyNuTau);
    ZeroVector(fMCEnergyOscNuE);
    ZeroVector(fMCEnergyBeamNuE);
    ZeroVector(fMCBeamWeight);
    ZeroVector(fMCBeamWeightBG);
    ZeroVector(fMCBeamWeightNuTau);
    ZeroVector(fMCBeamWeightOscNuE);
    ZeroVector(fMCBeamWeightBeamNuE);
    ZeroVector(fMCRecoShowerE);
    ZeroVector(fMCRecoMuonE);
    ZeroVector(fMCFlavor);
    ZeroVector(fMCRecoShowerEBG);
    ZeroVector(fMCRecoMuonEBG);
    ZeroVector(fMCFlavorBG);
    ZeroVector(fMCRecoShowerENuTau);
    ZeroVector(fMCRecoMuonENuTau);
    ZeroVector(fMCFlavorNuTau);
    ZeroVector(fMCRecoShowerEOscNuE);
    ZeroVector(fMCRecoMuonEOscNuE);
    ZeroVector(fMCFlavorOscNuE);
    ZeroVector(fMCRecoShowerEBeamNuE);
    ZeroVector(fMCRecoMuonEBeamNuE);
    ZeroVector(fMCFlavorBeamNuE);
  }
}

---