NCExtrapolationModule.cxx

Go to the documentation of this file.

//$Id: NCExtrapolationModule.cxx,v 1.234 2010/01/19 17:02:22 rodriges Exp $
//
//NCExtrapolationModule.cxx
//
//Module for analyzing beam data for nc analysis
//
//B. Rebel 10/2004

#include "NCUtils/Extrapolation/NCExtrapolationModule.h"

#include "NCUtils/Extrapolation/NCCoordinateConverter.h"
#include "NCUtils/Extrapolation/NCEventAdder.h"
#include "NCUtils/Extrapolation/NCExtrapolation.h"
#include "NCUtils/Extrapolation/NCFitMaster.h"
#include "NCUtils/Extrapolation/NCPOTCounter.h"
#include "NCUtils/NCOscProb.h"
#include "NCUtils/NCRunUtil.h"

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

#include "MessageService/MsgService.h"

#include "TCanvas.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TLegend.h"
#include "TRandom3.h"
#include "TSystem.h"
#include "TStopwatch.h"

#include <cassert>
#include <algorithm>

using namespace BeamType;
using namespace NCType;

CVSID("$Id: NCExtrapolationModule.cxx,v 1.234 2010/01/19 17:02:22 rodriges Exp $");

//......................................................................
NCExtrapolationModule::NCExtrapolationModule() :
  fTrueOscPars(0),
  fNumMockExperiments(0),
  fPOTCount(0)
{
  MSG("JobC", Msg::kDebug) << "NCExtrapolationModule::Constructor" << endl;

  // "Special" means studying systematic shifts.
  // These are the TOm selection variables.
  for(int k = 3; k < 6; ++k){
    fDQNearMCTotalSpecial.push_back(CreateSpecialPlot(k, "NearMCTotalSpecial"));
  }

  fEventInfo.header = new ANtpHeaderInfo;
  fEventInfo.beam = new ANtpBeamInfo;
  fEventInfo.reco = new ANtpRecoInfo;
  fEventInfo.analysis = new ANtpAnalysisInfo;
  fEventInfo.truth = new ANtpTruthInfoBeam;

  fPOTCount = new NCPOTCounter;
}

enum{
  kDQEventLength      = 3,
  kDQNumTracks        = 4,
  kDQTrackExtension   = 5
};

//----------------------------------------------------------------------
TH1* NCExtrapolationModule::CreateSpecialPlot(int k, TString type) const
{
  assert(k >= 3 && k < 6);

  // Definition of a data quality plot: name, title, binning
  struct DQDef
  {
    // Order of elements is important, this is how we are initialised
    TString name;
    TString title;
    int bins;
    double start;
    double end;
  };

  const DQDef kDQVars[6] = {
    {"?", "?", -1, -1, -1},
    {"?", "?", -1, -1, -1},
    {"?", "?", -1, -1, -1},
    {"eventLength",    "Event Length (planes)",    50,   0,  500},
    {"numTracks",      "Number of Tracks",         10,   0,   10},
    {"trackExtension", "Track Extension (planes)", 50, -30,   20}
  };

  TH1* ret = new TH1D(kDQVars[k].name+type, "", kDQVars[k].bins, kDQVars[k].start, kDQVars[k].end);

  ret->GetXaxis()->SetTitle(kDQVars[k].title);
  ret->GetYaxis()->SetTitle("Events");

  ret->GetXaxis()->CenterTitle();
  ret->GetYaxis()->CenterTitle();

  return ret;
}

//----------------------------------------------------------------------
NCExtrapolationModule::~NCExtrapolationModule()
{
  MSG("JobC", Msg::kDebug) << "NCExtrapolationModule::Destructor" << endl;

  for(unsigned int n = 0; n < fExtrapolations.size(); ++n) delete fExtrapolations[n];
  for(unsigned int n = 0; n < fFitMasters.size(); ++n) delete fFitMasters[n];

  if(fEventInfo.header) delete fEventInfo.header;
  if(fEventInfo.beam) delete fEventInfo.beam;
  if(fEventInfo.reco) delete fEventInfo.reco;
  if(fEventInfo.reco_nom) delete fEventInfo.reco_nom;
  if(fEventInfo.truth) delete fEventInfo.truth;
  if(fEventInfo.analysis) delete fEventInfo.analysis;

  if(fTrueOscPars) delete fTrueOscPars;

  if(fPOTCount) delete fPOTCount;

  for(unsigned int n = 0; n < fFDMCSample.size(); ++n){
    delete fFDMCSample[n]->header;
    delete fFDMCSample[n]->beam;
    delete fFDMCSample[n]->reco;
    delete fFDMCSample[n]->reco_nom;
    delete fFDMCSample[n]->truth;
    delete fFDMCSample[n]->analysis;
    delete fFDMCSample[n];
  }

  for(unsigned int n = 0; n < fMockFits.size(); ++n) delete fMockFits[n];
}

//----------------------------------------------------------------------
void NCExtrapolationModule::DrawDataQualityPlots()
{
  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* saveDir = gDirectory;

  //loop over the data quality histograms
  //make folders for data/mc and subfolders for near/far
  TDirectory* dir = saveDir->mkdir("DataQuality", "DataQuality");
  dir->cd();

  TDirectory* near = dir->mkdir("NearDetector", "NearDetector");

  near->cd();

  if(fMakeDQPlotsSpecial) DrawDataQualityPlotsSpecial();

  saveDir->cd();
}

//----------------------------------------------------------------------
void NCExtrapolationModule::DrawDataQualityPlotsSpecial()
{
  for(int dist = 0; dist < 3; ++dist){

    fDQNearMCTotalSpecial[dist]->Write();


    TString canvName = fDQNearMCTotalSpecial[dist]->GetName();
    //    int index = canvName.Index("MC");
    canvName += "Canv";
    TCanvas *canv = new TCanvas(canvName, canvName, 150, 10, 900, 600);
    TLegend *leg = new TLegend(0.75, 0.75, 0.95, 0.95);
    leg->SetBorderSize(0);

    fDQNearMCTotalSpecial[dist]->Draw("");
    leg->AddEntry(fDQNearMCTotalSpecial[dist], "MC Total", "l");
    leg->Draw();

    gDirectory->Append(canv);

  }//end loop over dists
}

//----------------------------------------------------------------------
void NCExtrapolationModule::Run()
{
  fConfig.Print(); // Print out the configuration in use.

  assert(fTrueOscPars);

  //get the pot for the different beams
  const NCType::EFileType fileType = fUseMockData ? NCType::kMockFile : NCType::kBeamFile;

  //use following map for FD MC
  fPOTCount->SetPOTValues(fDataPOTNear,          Detector::kNear,
                          NCType::kBeamFile,     NCType::kData);

  fPOTCount->SetPOTValues(fMCPOTNear,            Detector::kNear,
                          NCType::kBeamFile,     NCType::kMC);

  fPOTCount->SetPOTValues(fDataPOTFar,           Detector::kFar,
                          fileType,              NCType::kData);

  fPOTCount->SetPOTValues(fMCPOTFar,             Detector::kFar,
                          NCType::kBeamFile,     NCType::kMC);

  fPOTCount->SetPOTValues(fMCPOTFarTau,          Detector::kFar,
                          NCType::kTauFile,      NCType::kMC);

  fPOTCount->SetPOTValues(fMCPOTFarElectron,     Detector::kFar,
                          NCType::kElectronFile, NCType::kMC);

  MSG("NCExtrapolationModule", Msg::kInfo) << "Done setting POT values" << endl;


  ChainMap dataND;
  ChainMap mcND;
  ChainMap dataFD;
  ChainMap mcFD;
  ChainMap mcTauFD;
  ChainMap mcElectronFD;

  // The beamtype/runperiod pairs we're using
  vector<NCBeam::Info> infos;

  for(unsigned int bidx = 0; bidx < fBeamIndex.size(); ++bidx){
    const BeamType::BeamType_t bt = fBeamIndex[bidx];
    for(unsigned int runPidx=0; runPidx < fRunsToUse.size(); ++runPidx){
      const NC::RunUtil::ERunType runPeriod = fRunsToUse[runPidx];
      const NCBeam::Info info(bt, runPeriod);

      infos.push_back(info);
    }
  }

  for(unsigned int i=0; i<infos.size(); ++i){
    NCBeam::Info info=infos[i];
    dataND[info] = new TChain("uDST");
    mcND[info] = new TChain("uDST");
    dataFD[info] = new TChain("uDST");
    mcFD[info] = new TChain("uDST");
    mcTauFD[info] = new TChain("uDST");
    mcElectronFD[info] = new TChain("uDST");
  }

  AddFilesToChain(dataND,  Detector::kNear, NCType::kBeamFile, NCType::kData);
  AddFilesToChain(mcND,    Detector::kNear, NCType::kBeamFile, NCType::kMC);
  AddFilesToChain(dataFD,  Detector::kFar,  fileType,          NCType::kData);
  AddFilesToChain(mcFD,    Detector::kFar,  NCType::kBeamFile, NCType::kMC);
  AddFilesToChain(mcTauFD, Detector::kFar,  NCType::kTauFile,  NCType::kMC);
  AddFilesToChain(mcElectronFD, Detector::kFar, NCType::kElectronFile, NCType::kMC);

  const TString ext = NCType::kExtractionNames[fExtractionType];

  for(unsigned int i=0; i<infos.size(); ++i){
    NCBeam::Info info=infos[i];
    fEventInfo.SetChainBranches(dataND[info],       ext, fUseMCAsData);
    fEventInfo.SetChainBranches(dataFD[info],       ext, fUseMCAsData);
    fEventInfo.SetChainBranches(mcND[info],         ext, true);
    fEventInfo.SetChainBranches(mcFD[info],         ext, true);
    fEventInfo.SetChainBranches(mcTauFD[info],      ext, true);
    fEventInfo.SetChainBranches(mcElectronFD[info], ext, true);
  } 


  CreateExtrapolations(fExtrapolationsList);

  PrepareForExtrapolation();

  NC::IEventAdder* eventAdder = NC::EventAdderBuilder(fConfig);

  for(unsigned int i=0; i<infos.size(); ++i){
    NCBeam::Info info=infos[i];
    eventAdder->AddEvents(this, &fEventInfo,
                          dataND[info], mcND[info],
                          dataFD[info], mcFD[info], mcTauFD[info], mcElectronFD[info]);
  }
  delete eventAdder;

  ExtrapolateNearToFar();

  if(fNumMockExperiments > 0) DoMockExperiments();

  //make a file and the trees to go into the file

  MSG("NCExtrapolationModule", Msg::kInfo) << "fFileName = " << fFileName << endl;

  TFile ntpFile(fFileName, "RECREATE");
  ntpFile.cd();

  for(unsigned int k = 0; k < fExtrapolations.size(); ++k){
    TDirectory* saveDir = gDirectory;
    TString plotName = "plots"+fExtrapolations[k]->GetShortName();
    TDirectory* newDir = gDirectory->mkdir(plotName, plotName);

    // In case we failed to make the new directory then this way
    // probably gives us the best chance of rescuing the plots
    if(newDir) newDir->cd();
    fFitMasters[k]->WriteResources();
    if(newDir) newDir->cd();
    fExtrapolations[k]->WriteResources(fTrueOscPars);

    saveDir->cd();
  }

  if(!fUseMCAsData && !fUseMockData) DrawDataQualityPlots();

  ntpFile.cd();
  fConfig.SetName("config_registry");
  fConfig.Write();

  if(!fMockFits.empty()){
    ntpFile.cd();
    gDirectory->mkdir("mock_expts")->cd();
    for(unsigned int n = 0; n < fMockFits.size(); ++n){
      Registry* r = fMockFits[n];
      r->SetName(TString::Format("expt_%d", n));
      r->Write();
    }
    ntpFile.cd();
  }

  MSG("NCExtrapolationModule", Msg::kInfo) << "Writing Final File!" << endl;

  ntpFile.Write();
  ntpFile.Close();

  if(fSaveExtrapolationsFile != ""){
    TFile saveFile(fSaveExtrapolationsFile, "RECREATE");
    for(unsigned int k = 0; k < fExtrapolations.size(); ++k)
      fExtrapolations[k]->Write();
    saveFile.Close();
  }
}

//......................................................................
void NCExtrapolationModule::AddExtrapolation(NCExtrapolationFactory* e)
{
  std::vector<NCExtrapolationFactory*>& facts = GetExtrapolationFactories();
  for(unsigned int n = 0; n < facts.size(); ++n){
    // If this happens it means that extrapolations have got registered
    // twice somehow. Check your environment, and rebuild stuff, including
    // the macro.
    assert(e->GetCodeName() != facts[n]->GetCodeName());
  }
  facts.push_back(e);
}

//......................................................................
void NCExtrapolationModule::SetMCExposureForBeam(NCBeam::Info beam, double pot)
{
  fPOTCount->SetMCExposureForBeam(beam, pot);
}

//......................................................................
const Registry& NCExtrapolationModule::DefaultConfig() const
{
  static Registry r;

  r.UnLockValues();

  r.Set("DataMCPath",             "dataFiles/*.root");
  r.Set("FarDataPath",            "dataFiles/*.root");
  r.Set("FileName",               "extrapolationFile.root");

  r.Set("OscillationModel",       NCType::kSterileFraction);
  r.Set("UseMCAsData",            false);
  r.Set("SplitMC",                false);
  r.Set("SplitMCSeed",            -1);
  r.Set("ExtractionType",         NCType::kNCCCExtractionCuts);
  r.Set("SingleSpectrum",         false);
  r.Set("ExtrapolationsList",     "");
  r.Set("BeamType",               "L010z185i");
  r.Set("UseMockData",            false);
  r.Set("MockDataSet",            "F21910001");

  TString runsToUse;
  for(int n = 0; n <= NC::RunUtil::kMaxRun; ++n)
    runsToUse += TString::Format("%d ", n);
  r.Set("RunsToUse",              runsToUse);

  r.Set("OscPars",                0);
  r.Set("NumMockExperiments",     0);
  r.Set("MockExperimentsSeed",    0);
  r.Set("LoadExtrapolationsFile", "");
  r.Set("IncludeDataFromLoadedExtrapolations", false);
  r.Set("SaveExtrapolationsFile", "");

  r.Set("MakeDQPlotsSpecial",     false);

  // Get the default configurations for all the registered extrapolations
  vector<NCExtrapolationFactory*> facts =  GetExtrapolationFactories();
  for(unsigned int n = 0; n < facts.size(); ++n)
    r.Merge(facts[n]->DefaultConfig());

  //PID sensitivity study
  r.Set("UsePIDCustomCut",        false);
  r.Set("PIDCustomCut",           0.0);

  r.Set("MakeShiftedBeams",       false);
  r.Set("ShiftedBeamsSimpleOnly", false);
  r.Set("SingleShiftedBeam",      -1);

  //set the defaults for the systematic parameters
  //set them all off for fitting and changing mc by default
  TString adjust   = "Adjust";
  TString mcAsData = "ChangeMCAsData";
  TString fit      = "Fit";
  for(int i = 0; i < NCType::kNumSystematicParameters; ++i){
    r.Set(NCType::kParams[i].name,            false);
    r.Set(mcAsData + NCType::kParams[i].name, false);
    r.Set(adjust + NCType::kParams[i].name,   0); // measured in sigma
    r.Set(fit + NCType::kParams[i].name,      false);
  }

  // Add all the configuration that NCExtrapolation knows how to read
  r.Merge(NCExtrapolation::DefaultConfig());

  r.Merge(NC::FitMaster::DefaultConfig());

  r.Merge(NC::EventAdderDefaultConfig());

  r.LockValues();
  return r;
}

//......................................................................
void NCExtrapolationModule::Config(const Registry& r)
{
  fConfig = r;

  int         tmpb;  // a temp bool.
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string

  if (r.Get("DataMCPath",            tmps)){
    //for some reason the TSystemDirectory constructor can't resolve
    //environmental variables, so just take care of it here
    TString temp = tmps;
    // Assume the form $ENVVAR/path
    if(temp.Contains("$")){
      int pos = temp.Index("/"); // find the first /
      temp.Resize(pos); // and truncate at it
      temp.Remove(TString::kLeading, '$'); // drop the $ off the front
      fDataMCPath = tmps; // Copy the original string
      // and then replace the part that we think is the environment variable
      fDataMCPath.Replace(0, pos, gSystem->Getenv(temp));
    }
    else fDataMCPath = temp;
    // Add a trailing slash so that concatenating the path later will work
    if(!fDataMCPath.EndsWith("/")) fDataMCPath += "/";
  }
  if(r.Get("FileName",               tmps)) fFileName           = tmps;

  if(r.Get("OscillationModel",       tmpi)) fOscillationModel   = NCType::EOscModel(tmpi);
  if(r.Get("UseMCAsData",            tmpb)) fUseMCAsData        = tmpb;
  if(r.Get("ExtractionType",         tmpi)) fExtractionType     = NCType::EExtraction(tmpi);
  if(r.Get("SingleSpectrum",         tmpb)) fSingleSpectrum     = tmpb;
  if(r.Get("ExtrapolationsList",     tmps)) fExtrapolationsList = tmps;
  if(r.Get("UseMockData",            tmpb)) fUseMockData        = tmpb;
  if(r.Get("UsePIDCustomCut",        tmpb)) fUsePIDCustomCut    = tmpb;
  if(r.Get("PIDCustomCut",           tmpd)) fPIDCustomCut       = tmpd;

  if(r.Get("MakeShiftedBeams",       tmpb)) fShiftedBeams       = tmpb;
  if(r.Get("ShiftedBeamsSimpleOnly", tmpb)) fSimpleShiftedBeams = tmpb;
  if(r.Get("SingleShiftedBeam",      tmpi)) fSingleShiftedBeam  = tmpi;

  if(r.Get("RunsToUse",              tmps)){
    fRunsToUse.clear();
    vector<int> runs = NC::Utility::ParseNumberList(tmps);
    for(unsigned int n = 0; n < runs.size(); ++n)
      fRunsToUse.push_back(NC::RunUtil::ERunType(runs[n]));
  }

  if(r.Get("NumMockExperiments",     tmpi)) fNumMockExperiments  = tmpi;
  if(r.Get("MockExperimentsSeed",    tmpi)) fMockExperimentsSeed = tmpi;
  if(r.Get("LoadExtrapolationsFile", tmps)) fLoadExtrapolationsFile = tmps;
  if(r.Get("IncludeDataFromLoadedExtrapolations", tmpb))
     fIncludeDataFromLoadedExtrapolations = tmpb;
  if(r.Get("SaveExtrapolationsFile", tmps)) fSaveExtrapolationsFile = tmps;

  if(r.Get("MakeDQPlotsSpecial",     tmpb)) fMakeDQPlotsSpecial  = tmpb;

  assert(sizeof(NC::OscProb::OscPars*) == sizeof(int));
  if(r.Get("OscPars", tmpi)) fTrueOscPars = (NC::OscProb::OscPars*)(tmpi);

  if (r.Get("BeamType",              tmps))
    fBeamIndex = NCType::BeamListFromString(tmps);

  fPOTCount->Config(r, fBeamIndex, fRunsToUse);
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::AddFilesToChain(ChainMap& chainMap,
                                            Detector::Detector_t det,
                                            NCType::EFileType fileType,
                                            NCType::EDataMC dataMC)
{
  MSG("NCExtrapolationModule", Msg::kInfo) << "fDataMCPath = "
                                           << fDataMCPath << endl;

  for(ChainMap::iterator it=chainMap.begin();
      it != chainMap.end();
      ++it){

    const NCBeam::Info info=it->first;
    const TString beamType=BeamType::AsString(info.GetBeamType());

    vector<TString> files = fPOTCount->GetListOfFiles(info, fileType,
                                                      det, dataMC);
    
    for(unsigned int de = 0; de < files.size(); ++de) {
      it->second->Add(fDataMCPath+beamType+"/"+files[de]);
      MSG("NCExtrapolationModule", Msg::kInfo) 
        << "Adding "
        << fDataMCPath+beamType
        << "/"+files[de] << " to "
        << NCType::FileTypeAsString(fileType) << " chain" << endl;

    }//end loop over files in directory
  }//end loop over map
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::FillDataQualityPlotsSpecial(std::vector<TH1*>& dq, bool osc)
{

  osc=false;

  //dont bother with events outside the fiducial volume
  if(fEventInfo.reco->inFiducialVolume < 1) return;

  double weight;
  switch(NC::RunUtil::FindRunType(fEventInfo.header)){
  case NC::RunUtil::kRunI: 
    weight=fEventInfo.reco->weight;
    break;
  case NC::RunUtil::kRunII: 
    weight=fEventInfo.reco->weightRunII;
    break;
  case NC::RunUtil::kRunIII: 
    weight=fEventInfo.reco->weightRunIII;
    break;
  default:
    assert(0 && "Unknown runtype");
  }

  double trackExtension = -1.*fEventInfo.reco->trackExtension;
  double showerEnergy = fEventInfo.reco->showerEnergyNC;

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

  dq[kDQEventLength-3]->Fill(fEventInfo.reco->eventLength, weight);
  dq[kDQNumTracks-3]->Fill(fEventInfo.reco->numberTracks, weight);

  if(fEventInfo.reco->numberTracks > 0
     && fEventInfo.reco->trackExtension > -9000.
     && showerEnergy > 0.0)
    dq[kDQTrackExtension-3]->Fill(trackExtension, weight);


  //  cout << "Weight: " << weight << endl
  //    << "eventLength: " << fEventInfo.reco->eventLength
  //    << " histo: " << dq[0]->GetBinContent(15) << endl
  //    << "NumTracks: " << fEventInfo.reco->numberTracks
  //    << " histo: " << dq[1]->GetBinContent(2) << endl
  //       << "TrackExtension: " << trackExtension  << "  " << showerEnergy
  //    << " histo: " << dq[2]->GetBinContent(30) << endl
  //    << endl;
}

//-----------------------------------------------------------------------
//preparation for extrapolating Near to Far
void NCExtrapolationModule::PrepareForExtrapolation()
{
  MSG("NCExtrapolationModule", Msg::kInfo) << "PrepareForExtrapolation" << endl;
  //pass the module registry into the extrapolations to prepare them
  for(unsigned int i = 0; i < fExtrapolations.size(); ++i)
    fExtrapolations[i]->Prepare(fConfig);
}


//-----------------------------------------------------------------------
//extrapolation step does all extrapolation algorithms at once
//the farMC chain is the monte carlo for the far detector.
//some extrapolations need to do stuff with the near first,
//then use those results to apply to the far monte carlo events
//for the prediction
void NCExtrapolationModule::ExtrapolateNearToFar()
{
  /*
  const bool loadingExtrapolations = fLoadExtrapolationsFile != "";

  // Ugh - should instead break this function down so we can select parts of it
  if(fLoadExtrapolationsFile != ""){
    for(unsigned int i = 0; i < fExtrapolations.size(); ++i){
      NC::FitMaster* fm = new NC::FitMaster(fExtrapolations[i]);
      fm->Prepare(fConfig);
      fm->Run(fTrueOscPars);
      fFitMasters.push_back(fm);
    }
    return;
  }
  */

  MSG("NCExtrapolationModule", Msg::kInfo) << "ExtrapolateNearToFar" << endl;

  //predict the far detector and make sure to combine any running conditions
  for(unsigned int i = 0; i < fExtrapolations.size(); ++i){
    fExtrapolations[i]->DoneFilling();
    NC::FitMaster* fm = new NC::FitMaster(fExtrapolations[i]);
    fm->Prepare(fConfig);
    fm->Run(fTrueOscPars);
    fFitMasters.push_back(fm);
  }
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::DoMockExperiments()
{
  if(fMockExperimentsSeed == 0)
    MSG("NCExtrapolationModule", Msg::kWarning) << "DoMockExperiments: random seed not set, using default (0)" << endl;


  TRandom3 r;

  for(int n = 0; n < fNumMockExperiments; ++n){
    // Reset everyone's FD data
    for(unsigned int i = 0; i < fExtrapolations.size(); ++i)
      fExtrapolations[i]->Reset(false, true, false, false);

    r.SetSeed(fMockExperimentsSeed);
    r.SetSeed(n+r.Integer(UInt_t(-1)));

    for(unsigned int n = 0; n < fFDMCSample.size(); ++n){
      NCEventInfo* s = fFDMCSample[n];
      for(unsigned int i = 0; i < fExtrapolations.size(); ++i){
        double tmp_weight = s->reco->weight;
        s->reco->weight = r.Poisson(s->reco->weight);
        // Is this way faster?
        //      s->reco->weight = (s->reco->weight > r.Uniform()) ? 1 : 0;

        BeamType::BeamType_t beamType = (BeamType::BeamType_t)(fEventInfo.beam->beamType);
        NCBeam::Info beamInfo(beamType, NC::RunUtil::kRunI);

        if(s->reco->weight){
          fExtrapolations[i]->AddEvent(*s, true,
                                       NCType::kBeamFile, // TODO is kBeamFile right?
                                       beamInfo);
        }
        s->reco->weight = tmp_weight;
      }//end loop over extrapolation objects
    }

    for(unsigned int i = 0; i < fExtrapolations.size(); ++i){
      NCExtrapolation* extrap = fExtrapolations[i];
      //      extrap->CombineRuns(Detector::kFar);
      NC::FitMaster fm(extrap);
      Registry r = fConfig;
      // Counteract the "preciseParams" scaling normally done in the
      // initial fit. This is a bit ugly, would be better to disable the
      // scaling instead of fighting it.

      // Changed my mind, we need to be fair to the real fit

      //      r.UnLockValues();
      //      r.Set("PrecScale", 0.01);

      fm.Prepare(r);

      fm.SetFixedValuesForMarginalization(fTrueOscPars);
      fm.Run(0);

      const double trueChi = extrap->FindChiSqrForPars(fTrueOscPars,
                                                       NC::SystPars());

      // Slight problem here in that we don't distinguish which extrapolation
      // this result came from, so there's no way to tell in the output file.
      // Isn't a problem if you only run one extrapolation, which is the common
      // use case.
      Registry* reg = fm.GetBestFitPointAsRegistry();
      reg->Set("true_chisq", trueChi);
      fMockFits.push_back(reg);
    }
  }
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::
SystematicsFromRegistry(bool* systematicsToUse,
                        vector<double>& systematicsAdjust)
{
  //get the parameters to be altered
  int         tmpb;  // a temp bool. See comment under DefaultConfig...
  double      tmpd;  // a temp double.

  //set the defaults for the systematic parameters from the config
  //set the adjustment to be an actual value, the macro should have
  //passed in a number of sigma

  for(int i = 0; i < NCType::kNumSystematicParameters; ++i){
    if(fConfig.Get("ChangeMCAsData" + NCType::kParams[i].name, tmpb))
      systematicsToUse[i] = tmpb;

    if(fConfig.Get("Adjust" + NCType::kParams[i].name, tmpd))
      systematicsAdjust.push_back(NCType::kParams[i].defaultVal +
                                  tmpd*NCType::kParams[i].sigma);

    if(systematicsToUse[i]){
      MAXMSG("NCExtrapolationModule", Msg::kInfo, NCType::kNumSystematicParameters)
        << "Using MC As Data: "
        << NCType::kParams[i].name << " adjusted to "
        << systematicsAdjust[i] << endl;
    }
  }
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::ApplySystematicShifts(NC::RunUtil::ERunType runType)
{
  static bool systematicsToUse[NCType::kNumSystematicParameters];
  static vector<double> systematicsAdjust;

  static bool once = true;
  if(once){
    once = false;
    SystematicsFromRegistry(systematicsToUse, systematicsAdjust);
  }

  // always use the nue values when finding the event weight because
  // the size of the nue contrib is passed in by the macro
  fEventInfo.SetEventWeight(systematicsToUse, systematicsAdjust,
                            fTrueOscPars, 0, runType, true, true);

  if(TMath::Abs(fEventInfo.reco->weight) < 5.e-2 &&
     fEventInfo.reco->weight < 0.)
    fEventInfo.reco->weight = 0.;
}

//-----------------------------------------------------------------------
void NCExtrapolationModule::SetEnergies()
{
  if(fEventInfo.analysis->isCC > 0){
    // The RHSs here *should* be reco and not reco_nom, because we
    // want the shifted value of (nu|shower)EnergyCC if it has been
    // shifted
    fEventInfo.reco->nuEnergy = fEventInfo.reco->nuEnergyCC;
    fEventInfo.reco->showerEnergy = fEventInfo.reco->showerEnergyCC;

    // We need to do this because (inter alia) the calculation of the
    // ND cleaning syst relies on showerEnergy being set to the right
    // one. We're safe from considerations of whether things have been
    // systematically shifted because reco_nom is never systematically
    // shifted
    fEventInfo.reco_nom->nuEnergy = fEventInfo.reco_nom->nuEnergyCC;
    fEventInfo.reco_nom->showerEnergy = fEventInfo.reco_nom->showerEnergyCC;
  }
  else if(fEventInfo.analysis->isNC > 0){
    fEventInfo.reco->nuEnergy = fEventInfo.reco->showerEnergyNC;
    fEventInfo.reco->showerEnergy = fEventInfo.reco->showerEnergyNC;

    // See comment above
    fEventInfo.reco_nom->nuEnergy = fEventInfo.reco_nom->nuEnergyNC;
    fEventInfo.reco_nom->showerEnergy = fEventInfo.reco_nom->showerEnergyNC;

    // (PAR) This used to be set here, but that breaks FindEventWeight
    // for NC events if we're using CC energy for those. But
    // NCExtrapolation::DrawBestFitSpectra relies on (shower+track)
    // giving the NC or CC energy appropriate to the classification,
    // so setting track energy to zero is now done in NCEnergyBin::AddEventToBin
    //fEventInfo.reco->muEnergy = 0.;
  }

  // Energies < 0 occur when there's no shower and we get a default value
  if(fEventInfo.reco->showerEnergy < 0.) fEventInfo.reco->showerEnergy = 0.;
}

//-----------------------------------------------------------------------
double NCExtrapolationModule::GetMCScaleFactor(Detector::Detector_t det,
                                               NCType::EFileType fileType,
                                               NCBeam::Info info)
{
  // Choose the appropriate pot maps for this event and scale
  // down the event to match the data POT.
  POTMap* potData=0;
  POTMap* potMC=0;
  switch(det){
  case Detector::kNear:
    potData=&fDataPOTNear;
    potMC=&fMCPOTNear;
    break;
  case Detector::kFar:
    potData=&fDataPOTFar;
    switch(fileType){
    case NCType::kTauFile:
      potMC=&fMCPOTFarTau;
      break;
    case NCType::kElectronFile:
      potMC=&fMCPOTFarElectron;
      break;
    default:
      potMC=&fMCPOTFar;
    }
    break;
  default:
    assert(0 && "Unknown detector. This shouldn't happen");
  }

  double tmp=GetScaleFactorFromMaps(*potData, *potMC, info);
  if(tmp>0) return tmp;

  // ...otherwise try to find the same beam with no syst shifts
  const vector<NCType::EFitParam> adj;
  const vector<double> vals;
  info.SetSystShifts(adj, vals);

  return GetScaleFactorFromMaps(*potData, *potMC, info);
}

//-----------------------------------------------------------------------
double NCExtrapolationModule::GetScaleFactorFromMaps(POTMap& potData,
                                                     POTMap& potMC,
                                                     NCBeam::Info info)
{
  // See if we actually found the beam in both...
  const bool foundData = potData.find(info) != potData.end();
  const bool foundMC = potMC.find(info) != potMC.end();

  // If we have data POTs for this beam, we'd better have some MC to go with it.
  if(foundData && !foundMC){
    MSG("NCExtrapolationModule", Msg::kError)
      << "Found data POT but not MC for " 
      << info << endl;
    
    MSG("NCExtrapolationModule", Msg::kError) << "Data POTs: " << endl;
    for(POTMap::iterator it=potData.begin(); it!=potData.end(); ++it)
      MSG("NCExtrapolationModule", Msg::kError)
        << it->first << ": " << it->second << endl;
    
    MSG("NCExtrapolationModule", Msg::kError) << "MC POTs: " << endl;
    for(POTMap::iterator it=potMC.begin(); it!=potMC.end(); ++it)
      MSG("NCExtrapolationModule", Msg::kError)
        << it->first << ": " << it->second << endl;
    exit(1);
  }

  // If we have no data, return 0;
  if(!foundData) return 0;

  // Otherwise return the actual scale factor
  return potData.find(info)->second/potMC.find(info)->second;
}

//-----------------------------------------------------------------------
void NCExtrapolationModule::SetEventWeight(BeamType::BeamType_t beamType,
                                           NC::RunUtil::ERunType runType,
                                           bool fakeDataSet)
{
  Detector::Detector_t det = Detector::Detector_t(fEventInfo.header->detector);
  SimFlag::SimFlag_t sim = SimFlag::SimFlag_t(fEventInfo.header->dataType);
  NCType::EFileType fileType = NCType::FindFileType(fEventInfo.header);

  // Data events don't get scaled or systematically shifted
  if(sim==SimFlag::kData){
    SetEnergies();
    return;
  }
  else if(sim==SimFlag::kMC){
    if(fakeDataSet){
      // This includes SKZP via NCEventInfo::GetEventWeight
      ApplySystematicShifts(runType);
    }
    else{
      // Just apply the SKZP weight
      fEventInfo.reco->weight=fEventInfo.GetSKZPWeight(runType);
    }

    const NCBeam::Info info=NCBeam::Info(beamType, runType);
    fEventInfo.reco->weight *= GetMCScaleFactor(det, fileType, info);
    // ApplySystematicShifts will (indirectly) have altered nuEnergy and
    // showerEnergy, so make sure they're set right
    SetEnergies();
  }
  else{
    assert(0 && "This shouldn't happen");
  }
}

//-----------------------------------------------------------------------
bool NCExtrapolationModule::FinalEventCheck()
{
  //check any final cuts here
  if(!fEventInfo.FinalEventCheck()
     && !fSingleSpectrum) return false;

  return fEventInfo.header->isGoodData;
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::CreateExtrapolations(TString extraps)
{
  MSG("NCExtrapolationModule", Msg::kInfo) << "Registered extrapolations:\n";

  const vector<TString> extrapArgNames = NC::Utility::ParseStringList(extraps);

  for(unsigned int n = 0; n < GetExtrapolationFactories().size(); ++n){
    NCExtrapolationFactory* ex = GetExtrapolationFactories()[n];
    const TString codeName = ex->GetCodeName();

    MSG("NCExtrapolationModule", Msg::kInfo) << "  " << codeName << " - ";

    bool inExtraps = false;
    for(unsigned int i = 0; i < extrapArgNames.size(); ++i)
      if(extrapArgNames[i] == codeName) inExtraps = true;

    if(inExtraps){
      if(fLoadExtrapolationsFile != ""){
        MSG("NCExtrapolationModule", Msg::kInfo) << "Loading from file\n";
        TDirectory* curDir = gDirectory;
        TFile f(fLoadExtrapolationsFile);
        NCExtrapolation* e = (NCExtrapolation*)f.Get("NCExtrapolation"+codeName);
        assert(e);
        fExtrapolations.push_back(e);
        f.Close();
        curDir->cd();
        // Restoring data spectra from the file is not generally useful.
        if(!fIncludeDataFromLoadedExtrapolations)
          e->Reset(true, true, false, false);
      }
      else{
        MSG("NCExtrapolationModule", Msg::kInfo) << "Creating\n";
        fExtrapolations.push_back(ex->Create());
      } // end if fLoadExtrapolationsFile
    }
    else{
      MSG("NCExtrapolationModule", Msg::kInfo) << "Not creating\n";
    } // end if inExtraps
  } // end for n
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::AddEventToExtrapolations(bool fakeDataSet)
{
  // This is pretty much immediately after the event comes off disk. Call
  // SetEnergies here, as early as possible, so showerEnergy and trackEnergy
  // are set to the right NC/CC variants. This is certainly needed for getting
  // systematic shifts right, and maybe other reasons.
  SetEnergies();

  SimFlag::SimFlag_t sim = SimFlag::SimFlag_t(fEventInfo.header->dataType);

  //Use PID custom selection cut if so desired
  if(fUsePIDCustomCut &&
     TMath::Abs(fPIDCustomCut-fEventInfo.analysis->separationParameterCut) > 1e-3){
    fEventInfo.analysis->separationParameterCut = fPIDCustomCut;

    fEventInfo.analysis->isNC = int(fEventInfo.analysis->separationParameter
                                    < fPIDCustomCut);
    fEventInfo.analysis->isCC = int(fEventInfo.analysis->separationParameter
                                    >= fPIDCustomCut);
  }

  //check for any extra cuts you may want to apply, eg good shower cut
  if(!FinalEventCheck()) return;

  //check if the event is from a tau file.
  const NCType::EFileType fileType = NCType::FindFileType(fEventInfo.header);
  const NC::RunUtil::ERunType runType  = NC::RunUtil::FindRunType(fEventInfo.header,
                                                                  fEventInfo.reco);

  if(fileType == NCType::kTauFile ||
     fileType == NCType::kElectronFile){
    //dont add in NC tau or nue appearance
    //events as you cant distinguish flavor
    //with NC events.  NCEnergyBin wouldnt
    //accept the event anyway.  also reject
    //events from these files that start off as
    //nue...just say they dont oscillate as
    //they dont contribute much anyway
    if(fEventInfo.truth->interactionType == NCType::kNC ||
       TMath::Abs(fEventInfo.truth->nonOscNuFlavor) == 12) return;
  }
  else if(fileType == NCType::kUnknownFile) return;

  // Single spectrum fitting
  if(fSingleSpectrum){
    fEventInfo.analysis->isNC = 0;
    fEventInfo.analysis->isCC = 1;
  }


  if(sim == SimFlag::kData || fileType == NCType::kMockFile)
    fEventInfo.truth = 0;

  // Mock data should have simflag set to data
  if (fileType == NCType::kMockFile)
    fEventInfo.header->dataType=SimFlag::kData;

  //--------------------------------------------------------------------------------
  //put in for testing purposes, should be commented out for real fits
  //if(TMath::Abs(fEventInfo.truth->nuFlavor) == 12) continue;
  //     fEventInfo.analysis->isNC = 0;
  //     fEventInfo.analysis->isCC = 0;
  //     if(fEventInfo.truth->interactionType < 1) fEventInfo.analysis->isNC = 1;
  //     else if(fEventInfo.truth->interactionType > 0) fEventInfo.analysis->isCC = 1;
  //     fEventInfo.reco->nuEnergy = truthInfo->nuEnergy;
  //--------------------------------------------------------------------------------

  // TODO - correct if statement?
  if(sim == SimFlag::kMC && fakeDataSet &&
     fEventInfo.header->detector == Detector::kFar &&
     fNumMockExperiments > 0){
    NCEventInfo* evt = new NCEventInfo;
    evt->DeepCopy(&fEventInfo);

    fFDMCSample.push_back(evt);
  }

  const BeamType::BeamType_t beamType =
    (BeamType::BeamType_t)(fEventInfo.beam->beamType);
  const NCBeam::Info beamInfo(beamType, runType);
  const NCBeam::Info infoAll(beamType, NC::RunUtil::kRunAll);


  //add event to each extrapolation
  for(unsigned int i = 0; i < fExtrapolations.size(); ++i){
    if(sim==SimFlag::kMC){
        SetEventWeight(beamType, runType, fakeDataSet);

        // If we have MC for this beam but not data, SetEventWeight
        // will have set the weight to zero, so don't add it at
        // all. This ensures that the MC beams that have events added
        // to them are the ones that were selected in the RunsToUse
        // registry key *and* have data available
        if(fEventInfo.reco->weight!=0){
          
          //NCBeam::Info info=NCBeam::Info(beamType, thisRunType);
          fExtrapolations[i]->AddEvent(fEventInfo, fakeDataSet,
                                       fileType, beamInfo);
          // Also add it to the RunAll beam
          fExtrapolations[i]->AddEvent(fEventInfo, fakeDataSet,
                                       fileType, infoAll);
          
          if(fShiftedBeams){
            AddShiftedEventToExtrapolations(fEventInfo, fakeDataSet, fileType,
                                            runType, beamInfo);
            AddShiftedEventToExtrapolations(fEventInfo, fakeDataSet, fileType,
                                            runType, infoAll);
          }

        }// end if non-zero weight
    }// end if MC
    else{
      SetEnergies();
      // Data/mock data goes in whatever run it really is...
      fExtrapolations[i]->AddEvent(fEventInfo, fakeDataSet,
                                   fileType, beamInfo);
      // ...and RunAll
      fExtrapolations[i]->AddEvent(fEventInfo, fakeDataSet,
                                   fileType, infoAll);
    }
  }//end loop over extrapolation objects


  if(fMakeDQPlotsSpecial){
    if(fEventInfo.header->dataType == int(SimFlag::kMC)){
      if(fEventInfo.header->detector == int(Detector::kNear)){
        FillDataQualityPlotsSpecial(fDQNearMCTotalSpecial, false);
      }
    }
  }
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::
GetListOfShifts(vector<bool*>& useParameters,
                vector<vector<double> >& adjustedValues) const
{
  assert(useParameters.empty() && adjustedValues.empty());

  // This allows us to find out if a particular variable is in the fit
  NC::CoordinateConverter cc;
  cc.Prepare(fConfig, true);

  // First work out all the combinations of systematics we need. We use
  // every combination of nominal and plus/minus one sigma. ie we should
  // add each event 3^n-1 times (we don't do the completely nominal case)

  // We start with the nominal arrangement - ie no shifts enabled, and all
  // set to their default values
  useParameters.resize(1);
  adjustedValues.resize(1);

  useParameters[0] = new bool[NCType::kNumSystematicParameters];
  adjustedValues[0].resize(NCType::kNumSystematicParameters);
  for(int n = 0; n < NCType::kNumSystematicParameters; ++n){
    useParameters[0][n] = false;
    adjustedValues[0][n] = NCType::kParams[n].defaultVal;
  }

  // For every possible systematic
  for(unsigned int syst = 0; syst < NCType::kNumSystematicParameters; ++syst){
    // We check if we're fitting for it
    if(!cc.IsFit(NCType::EFitParam(syst))) continue;

    // and if so, duplicate all the configurations we have so far in two
    // more modes, for plus/minus one sigma in this new parameter.

    // If the option for simple shifted beams is on then we only base our
    // shifting on the nominal beam -> 2n shifted beams
    const unsigned int N = fSimpleShiftedBeams ? 1 : useParameters.size();

    for(unsigned int n = 0; n < N; ++n){
      bool* up = new bool[NCType::kNumSystematicParameters];
      for(int i = 0; i < NCType::kNumSystematicParameters; ++i)
        up[i] = useParameters[n][i];

      vector<double> av = adjustedValues[n];
      up[syst] = true;
      av[syst] -= NCType::kParams[syst].sigma;
      useParameters.push_back(up);
      adjustedValues.push_back(av);

      // Add twice the shift to cancel out the previous shift the other way

      av[syst] += 2*NCType::kParams[syst].sigma;
      useParameters.push_back(up);
      adjustedValues.push_back(av);
    }
  }

  // Remove the nominal case that we started with
  delete[] useParameters[0];
  useParameters.erase(useParameters.begin());
  adjustedValues.erase(adjustedValues.begin());

  // Check that we didn't lose track somewhere
  assert(useParameters.size() == adjustedValues.size());
}


//-----------------------------------------------------------------------
void NCExtrapolationModule::
AddShiftedEventToExtrapolations(NCEventInfo eventInfo,
                                bool useMCAsData,
                                NCType::EFileType fileType,
                                NC::RunUtil::ERunType runType,
                                NCBeam::Info beamInfo)
{
  // No point adding/shifting data events
  if(eventInfo.header->dataType == NCType::kData || useMCAsData) return;

  Detector::Detector_t det = Detector::Detector_t(fEventInfo.header->detector);
  // On the assumption that the contents of the registry don't change between
  // successive calls to us, then save time by only working out what shifts
  // to do once.
  static vector<bool*> useParameters;
  static vector<vector<double> > adjustedValues;
  static bool once = true;
  if(once){
    once = false;
    GetListOfShifts(useParameters, adjustedValues);
  }


  NC::OscProb::OscPars* noosc = new NC::OscProb::NoOscillations;

  for(unsigned int i = 0; i < useParameters.size(); ++i){
    if(fSingleShiftedBeam >= 0 && int(i) != fSingleShiftedBeam) continue;

    static NCEventInfo backup; // Suspect that construction is expensive

    backup.DeepCopy(&eventInfo);

    fEventInfo = eventInfo;

    beamInfo.SetSystShifts(useParameters[i], &adjustedValues[i][0]);

    if(FinalEventCheck()){
      for(unsigned int j = 0; j < fExtrapolations.size(); ++j){

        // This function applies weight and syst shifts too to the
        // object (that it happens to ~ own)
        fEventInfo.SetEventWeight(useParameters[i], adjustedValues[i],
                                  noosc, 0, runType,
                                  true, false);
        fEventInfo.reco->weight *= GetMCScaleFactor(det, fileType, beamInfo);
        // SetEventWeight will (indirectly) have altered nuEnergy and
        // showerEnergy, so make sure they're set right
        SetEnergies();
        fExtrapolations[j]->AddEvent(eventInfo, useMCAsData, fileType, beamInfo);
      }
    }

    eventInfo.DeepCopy(&backup);
  } // end for i

  delete noosc;
}


//-----------------------------------------------------------------------
std::vector<NCExtrapolationFactory*>& NCExtrapolationModule::
GetExtrapolationFactories()
{
  // Looks like we were falling foul of some variant of the
  // "static initialization order fiasco". Apply the workaround from
  // http://www.parashift.com/c++-faq-lite/ctors.html#faq-10.13
  static std::vector<NCExtrapolationFactory*> facts;
  return facts;
}

---