NuShiftableUnbinnedSpectrum.cxx

Go to the documentation of this file.

#include "TChain.h"
#include "TChainElement.h"
#include "TFile.h"
#include "TObjArray.h"

#include "MessageService/MsgService.h"
#include "NtupleUtils/NuFCEvent.h"
#include "NtupleUtils/NuMMParameters.h"
#include "NtupleUtils/NuShiftableUnbinnedSpectrum.h"
#include "NtupleUtils/NuTreeWrapper.h"

ClassImp(NuShiftableUnbinnedSpectrum)

CVSID("$Id: NuShiftableUnbinnedSpectrum.cxx,v 1.1 2008/10/17 11:47:03 evansj Exp $");

//____________________________________________________________________72
NuShiftableUnbinnedSpectrum::NuShiftableUnbinnedSpectrum()
  : NuShiftableSpectrum(),
    fCharge(1),
    fBinningScheme(NuBinningScheme::kNuMuBar0325Std2)
{
}

//____________________________________________________________________72
NuShiftableUnbinnedSpectrum
::NuShiftableUnbinnedSpectrum(const std::string filename)
  : NuShiftableSpectrum(),
    fCharge(1)
{
  /*
  const NuTreeWrapper dst(filename);
  this->PoT(dst.GetPoT());
  Int_t numEntries = dst.GetEntries();
  for (Int_t entry=0; entry<numEntries; ++entry){
    NuEvent event = dst.GetInfoObject(entry);
    if (event.charge == fCharge || 0==fCharge){
      NuMiniEvent miniEvent;
      miniEvent.recoTrackEnergy = event.trkEn;
      miniEvent.recoNuEnergy = event.energy;
      fEventList.push_back(miniEvent);
    }
  }
  */
  TChain fcDST("FCTree");
  fcDST.Add(filename.c_str());
  NuFCEvent* fcEvent = 0;
  fcDST.SetBranchAddress("NuFCEvent",&fcEvent);
  Int_t numEntries = fcDST.GetEntries();
  for (Int_t entry=0; entry<numEntries; ++entry){
    fcDST.GetEvent(entry);
    if (this->Charge() == fcEvent->charge || 0==this->Charge()){
      NuMiniEvent miniEvent;
      miniEvent.recoTrackEnergy = fcEvent->trkEn;
      miniEvent.recoNuEnergy = fcEvent->energy;
      fEventList.push_back(miniEvent);
    }
  }
  
  TObjArray* fileList = fcDST.GetListOfFiles();
  Double_t totPoT = 0.0;
  for (Int_t counter=0; counter<fileList->GetEntries(); ++counter){
    TChainElement* chainElement =
      dynamic_cast<TChainElement*> (fileList->At(counter));
    if (chainElement){
      string fileName = chainElement->GetTitle();
      TFile f(fileName.c_str(),"READ");
      TH1F* hTotalPot = (TH1F*) f.Get("hTotalPot");
      if (hTotalPot){
        totPoT += hTotalPot->Integral();
      }
      f.Close();
    }
  }
  this->PoT(totPoT);
  MSG("NuTreeWrapper.cxx",Msg::kInfo)
    << "PoT: " << totPoT << endl;
}

//____________________________________________________________________72
const TH1D NuShiftableUnbinnedSpectrum::Spectrum(const NuMMParameters& pars) const
{
  return this->Spectrum(pars.ShwEnScale());
}

//____________________________________________________________________72
const TH1D NuShiftableUnbinnedSpectrum::Spectrum(const Double_t shift) const
{
  const NuUtilities utils;
  const vector<Double_t> recoBins = utils.RecoBins(fBinningScheme);
  const Int_t numRecoBins = recoBins.size()-1;
  Float_t *recoBinsArray;
  recoBinsArray=new Float_t[numRecoBins+1];
  {
    Int_t i=0;
    for (vector<Double_t>::const_iterator itBin = recoBins.begin();
         itBin != recoBins.end();
         ++itBin, ++i){
      recoBinsArray[i] = *itBin;
    }
  }
  TH1D outSpectrum("RecoEnergy_FD","",
                   numRecoBins,recoBinsArray);

  for (vector<NuMiniEvent>::const_iterator it = fEventList.begin();
       it != fEventList.end();
       ++it){
    Double_t energy = it->recoNuEnergy;
    Double_t trkEn = it->recoTrackEnergy;
    Double_t shiftedEnergy = energy + shift*trkEn;
    outSpectrum.Fill(shiftedEnergy);
  }

  return outSpectrum;
}

---