NCExtrapolationNone.cxx

Go to the documentation of this file.

//
// NCExtrapolationNone.cxx
//
// C. Backhouse 8/2008
//
// $Id: NCExtrapolationNone.cxx,v 1.27 2009/09/14 13:49:53 tinti Exp $

#include "NCUtils/Extrapolation/NCExtrapolationNone.h"
#include "NCUtils/Extrapolation/NCEnergyBin.h"

#include "MessageService/MsgService.h"

#include "AnalysisNtuples/ANtpAnalysisInfo.h"
#include "AnalysisNtuples/ANtpBeamInfo.h"
#include "AnalysisNtuples/ANtpRecoInfo.h"
#include "AnalysisNtuples/ANtpHeaderInfo.h"
#include "AnalysisNtuples/ANtpTruthInfoBeam.h"

#include "TDirectory.h"
#include "TFile.h"

#include <cassert>

CVSID("$Id: NCExtrapolationNone.cxx,v 1.27 2009/09/14 13:49:53 tinti Exp $");

#include "NCExtrapolationModule.h"

REGISTER_NCEXTRAPOLATION(NCExtrapolationNone, None)

//......................................................................
NCExtrapolationNone::NCExtrapolationNone()
{
  nCalls = 0;

  fEventsTree = new TTree("sel_evts", "sel_evts");

  fEventsTree->Branch("header.",   "ANtpHeaderInfo",   (void**)0);
  fEventsTree->Branch("beam.",     "ANtpBeamInfo",     (void**)0);
  fEventsTree->Branch("reco.",     "ANtpRecoInfo",     (void**)0);
  fEventsTree->Branch("analysis.", "ANtpAnalysisInfo", (void**)0);
}
//----------------------------------------------------------------------


//......................................................................
NCExtrapolationNone::~NCExtrapolationNone()
{
  MSG("NCExtrapolationNone", Msg::kInfo) << nCalls << " calls in total to NCExtrapolationNone" << endl;

  if(fEventsTree) delete fEventsTree;
}


const Registry& NCExtrapolationNone::DefaultConfig()
{
  static Registry r;

  r.UnLockValues();

  r.Set("SelectedEventsFilename", "/tmp/sel_evts.root");

  r.LockValues();
  return r;
}

void NCExtrapolationNone::Prepare(const Registry& r)
{
  NCExtrapolation::Prepare(r);

  const char* tmps;
  if(r.Get("SelectedEventsFilename", tmps)) fSelEvtsFname = tmps;
}
//----------------------------------------------------------------------



//----------------------------------------------------------------------
void NCExtrapolationNone::AddEvent(NCEventInfo info,
                                   bool useMCAsData,
                                   NCType::EFileType fileType,
                                   NCBeam::Info beamInfo)
{
  // Chain up to base-class's implementation so the NCBeam objects get filled
  NCExtrapolation::AddEvent(info, useMCAsData, fileType, beamInfo);

  // Then go ahead and put the events in the selected events tree too

  if(info.reco->inFiducialVolume < 1 ||
     (info.header->detector == int(Detector::kNear) &&
      info.reco->isSimpleCutsClean < 1))
     return;

  if(info.header->dataType != int(SimFlag::kData)) return;

  assert(fEventsTree);
  fEventsTree->SetBranchAddress("header.", &info.header);
  fEventsTree->SetBranchAddress("beam.", &info.beam);
  fEventsTree->SetBranchAddress("reco.", &info.reco);
  fEventsTree->SetBranchAddress("analysis.", &info.analysis);

  fEventsTree->Fill();
}

//----------------------------------------------------------------------
void NCExtrapolationNone::FillMCSpectrum(TH1* exp,
                                         const NCEnergyBin* bin,
                                         InfoFunction func,
                                         int size,
                                         const NC::OscProb::OscPars* pars,
                                         NCType::EEventType nccc,
                                         NCType::EOscMode oscmode)
{
  for(int e = 0; e < size; ++e){
    double trueE, showerE, trackE, weight;
    int flavour;

    // Call the member function specified by func
    (bin->*func)(trueE, showerE, trackE, weight, flavour, e);

    const double oscweight = pars->TransitionProbability(oscmode,
                                                         nccc,
                                                         NCType::kBaseLineFar,
                                                         trueE);

     exp->Fill(showerE+trackE, weight*oscweight);
  }
}


//----------------------------------------------------------------------
void NCExtrapolationNone::
FindSpectraForPars(const NC::OscProb::OscPars* oscPars,
                   const NC::SystPars& systPars,
                   vector<NCBeam::Info> /*beamsToUse*/,
                   vector<TH1*>& expvec,
                   vector<TH1*>& obsvec)
{
  ++nCalls;

  const double norm = systPars.NormScale();

  static TH1D* obsCache[2] = {0, 0};


  // For the moment just hardcode the beam here
  NCBeam* beam = GetBeam(Detector::kFar,
                         NCBeam::Info(BeamType::kL010z185i,
                                      NC::RunUtil::kRunAll));

  for(int nccci = 0; nccci < 2; ++nccci){
    NCType::EEventType nccc;
    if(nccci == 0) nccc = NCType::kNC; else nccc = NCType::kCC;

    if(nccc == NCType::kNC && !fUseNC) continue;
    if(nccc == NCType::kCC && !fUseCC) continue;

    if(!obsCache[nccci]){
      obsCache[nccci] = new TH1D("", "", kNumEnergyBinsFar,
                                 kEnergyBinsFar);
      for(int b = 0; b < beam->GetNumberEnergyBins(nccc); ++b){
        const NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);

        const int dataSize = bin->GetDataVectorSize();

        for(int e = 0; e < dataSize; ++e){
          double energy, weight;
          bin->GetDataInformation(energy, weight, e);

          obsCache[nccci]->Fill(energy, weight);
        }
      } // end for b
    } // end if no cache

    TH1D* obs = obsCache[nccci];

    TH1D* exp = new TH1D("", "", kNumEnergyBinsFar,
                         kEnergyBinsFar);

    for(int b = 0; b < beam->GetNumberEnergyBins(nccc); ++b){
      const NCEnergyBin* bin = beam->GetEnergyBin(b, nccc);

      FillMCSpectrum(exp, bin, &NCEnergyBin::GetMCInformation,
                     bin->GetMCSignalVectorSize(),
                     oscPars, nccc, NCType::kNuMuToNuMu);

      const NCType::EEventType bgtype = (nccc == NCType::kNC) ? NCType::kCC : NCType::kNC;

      FillMCSpectrum(exp, bin, &NCEnergyBin::GetMCBackgroundInformation,
                     bin->GetMCBackgroundVectorSize(),
                     oscPars, bgtype, NCType::kNuMuToNuMu);

      FillMCSpectrum(exp, bin, &NCEnergyBin::GetMCNuTauInformation,
                     bin->GetMCNuTauVectorSize(),
                     oscPars, nccc, NCType::kNuMuToNuTau);

      FillMCSpectrum(exp, bin, &NCEnergyBin::GetMCBeamNuEInformation,
                     bin->GetMCBeamNuEVectorSize(),
                     oscPars, nccc, NCType::kNuEToNuE);

      // Not handling nu_e appearance

      exp->Scale(norm);

    } // end for b

    expvec.push_back(exp);
    obsvec.push_back(obs);
  } // end for nccci

  assert(expvec.size() == obsvec.size());
}


//----------------------------------------------------------------------
void NCExtrapolationNone::CleanupSpectra(vector<TH1*> expvec,
                                         vector<TH1*> /*obsvec*/)
{
  for(unsigned int n = 0; n < expvec.size(); ++n) delete expvec[n];
}


//----------------------------------------------------------------------
void NCExtrapolationNone::
WriteResources(const NC::OscProb::OscPars* trueOscPars)
{
  NCExtrapolation::WriteResources(trueOscPars);

  TDirectory* fileDir = gDirectory;

  TFile f(fSelEvtsFname, "RECREATE");
  fEventsTree->Write();
  f.Close();

  fileDir->cd();
}

---