NCExtrapolationPID.cxx

Go to the documentation of this file.

#include "NCUtils/Extrapolation/NCExtrapolationPID.h"
#include "NCUtils/Extrapolation/PIDSpectrum.h"
#include "NCUtils/NCUtility.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 "Conventions/SimFlag.h"

#include "NCExtrapolationModule.h"
#include "NCSpectrumInterpolator.h"

#include "TH1.h"
#include "TFile.h"
#include "TList.h"
#include "TKey.h"
#include "TCanvas.h"

REGISTER_NCEXTRAPOLATION(NCExtrapolationPID, PID)

#include <iostream>

CVSID("$Id: NCExtrapolationPID.cxx,v 1.79 2009/12/18 01:30:21 rodriges Exp $");

ClassImp(NCExtrapolationPID)

//==============================================================================

NCExtrapolationPID::NCExtrapolationPID()
:  fUseCCEnergy(true), fUseSelectionEnergy(true), fTrueBinFactor(5), fNumPIDBins(20),
  fPIDmin(0.2), fPIDmax(1.1),
  fMCInFilename(""), fMCOutFilename(""),
  fDoFarNearCorrection(true), fEmulateFarNear(false)
{
}

//==============================================================================

NCExtrapolationPID::~NCExtrapolationPID()
{
  for (SpecMap::iterator it=fSpectra.begin(); it!=fSpectra.end(); ++it)
    delete it->second;
}


//==============================================================================

void NCExtrapolationPID::AddEvent(NCEventInfo         info,
                                  bool                useMCAsData,
                                  NCType::EFileType   fileType,
                                  NCBeam::Info        beamInfo)
{
  // This call fills the NCBeam objects with the event. We don't need it
  // for the fit, but it does make some pretty plots at the end.
  NCExtrapolation::AddEvent(info, useMCAsData, fileType, beamInfo);

  // Some events don't get into the NC or the CC spectrum - exclude
  // them if we're emulating the FarNear method
  if (fEmulateFarNear &&
      !(info.analysis->isNC || info.analysis->isCC)) return;

  // Add the necessary spectrum for this beam if it's not already there
  if (fSpectra.find(beamInfo)==fSpectra.end()) {
    // If this happened when we got our MC spectra from disk,
    // something is probably wrong
    assert(fMCInFilename=="");

    string name = string("PID_") + beamInfo.GetDescription().Data();
    fSpectra[beamInfo]=new PIDSpectrum(name, "title",
                                       fNumPIDBins, fPIDmin, fPIDmax, // PID bins
                                       // fNumPIDBins will be -1 if emulating FarNear
                                       fTrueBinFactor, beamInfo,
                                       info.analysis->separationParameterCut);
  }

  // If we loaded the spectra from file, we don't need to add any MC
  // events.  But add fake data events in the far detector - this
  // isn't the main use case, but it is useful for testing
  if (fMCInFilename!="" &&
      (info.header->detector==Detector::kNear ||
       (info.header->dataType==SimFlag::kMC && !useMCAsData))) return;

  if(!fUseNC && info.analysis->isNC) return;
  if(!fUseCC && info.analysis->isCC) return;

  if (!info.reco->inFiducialVolume) return;

  Detector::Detector_t det=Detector::Detector_t(info.header->detector);
  if (det==Detector::kNear && !info.reco->isSimpleCutsClean) return;

  double energy;
  if (fUseSelectionEnergy) {
    energy=info.reco->nuEnergy;
  } else {
    energy=fUseCCEnergy ? info.reco->nuEnergyCC : info.reco->nuEnergyNC;
  }

  // Don't use events with energy<0: this happens for NC energy for
  // events with no shower. Also obey energy threshold option.
  if (energy<0 || energy<fEnergyThreshold) return;

  if (info.header->dataType==SimFlag::kData || useMCAsData) {
    // For fake data, reco->weight contains all the scaling we need (including POT)
    fSpectra[beamInfo]->FillData(det,
                                info.truth, info.reco, info.analysis, energy,
                                info.reco->weight);
  } else {
    fSpectra[beamInfo]->FillMC(det,
                              info.truth, info.reco, info.analysis, energy,
                              info.reco->weight);
  }

}

//==============================================================================

void NCExtrapolationPID::FindSpectraForPars(const NC::OscProb::OscPars* pars,
                                            const NC::SystPars& s,
                                            vector<NCBeam::Info> beamsToUse,
                                            vector<TH1*>& exps,
                                            vector<TH1*>& obss)
{
  using namespace NC::Utility;

  if(fDoSystematicInterpolation) InitializeInterpolator(pars);

  if(!fInterpolator){
    exps = GetExpHists(beamsToUse, pars);
    obss = GetObsHists(beamsToUse);
    return;
  }


  const vector<double> shift = fCoordConv.VectorFromSystPars(s);

  vector<TH1*> interp = fInterpolator->GetInterpolatedSpectra(shift);

  vector<TH1*> nd_ratios;
  assert(interp.size()%2 == 0); // Every FD has a corresponding ND
  // The ND ratios are all the odd numbered ones
  for(unsigned int n = 1; n < interp.size(); n += 2){
    nd_ratios.push_back(interp[n]);
  }

  vector<const TH1*> nd_noms = GetNearMCSpectra(beamsToUse);

  assert(nd_ratios.size() == nd_noms.size());

  vector<TH1*> nd_spectra;
  for(unsigned int n = 0; n < nd_noms.size(); ++n){
    nd_spectra.push_back(CloneFast(nd_noms[n]));
    MultiplyFast(nd_spectra[n], nd_ratios[n]);
  }
  assert(nd_spectra.size() == nd_noms.size());

  // Pass nd_spectra to GetExpHists so the FD spectra take account of ND shifts
  exps = GetExpHists(beamsToUse, pars, nd_spectra);
  // Get observed spectra - no interpolation needs to be done on these
  obss = GetObsHists(beamsToUse);

  assert(interp.size() == 2*exps.size());
  // The FD ratios are the even numbered elements, so their indices correspond
  // to twice those in 'exps'
  for(unsigned int n = 0; n < exps.size(); ++n){
    MultiplyFast(exps[n], interp[2*n]);
  }

  // TODO - deleting spectra takes silly amounts of time
  for(unsigned int n = 0; n < interp.size(); ++n) delete interp[n];
  for(unsigned int n = 0; n < nd_spectra.size(); ++n) delete nd_spectra[n];
}

//==============================================================================

vector<const TH1*> NCExtrapolationPID::
GetNearMCSpectra(vector<NCBeam::Info> beamsToUse)
{
  vector<const TH1*> expHists;
  for (int i=0; i<(int)beamsToUse.size(); ++i) {
    NCBeam::Info beamInfo=beamsToUse[i];
    SpecMap::const_iterator spectrumForBeamInfo=fSpectra.find(beamInfo);
    assert(spectrumForBeamInfo!=fSpectra.end());
    assert(spectrumForBeamInfo->second);
    expHists.push_back(spectrumForBeamInfo->second->GetNearMC());
  }

  return expHists;
}

//==============================================================================

void NCExtrapolationPID::CleanupSpectra(vector<TH1*> exp, vector<TH1*> /*obs*/)
{
  // TODO - lots of time goes into deleting histograms.
  for (unsigned int n = 0; n < exp.size(); ++n) delete exp[n];
}

//==============================================================================

void NCExtrapolationPID::
WriteResources(const NC::OscProb::OscPars* trueOscPars)
{

  NC::OscProb::OscPars* bestFit=GetBestFitOscPars();
  NC::OscProb::OscPars* oscParsNoOsc = new NC::OscProb::NoOscillations;

  // Hack to fill NCBeam's histograms with something semi-reasonable
  // that includes the F/N ratio, for ease of comparison with other
  // methods
  
  // Loop over all of the beams
  for(fBeams_t::iterator it=fBeams.begin(); it!=fBeams.end(); ++it){
    NCBeam::Info beamInfo=it->first.first;
    // it->first.second is the detector. We'll just do the FD for now, so ignore it
    NCBeam* beam=it->second;
    assert(fSpectra.find(beamInfo)!=fSpectra.end());
    PIDSpectrum* spectrum=fSpectra[beamInfo];
    
    for(int inccc=NCType::kNC; inccc<=NCType::kCC; ++inccc){
      // Gah, I still don't understand int <-> enum implicit casting rules
      NCType::EEventType nccc=NCType::EEventType(inccc);
      TH1D* pseudoNoOsc=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kAll, 
                                                        oscParsNoOsc);
      TH1D* pseudoFit=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kAll, 
                                                      bestFit);
      TH1D* pseudoCC=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kNumu, 
                                                     bestFit);
      TH1D* pseudoBeamNue=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kBeamNue, 
                                                          bestFit);
      TH1D* pseudoOscNue=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kOscNue, 
                                                         bestFit);
      TH1D* pseudoNutau=spectrum->GetPseudoNCCCSpectrum(nccc, PIDSpectrum::kNutau, 
                                                        bestFit);
      
      TH1D* pseudoData=spectrum->GetPseudoNCCCComponent(spectrum->GetDataHist(), nccc);
      
      // The NCBeam plots want background for both spectra (ie, CC in
      // the NC spectrum, and NC in the CC spectrum), but we don't
      // have the NC component alone, so we have to form it here
      TH1D* pseudoBG;
      if(nccc==NCType::kNC) pseudoBG=pseudoCC;
      else{
        // Form the NC BG in the CC spectrum as the difference between
        // the total and the various CC parts
        pseudoBG=(TH1D*)pseudoFit->Clone("pseudoBG");
        pseudoBG->SetDirectory(0);
        
        pseudoBG->Add(pseudoCC, -1);
        // These are going to be tiny and not matter, but may as well
        // subtract them anyway
        pseudoBG->Add(pseudoBeamNue, -1);
        pseudoBG->Add(pseudoOscNue, -1);
        pseudoBG->Add(pseudoNutau, -1);
      }
      
#define COPYBINS(from, to) CopyBinContents(from,                        \
                                           beam->Get##to##Histogram(nccc)) 
      
      COPYBINS(pseudoData,    Data);
      COPYBINS(pseudoNoOsc,   MC);
      COPYBINS(pseudoFit,     MCFit);
      COPYBINS(pseudoBG,      MCBackground);
      COPYBINS(pseudoBeamNue, MCFitNuEToNuE);
      COPYBINS(pseudoOscNue,  MCFitNuMuToNuE);
      COPYBINS(pseudoNutau,   MCFitNuMuToNuTau);

#undef COPYBINS

    }      
    // We filled the histograms manually, and don't need NCBeam's automatic stuff
    // to happen. So inform it of that fact.
    beam->MarkHistogramsFilled();
  }

  delete oscParsNoOsc;

  NCExtrapolation::WriteResources(trueOscPars);

  // Done with the NCBeam stuff, so now we can get on with making my plots
    
  // Store my plots in a debug directory
  gDirectory->mkdir("debug", "debug")->cd();
    
  if(bestFit){
    NC::SystPars systs=GetBestFitSysts();
    
    WriteSpectra(bestFit, systs, "bestFit");
    
    if(trueOscPars) WriteSpectra(trueOscPars, systs, "truePars");
    
  }
    
  // If requested, save the PIDSpectrum objects out to file
  if (fMCOutFilename!="") WriteMCSpectra(fMCOutFilename.c_str());
}
  
//==============================================================================

void NCExtrapolationPID::PrintMap(POTMap& m, const char* name)
{
  cout << "Map " << name << endl;
  for (POTMap::const_iterator it=m.begin();
       it!=m.end(); ++it) {
    cout << it->first.GetDescription() << "\t" << it->second << endl;
  }
}

//==============================================================================

void NCExtrapolationPID::PrintSpectraMap()
{
  cout << "fSpectra:" << endl;
  for (SpecMap::const_iterator it=fSpectra.begin();
       it!=fSpectra.end(); ++it) {
    cout << it->first.GetDescription() << "\t" << hex << it->second << dec
         << "\t" << it->second->GetBeamInfo().GetTitle() << endl;
  }
}


//=================================================================================

vector<TH1*> NCExtrapolationPID::GetExpHists(vector<NCBeam::Info> beamsToUse,
                                             const NC::OscProb::OscPars* pars,
                                             vector<TH1*> ndSpectra) const
{
  vector<TH1*> expHists;
  for (int i=0; i<(int)beamsToUse.size(); ++i) {
    NCBeam::Info beamInfo=beamsToUse[i];
    SpecMap::const_iterator spectrumForBeamInfo=fSpectra.find(beamInfo);
    assert(spectrumForBeamInfo!=fSpectra.end());
    assert(spectrumForBeamInfo->second);
    expHists.push_back(spectrumForBeamInfo->second->
                       GetPredicted(pars, 
                                    "", 
                                    fDoFarNearCorrection,
                                    ndSpectra.empty() ? 0 : ndSpectra[i]));
  }

  return expHists;
}

//=================================================================================

vector<TH1*> NCExtrapolationPID::GetObsHists(vector<NCBeam::Info> beamsToUse) const

{
  vector<TH1*> obsHists;
  for (int i=0; i<(int)beamsToUse.size(); ++i) {
    NCBeam::Info beamInfo=beamsToUse[i];
    obsHists.push_back(fSpectra.find(beamInfo)->second->GetDataHist());
  }

  return obsHists;
}

//=================================================================================

void NCExtrapolationPID::WriteSpectra(const NC::OscProb::OscPars* coords, 
                                      const NC::SystPars& s,
                                      string suffix)
{
  if(suffix!="") gDirectory->mkdir(suffix.c_str(), suffix.c_str())->cd();

  vector<TH1*> fd_ratios;
  vector<TH1*> nd_ratios;

  if (fInterpolator) {
    const vector<double> shift = fCoordConv.VectorFromSystPars(s);
    vector<TH1*> ratios = fInterpolator->GetInterpolatedSpectra(shift);

    assert(ratios.size()%2 == 0); // near and far

    for(unsigned int n = 0; n < ratios.size(); n += 2)
      nd_ratios.push_back(ratios[n]);
    for(unsigned int n = 1; n < ratios.size(); n += 2)
      fd_ratios.push_back(ratios[n]);
  }

  int specIdx = 0;

  for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
    PIDSpectrum* spec=sit->second;

    TString beamName=spec->GetBeamInfo().GetDescription();
    assert(beamName!="");
    gDirectory->mkdir(beamName, beamName)->cd();

    TH1* nd_ratio = 0;
    TH1* fd_ratio = 0;

    // Need an interpolator obviously. Interpolating the shifted spectra makes
    // no sense. Interpolator doesn't return a kRunAll ratio which
    // unfortunately means that the output RunAll spectrum won't be
    // interpolated at all.
    const bool canInterpolate =
      fInterpolator &&
      spec->GetBeamInfo().IsNominal() &&
      spec->GetBeamInfo().GetRunType() != NC::RunUtil::kRunAll;

    if(canInterpolate){
      nd_ratio=nd_ratios[specIdx];
      fd_ratio=fd_ratios[specIdx];

      // These should match up, if not we've done the indexing wrong
      MSG("NCExtrapolationPID", Msg::kDebug)
        << "Wanted " << beamName << " found "
        << nd_ratio->GetName() << " and "
        << fd_ratio->GetName() << endl;

      ++specIdx;
    }

    vector<TCanvas*> v;
    v.push_back( spec->DrawSpectrum(coords, fd_ratio, suffix) );

    vector<TCanvas*> vNear=spec->DrawNearSpectra(nd_ratio, suffix);
    vector<TCanvas*> vE=spec->DrawEnergySlices(coords, fd_ratio, suffix);
    vector<TCanvas*> vPID=spec->DrawPIDSlices(coords, fd_ratio, suffix);

    v.insert(v.end(), vE.begin(),    vE.end());
    v.insert(v.end(), vPID.begin(),  vPID.end());
    v.insert(v.end(), vNear.begin(), vNear.end());

    for (vector<TCanvas*>::iterator it=v.begin();
         it!=v.end(); ++it) {
      (*it)->Write();
    }

    vector<PIDSpectrum::evtType> evtTypes=spec->EventTypeList();
    for (int j=0; j<(int)evtTypes.size(); ++j) {
      spec->GetOnePredicted(evtTypes[j], coords, suffix.c_str())->Write();
    }

    gDirectory->cd("..");
  }


  gDirectory->cd("..");
  // Also get the histograms of the tau component
  // (for the "can we see taus?" study)
  //for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit)
  //  sit->second->GetOnePredicted(PIDSpectrum::kNutau, coords, suffix.c_str())->Write();
}

//=================================================================================

void NCExtrapolationPID::DoneFilling()
{
  cout << "NCExtrapolationPID::DoneFilling()" << endl;
  // Only need to do this if we didn't load the spectra from file
  if(fMCInFilename == ""){
    for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
      DoneFillingOneBeam(sit->first);
    }
  }
}

//=================================================================================

void NCExtrapolationPID::DoneFillingOneBeam(NCBeam::Info beamInfo)
{
  cout << "NCExtrapolationPID::DoneFillingOneBeam("
       << beamInfo.GetDescription() << ")" << endl;

  cout << "Normalizing true-to-recos" << endl;
  fSpectra[beamInfo]->NormalizeTrueToRecos();

}

//=================================================================================

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

  r.UnLockValues();

  // Energy measure to use. Can be "CC", "NC" or "Selection" (which
  // uses NC or CC based on the event type given by the selection)
  r.Set("PIDFitEnergyType",   "CC");

  // Number of bins in PID to use
  r.Set("PIDFitPIDNbins",          20);
  // Minimum, maximum value of PID
  r.Set("PIDFitPIDmin",          0.2);
  r.Set("PIDFitPIDmax",          1.1);

  // How much finer to make the true energy binning than the reco
  // energy binning in the PIDSpectrum objects
  r.Set("PIDFitTrueBinFactor",   5);

  // The file to read the MC spectra from/write them to
  r.Set("PIDFitWriteMCToFile", "");
  r.Set("PIDFitReadMCFromFile", "");

  // Whether to do the far/near correction
  r.Set("PIDFitDoFarNearCorrection", true);

  // Emulate the FarNear method - overrides the number of bins setting
  // and the energy type if set to true
  r.Set("PIDFitEmulateFarNear", false);

  r.LockValues();
  return r;

}

//=================================================================================

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

  cout << "NCExtrapolationPID::Prepare" << endl;
  const char* tmps;
  int tmpi;
  double tmpd;

  if (r.Get("PIDFitEnergyType", tmps)) {
    string energyType(tmps);
    assert(energyType=="NC" || energyType=="CC" || energyType=="Selection");
    fUseSelectionEnergy = (energyType=="Selection");
    fUseCCEnergy = (energyType=="CC");
  }

  if (r.Get("PIDFitPIDNbins", tmpi)) fNumPIDBins=tmpi;

  // Minimum, maximum value of PID
  if (r.Get("PIDFitPIDmin", tmpd)) fPIDmin=tmpd;
  if (r.Get("PIDFitPIDmax", tmpd)) fPIDmax=tmpd;

  cout << "fNumPIDBins=" << fNumPIDBins
       << " PID range=[" << fPIDmin << ", " << fPIDmax << "]" << endl;

  if (r.Get("PIDFitTrueBinFactor", tmpi))       fTrueBinFactor=tmpi;
  if (r.Get("PIDFitDoFarNearCorrection", tmpi)) fDoFarNearCorrection=tmpi;

  if (r.Get("PIDFitEmulateFarNear", tmpi)) fEmulateFarNear=tmpi;

  if (fEmulateFarNear) {
    // If we're emulating the FarNear method, we just have two PID
    // bins, but we use fNumPIDBins=-1 to signal to PIDSpectrum that it
    // should put the separation point where the cut value is
    fNumPIDBins=-1;
    fUseSelectionEnergy=true;

    if (r.Get("PIDFitPIDNbins", tmpi) &&
        tmpi!=DefaultConfig().GetInt("PIDFitPIDNbins"))
      cout << "WARNING: PIDFitPIDNbins was set, but overridden by PIDFitEmulateFarNear"
           << endl;
  }

  // Files to read the MC spectra from/write them to
  if (r.Get("PIDFitReadMCFromFile", tmps)) fMCInFilename=tmps;
  if (r.Get("PIDFitWriteMCToFile", tmps))  fMCOutFilename=tmps;

  if (fMCInFilename!="") {
    ReadMCSpectra(fMCInFilename.c_str());
    // reset the observed spectra - if we got the PIDSpectrum's from file, we'll be
    // filling them with mock data presumably anyway
    for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
      TH2F* h=sit->second->GetDataHist();
      h->Reset();
    }
  }

}

//=================================================================================
void NCExtrapolationPID::Reset(bool nearData, bool farData,
                               bool nearMC, bool farMC)
{
  for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
    if(nearData) sit->second->ResetNearData();
    if(farData) sit->second->ResetData();
    if(nearMC) sit->second->ResetNearMC();
    if(farMC) sit->second->ResetMC();
  }
}

//=================================================================================
bool NCExtrapolationPID::EnableNearToFarExtrapolation(bool enable)
{
  const bool ret = fDoFarNearCorrection;
  fDoFarNearCorrection = enable;
  return ret;
}

//=================================================================================

void NCExtrapolationPID::WriteMCSpectra(const char* filename)
{
  TFile f(filename, "RECREATE");
  for (SpecMap::iterator sit=fSpectra.begin(); sit!=fSpectra.end(); ++sit) {
    sit->second->Write();
  }
}

//=================================================================================

void NCExtrapolationPID::ReadMCSpectra(const char* filename)
{

  cout << "Reading MC spectra from file " << filename << endl;
  TFile f(filename);

  TList *l = f.GetListOfKeys();
  TIter it(l);
  TObject *obj;
  TKey* k;
  while (( k=(TKey*)it() )) {
    obj=f.Get(k->GetName());
    cout << "Trying object " << obj->GetName() << endl;
    if (! (obj->IsA() == PIDSpectrum::Class()) ) {
      cerr << "Object with name " << obj->GetName() << " not a PIDSpectrum. Bailing" << endl;
      exit(1);
    }
    cout << "...object is a PIDSpectrum. Adding to list" << endl;
    PIDSpectrum* p=(PIDSpectrum*)obj;

    // fNumPIDBins is mangled when using fEmulateFarNear, so need this branch
    if (fEmulateFarNear) assert(p->GetNPIDBins()==2);
    else                 assert(p->GetNPIDBins()==fNumPIDBins);

    fSpectra[p->GetBeamInfo()]=p;
  }

}

//=================================================================================#

void NCExtrapolationPID::CopyBinContents(TH1D* hFrom, TH1F* hTo) const
{
  for(int i = 0; i < kNumEnergyBinsFar; ++i)
    hTo->SetBinContent(i+1, hFrom->GetBinContent(i+1));
}

---