NCAnalysisModule.cxx

Go to the documentation of this file.

//$Id: NCAnalysisModule.cxx,v 1.45 2008/04/02 20:32:51 brebel Exp $
//
//NCAnalysisModule.cxx
//
//Module for analyzing beam data for nc analysis
//
//B. Rebel 10/2004

#include "NCUtils/Extraction/NCAnalysisModule.h"
#include "NCUtils/NCOscProb.h"

#include "NCUtils/Extraction/NCExtractionADM.h"
#include "NCUtils/Extraction/NCExtractionDP.h"
#include "NCUtils/Extraction/NCExtractionAS.h"
#include "NCUtils/Extraction/NCExtractionTR.h"
#include "NCUtils/Extraction/NCExtractionTRann.h"
#include "NCUtils/Extraction/NCExtractionTO.h"
#include "NCUtils/Extraction/NCExtractionTOm.h"
#include "NCUtils/Extraction/NCExtractionNS.h"
#include "NCUtils/Extraction/NCExtractionKA.h"
#include "NCUtils/Extraction/NCExtractionKAD.h"
#include "NCUtils/Extraction/NCExtractionRO.h"
#include "NCUtils/Extraction/NCExtractionPL.h"
#include "NCUtils/Extraction/NCExtractionRPann.h"
#include "NCUtils/NCType.h"

#include "Validity/VldTimeStamp.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "DataUtil/EnergyCorrections.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "JobControl/JobCommand.h"
#include "Conventions/Detector.h"
#include "Conventions/ReleaseType.h"
#include "Conventions/SimFlag.h"
#include "Conventions/Munits.h"
#include "DataUtil/MCInfo.h"

#include "TSystemFile.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TLegend.h"

#include <cmath>
#include <cassert>
#include <algorithm>

using namespace EnergyCorrections;
//using namespace BeamType;
//using namespace NCType;

static const int kLT05GeV = 0;
static const int kLT1GeV = 1;
static const int kGT1GeV = 2;
static const int kAllGeV = 3;
static const int kNumEventsVsPOT = 4;

CVSID("$Id: NCAnalysisModule.cxx,v 1.45 2008/04/02 20:32:51 brebel Exp $");

// Declare this module to JobControl. Arguments are:
//  (1) The class name 
//  (2) The human-readable name 
//  (3) A short, human-readable description of what the module does
JOBMODULE(NCAnalysisModule, 
          "NCAnalysisModule",
          "A module used for analyzing beam data");

//----------------------------------------------------------------------
NCAnalysisModule::NCAnalysisModule() :
  fFileName("ncccSeparation.root"),
  fTreeName("PAN"),
  fMCPath("mcFiles/*.root"),
  fDataPath("dataFiles/*.root"),
  fNtpFile(0),
  fNtuple(0),
  fStrippedNtpFile(0)
{
  MSG("JobC", Msg::kDebug) << "NCAnalysisModule::Constructor" << endl;  

  fHeaderInfo = new ANtpHeaderInfo();
  fBeamInfo = new ANtpBeamInfo();
  fRecoInfo = new ANtpRecoInfo();
  fEventInfo = new ANtpEventInfoNC();
  fShowerInfo = new ANtpShowerInfoNC();
  fTrackInfo = new ANtpTrackInfoNC();
  fTruthInfo = new ANtpTruthInfoBeam();
  
  fUtils = new NCAnalysisUtils();

  NCType::ConstructBadRuns();

  VldTimeStamp start(2005,5,1,0,0,0);
  VldTimeStamp end(2007,5,1,0,0,0);
  fPOTVsTime = new TH1F("potVsTime", "", 365*2, start.GetSec(), end.GetSec());
  fPOTVsTime->GetXaxis()->SetTimeDisplay(1);
  fPOTVsTime->GetXaxis()->SetTimeFormat("%m/%d");
  fPOTVsTime->SetXTitle("Date");
  fPOTVsTime->SetYTitle("POT#times10^{12}");

  NCType::MakeDataQualityPlots(fDataQualityTotal, "DataTotal");
  NCType::MakeDataQualityPlots(fDataQualityNC,    "DataNC");
  NCType::MakeDataQualityPlots(fDataQualityCC,    "DataCC");

  NCType::MakeDataQualityPlots(fDataQualityTotalLowIntensity, "DataTotalLowIntensity");
  NCType::MakeDataQualityPlots(fDataQualityNCLowIntensity,    "DataNCLowIntensity");
  NCType::MakeDataQualityPlots(fDataQualityCCLowIntensity,    "DataCCLowIntensity");

  NCType::Make2DDataQualityPlots(fDataQuality2DTotal, "DataTotal");
  NCType::Make2DDataQualityPlots(fDataQuality2DNC,    "DataNC");
  NCType::Make2DDataQualityPlots(fDataQuality2DCC,    "DataCC");

  NCType::MakeDataQualityPlots(fDataQualityMCTotal, "MCTotal");
  NCType::MakeDataQualityPlots(fDataQualityMCNC,    "MCNC");
  NCType::MakeDataQualityPlots(fDataQualityMCCC,    "MCCC");

  NCType::MakeDataQualityPlots(fDataQualityMCTotalOsc, "MCTotalOsc");
  NCType::MakeDataQualityPlots(fDataQualityMCNCOsc,    "MCNCOsc");
  NCType::MakeDataQualityPlots(fDataQualityMCCCOsc,    "MCCCOsc");

}

//----------------------------------------------------------------------
NCAnalysisModule::~NCAnalysisModule()
{

  MSG("JobC", Msg::kDebug) << "NCAnalysisModule::Destructor" << endl;

  // free all the memory allocated dynamically
  if (fHeaderInfo) delete fHeaderInfo; 
  if (fBeamInfo) delete fBeamInfo; 
  if (fRecoInfo) delete fRecoInfo;
  if (fEventInfo) delete fEventInfo;
  if (fShowerInfo) delete fShowerInfo;
  if (fTruthInfo) delete fTruthInfo;
  if (fTrackInfo) delete fTrackInfo;
  
  if (fCuts) delete fCuts;
  
  if (fUtils) delete fUtils;
  
}

//----------------------------------------------------------------------
void NCAnalysisModule::BeginJob()
{
  MSG("NCAnalysisModule", Msg::kDebug) << "start BeginJob" << endl;
  
  return;
}

//----------------------------------------------------------------------
//here is where you would loop over the events in the snarl to fill the 
//analysis tree
JobCResult NCAnalysisModule::Ana(const MomNavigator *mom)
{
  MSG("NCAnalysisModule", Msg::kDebug) << "start ana" << endl;
  
  JobCResult result(JobCResult::kPassed);
  
  //get the records from MOM
  assert(mom);
  
  return result;
}

//----------------------------------------------------------------------
void NCAnalysisModule::Help() 
{
  MSG("JobC", Msg::kInfo) 
    << "NCAnalysisModule::Help\n"
    <<"NCAnalysisModule is a module which analyzes beam data."
    << endl;
}

//----------------------------------------------------------------------
void NCAnalysisModule::EndJob() 
{

  MakeuDST(NCType::kMC);
  MakeuDST(NCType::kData);

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::MakeuDST(int dataMC)
{
  TString name = "data_" + fDataSet;
  bool fakeData = false;
  if(dataMC == NCType::kMC) name = "mc";
  else if(dataMC == NCType::kData && fUseMCAsData) fakeData = true;

  //make the chain for the data you want to use
  TChain *chain = new TChain(fTreeName);

  TString path = "";

  //add files to the chain
  if(dataMC == NCType::kMC)
    path = fMCPath;
  else{
    path = fDataPath;
  }//end if doing "data"

  chain->Add(path);

  MSG("NCAnalysisModule", Msg::kInfo) << "adding path " << path 
                                      << " to chain" << endl;

  MSG("NCAnalysisModule", Msg::kInfo) << "Tree Name = " 
                                      << fTreeName << endl;

  chain->SetBranchAddress("header.", &fHeaderInfo);
  chain->SetBranchAddress("beam.",   &fBeamInfo);
  chain->SetBranchAddress("event.",  &fEventInfo);
  chain->SetBranchAddress("shower.", &fShowerInfo);
  chain->SetBranchAddress("track.",  &fTrackInfo);
  if(fUseMCAsData || dataMC == NCType::kMC) 
    chain->SetBranchAddress("truth.",  &fTruthInfo);

  chain->GetEntry(0);

  //figure out which kind of cuts object to use
  MSG("NCAnalysisModule", Msg::kInfo) << "fCutSuite = " << fCutSuite << endl;
  if(fCutSuite == NCType::kCCCuts)
    fCuts = new NCAnalysisCutsCC();
  else if(fCutSuite == NCType::kNCCuts)
    fCuts = new NCAnalysisCutsNC();
  else if(fCutSuite == NCType::kNCCCFidCuts)
    fCuts = new NCAnalysisCutsNCCCFid();
  else if(fCutSuite == NCType::kTOCuts)
    fCuts = new NCAnalysisCutsTO();
  else if(fCutSuite == NCType::kOxCuts)
    fCuts = new NCAnalysisCutsOx();

  fCuts->SetInfoObjects(fHeaderInfo, fBeamInfo, fEventInfo, 
                        fTrackInfo, fShowerInfo, fTruthInfo);
  fUtils->SetInfoObjects(fHeaderInfo, fBeamInfo, fEventInfo, 
                         fTrackInfo, fShowerInfo, 0, 0, fTruthInfo);
  fUtils->SetCutsObject(fCuts);


  //Setting type of hadron production reweighting.
  fUtils->SetMEGAWeightConfig(fMEGAWeightConfig);
  
  // beam type from the BeamAna stuff
  fCuts->SetBeamType(fBeamType);

  //check if should be using all L010z or 185i data
  fCuts->SetUseAllBeams(fUseAllBeams);
  fCuts->SetUseAllL010z(fUseAllL010z);
  fCuts->SetUseAll185i(fUseAll185i);
  
  //fill the vectors holding the extraction and analysis objects
  FillExtractionAndAnalysisVectors(dataMC);

  //save the current directory so as to not attach 
  //files to each other
  TDirectory* saveDir = gDirectory;

  //make a file and the trees to go into the file
  TString fileName = fFileName;
  //data set name is inserted into base file name
  int pos = fileName.Index(".uDST");
  fileName.Insert(pos, "_"+name);

  //indicate you have a stripped file
  TString stripFileName = fileName;
  pos = stripFileName.Index(".root");
  stripFileName.Insert(pos, "_strip");
  TFile *strippedNtpFile = new TFile(stripFileName,"RECREATE");      
  TTree *strippedNtuple = new TTree("uDST", "Stripped AnalysisTree");
  strippedNtuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 
                          64000, 2);
  if(fUseMCAsData || dataMC == NCType::kMC) 
    strippedNtuple->Branch("truth.", "ANtpTruthInfoBeam", &fTruthInfo, 
                            64000, 2);
  strippedNtuple->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
  strippedNtuple->Branch("reco.", "ANtpRecoInfo", &fRecoInfo, 64000, 2);
  strippedNtuple->SetDirectory(strippedNtpFile);
  // construct branch names from class names
  // and add the relevant branches to ntuple files
  for(UInt_t i = 0; i < fExtractions.size(); ++i) {
    string branchName = fExtractions[i]->ClassName();
    branchName = branchName.erase(0,12); // this gives the ending, 
    // "DP", "ADM", etc.
    branchName = "analysis" + branchName + ".";
    
    // add the analysis branches
    MSG("NCAnalysisModule", Msg::kInfo) 
      << "Setting " << branchName << " branch..." << endl;
    strippedNtuple->Branch(branchName.c_str(), "ANtpAnalysisInfo", 
                           &fAnalysisInfos[i], 64000, 2);
  }
  TTree *beamntps = new TTree("beam", "AnalysisTree");
  beamntps->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  beamntps->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
  beamntps->SetDirectory(strippedNtpFile);

  //return to original directory
  saveDir->cd();

  TFile *ntpFile = new TFile(fileName,"RECREATE");      
  TTree *ntuple = new TTree("uDST", "AnalysisTree");
  ntuple->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  if(fUseMCAsData || dataMC == NCType::kMC) 
    ntuple->Branch("truth.", "ANtpTruthInfoBeam", &fTruthInfo, 64000, 2);
  ntuple->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
  ntuple->Branch("reco.", "ANtpRecoInfo", &fRecoInfo, 64000, 2);
  for(UInt_t i = 0; i < fExtractions.size(); ++i) {
    string branchName = fExtractions[i]->ClassName();
    branchName = branchName.erase(0,12); // this gives the ending, 
    // "DP", "ADM", etc.
    branchName = "analysis" + branchName + ".";
    
    // add the analysis branches
    MSG("NCAnalysisModule", Msg::kInfo) 
      << "Setting " << branchName << " branch..." << endl;
    ntuple->Branch(branchName.c_str(), "ANtpAnalysisInfo", 
                   &fAnalysisInfos[i], 64000, 2);
  }
  ntuple->SetDirectory(ntpFile);
  TTree *beamntp = new TTree("beam", "AnalysisTree");
  beamntp->Branch("header.", "ANtpHeaderInfo", &fHeaderInfo, 64000, 2);
  beamntp->Branch("beam.", "ANtpBeamInfo", &fBeamInfo, 64000, 2);
  beamntp->SetDirectory(ntpFile);

  //return to original directory
  saveDir->cd();

  fStrippedNtpFile = strippedNtpFile;
  fNtpFile = ntpFile;

  MSG("NCAnalysisModule", Msg::kInfo) << "filenames " << fileName 
                                      << " " << stripFileName << endl;

  fNtuple = ntuple;
  fStrippedNtuple = strippedNtuple;
  fBeamNtuple = beamntp;
  fBeamNtupleStp = beamntps;

  //clear the data POT map
  fDataPOTs.clear();

  //next lines do all the work - thank goodness for methods
  //only prepare if dealing with MC - data pdfs get filled 
  //in the extraction stage because that is the only time you 
  //can say if they are NC or CC (ie no truth information for data)
  if(dataMC == NCType::kMC) PrepareForExtraction(chain);
  ExtractNCCC(chain, fakeData);

  int ntupleEntries = fNtuple->GetEntries();

  //write out the histograms to a file.  set the normalization 
  //if data and mc pot total is > 0.
  double mcPOT = GetTotalPOTCount(NCType::kMC);
  double dataPOT = GetTotalPOTCount(NCType::kData);
  for(int k = 0; k < static_cast<int>(fExtractions.size()); ++k){
    if(dataPOT > 0. && mcPOT > 0.) 
      fExtractions[k]->SetNormalization(dataPOT/mcPOT);
  }
  
  WriteFile(fNtpFile, fNtuple, fBeamNtuple, dataMC);
  WriteFile(fStrippedNtpFile, fStrippedNtuple, fBeamNtupleStp, dataMC);

  fNtuple = 0;
  fStrippedNtuple = 0;
  fNtpFile = 0;
  fStrippedNtpFile = 0;

  //check to see if there was anything in the ntuples, 
  //if not remove the file from the system
  if(ntupleEntries < 1){
    TSystemFile ntp(fileName, "./");
    TSystemFile strip(stripFileName, "./");

    ntp.Delete();
    strip.Delete();
  }

  //if(strippedNtpFile) delete strippedNtpFile;
  //if(ntpFile) delete ntpFile;

  delete chain;

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::FillExtractionAndAnalysisVectors(int dataMC)
{
  MSG("NCAnalysisModule", Msg::kInfo) << "SimFlag =  " 
                                      << fHeaderInfo->dataType << endl;
  Detector::Detector_t det = Detector::kNear;
  if(fHeaderInfo->detector == (int)Detector::kFar) det = Detector::kFar;

  if(dataMC == NCType::kMC){
    fExtractions.clear();
    fTOExtraction = -1;
    fTOmExtraction = -1;

    if(fDoDPExtraction) 
      fExtractions.push_back(new NCExtractionDP(fCuts, fUtils, this->GetConfig() ));
    if(fDoADMExtraction){ 
      fExtractions.push_back(new NCExtractionADM(fCuts, fUtils, this->GetConfig() ));
      fResultsADM = new NCExtractionResultsADM(det);
    }
    if(fDoTRExtraction) 
      fExtractions.push_back(new NCExtractionTR(fCuts, fUtils, this->GetConfig() ));
    if(fDoTOExtraction){
      fTOExtraction = fExtractions.size();
      fExtractions.push_back(new NCExtractionTO(fCuts, fUtils, this->GetConfig() ));
    }
    if(fDoTOmExtraction){
      fTOmExtraction = fExtractions.size();
      fExtractions.push_back(new NCExtractionTOm(fCuts, fUtils, this->GetConfig() ));
    }
    if(fDoTRannExtraction) 
      fExtractions.push_back(new NCExtractionTRann(fCuts, fUtils, this->GetConfig() ));
    if(fDoRPannExtraction) 
      fExtractions.push_back(new NCExtractionRPann(fCuts, fUtils, this->GetConfig() ));
    if(fDoASExtraction)
        fExtractions.push_back(new NCExtractionAS(fCuts, fUtils, this->GetConfig() ));
    if(fDoNSExtraction) 
      fExtractions.push_back(new NCExtractionNS(fCuts, fUtils, this->GetConfig() ));
    if(fDoKAExtraction) 
      fExtractions.push_back(new NCExtractionKA(fCuts, fUtils, this->GetConfig() ));
    if(fDoKADExtraction) 
      fExtractions.push_back(new NCExtractionKAD(fCuts, fUtils, this->GetConfig() ));
    if(fDoROExtraction){
      fExtractions.push_back(new NCExtractionRO(fCuts, fUtils, this->GetConfig() ));
    }
    if(fDoPLExtraction) 
      fExtractions.push_back(new NCExtractionPL(fCuts, fUtils, this->GetConfig() ));
      
    //fill the vector of analysis info objects (one per extraction) 
    MSG("NCAnalysisModule", Msg::kInfo) 
      << "Using " << fExtractions.size() << " extraction  methods."
      << " Now adding the corresponding ANtpAnalysisInfo objects..."
      << endl;
    
    fAnalysisInfos.clear();
    for (uint i = 0; i < fExtractions.size(); ++i) 
      fAnalysisInfos.push_back(new ANtpAnalysisInfo());  


    int numEx = (int)fAnalysisInfos.size();
    TString dataType[2] = {"MC", "Data"};
    TString label =  "";
    for(int i = 0; i < 2; ++i){
      fNCLikeNumerator.push_back(new TH2D("ncLikeNumerator"+dataType[i], "", 
                                          numEx, 0., 1.*numEx, numEx, 0., 1.*numEx));
      fNCLikeDenominator.push_back(new TH2D("ncLikeDenominator"+dataType[i], "", 
                                            numEx, 0., 1.*numEx, numEx, 0., 1.*numEx));
      fNCLikeCorrelation.push_back(new TH2D("ncLikeCorrelation"+dataType[i], "", 
                                            numEx, 0., 1.*numEx, numEx, 0., 1.*numEx));
      fNCLikeHist.push_back(new TH1D("ncLike"+dataType[i], "", numEx+1, 0., 1.*(numEx+1)));
      fNCLikeHist[i]->SetXTitle("# Methods");
      fNCLikeHist[i]->SetYTitle("# Events");

      fCCLikeNumerator.push_back(new TH2D("ccLikeNumerator"+dataType[i], "", 
                                          numEx, 0., 1.*numEx, numEx, 0., 1.*numEx));
      fCCLikeDenominator.push_back(new TH2D("ccLikeDenominator"+dataType[i], "", 
                                            numEx, 0., 1.*numEx, numEx, 0., 1.*numEx));
      fCCLikeCorrelation.push_back(new TH2D("ccLikeCorrelation"+dataType[i], "", 
                                            numEx, 0., 1.*numEx, numEx, 0., 1.*numEx));
      fCCLikeHist.push_back(new TH1D("ccLike"+dataType[i], "", numEx+1, 0., 1.*(numEx+1)));
      fCCLikeHist[i]->SetXTitle("# Methods");
      fCCLikeHist[i]->SetYTitle("# Events");

      if(i == 0){
        fNCLikeNumeratorTrue = new TH2D("ncLikeNumeratorTrue", "", 
                                        numEx, 0., 1.*numEx, numEx, 0., 1.*numEx);
        fNCLikeDenominatorTrue = new TH2D("ncLikeDenominatorTrue", "", 
                                          numEx, 0., 1.*numEx, numEx, 0., 1.*numEx);
        fNCLikeCorrelationTrue  = new TH2D("ncLikeCorrelationTrue"+dataType[i], "", 
                                           numEx, 0., 1.*numEx, numEx, 0., 1.*numEx);

        fCCLikeNumeratorTrue = new TH2D("ccLikeNumeratorTrue", "", 
                                        numEx, 0., 1.*numEx, numEx, 0., 1.*numEx);
        fCCLikeDenominatorTrue = new TH2D("ccLikeDenominatorTrue", "", 
                                          numEx, 0., 1.*numEx, numEx, 0., 1.*numEx);
        fCCLikeCorrelationTrue  = new TH2D("ccLikeCorrelationTrue"+dataType[i], "", 
                                           numEx, 0., 1.*numEx, numEx, 0., 1.*numEx);
      }

      for(int x = 0; x < numEx; ++x){
        //set the axis labels
        label = fExtractions[x]->ClassName();
        label.Remove(0,12);
        fNCLikeDenominator[i]->GetXaxis()->SetBinLabel(x+1, label);
        fNCLikeNumerator[i]->GetXaxis()->SetBinLabel(x+1, label);
        fNCLikeCorrelation[i]->GetXaxis()->SetBinLabel(x+1, label);
        fNCLikeDenominator[i]->GetYaxis()->SetBinLabel(x+1, label);
        fNCLikeNumerator[i]->GetYaxis()->SetBinLabel(x+1, label);
        fNCLikeCorrelation[i]->GetYaxis()->SetBinLabel(x+1, label);

        fCCLikeDenominator[i]->GetXaxis()->SetBinLabel(x+1, label);
        fCCLikeNumerator[i]->GetXaxis()->SetBinLabel(x+1, label);
        fCCLikeCorrelation[i]->GetXaxis()->SetBinLabel(x+1, label);
        fCCLikeDenominator[i]->GetYaxis()->SetBinLabel(x+1, label);
        fCCLikeNumerator[i]->GetYaxis()->SetBinLabel(x+1, label);
        fCCLikeCorrelation[i]->GetYaxis()->SetBinLabel(x+1, label);
        if(i == 0){
          fNCLikeDenominatorTrue->GetXaxis()->SetBinLabel(x+1, label);
          fNCLikeNumeratorTrue->GetXaxis()->SetBinLabel(x+1, label);
          fNCLikeCorrelationTrue->GetXaxis()->SetBinLabel(x+1, label);
          fNCLikeDenominatorTrue->GetYaxis()->SetBinLabel(x+1, label);
          fNCLikeNumeratorTrue->GetYaxis()->SetBinLabel(x+1, label);
          fNCLikeCorrelationTrue->GetYaxis()->SetBinLabel(x+1, label);

          fCCLikeDenominatorTrue->GetXaxis()->SetBinLabel(x+1, label);
          fCCLikeNumeratorTrue->GetXaxis()->SetBinLabel(x+1, label);
          fCCLikeCorrelationTrue->GetXaxis()->SetBinLabel(x+1, label);
          fCCLikeDenominatorTrue->GetYaxis()->SetBinLabel(x+1, label);
          fCCLikeNumeratorTrue->GetYaxis()->SetBinLabel(x+1, label);
          fCCLikeCorrelationTrue->GetYaxis()->SetBinLabel(x+1, label);
        }
      }//end loop over x axis labels


    }//end initialization of correlation matricies

  }
  else{

    if(fExtractions.size() < 1) 
      MSG("NCAnalysisModule", Msg::kWarning) 
        << "trying to reset extraction objects"
        << " in vector with no entries." 
        << " must fill mc first" << endl;

    //reset the data histograms in each of the extraction objects
    //because you are changing data sets
    for(uint i = 0; i < fExtractions.size(); ++i) 
      fExtractions[i]->ResetHistograms(NCType::kData);

    //reset the correlation histograms too
    fNCLikeHist[NCType::kData]->Reset();
    fNCLikeNumerator[NCType::kData]->Reset();
    fNCLikeDenominator[NCType::kData]->Reset();
    fNCLikeCorrelation[NCType::kData]->Reset();
    fCCLikeHist[NCType::kData]->Reset();
    fCCLikeNumerator[NCType::kData]->Reset();
    fCCLikeDenominator[NCType::kData]->Reset();
    fCCLikeCorrelation[NCType::kData]->Reset();
  }

  return;
}

//----------------------------------------------------------------------
const Registry& NCAnalysisModule::DefaultConfig() const
{
        
  int itrue = 1;  // work around for lack of bool in registry
  int ifalse = 0; // work around for lack of bool in registry
  
  static Registry r;
  
  r.UnLockValues();
  
  r.Set("FileName",               "ncccSeparation.root");
  r.Set("TreeName",               "antp");
  r.Set("MCPath",                 "mcFiles/*.root");
  r.Set("DataPath",               "dataFiles/*.root");
  r.Set("PDFPath",                "DP_histos_far.root");
  r.Set("FillNtuples",            ifalse);
  r.Set("ReadPDFs",               ifalse);
  r.Set("UseMCAsData",            itrue);
  r.Set("BeamType",               "L010z185i");
  r.Set("UseAllBeams",            ifalse);
  r.Set("UseAllL010z",            ifalse);
  r.Set("UseAll185i",             ifalse);
  r.Set("DoDPExtraction",         itrue); 
  r.Set("DoADMExtraction",        ifalse); 
  r.Set("DoTOExtraction",         ifalse);
  r.Set("DoTOmExtraction",         ifalse);
  r.Set("DoTRExtraction",         ifalse);
  r.Set("DoTRannExtraction",      ifalse);
  r.Set("DoASExtraction",         ifalse);
  r.Set("DoNSExtraction",         ifalse);
  r.Set("DoKAExtraction",         ifalse);
  r.Set("DoKADExtraction",        ifalse);
  r.Set("DoROExtraction",         ifalse);
  r.Set("DoRPannExtraction",      ifalse);
  r.Set("DoPLExtraction",         ifalse);
  r.Set("RPAnnUseLowETrain",      ifalse);
  r.Set("NDSubRunToUse",          6);
  r.Set("FileCountLimit",        1000000000);
  r.Set("DataSet",                "/");
  r.Set("CutSuite",               NCType::kCCCuts);
  r.Set("MEGAWeightConfig",       "Boston");
  r.Set("AnnNSPrior",              1);
  r.Set("AnnNSFiducial",           1);
  r.Set("CountChainEntries",       itrue);
  r.LockValues();
  return r;
}

//----------------------------------------------------------------------
void NCAnalysisModule::Config(const Registry& r)
{
  int         tmpb;  // a temp bool. See comment under DefaultConfig...
  int         tmpi;  // a temp int.
  double      tmpd;  // a temp double.
  const char* tmps;  // a temp string
  
  if (r.Get("FileName",         tmps)) fFileName          = tmps;
  if (r.Get("TreeName",         tmps)) fTreeName          = tmps;
  if (r.Get("MCPath",           tmps)) fMCPath            = tmps;
  if (r.Get("DataPath",         tmps)) fDataPath          = tmps;
  if (r.Get("FillNtuples",      tmpb)) fFillNtuples       = tmpb;
  if (r.Get("UseMCAsData",      tmpb)) fUseMCAsData       = tmpb;
  if (r.Get("BeamType",         tmps)) fBeamType          = BeamType::TagToEnum(tmps);
  if (r.Get("UseAllBeams",      tmpb)) fUseAllBeams       = tmpb;
  if (r.Get("UseAllL010z",      tmpb)) fUseAllL010z       = tmpb;
  if (r.Get("UseAll185i",       tmpb)) fUseAll185i        = tmpb;
  if (r.Get("DoDPExtraction",   tmpb)) fDoDPExtraction    = tmpb;
  if (r.Get("DoADMExtraction",  tmpb)) fDoADMExtraction   = tmpb;
  if (r.Get("DoTOExtraction"   ,tmpb)) fDoTOExtraction    = tmpb;
  if (r.Get("DoTOmExtraction"  ,tmpb)) fDoTOmExtraction   = tmpb;
  if (r.Get("DoTRExtraction"   ,tmpb)) fDoTRExtraction    = tmpb;
  if (r.Get("DoTRannExtraction",tmpb)) fDoTRannExtraction = tmpb;
  if (r.Get("DoASExtraction",   tmpb)) fDoASExtraction    = tmpb;
  if (r.Get("DoNSExtraction",   tmpb)) fDoNSExtraction    = tmpb;
  if (r.Get("DoKAExtraction",   tmpb)) fDoKAExtraction    = tmpb;
  if (r.Get("DoKADExtraction",  tmpb)) fDoKADExtraction   = tmpb;
  if (r.Get("DoROExtraction",   tmpb)) fDoROExtraction    = tmpb;
  if (r.Get("DoPLExtraction",   tmpb)) fDoPLExtraction    = tmpb;
  if (r.Get("DoRPannExtraction",tmpb)) fDoRPannExtraction = tmpb;
  if (r.Get("CutSuite",         tmpi)) fCutSuite          = tmpi;
  if (r.Get("MEGAWeightConfig", tmps)) fMEGAWeightConfig  = tmps;
  if (r.Get("NDSubRunToUse",    tmpi)) fNDSubRunToUse     = tmpi;
  if (r.Get("FileCountLimit",  tmpi)) fFileCountLimit   = tmpi;
  if (r.Get("CountChainEntries",tmpb)) fCountChainEntries = tmpb;
  if (r.Get("DataSet",          tmps)){
    fDataSet           = tmps;
    fDataSet.Remove(TString::kLeading, '/');
  }

  tmpd = 0.;

  return;
}

//-----------------------------------------------------------------------
//preparation for extracting NC/CC - can be done with either data or MC
//useful for comparing data and MC pdfs later
//preparation step prepares all extraction objects at once
void NCAnalysisModule::PrepareForExtraction(TChain *chain)
{
  MSG("NCAnalysisModule", Msg::kInfo) << "PrepareForExtraction " 
                                      << chain->GetName() << endl;      

  //do whatever needs to be done for each extraction object     
  for(uint i = 0; i < fExtractions.size(); ++i)
    fExtractions[i]->Prepare();

  //*******************
  //the following lines of code are deprecated, but could be 
  //resurrected if any one decides to put training steps into 
  //the extraction objects
  //*******************

//   int ctr = 0;
  
//   int tot = -999;
//   if(fCountChainEntries) tot = chain->GetEntries();

//   //loop over the entries in the tree
//   while( chain->GetEntry(ctr) && ctr < fEventCountLimit){
//     if(ctr % 20000 == 0)
//       MSG("NCAnalysisModule", Msg::kInfo) << "On " << ctr << " / " 
//                                           << tot << endl;
//     ++ctr;
    
//     fBeamInfo->protonIntensity = fBeamInfo->tortgt;

//     //make sure the mc beam type is the correct one
//     if(fHeaderInfo->dataType == (int)SimFlag::kMC)
//       fBeamInfo->beamType = BeamType::ToZarko(fBeamType);
    
//     //set the info objects to use
//     fCuts->SetInfoObjects(fHeaderInfo, fBeamInfo, fEventInfo, 
//                        fTrackInfo, fShowerInfo, fTruthInfo);
//     fUtils->SetInfoObjects(fHeaderInfo, fBeamInfo, fEventInfo, 
//                         fTrackInfo, fShowerInfo, 0, 0, fTruthInfo);

// //     MAXMSG("NCAnalysisModule", Msg::kInfo, 200) 
// //       << "run = " 
// //       << fHeaderInfo->run 
// //       << " snarl = " 
// //       << fHeaderInfo->snarl 
// //       << " good snarl " 
// //       << fCuts->IsGoodBeamSnarl() 
// //       << " good event " 
// //       << fCuts->IsGoodBeamEvent()
// //       << " in fiducial volume " 
// //       << fCuts->InBeamFiducialVolume()
// //       << endl;    

//     //check that the snarl for this event is a good one 
//     //and that the event is also good
//     //throw out events outside the fiducial volume    
//     if(!fCuts->IsGoodBeamSnarl() 
//        || !fCuts->IsGoodBeamEvent()
//        || !fCuts->InBeamFiducialVolume()
//       )continue;

//     //*******************
//     //put call to do training steps or fill histograms for training steps here
//     //*******************

//   }//end loop over chain


  return;
}

//-----------------------------------------------------------------------
//extraction can be done with either data or MC
//extraction step does all extraction algorithms at once
void NCAnalysisModule::ExtractNCCC(TChain *chain, bool fakeData)
{
  
  MSG("NCAnalysisModule", Msg::kInfo) << "ExtractNCCC" << endl; 
  
  int ctr = 0;
  int passSnarl = 0;
  int passEvt = 0;
  int passFid = 0;
  bool passCut1 = false;
  bool passCut2 = false;
  bool passCut3 = false;
  bool newSnarl = false;

  int fileCtr = 0;
  int prevRun = -1;
  int prevSubRun = -1;
  int tot = fFileCountLimit;
  if(fCountChainEntries) tot = chain->GetEntries();
  
  chain->GetEntry(0);
  //Get The proper POT/Snarl value
  Detector::Detector_t  detType = Detector::kUnknown;
  if(fHeaderInfo->detector == static_cast<int>(Detector::kFar)) 
    detType = Detector::kFar;  
  if(fHeaderInfo->detector == static_cast<int>(Detector::kNear)) 
    detType = Detector::kNear;  

  int dataMC = NCType::kMC;
  if(fHeaderInfo->dataType == (int)SimFlag::kMC){
    const char* mcString = (fHeaderInfo->softVersion).Data();
    ReleaseType::Release_t mcType=ReleaseType::StringToType(mcString);
    if (mcType==ReleaseType::kUnknown) {
        MSG("NCAnalysisModule",Msg::kWarning) 
            << "Can't figure out MC version, defaulting to carrot.\n";
        mcType =  ReleaseType::kCarrot;
    }
    
    fMCPOTPerSnarl = MCInfo::GetMCPoT(detType,fBeamType,mcType);
    if(fHeaderInfo->dataType == (int)SimFlag::kMC)  
      MSG("NCAnalysisModule",Msg::kInfo) << "MC POT per snarl: "
                                         << fMCPOTPerSnarl << endl;
  }
  if(fakeData || fHeaderInfo->dataType == (int)SimFlag::kData){
    dataMC = NCType::kData;
    tot = 100000;
  }

  double multiplicityNC[kNumEventsVsPOT] = {0.};
  double multiplicityCC[kNumEventsVsPOT] = {0.};
  double prevPOT = 0.;
  //loop over the entries in the tree
  while( chain->GetEntry(ctr) && fileCtr < tot){
    if(ctr % 20000 == 0) 
      MSG("NCAnalysisModule", Msg::kInfo)
        << "On " << ctr << " / " << tot << " file " << fileCtr << endl;
    
    ++ctr;

    if(fHeaderInfo->run != prevRun || fHeaderInfo->subRun != prevSubRun){
      prevRun = fHeaderInfo->run;
      prevSubRun = fHeaderInfo->subRun;
      ++fileCtr;
    }

    if(fHeaderInfo->newSnarl > 0 
       && fHeaderInfo->dataType == (int)SimFlag::kData 
       && detType == Detector::kNear){
      for(uint i = 0; i < fDataQuality2DTotal.size(); ++i){
        fDataQuality2DTotal[i]->Fill(prevPOT, multiplicityCC[i]+multiplicityNC[i]);
        fDataQuality2DNC[i]->Fill(prevPOT, multiplicityNC[i]);
        multiplicityNC[i] = 0.;
        fDataQuality2DCC[i]->Fill(prevPOT, multiplicityCC[i]);
        multiplicityCC[i] = 0.;
      }
      prevPOT = 0.;
    }

    if(fHeaderInfo->dataType == (int)SimFlag::kData){
      fTruthInfo = 0;
      //cascade which pot value you use for accounting
//       fBeamInfo->protonIntensity = fBeamInfo->trtgtd;
//       if(fBeamInfo->trtgtd <= 0.0) fBeamInfo->protonIntensity = fBeamInfo->tortgt;
//       else if(fBeamInfo->tortgt <= 0.0) fBeamInfo->protonIntensity = fBeamInfo->tor101;
      fBeamInfo->protonIntensity = fBeamInfo->tortgt;
      if(fBeamInfo->tortgt <= 0.0) fBeamInfo->protonIntensity = fBeamInfo->trtgtd;
      else if(fBeamInfo->trtgtd <= 0.0) fBeamInfo->protonIntensity = fBeamInfo->tor101;
      //dont want to count POT or events from bad runs
      if(!NCType::IsNDRunGood(fHeaderInfo->run, fHeaderInfo->subRun)) continue;
    }
    //make sure the mc beam type is the correct one
    else if(fHeaderInfo->dataType == (int)SimFlag::kMC){
      fBeamInfo->beamType = BeamType::ToZarko(fBeamType);
      fBeamInfo->protonIntensity = fMCPOTPerSnarl;
    }
    prevPOT = fBeamInfo->protonIntensity;

    //set the info objects to use
    fCuts->SetInfoObjects(fHeaderInfo, fBeamInfo, fEventInfo, 
                          fTrackInfo, fShowerInfo, fTruthInfo);
    fUtils->SetInfoObjects(fHeaderInfo, fBeamInfo, fEventInfo, 
                           fTrackInfo, fShowerInfo, 0, 0, fTruthInfo);

    newSnarl = fCuts->IsNewSnarl();

    if (passCut1 = (fCuts->IsGoodBeamSnarl()) ) passSnarl++;
    if (passCut2 = (fCuts->IsGoodBeamEvent() && passCut1) ) passEvt++;
    if (passCut3 = (fCuts->InBeamFiducialVolume() && passCut2)) passFid++;

    if( newSnarl && passCut1 ){
      CountPOTs(fakeData);
      AddToNtuple(newSnarl);
    }

    if (!passCut2) continue;

    FillRecoInfo(fRecoInfo, fHeaderInfo, 
                 fEventInfo, fShowerInfo, 
                 fTrackInfo, fBeamInfo, fCuts,
                 fUtils, fBeamInfo->beamType);

    for(uint i = 0; i < fExtractions.size(); ++i)
      fExtractions[i]->DoExtraction(fAnalysisInfos[i], fRecoInfo,
                                    fHeaderInfo, fEventInfo, 
                                    fShowerInfo, fTrackInfo, 
                                    fTruthInfo, fBeamInfo->beamType);

    AddToNtuple();

    if(fHeaderInfo->dataType == (int)SimFlag::kData 
       && detType == Detector::kNear 
       && fRecoInfo->inFiducialVolume > 0
       && fTOmExtraction >= 0){
      if(fAnalysisInfos[fTOmExtraction]->isNC > 0 
         && fRecoInfo->isCleanHighMultSnarl > 0){
        if(fRecoInfo->nuEnergy < 0.5) multiplicityNC[kLT05GeV] += 1.; 
        else if(fRecoInfo->nuEnergy < 1.0) multiplicityNC[kLT1GeV] += 1.; 
        else if(fRecoInfo->nuEnergy > 1.0) multiplicityNC[kGT1GeV] += 1.; 
        multiplicityNC[kAllGeV] += 1.;
      }
      if(fAnalysisInfos[fTOmExtraction]->isCC > 0){
        if(fRecoInfo->nuEnergy < 0.5) multiplicityCC[kLT05GeV] += 1.; 
        else if(fRecoInfo->nuEnergy < 1.0) multiplicityCC[kLT1GeV] += 1.; 
        else if(fRecoInfo->nuEnergy > 1.0) multiplicityCC[kGT1GeV] += 1.; 
        multiplicityCC[kAllGeV] += 1.;
      }
    }

    //find correlations between the methods
    if( fCuts->PassesFinalSelection(fRecoInfo) ) FindCorrelations(dataMC);

  }//end loop over chain

  //looped over all the data, now do the fit for the ADM extraction
  //DoADMFit();

  // print some details about cuts
  MSG("NCAnalysisModule", Msg::kInfo) 
    << passSnarl << "/" << ctr << " passed IsGoodBeamSnarl()\n" 
    << passEvt << "/" << passSnarl << " passed IsGoodBeamEvent()\n"
    << passFid << "/" << passEvt << " passed fiducial cuts" 
    << endl;    
  
  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::AddToNtuple(bool beam)
{

  //fill the appropriate uDST ntpule

  if(beam){
    if(fNtpFile->GetSize() < 1850400000)
      fBeamNtuple->Fill();
    if(fStrippedNtpFile->GetSize() < 1850400000)
      fBeamNtupleStp->Fill();
  }

  else{
    //keep filling the ntuple if the total file size < 1.8 Gb
    if(fNtpFile->GetSize() < 1850400000)
      fNtuple->Fill();
    else{ 
      MAXMSG("NCAnalysisModule", Msg::kWarning,1) 
        << "Size limit exceeded no longer filling full ntuple" << endl;
    }

    if( fCuts->PassesFinalSelection(fRecoInfo) ){
      if(fStrippedNtpFile->GetSize() < 1850400000)
        fStrippedNtuple->Fill();
      else{
        MAXMSG("NCAnalysisModule", Msg::kWarning,1) 
          << "Size limit exceeded no longer filling stripped ntuple"
          << endl;
      }
      //fill up the data quality histograms now
      ANtpAnalysisInfo *ana = 0;
      if(fTOmExtraction >= 0){
        ana = fAnalysisInfos[fTOmExtraction];
        if( fUtils->FinalEventCheck(ana, fRecoInfo) ){
          if(fHeaderInfo->dataType == (int)SimFlag::kData){
            FillDataQualityPlots(fDataQualityTotal, NCType::kUnknown, false);
            if(ana->isNC > 0)
              FillDataQualityPlots(fDataQualityNC, NCType::kNC, false);
            else if(ana->isCC > 0)
              FillDataQualityPlots(fDataQualityCC, NCType::kCC, false);
            if(fBeamInfo->protonIntensity < 4.){
              FillDataQualityPlots(fDataQualityTotalLowIntensity, NCType::kUnknown, false);
              if(ana->isNC > 0)
                FillDataQualityPlots(fDataQualityNCLowIntensity, NCType::kNC, false);
              else if(ana->isCC > 0)
                FillDataQualityPlots(fDataQualityCCLowIntensity, NCType::kCC, false);
            }//end if low intensity
          }//end if data
          else if(fHeaderInfo->dataType == (int)SimFlag::kMC){
            FillDataQualityPlots(fDataQualityMCTotal, NCType::kUnknown, false);
            FillDataQualityPlots(fDataQualityMCTotalOsc, NCType::kUnknown, true);
            if(fHeaderInfo->detector == (int)Detector::kFar){
              if(ana->isNC > 0){
                FillDataQualityPlots(fDataQualityMCNC, NCType::kNC, false);
                FillDataQualityPlots(fDataQualityMCNCOsc, NCType::kNC, true);
              }
              else if(ana->isCC > 0){
                FillDataQualityPlots(fDataQualityMCCC, NCType::kCC, false);
                FillDataQualityPlots(fDataQualityMCCCOsc, NCType::kCC, true);
              }
            }//end if far detector
            else{
              FillDataQualityPlots(fDataQualityMCNC, NCType::kNC, false);
              FillDataQualityPlots(fDataQualityMCNCOsc, NCType::kNC, true);
              FillDataQualityPlots(fDataQualityMCCC, NCType::kCC, false);
              FillDataQualityPlots(fDataQualityMCCCOsc, NCType::kCC, true);
            }//end if near detector 
          }//end if MC
        }//end if passes final event checks
      }//end if TO extraction is there
    }//end if new snarl or passes final selection cuts
  } // else of [if(beam)]
  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::FillDataQualityPlots(std::vector<TH1D *> &dq, 
                                            int nccc,
                                            bool osc)
{
  //dont bother with events outside the fiducial volume
  if(fRecoInfo->inFiducialVolume < 1) return;
  
  double weight = fRecoInfo->weight;
  if(NCType::FindRunType(fHeaderInfo) == NCType::kRunII){
    weight = fRecoInfo->weightRunII;
  }

  double trackExtension = -1.*fRecoInfo->trackExtension;
  double month = (1.*fHeaderInfo->month-4.5)+12.*(fHeaderInfo->year-2005);
  double showerEnergy = fRecoInfo->showerEnergyNC;
  double energy = fRecoInfo->nuEnergyNC;
  if(nccc == NCType::kNC){

    //in the ND the NC and CC histograms for the MC are 
    //filled based on truth, not reco
    if(fHeaderInfo->detector == (int)Detector::kNear
       && fHeaderInfo->dataType == (int)SimFlag::kMC
       && fTruthInfo->interactionType == NCType::kCC) return;
  }
  else if(nccc == NCType::kCC){
    showerEnergy = fRecoInfo->showerEnergyCC;
    energy = fRecoInfo->nuEnergyCC;

    //in the ND the NC and CC histograms for the MC are 
    //filled based on truth, not reco
    if(fHeaderInfo->detector == (int)Detector::kNear
       && fHeaderInfo->dataType == (int)SimFlag::kMC
       && fTruthInfo->interactionType == NCType::kNC) return;
  }

  VldTimeStamp ts(fHeaderInfo->year, fHeaderInfo->month, fHeaderInfo->day,
                  fHeaderInfo->hour, fHeaderInfo->minute, (int)fHeaderInfo->second);

  double deltaT = 0.;
  if(fHeaderInfo->detector == (int)Detector::kFar){
    if(fHeaderInfo->dataType == (int)SimFlag::kMC && osc){
      vector<double> oscPars;
      oscPars.push_back(0.25*TMath::Pi());
      oscPars.push_back(2.38e-3);
      oscPars.push_back(0.);
      oscPars.push_back(0.61);//SNO+KamLAND
      oscPars.push_back(7.59e-5);//SNO+KamLAND
      oscPars.push_back(2.7);//g/cc standard rock
      oscPars.push_back(0.);//why not?
      oscPars.push_back(1.);//again, why not?
      
      weight *= fUtils->Find3FlavorProbability(NCType::kNuMuToNuMu,
                                               fTruthInfo->interactionType,
                                               NCType::kBaseLineFar,
                                               fTruthInfo->nuEnergy,
                                               oscPars);
    }//end if MC
    else if(fHeaderInfo->dataType == (int)SimFlag::kData){
      deltaT = 1.e6*(fEventInfo->stripTime1st-1.e-9*fBeamInfo->nearestNSToSpill);
      fVtxX.push_back(fRecoInfo->vtxX);
      fVtxY.push_back(fRecoInfo->vtxY);
    }
  }//end if far 

  if(fHeaderInfo->detector == (int)Detector::kNear
     && fRecoInfo->isCleanHighMultSnarl < 1 ) return;

  dq[NCType::kDQVtxX]->Fill(fRecoInfo->vtxX, weight);
  dq[NCType::kDQVtxY]->Fill(fRecoInfo->vtxY, weight);
  dq[NCType::kDQVtxZ]->Fill(fRecoInfo->vtxZ, weight);
  dq[NCType::kDQEventLength]->Fill(fRecoInfo->eventLength, weight);
  dq[NCType::kDQNumTracks]->Fill(fRecoInfo->numberTracks, weight);
  dq[NCType::kDQShowerE]->Fill(showerEnergy, weight);
  dq[NCType::kDQTrackMomentum]->Fill(fRecoInfo->trackMomentum, weight);
  dq[NCType::kDQEventsVsTime]->Fill(month);
  dq[NCType::kDQRadialVtx]->Fill(TMath::Sqrt(TMath::Power(fRecoInfo->vtxX, 2.)
                                             + TMath::Power(fRecoInfo->vtxY, 2.)), weight);
  dq[NCType::kDQRadialSqrVtx]->Fill(TMath::Power(fRecoInfo->vtxX, 2.)
                                    + TMath::Power(fRecoInfo->vtxY, 2.), weight);
  if(fHeaderInfo->newSnarl > 0)
    dq[NCType::kDQEventsPerSnarl]->Fill(fHeaderInfo->events);
  dq[NCType::kDQPHFraction]->Fill(fEventInfo->pulseHeight/fHeaderInfo->snarlPulseHeight, weight);
  dq[NCType::kDQDeltaTSpill]->Fill(deltaT, weight);
  dq[NCType::kDQAvStripsPerPlane]->Fill(fEventInfo->totalStrips/TMath::Power(fEventInfo->planes,2.), weight);
  if(fEventInfo->showers > 0.)
    dq[NCType::kDQTransverseRMS]->Fill(TMath::Sqrt(TMath::Power(fShowerInfo->transverseRMSU,2.)
                                                   +TMath::Power(fShowerInfo->transverseRMSV,2.)), 
                                       weight);

  if(fEventInfo->tracks > 0.){
    dq[NCType::kDQTrackDCosZVtx]->Fill(fTrackInfo->dcosZVtx, weight);
    dq[NCType::kDQTrackDCosYVtx]->Fill(fTrackInfo->dcosYVtx, weight);
    dq[NCType::kDQTrackEndZ]->Fill(fTrackInfo->endZ, weight);
    dq[NCType::kDQTrackEndY]->Fill(fTrackInfo->endY, weight);
    dq[NCType::kDQEndMetersToEdge]->Fill(fTrackInfo->endMetersToCloseEdge, weight);
  }
  if(fRecoInfo->numberTracks > 0 && fEventInfo->showers > 0.)
    dq[NCType::kDQTrackExtension]->Fill(trackExtension, weight);

  dq[NCType::kDQTotalStrips]->Fill(fEventInfo->totalStrips, weight);
  dq[NCType::kDQPulseHeight]->Fill(1.e-3*fEventInfo->pulseHeight, weight);

  //make the stability plots that Karol loves so much
  dq[NCType::kDQEnergy]->Fill(energy, weight);

  //different quadrants in the detector - beam center is basically at (1.5, 0.)
  if(fRecoInfo->vtxX > 1.5 && fRecoInfo->vtxY > 0.)
    dq[NCType::kDQEnergyQ1]->Fill(energy, weight);
  else if(fRecoInfo->vtxX < 1.5 && fRecoInfo->vtxY > 0.)
    dq[NCType::kDQEnergyQ2]->Fill(energy, weight);
  else if(fRecoInfo->vtxX < 1.5 && fRecoInfo->vtxY < 0.)
    dq[NCType::kDQEnergyQ3]->Fill(energy, weight);
  else if(fRecoInfo->vtxX > 1.5 && fRecoInfo->vtxY < 0.)
    dq[NCType::kDQEnergyQ4]->Fill(energy, weight);

  //different ranges in z - fiducial volume is 1.73 < z < 4.73 m
  if(fRecoInfo->vtxZ > 1.73 && fRecoInfo->vtxZ < 2.73)
    dq[NCType::kDQEnergyZ1]->Fill(energy, weight);
  else if(fRecoInfo->vtxZ > 2.73 && fRecoInfo->vtxZ < 3.73)
    dq[NCType::kDQEnergyZ2]->Fill(energy, weight);
  else if(fRecoInfo->vtxZ > 3.73 && fRecoInfo->vtxZ < 4.73)
    dq[NCType::kDQEnergyZ3]->Fill(energy, weight);

  //different annuli - just do 2, the first with radius < 0.5, the second beyond that
  if(TMath::Sqrt(TMath::Power(fRecoInfo->vtxX-1.4885, 2.)
                 + TMath::Power(fRecoInfo->vtxY-0.1397, 2.)) < 0.5)
    dq[NCType::kDQEnergyA1]->Fill(energy, weight);
  else if(TMath::Sqrt(TMath::Power(fRecoInfo->vtxX-1.4885, 2.)
                      + TMath::Power(fRecoInfo->vtxY-0.1397, 2.)) > 0.5)
    dq[NCType::kDQEnergyA2]->Fill(energy, weight);

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::FindCorrelations(int dataMC)
{

  //count up the number of extractions that categorize an
  //event as nc
  int ncLike = 0;
  int ccLike = 0;
  vector<bool> isNCLike;
  vector<bool> isCCLike;
  for(uint i = 0; i < fAnalysisInfos.size(); ++i){
    isNCLike.push_back((bool)fAnalysisInfos[i]->isNC);
    isCCLike.push_back((bool)fAnalysisInfos[i]->isCC);
    ncLike += fAnalysisInfos[i]->isNC;
    ccLike += fAnalysisInfos[i]->isCC;
  }

  fNCLikeHist[dataMC]->Fill(ncLike*1. + 0.5);
  fCCLikeHist[dataMC]->Fill(ccLike*1. + 0.5);

  //find the correlations between each method
  //count up total nclike for each method and 
  //total that are common between each method and the others
  int numX = fNCLikeDenominator[dataMC]->GetXaxis()->GetNbins();
  int numY = fNCLikeDenominator[dataMC]->GetYaxis()->GetNbins();

  //check that numX = numY = fAnalysisInfos.size()
  if(numX != numY || numX != (int)isNCLike.size()){
    MSG("NCAnalysisModule", Msg::kWarning) << "cant do correlations"
                                           << " " << numX << " " 
                                           << numY << " " << isNCLike.size()
                                           << " - bail" << endl;
    return;
  }

  for(int x = 0; x < numX; ++x){
    for(int y = 0; y < numY; ++y){
      
      if(isNCLike[x]){
        fNCLikeDenominator[dataMC]->Fill(x*1.+0.5, y*1.+0.5);
        if(isNCLike[y]) fNCLikeNumerator[dataMC]->Fill(x*1.+0.5, y*1.+0.5);
        
        if(dataMC == NCType::kMC 
           && fTruthInfo->interactionType == NCType::kNC){
          fNCLikeDenominatorTrue->Fill(x*1.+0.5, y*1.+0.5);
          if(isNCLike[y]) fNCLikeNumeratorTrue->Fill(x*1.+0.5, y*1.+0.5);
        }//end if truth nc
        if(dataMC == NCType::kMC 
           && fTruthInfo->interactionType == NCType::kCC){
          fCCLikeDenominatorTrue->Fill(x*1.+0.5, y*1.+0.5);
          if(isNCLike[y]) fCCLikeNumeratorTrue->Fill(x*1.+0.5, y*1.+0.5);
        }//end if truth nc
      }//end if nclike for method x

      if(isCCLike[x]){
        fCCLikeDenominator[dataMC]->Fill(x*1.+0.5, y*1.+0.5);
        if(isCCLike[y]) fCCLikeNumerator[dataMC]->Fill(x*1.+0.5, y*1.+0.5);
        
      }//end if cclike for method x
      
    }//end loop over y
  }//end loop over x

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::CountPOTs(bool fakeData)
{
  map<int, double>::iterator dataItr = fDataPOTs.begin();
  map<int, double>::iterator mcItr = fMCPOTs.begin();
  int run = 100*fHeaderInfo->run + fHeaderInfo->subRun;
  if(fHeaderInfo->dataType == (int)SimFlag::kData){
    dataItr = fDataPOTs.find(run);
    if( dataItr != fDataPOTs.end() ) 
      fDataPOTs[run] += fBeamInfo->protonIntensity;
    else 
      fDataPOTs[run] = fBeamInfo->protonIntensity;

    double month = (1.*fHeaderInfo->month-4.5)+12.*(fHeaderInfo->year-2005);
    fDataQualityTotal[NCType::kDQPOTVsTime]->Fill(month, fBeamInfo->protonIntensity);
    fDataQualityCC[NCType::kDQPOTVsTime]->Fill(month, fBeamInfo->protonIntensity);
    fDataQualityNC[NCType::kDQPOTVsTime]->Fill(month, fBeamInfo->protonIntensity);
  }
  else if(fakeData){
    
    dataItr = fDataPOTs.find(run);
    
    //ND files have POT per snarl, FD files have POT per Run
    if( dataItr != fDataPOTs.end() 
        && fHeaderInfo->detector == (int)Detector::kNear){
      MAXMSG("NCAnalysisModule", Msg::kDebug, 100) 
        << "data " << fMCPOTPerSnarl << " " 
        << fHeaderInfo->run << " " << fHeaderInfo->subRun 
        << " " << run << " " << fDataPOTs[run] << endl;
      fDataPOTs[run] += fMCPOTPerSnarl;
      MAXMSG("NCAnalysisModule", Msg::kDebug, 100) 
        << fDataPOTs[run] << endl;
    }
    else 
      fDataPOTs[run] = fMCPOTPerSnarl;
  }
  else if(fHeaderInfo->dataType == (int)SimFlag::kMC){ 

    mcItr = fMCPOTs.find(run);
    if( mcItr != fMCPOTs.end() 
        && fHeaderInfo->detector == (int)Detector::kNear){
      MAXMSG("NCAnalysisModule", Msg::kDebug,100) 
        << "mc " << fMCPOTPerSnarl << " " 
        << fHeaderInfo->run << " " 
        << fMCPOTs[run] << endl;
      fMCPOTs[run] += fMCPOTPerSnarl;
      MAXMSG("NCAnalysisModule", Msg::kDebug,100) 
        << fMCPOTs[run] << endl;
    }
    else 
      fMCPOTs[run] = fMCPOTPerSnarl;
    MAXMSG("NCAnalysisModule", Msg::kDebug,100) 
      << fMCPOTs[run] << endl;
  }

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::WritePOTInfo(int dataMC)
{
  if(fDataPOTs.size() != 0 && dataMC == NCType::kData){
    int firstDataRun = fDataPOTs.begin()->first;
    //last run is given by the first entry in the reversed map
    int lastDataRun = fDataPOTs.rbegin()->first;
    TH1F *dataPOT = new TH1F("dataPOT", "", lastDataRun-firstDataRun+1, 
                             1.*firstDataRun, 1.*(lastDataRun+1));
    map<int, double>::iterator dataItr = fDataPOTs.begin();
    while( dataItr != fDataPOTs.end() ){
      dataPOT->Fill(1.*dataItr->first + 0.5, dataItr->second);
      ++dataItr;
    }
    dataPOT->Write();
    MSG("NCAnalysisModule", Msg::kInfo) << "total data POT = " 
                                        << dataPOT->Integral()
                                        << " x 10^{12}" << endl;

    fPOTVsTime->Write();
    delete dataPOT;
  }

  if(fMCPOTs.size() != 0 && dataMC == NCType::kMC){
    int firstMCRun = fMCPOTs.begin()->first;
    //last run is given by the first entry in the reversed map
    int lastMCRun = fMCPOTs.rbegin()->first;
    TH1F *mcPOT = new TH1F("mcPOT", "", lastMCRun-firstMCRun+1, 
                           1.*firstMCRun, 1.*(lastMCRun+1));
    map<int, double>::iterator mcItr = fMCPOTs.begin();
    while( mcItr != fMCPOTs.end() ){
      mcPOT->Fill(1.*mcItr->first + 0.5, mcItr->second);
      ++mcItr;
    }
    mcPOT->Write();
    MSG("NCAnalysisModule", Msg::kInfo) << "total MC POT = " 
                                        << mcPOT->Integral()
                                        << " x 10^{12}" << endl;
    delete mcPOT;
  }
                                     
  return;
}

//-----------------------------------------------------------------------
double NCAnalysisModule::GetTotalPOTCount(int dataMC)
{
  double total = 0.;

  if(fDataPOTs.size() != 0 && dataMC == NCType::kData){
    map<int, double>::iterator dataItr = fDataPOTs.begin();
    while( dataItr != fDataPOTs.end() ){
      total += dataItr->second;
      ++dataItr;
    }
  }

  if(fMCPOTs.size() != 0 && dataMC == NCType::kMC){
    map<int, double>::iterator mcItr = fMCPOTs.begin();
    while( mcItr != fMCPOTs.end() ){
      total += mcItr->second;
      ++mcItr;
    }
  }

  return total;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::WriteFile(TFile *file, TTree *events, TTree *beam, int dataMC)
{

  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* saveDir = gDirectory;
  
  file->cd();

  //normalization already set above
  for(int k = 0; k < static_cast<int>(fExtractions.size()); ++k){
    fExtractions[k]->WriteResources();
  }
  
  events->Write();
  beam->Write();
  WritePOTInfo(dataMC);
  for(int i = 0; i < NCType::kNumEfficiencyAndPurityBaseNames; ++i){
    for(int j = 0; j < 2; ++j){
      WriteComparisonCanvas(i, j);
    }
  }

  //set the names of the histogram bins and set up the correlation histogram
  fNCLikeCorrelation[dataMC]->Divide(fNCLikeNumerator[dataMC],
                                     fNCLikeDenominator[dataMC]);

  fNCLikeHist[dataMC]->Write();
  fNCLikeNumerator[dataMC]->Write();
  fNCLikeDenominator[dataMC]->Write();
  fNCLikeCorrelation[dataMC]->Write();

  fCCLikeCorrelation[dataMC]->Divide(fCCLikeNumerator[dataMC],
                                     fCCLikeDenominator[dataMC]);

  fCCLikeHist[dataMC]->Write();
  fCCLikeNumerator[dataMC]->Write();
  fCCLikeDenominator[dataMC]->Write();
  fCCLikeCorrelation[dataMC]->Write();

  if(dataMC == NCType::kMC){
    fNCLikeCorrelationTrue->Divide(fNCLikeNumeratorTrue,
                                   fNCLikeDenominatorTrue);
    fNCLikeNumeratorTrue->Write();
    fNCLikeDenominatorTrue->Write();
    fNCLikeCorrelationTrue->Write();

    fCCLikeCorrelationTrue->Divide(fCCLikeNumeratorTrue,
                                   fCCLikeDenominatorTrue);
    fCCLikeNumeratorTrue->Write();
    fCCLikeDenominatorTrue->Write();
    fCCLikeCorrelationTrue->Write();
  }

  TDirectory *dir = file->mkdir("DataQuality", "DataQuality");
  dir->cd();
  
  if(dataMC == NCType::kData){
    for(int i = 0; i < NCType::kDQNumDists; ++i){
      fDataQualityTotal[i]->Write();
      fDataQualityNC[i]->Write();
      fDataQualityCC[i]->Write();
    }

    for(int i = 0; i < NCType::kDQ2DNumDists; ++i){
      fDataQuality2DTotal[i]->Write();
      fDataQuality2DNC[i]->Write();
      fDataQuality2DCC[i]->Write();
    }

    TGraph *vtxXY = new TGraph(fVtxX.size(), &fVtxX[0], &fVtxY[0]);
    vtxXY->SetTitle("");
    vtxXY->SetName("vertexXY");
    vtxXY->GetHistogram()->SetXTitle("x_{vtx} (m)");
    vtxXY->GetHistogram()->SetYTitle("y_{vtx} (m)");
    vtxXY->GetHistogram()->GetXaxis()->CenterTitle();
    vtxXY->GetHistogram()->GetYaxis()->CenterTitle();
    vtxXY->Write();

    TDirectory *dirli = file->mkdir("LowIntensity", "LowIntensity");
    dirli->cd();

    for(int i = 0; i < NCType::kDQNumDists; ++i){
      fDataQualityTotalLowIntensity[i]->Write();
      fDataQualityNCLowIntensity[i]->Write();
      fDataQualityCCLowIntensity[i]->Write();
    }

  }
  else if(dataMC == NCType::kMC){
    //scale all the MC histograms to 2.5e20 POT
    double mcPOT = 1.e12*GetTotalPOTCount(NCType::kMC);
    MSG("NCAnalysisModule", Msg::kInfo) << "mcPOT = " << mcPOT << endl;

    for(int i = 0; i < NCType::kDQNumDists; ++i){
      fDataQualityMCTotal[i]->Scale(2.5e20/mcPOT);
      fDataQualityMCTotal[i]->Write();
      fDataQualityMCTotal[i]->Scale(mcPOT/2.5e20);

      fDataQualityMCNC[i]->Scale(2.5e20/mcPOT);
      fDataQualityMCNC[i]->Write();
      fDataQualityMCNC[i]->Scale(mcPOT/2.5e20);

      fDataQualityMCCC[i]->Scale(2.5e20/mcPOT);
      fDataQualityMCCC[i]->Write();
      fDataQualityMCCC[i]->Scale(mcPOT/2.5e20);

      fDataQualityMCTotalOsc[i]->Scale(2.5e20/mcPOT);
      fDataQualityMCTotalOsc[i]->Write();
      fDataQualityMCTotalOsc[i]->Scale(mcPOT/2.5e20);

      fDataQualityMCNCOsc[i]->Scale(2.5e20/mcPOT);
      fDataQualityMCNCOsc[i]->Write();
      fDataQualityMCNCOsc[i]->Scale(mcPOT/2.5e20);

      fDataQualityMCCCOsc[i]->Scale(2.5e20/mcPOT);
      fDataQualityMCCCOsc[i]->Write();
      fDataQualityMCCCOsc[i]->Scale(mcPOT/2.5e20);
    }
  }
  
  file->cd();

  file->Write("",TObject::kOverwrite);

  MSG("NCAnalysisModule", Msg::kInfo) << "file " << file->GetName() 
                                      << " written" << endl;
  //file->Close();

  saveDir->cd();

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::WriteComparisonCanvas(int hist, int nccc)
{
  if(hist > NCType::kNumEfficiencyAndPurityBaseNames){
    MSG("NCAnalysisModule", Msg::kWarning) << "trying to compare non-existant"
                                           << " histogram in WriteComparisonCanvas"
                                           << endl;
    return;
  }
  if(fExtractions.size() == 0){
    MSG("NCAnalysisModule", Msg::kWarning) << "no extractions to compare"
                                           << " in WriteComparisonCanvas"
                                           << endl;
    return;
  }

  TString name = fExtractions[0]->GetIDHistogram(hist, nccc)->GetName();
  int pos = name.Index(fExtractions[0]->GetDataType());
  name.Remove(pos);
  if(nccc == NCType::kNC) name += "NC";
  else if(nccc == NCType::kCC) name += "CC";
  name += "Comparisons";

  TCanvas *canv = new TCanvas(name, name, 150, 10, 600, 600);

  TLegend *leg = new TLegend(0.75, 0.75, 0.9, 0.9);
  leg->SetBorderSize(0);

  //make the pads
  int numX = 2;
  int numY = 2;
  double xSize = 0.5;
  double ySize = 0.5;

  TString xNames[2] = {"1", "2"};
  TString yNames[2] = {"1", "2"};
  TString padNames[4] = {"Efficiency", "Purity", "FoM1", "FoM2"};
  TString drawOption = "ace1";

  //make a vector of TPads
  vector<TPad *> pads;
  int ctr = 0;
  
  //loop over and make the pads, position them from left to right, top to bottom
  for(int y = 0; y < numY; ++y){
    for(int x = 0; x < numX; ++x){

      pads.push_back(new TPad(padNames[ctr]+xNames[x]+yNames[y], "", 
                               0.0+x*xSize, 1.0-(y+1)*ySize, 
                               (x+1)*xSize, 1.0-y*ySize));

      canv->cd();
      pads[ctr]->Draw();
      pads[ctr]->SetGridx();
      pads[ctr]->SetGridy();
      pads[ctr]->cd();

      for(uint i = 0; i < fExtractions.size(); ++i){
        if(i == 0) drawOption = "ace1";
        else drawOption = "ce1";

        if(ctr == 0){
          fExtractions[i]->GetEfficiency(hist, nccc)->Draw(drawOption);
          leg->AddEntry(fExtractions[i]->GetEfficiency(hist, nccc), 
                        fExtractions[i]->GetDataType(), "l");
        }
        else if(ctr == 1)
          fExtractions[i]->GetPurity(hist, nccc)->Draw(drawOption);
        else if(ctr == 2)
          fExtractions[i]->GetFoM1(hist, nccc)->Draw(drawOption);
        else if(ctr == 3)
          fExtractions[i]->GetFoM2(hist, nccc)->Draw(drawOption);
      }//end loop over extractions

      ++ctr;
    }//end x loop
  }//end y loop

  leg->Draw();  
  canv->Update();
  canv->Write();

  delete leg;
  delete canv;

  return;
}

//-----------------------------------------------------------------------
void NCAnalysisModule::DoADMFit()
{
  NCExtractionADM* extADM
    = dynamic_cast<NCExtractionADM *>(fExtractions[1]); 
  
  if (extADM!=0x0) extADM->FitForNCCC();
 
  double nc = 0.;
  double cc = 0.;
  double ncError = 0.;
  double ccError = 0.;
  double ncTrue = 0.;
  double ccTrue = 0.;
  double binCenter = 0.;  
  vector<double> nctrueVec;
  vector<double> cctrueVec;

  for(int i = 0; i < NCType::kNumRangeLowerLimitsADM; ++i){
//     fMCADM->GetEnergyBinFitParameters(i, nc, cc, ncError, ccError, 
//                                    ncTrue, ccTrue, binCenter);
    MSG("NCAnalysisModule", Msg::kInfo) << binCenter << " nc - " << nc 
                                        << " +/- " << ncError
                                        << "/" << ncTrue 
                                        << " cc - " << cc 
                                        << " +/- " << ccError
                                        << "/" << ccTrue << endl;

    cctrueVec.push_back(ccTrue);
    nctrueVec.push_back(ncTrue);
  }//end loop over energy bins

  vector<double> ncVec;
  vector<double> ncErrorVec;
  vector<double> ccVec;
  vector<double> ccErrorVec;

  if (extADM!=0x0) extADM->GetEnergyBinFitParameters(ccVec, ccErrorVec, 
                                                     ncVec, ncErrorVec);

                                    
  //setting total pot to 0 for now - need to change the object later 
  //as the total pot are now stored in a histogram written to the root 
  //file
  fResultsADM->SetResults(ccVec, ccErrorVec, ncVec, 
                          ncErrorVec, 0., 
                          static_cast<int>(ncVec.size()));


  return;
}


//----------------------------------------------------------------------
void NCAnalysisModule::FillRecoInfo(ANtpRecoInfo *recoInfo,
                                    ANtpHeaderInfo *headerInfo, 
                                    ANtpEventInfoNC *eventInfo,
                                    ANtpShowerInfoNC *showerInfo, 
                                    ANtpTrackInfoNC *trackInfo,
                                    ANtpBeamInfo *beamInfo,
                                    NCAnalysisCuts *cuts,
                                    NCAnalysisUtils *utils,
                                    int beamType)
{
  void *junk = showerInfo;
  junk = 0;

  //reset the values
  recoInfo->Reset();
    
  recoInfo->inFiducialVolume = (int)cuts->InBeamFiducialVolume();
  recoInfo->isFullyContained = (int)(cuts->InBeamFiducialVolume()
                                     && cuts->IsStoppingBeamMuon());
  
  recoInfo->passesCuts = 1;
  recoInfo->pass = 1;
 
  //GetEventEnergy returns sum of track and shower energy - NC events
  //dont count the track energy so just make the reco'd nu energy
  //for NC events the shower energy only
  recoInfo->nuEnergyCC = utils->GetEventEnergy(CandShowerHandle::kCC);
  recoInfo->nuEnergyNC = utils->GetShowerEnergy(CandShowerHandle::kNC);
  recoInfo->nuEnergy   = utils->GetEventEnergy();
  
  if(eventInfo->tracks > 0)
    recoInfo->muEnergy = utils->GetTrackEnergy(cuts->IsStoppingBeamMuon());
  else
    recoInfo->muEnergy = 0.;

  if(eventInfo->showers > 0){
    recoInfo->showerEnergy = utils->GetShowerEnergy(CandShowerHandle::kCC);
    recoInfo->showerEnergyCC = utils->GetShowerEnergy(CandShowerHandle::kCC);
    recoInfo->showerEnergyNC = utils->GetShowerEnergy(CandShowerHandle::kNC);
  }
  else{
    recoInfo->showerEnergy = 0.;
    recoInfo->showerEnergyCC = 0.;
    recoInfo->showerEnergyNC = 0.;
  }

  //check that this is a reasonable shower
  recoInfo->isGoodShower = (int)cuts->IsGoodShower();

  //these arent filled yet 
//   recoInfo->recoQENuEnergy = 1.;
//   recoInfo->recoQEQ2 = 1.;
//   recoInfo->recoNuDCos = 1.;

  if(recoInfo->nuEnergy > 0.){ 
    recoInfo->hadronicY = recoInfo->showerEnergy;
    recoInfo->hadronicY /= recoInfo->nuEnergy;
  }

  recoInfo->muDCosZVtx = trackInfo->dcosZVtx;
  double vtxX = 0., vtxY = 0., vtxZ = 0.;
  recoInfo->vtxR = utils->GetEventVertex(vtxX, vtxY, vtxZ);
  recoInfo->vtxX = vtxX;
  recoInfo->vtxY = vtxY;
  recoInfo->vtxZ = vtxZ;
  recoInfo->eventLength = eventInfo->lengthInPlanes;
  recoInfo->trackExtension = eventInfo->trackExtension;
  recoInfo->numberTracks = eventInfo->tracks;
  recoInfo->trackLength = trackInfo->length;

  bool isdata = true;
  if(headerInfo->dataType == (int)SimFlag::kMC) 
    isdata = false;

  Detector::Detector_t  detType = Detector::kNear;
  if(headerInfo->detector == (int)Detector::kFar) 
    detType = Detector::kFar;

  recoInfo->trackMomentum = CorrectSignedMomentumFromCurvature(trackInfo->fitMomentum, 
                                                               isdata, detType );
  
  //RPL 17/08/04: protect against passing -ive momenta to correction fn
  if( 0 < trackInfo->rangeMomentum ){
    recoInfo->trackRange = CorrectMomentumFromRange(trackInfo->rangeMomentum,  
                                                    isdata, detType);
  }
  else{
    recoInfo->trackRange = trackInfo->rangeMomentum;
  }
  
  recoInfo->sigmaQoverP = trackInfo->sigmaQoverP;
  recoInfo->eventInSnarl = eventInfo->index;
   
  recoInfo->weight = 1.0;
  recoInfo->weightRunII = 1.0;
  // calculate the MEGA weight for MC events in fiducial for MC set
  if(headerInfo->dataType == (int)SimFlag::kMC){

    // if the second argument is false, it uses the mock data weight
    // if it is true, it uses the weight from the best fit to the data
    recoInfo->weight = utils->FindMEGAFitWeight(beamType, NCType::kRunI);
    recoInfo->weightRunII = utils->FindMEGAFitWeight(beamType, NCType::kRunII);

    //SetWeight(recoInfo, headerInfo, beamType, utils);

  }//end if in fiducial and MC

  // data cleaning variables
  recoInfo->minTimeSeparation = eventInfo->minTimeSeparation;
  recoInfo->closeTimeEvent = eventInfo->closeTimeEvent;

  recoInfo->totalStrips = eventInfo->totalStrips;
  recoInfo->planes = eventInfo->planes;
  recoInfo->trackPlanes = trackInfo->planes;
  recoInfo->showerPlanes = showerInfo->planes;
  recoInfo->showerTotalStrips = showerInfo->totalStrips;
  recoInfo->evtEnergyGeV = eventInfo->energyGeV;

  if(detType == Detector::kNear){
    recoInfo->closeTimeDeltaZ = eventInfo->closeTimeDeltaZ;
    recoInfo->edgeActivityPH = eventInfo->edgeActivityPH;
    recoInfo->edgeActivityStrips = eventInfo->edgeActivityStrips;
    recoInfo->oppEdgePH = eventInfo->oppEdgePH;
    recoInfo->oppEdgeStrips = eventInfo->oppEdgeStrips;
    recoInfo->isCleanLowMultSnarl = (int)cuts->IsCleanLowMultSnarl();
    recoInfo->isCleanHighMultSnarl = (int)cuts->IsCleanHighMultSnarl();
  }
  
  recoInfo->deltaTimeSpill = (eventInfo->stripTime1st - 1.e-9*beamInfo->nearestNSToSpill)*1.e6;//in us
  
  return;
}

//----------------------------------------------------------------------
void NCAnalysisModule::SetWeight(ANtpRecoInfo *recoInfo,
                                 ANtpHeaderInfo *headerInfo,
                                 int beamType,
                                 NCAnalysisUtils *utils)
{
  double xfptWeight = recoInfo->weight;

  //these numbers are from using TRandom2::Gaus() with the appropriate means
  //and sigmas.  the means for the Delta m^2 and sin^2 2 theta are the 
  //best fit values of the CC analysis, the f_s mean is 0.
  double showerEnergyShift = 0.1098;
  double trackEnergyShift = -0.0176;
  double normalizationFar = 0.0014;
  double skzpShift = 1.3822;  
  double deltaMSqr = 0.002867;
  double sinSqr2Theta = 0.8472;
  double fSterile = 0.1452;

  //find the weight for the skzp weighting.  the uncertainty method returns
  //an relative uncertainty, so to apply it have to multiply by the actual 
  //weight.  for example if the uncertainty is 2% then the method returns 1.02
  xfptWeight *= utils->FindMEGAFitWeightUncertainty(beamType, 
                                                    NCType::kRunI,
                                                    skzpShift);

  //get the survival probability for this neutrino
  double survivalProb = 1.;
  std::vector<double> oscPars;
  oscPars.push_back(sinSqr2Theta);
  oscPars.push_back(deltaMSqr);
  oscPars.push_back(fSterile);

  survivalProb = utils->FindSurvivalProbability(oscPars);

  //near detector only has the xfpt weight applied
  recoInfo->weight = xfptWeight;

  //far detector gets the relative normalization and survival probability
  if(headerInfo->detector == (int)Detector::kFar)
    recoInfo->weight *= normalizationFar*survivalProb;

  //apply the energy scale shifts - do uniform shower energy shifts 
  //no matter what type of event you are looking at.
  recoInfo->muEnergy *= 1. + trackEnergyShift;
  
  recoInfo->showerEnergy *= 1. + showerEnergyShift;
  recoInfo->showerEnergyCC *= 1. + showerEnergyShift;
  recoInfo->showerEnergyNC *= 1. + showerEnergyShift;

  recoInfo->nuEnergyCC = recoInfo->muEnergy + recoInfo->showerEnergyCC;
  recoInfo->nuEnergyNC = recoInfo->showerEnergyNC;
  recoInfo->nuEnergy   = recoInfo->muEnergy + recoInfo->showerEnergy;
  
  if(headerInfo->detector == (int)Detector::kFar)
    headerInfo->run += 899000;
  else if(headerInfo->detector == (int)Detector::kNear)
    headerInfo->run += 500000;

  return;
}

---