NueModule.cxx

Go to the documentation of this file.

// $Id: NueModule.cxx,v 1.30 2008/03/03 17:52:25 boehm Exp $
//
// FILL_IN: [Document your code!!]
//
// vahle@hep.ucl.ac.uk
#include <dirent.h>
#include "TFile.h"
#include "TChain.h"
#include "TClonesArray.h"
#include "TProfile2D.h"
#include "StandardNtuple/NtpStRecord.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro
#include "Calibrator/CalMIPCalibration.h"
#include "MCNtuple/NtpMCRecord.h"
#include "MCNtuple/NtpMCTruth.h"
#include "TruthHelperNtuple/NtpTHRecord.h"
#include "CandNtupleSR/NtpSRRecord.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "NueAna/NueAnaTools/SntpHelpers.h"
#include "NueAna/Module/NueModule.h"
#include "NueAna/NueRecord.h"
#include "NueAna/NueRecordAna.h"
#include "BeamData/ana/Summary/BeamSummary.h"
#include "MCReweight/MuParentHelper.h"
#include "MCReweight/Zbeam.h"
#include "MCReweight/Zfluk.h"
#include "MCReweight/Kfluk.h"
#include "MCReweight/NeugenWeightCalculator.h"
#include "MCReweight/MCReweight.h"
#include "MuonRemoval/NtpMRRecord.h"
#include "CalDetDST/UberRecord.h"
#include "MCReweight/SKZPWeightCalculator.h"

#include "NueAna/NueStandard.h"
#include "NueAna/NueAnaTools/Selection.h"
#include "DataUtil/GetTempTags.h"   
#include "TSystem.h"


CVSID("$Id: NueModule.cxx,v 1.30 2008/03/03 17:52:25 boehm Exp $");

#include "DatabaseInterface/DbiResultPtr.tpl"


JOBMODULE(NueModule, "NueModule",
          "");
//......................................................................

NueModule::NueModule():
   tmpltfile(""),
   kDPlaneCut(-1),
   kLoPhNStripCut(-1),
   kLoPhNPlaneCut(-1),
   kPhStripCut(-1),
   kPhPlaneCut(-1),
   kBeamType(BeamType::kUnknown),
   kMuPiDir(""),
   kOpenedMuPiFile(false),
   counter(0),
   passcounter(0),
   passcutcounter(0),
   failcounter(0),
   writecounter(0),
   foundmeu(false),
   MSTminsig(0.),
   MSTmaxmetric(0.),
   MSTminfarsig(0.),
   MSTmaxmetriclowz(0.),
   SIGCORRMEU(1.),
   MSTetemplate(0),
   MSTbtemplate(0),
   MSTemtemplate(0),
   MSTbmtemplate(0),
   threshCut(0.),
   sasFileName(""),
   pot(0.),
   MEUPERGEV(25),   //not correct but close, correctly set by NueCalibration
                    //  adjusted on 01-08-2008 by JAB
   foundRelease(false),
   release(0),
   skzpcfg("PiMinus_CedarDaikon")
{
  mupar = new MuParentHelper();
//  zbeam = new Zbeam();
//  zfluk = new Zfluk();
//  zfluk->UseParameterization("SKZP"); //this is the default
  kfluk = new Kfluk();

  skzpCalc = new SKZPWeightCalculator(skzpcfg, true);
  mcr = &MCReweight::Instance();
  NeugenWeightCalculator *n=new NeugenWeightCalculator();
  mcr->AddWeightCalculator(n);

  xsecreweightreg = new Registry();
  xsecreweightreg->UnLockValues();
  xsecreweightreg->UnLockKeys();
  xsecreweightreg->Clear();
  xsecreweightreg->Set("neugen:config_name","MODBYRS");
  xsecreweightreg->Set("neugen:config_no",3);
  xsecreweightreg->LockValues();
  xsecreweightreg->LockKeys();

  mrccRun = false;
  
  kPidName = "None";
  kPIDHighCut = ANtpDefVal::kDouble;
  kPIDLowCut = ANtpDefVal::kDouble;
  kCCPIDCut = ANtpDefVal::kDouble;
  kHighECut = ANtpDefVal::kDouble;
  kLowECut = ANtpDefVal::kDouble;
}

//......................................................................

NueModule::~NueModule()
{
  if(mupar){
    delete mupar;
    mupar=0;
  }
  if(zbeam){
    delete zbeam;
    zbeam=0;
  }

  if(skzpCalc){
    delete skzpCalc;
    skzpCalc=0;
  }
  if(zfluk){
    delete zfluk;
    zfluk=0;
  }
  if(kfluk){
    delete kfluk;
    kfluk=0;
  }
}

//......................................................................

JobCResult NueModule::Reco(MomNavigator* mom)
{
//======================================================================
// Reads in sue-style ntuples from mom, calculates a bunch of useful
// nue quantites
//======================================================================
   MSG("NueModule",Msg::kDebug)<<"In NueModule::Reco"<<endl;

   if(counter%1000==0){      
      MSG("NueModule",Msg::kInfo)<<"On entry "<<counter<<endl;
   }
   counter++;
   bool foundSR=false;
   bool foundMC=false;
   bool foundTH=false;
   bool foundMR=false;
   bool foundUR=false;

   bool foundST=false;
   bool passesCuts = false;

   NtpStRecord *str = 0;
   NtpStRecord *oldst = 0;
   NtpMRRecord *mr = 0;

   //if we didn't find a NtpStRecord, it could be because we're reading an 
   //old tree, check for the other records
   NtpMCRecord *mc=0;
   NtpSRRecord *sr=0;
   NtpTHRecord *th=0;

   VldContext vc;  
   str = dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord",
                                                      "Primary"));
   if(str){
     foundST=true;
     vc=str->GetHeader().GetVldContext();
     if(!foundRelease){
        release = str->GetRelease();
        foundRelease = true; 
        title = ReleaseType::AsString(release);

        string file = DataUtil::GetTempTagString(str, "file");
        filename = gSystem->BaseName(file.c_str());
     }
     //check for MR:
     mr=dynamic_cast<NtpMRRecord *>(mom->GetFragment("NtpMRRecord"));
     if(mr){
       if(ReleaseType::IsBirch(release)) {
         oldst=dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord",
                                                            "NtpStRecordOld")); 
       }
       else if(ReleaseType::IsCedar(release)){
         oldst=str;
         str=dynamic_cast<NtpStRecord *>(mom->GetFragment("NtpStRecord",
                                                          "MuonRemoved"));
       }
       if(oldst) foundMR = true;
     }
   }
   else {
     release = ReleaseType::kUnknown;
     title = ReleaseType::AsString(release);
     sr = dynamic_cast<NtpSRRecord *>(mom->GetFragment("NtpSRRecord"));
     if (sr) {
       foundSR = true;
       vc=sr->GetHeader().GetVldContext();        
     }
     mc = static_cast<NtpMCRecord *>(mom->GetFragment("NtpMCRecord"));
     if (mc) {
       foundMC = true;
       th = static_cast<NtpTHRecord *>(mom->GetFragment("NtpTHRecord"));
       if (th) foundTH = true;
     }
   }

   if(!foundSR&&!foundMC&&!foundST&&!foundMR){
      //got nothing useful from mom
      MSG("NueModule",Msg::kWarning)<<"Got Nothing from mom"<<endl;
      failcounter++;
      return JobCResult::kFailed;
   }
   if(mrccRun&&!foundMR&&foundST){
      return JobCResult::kFailed;
   }

   UberRecord *ur = dynamic_cast<UberRecord *>(mom->GetFragment("UberRecord"));
   if(ur) foundUR=true;
   
   //UGLY!!!!
   //need to be able to compare single strip 
   //pulseheight between ND and FD, but we can't
   //in the ntuple.  Will ask the db for the
   //conversion from sigcor to MEU 
   //on the first event, and use throughout.  
   //Since sigcorrs are already corrected for all strip variations, 
   //it's ok to ask for 1 number.  
   //Should be ok for MC, maybe ok for 1 run of data at a time, but will
   //definitely be wrong if many weeks worth of data are strung together
   //we should revisit this and make it better in the future!!!!!!!
   if(!foundmeu){
     DbiResultPtr<CalMIPCalibration> dbp(vc);
     if(dbp.GetNumRows()>0){
       const CalMIPCalibration *m = dbp.GetRow(0);
       float mip=m->GetMIP(1.); 
       if(mip>0){
         SIGCORRMEU=1./mip;
         foundmeu=true;
       }
     }
   }

  int evtn = 0;
  const RecCandHeader *header = 0;

   if(foundST){
     evtn=str->evthdr.nevent;
     header = &(str->GetHeader());
   }

   if(!foundSR&&foundMC){
      vc = mc->GetHeader().GetVldContext();
      header = &(mc->GetHeader());
      evtn = 0;
   }

   if(foundSR){
     MSG("NueModule",Msg::kDebug)<<"Got SR from mom"<<endl;
     evtn=sr->evthdr.nevent;
     header = &(sr->GetHeader());
   }

   //if first event and we want to fix the Mu parents,
   // set the mupar file directory
   if(header->GetVldContext().GetSimFlag() == SimFlag::kData)
     kFixMuParents = 0;

   if(ReleaseType::IsData(release))
      kBeamType = DetermineBeamType(vc);
   else
      kBeamType = DetermineBeamType(filename);

   if(!kOpenedMuPiFile&&kFixMuParents && ReleaseType::IsCarrot(release)){
     string flxsfx="le010z185i";
     if(kBeamType==BeamType::kL010z185i){     flxsfx="le010z185i";  }
     else if(kBeamType==BeamType::kL100z200i){flxsfx="le100z200i";  }
     else if(kBeamType==BeamType::kL250z200i){flxsfx="le250z200i";  }
     else if(kBeamType==BeamType::kL150z200i){flxsfx="le150z200i";  }
     else if(kBeamType==BeamType::kL010z170i){flxsfx="le010z170i";  }
     else if(kBeamType==BeamType::kL010z200i){flxsfx="le010z200i";  }
     else if(kBeamType==BeamType::kL010z000i){flxsfx="le010z000i";  }
        
     string fullPath = kMuPiDir + flxsfx;
     mupar->SetFileDir(kMuPiDir,true);
     kOpenedMuPiFile=true;
   }

   MSG("NueModule",Msg::kDebug)<<"Events in this snarl: "<<evtn<<endl;

   if(ReleaseType::IsData(release))
      kBeamType = DetermineBeamType(vc);
               
   //int kzarkoBeamType = BeamType::ToZarko(kBeamType);
                                             
   if(evtn == 0)
   {
      NueHeader h(vc);
      h.SetSnarl(header->GetSnarl());
      h.SetRun(header->GetRun());
      h.SetSubRun(header->GetSubRun());
      h.SetEventNo(-1);
      h.SetEvents(0);
      h.SetTrackLength(0);
      h.SetRelease(release);
      h.SetBeamType(kBeamType);
      h.SetFoundBits(foundSR, foundMC || foundST, foundTH,foundMR);

     //make an ana_nue object
     NueRecord *nue = new NueRecord(h);
     //make a nuerecordana object
     NueRecordAna ana_nue(*nue);
     ana_nue.SetRelease(release);
     ana_nue.SetBeamType(kBeamType);

     nue->SetTitle(title.c_str());
     //fill the truth info
     //if it's data, we can also do beam mon analyze
     //  ana_nue.bmona.SetBeamSummary(bsum);
                                                                                
     //must set the beam and parent helper to fill flux info
//     ana_nue.anaia.SetBeam(kzarkoBeamType);
//     ana_nue.nuefwa.SetBeam(kzarkoBeamType);
//     ana_nue.anaia.SetBeam(kBeamType);
//     ana_nue.nuefwa.SetBeam(kBeamType);
     ana_nue.nuefwa.SetSKZPCalc(skzpCalc, skzpcfg);
//     ana_nue.nuefwa.SetZfluk(zfluk);
     ana_nue.nuefwa.SetKfluk(kfluk);
     ana_nue.fiana.SetMuParentHelper(mupar);
     ana_nue.fiana.FixMuParents(kFixMuParents);
               
     //set xsec registry and mcreweight to do xsec reweighting
     ana_nue.nuexsa.SetMCReweight(mcr);
     ana_nue.nuexsa.SetRegistry(xsecreweightreg);
     
     if(foundST)
       ana_nue.FillTrue(0,str);
     if(foundMC || foundSR)
       ana_nue.FillTrue(0,sr,mc,th);

     ana_nue.aneia.Analyze(-10,str);  //need to fill srevent to get quality
                                                                                
     //give the ana_nue object to mom, she'll write it to the file
     mom->AdoptFragment(nue);
     writecounter++;
     failcounter++;
     return JobCResult::kPassed;
   }

   bool anypassed = false;
  for(int i=0;i<evtn;i++){
     passesCuts = false;
    
     const int SIZE = 5500;
     float ph0[SIZE] = {-1};
     float ph1[SIZE] = {-1};
     for(int k=0; k < SIZE; k++) { ph0[k] = ph1[k] = -1;}

     SntpHelpers::FillEventEnergy(ph0, ph1, i, str, SIZE);
 
     //make a header
     NueHeader h(vc);
     
     //fill header info
     h.SetSnarl(header->GetSnarl());
     h.SetRun(header->GetRun());
     h.SetSubRun(header->GetSubRun());
     h.SetEventNo(i);
     h.SetEvents(evtn);
     h.SetRelease(release);
     h.SetBeamType(kBeamType);
     h.SetFoundBits((foundSR || foundST), foundMC || foundST,
                    foundTH || foundST, foundMR);

     MSG("NueModule",Msg::kDebug)<<"Getting event"<<endl;
     
     NtpSREvent *event = 0;
     if(foundST) event = SntpHelpers::GetEvent(i,str);
     if(foundSR) event = SntpHelpers::GetEvent(i,sr);
     
     //loop over tracks in this event, find longest
     int longtrack=0;
     if(event){
       for(int j=0;j<event->ntrack;j++){
           int tindex = SntpHelpers::GetTrackIndex(j,event);
           NtpSRTrack *track = 0;
           if(foundST) track = SntpHelpers::GetTrack(tindex,str);
           if(foundSR) track = SntpHelpers::GetTrack(tindex,sr);
           if(track && longtrack<track->plane.n){
              longtrack = track->plane.n;
           }
        }
     }
     
     h.SetTrackLength(longtrack);
     
     NueRecord *nue = new NueRecord(h);
     //make a nuerecordana object
     NueRecordAna ana_nue(*nue);
     ana_nue.SetRelease(release);
     ana_nue.SetBeamType(kBeamType);
     ana_nue.SetEventEnergyArray(ph0, ph1);

     nue->SetTitle(title.c_str());

     
     if(foundST)
     {                   
        fCut.SetInfoObject(i, str);
        passesCuts = fCut.PassesAllCuts();
     }

     if(foundSR)
     {
       passesCuts = true;
       MSG("NueModule",Msg::kWarning)
          <<"Unable to cut on pre-Birch files - why are you running these files?"<<endl;
     }
     
     ana_nue.aneia.SetParams(SIGCORRMEU, MEUPERGEV);
     ana_nue.ansia.SetParams(SIGCORRMEU, MEUPERGEV);
     ana_nue.antia.SetParams(SIGCORRMEU, MEUPERGEV);
     ana_nue.hca.SetParams(SIGCORRMEU, MEUPERGEV);

     //must set the beam and parent helper to fill flux info
//     ana_nue.nuefwa.SetBeam(kzarkoBeamType);
//     ana_nue.nuefwa.SetBeam(kBeamType);
     ana_nue.nuefwa.SetSKZPCalc(skzpCalc, skzpcfg);

//     ana_nue.nuefwa.SetZbeam(zbeam);
//     ana_nue.nuefwa.SetZfluk(zfluk);
     ana_nue.nuefwa.SetKfluk(kfluk);
     ana_nue.fiana.SetMuParentHelper(mupar);
     ana_nue.fiana.FixMuParents(kFixMuParents);

     //set xsec registry and mcreweight to do xsec reweighting
     ana_nue.nuexsa.SetMCReweight(mcr);
     ana_nue.nuexsa.SetRegistry(xsecreweightreg);
  

     if(foundST){
         ana_nue.FillTrue(i,str);
         ana_nue.FillReco(i,str);
     }
     if(foundSR){
         ana_nue.FillTrue(i,sr,mc,th);
         ana_nue.FillReco(i,sr);
     }
                                                                              
     if(passesCuts)
     {
         MSG("NueModule",Msg::kDebug)<<"Tracklength: "<<longtrack
            <<"Event Length: "<<event->plane.n
            <<"Event Energy: "<<event->ph.sigcor
            <<"Event Energy : "<<event->ph.sigcor<<endl;

         passcutcounter++;
         //analyze calculates everybodies variables, then calls the Truth and Reco Object Fillers
         //configure MSTCalcAna routine
         MSG("NueModule",Msg::kDebug)<<"Setting templates: "
                                     <<MSTetemplate->GetName()<<" "
                                     <<MSTbtemplate->GetName()<<" "
                                     <<MSTemtemplate->GetName()<<" "
                                     <<MSTbmtemplate->GetName()<<endl;
         
         ana_nue.msta.SetSigTemplate(MSTetemplate);
         ana_nue.msta.SetBGTemplate(MSTbtemplate);
         ana_nue.msta.SetSigMIPTemplate(MSTemtemplate);
         ana_nue.msta.SetBGMIPTemplate(MSTemtemplate);
         ana_nue.msta.SetMSTParams(MSTminsig,MSTmaxmetric,
                                   MSTminfarsig,MSTmaxmetriclowz,SIGCORRMEU);
         ana_nue.sfa.SetParams(SIGCORRMEU, MEUPERGEV);
         ana_nue.sfa.SetCutParams(kDPlaneCut,kPhStripCut,kPhPlaneCut);
         ana_nue.fva.SetParams(SIGCORRMEU, MEUPERGEV);
         ana_nue.mda.SetMdaParams(threshCut, sasFileName);
                               
         if(foundST)  ana_nue.Analyze(i,str);
         if(foundSR)  ana_nue.Analyze(i,sr);
         if(foundMR)  ana_nue.Analyze(i,mr,oldst);
         if(foundUR)  ana_nue.Analyze(ur);
     }// end of If(passcuts)
     else{
       cout<<"Rejecting based on cuts"<<endl;
     }

     if(PassesBlindingCuts(nue)){  
       anypassed = true;
       //give the ana_nue object to mom, she'll write it to the file
       MSG("NueModule",Msg::kDebug)<<"Giving Fragment to mom"<<endl;
       writecounter++;
       mom->AdoptFragment(nue);
       MSG("NueModule",Msg::kDebug)<<"Mom took it"<<endl;
     }
 
     if(i+1 == evtn && !anypassed){ 
       //push through something so the POT are still counted
       h.SetEvents(evtn);
       NueRecord *nue = new NueRecord(h);
       NueRecordAna ana_nue(*nue);
       ana_nue.SetRelease(release);
       ana_nue.SetBeamType(kBeamType);
       nue->SetTitle(title.c_str());
       ana_nue.aneia.Analyze(-10,str);
       writecounter++;
       mom->AdoptFragment(nue);
     }
   }
   MSG("NueModule",Msg::kDebug)<<"Done with snarl"<<endl;
   passcounter++;
   return JobCResult::kPassed; // kNoDecision, kFailed, etc.
}

//......................................................................

const Registry& NueModule::DefaultConfig() const
{
//======================================================================
// Supply the default configuration for the module
//======================================================================
   MSG("NueModule",Msg::kDebug)<<"In NueModule::DefaultConfig"<<endl;

  static Registry r  = fCut.DefaultConfig(); // Default configuration for module

  // Set name of config
  std::string name = this->GetName();
  name += ".config.default";
  r.SetName(name.c_str());

  // Set values in configuration
  r.UnLockValues();
  r.Set("MSTTmpltFile","templates.root");
  r.Set("DPlaneCut",-1);
  r.Set("LoPhNStripCut",-1);
  r.Set("LoPhNPlaneCut",-1);
  r.Set("PhStripCut",-1);
  r.Set("PhPlaneCut",-1);

  r.Set("MSTminsig",0.5);
  r.Set("MSTmaxmetric",1000.);
  r.Set("MSTminfarsig",1.5);
  r.Set("MSTmaxmetriclowz",20.);
  
  r.Set("MdaThreshCut",0.0);
  r.Set("MdaSASFile","");
  r.Set("BeamString", "");
  r.Set("MuPiDir","");
  r.Set("BeamType",2);
  r.Set("FixMuParents",0);
  r.Set("MRCC", 0);

  r.Set("HighEnergyCut", ANtpDefVal::kDouble);
  r.Set("LowEnergyCut", ANtpDefVal::kDouble);
  r.Set("HighPIDCut", ANtpDefVal::kDouble);
  r.Set("LowPIDCut", ANtpDefVal::kDouble);
  r.Set("ANTIPID", "None");
  r.Set("CCPID", "None");

  r.LockValues();

  return r;
}

//......................................................................

void NueModule::Config(const Registry& r)
{
//======================================================================
// Configure the module given the Registry r
//======================================================================
  MSG("NueModule",Msg::kDebug)<<"In NueModule::Config"<<endl;

  fCut.Config(r);

  const char* tmps;

  if(r.Get("ANTIPID", tmps)) {kPidName = tmps;}
  if(r.Get("CCPID", tmps))      { kCCPidName = tmps;}

  if(r.Get("MSTTmpltFile",tmps)) { tmpltfile = tmps;}
  int imps;
  if(r.Get("DPlaneCut",imps)) { kDPlaneCut=imps;}
  if(r.Get("LoPhNStripCut",imps)) { kLoPhNStripCut=imps;}
  if(r.Get("LoPhNPlaneCut",imps)) { kLoPhNPlaneCut=imps;}

  if(r.Get("MdaSASFile",tmps)) { sasFileName = tmps;}

  double fmps;
  if(r.Get("HighPIDCut", fmps)) { kPIDHighCut = fmps; }
  if(r.Get("LowPIDCut", fmps))  { kPIDLowCut = fmps; }
  if(r.Get("CCPIDCut", fmps))   { kCCPIDCut = fmps;  }
  if(r.Get("HighEnergyCut", fmps)) { kHighECut = fmps; }
  if(r.Get("LowEnergyCut", fmps))  { kLowECut = fmps; }

  if(r.Get("PhStripCut",fmps)) { kPhStripCut=fmps;}
  if(r.Get("PhPlaneCut",fmps)) { kPhPlaneCut=fmps;}
  if(r.Get("MSTminsig",fmps)){ MSTminsig=fmps; }
  if(r.Get("MSTmaxmetric",fmps)){ MSTmaxmetric=fmps; }
  if(r.Get("MSTminfarsig",fmps)){ MSTminfarsig=fmps; }
  if(r.Get("MSTmaxmetriclowz",fmps)){ MSTmaxmetriclowz=fmps; }
  if(r.Get("MdaThreshCut",fmps)) { threshCut = fmps;}

  if(r.Get("MuPiDir",tmps)){kMuPiDir=(string)(tmps);}

//  if(r.Get("BeamType",imps)){kBeamType=imps;}

  if(r.Get("BeamString", tmps)){
     beamstring = tmps;  
     BeamType::BeamType_t beam = BeamType::kUnknown;
     if(beamstring.find("le010z185i")!=string::npos){ beam=BeamType::kL010z185i; }
     if(beamstring.find("le100z200i")!=string::npos){ beam=BeamType::kL100z200i; }
     if(beamstring.find("le250z200i")!=string::npos){ beam=BeamType::kL250z200i; }
     if(beamstring.find("le150z200i")!=string::npos){ beam=BeamType::kL150z200i; }
     if(beamstring.find("le010z200i")!=string::npos){ beam=BeamType::kL010z200i; }
     if(beamstring.find("le010z170i")!=string::npos){ beam=BeamType::kL010z170i; }
     if(beamstring.find("le010z000i")!=string::npos){ beam=BeamType::kL010z200i; }
     kBeamType = beam;
  }

  if(r.Get("FixMuParents",imps)){ 
    if(imps) kFixMuParents=true; 
    else kFixMuParents=false;
  }

  if(r.Get("MRCC",imps)){ SetMRCCRunning(imps);}
}

void NueModule::BeginJob()
{
  MSG("NueModule",Msg::kDebug)<<"In NueModule::BeginJob"<<endl;
   
  if(tmpltfile!=""){
    MSG("NueModule",Msg::kDebug)<<"opening MST template file "
                                <<tmpltfile<<endl;
    TFile f(tmpltfile.c_str());
    if(f.IsOpen()){
      MSTetemplate = (TProfile2D *)f.Get("nlambdanele");
      MSTbtemplate = (TProfile2D *)f.Get("nlambdanother");
      MSTemtemplate = (TProfile2D *)f.Get("mipdistele");
      MSTbmtemplate = (TProfile2D *)f.Get("mipdistother");
      MSG("NueModule",Msg::kDebug)<<"Did I get the histos? "
                                  <<MSTetemplate<<" "<<MSTbtemplate<<" "
                                  <<MSTemtemplate<<" "<<MSTbmtemplate<<endl;
      MSTetemplate->SetDirectory(0);
      MSTbtemplate->SetDirectory(0);
      MSTemtemplate->SetDirectory(0);
      MSTbmtemplate->SetDirectory(0);
      f.Close();
      MSG("NueModule",Msg::kDebug)<<"Do I still have them? "
                                  <<MSTetemplate->GetName()<<endl;

    }
  }
// Beam Summary ntuples are obsolete!
//  LoadBeamSummaryData(bmondir.c_str());
}

void NueModule::EndJob()
{

   MSG("NueModule",Msg::kInfo)<<"Counter "<<counter
                              <<" passcutcounter "<<passcutcounter
                              <<" passcounter "<<passcounter
                              <<" failcounter "<<failcounter
                              <<" writecounter "<<writecounter<<endl;
//   MSG("NueModule",Msg::kInfo)<<"Number of POT in this run: "<<pot<<endl;

   if(MSTetemplate!=0){
     delete MSTetemplate;
     MSTetemplate=0;
   }
   if(MSTbtemplate!=0){
     delete MSTbtemplate;
     MSTbtemplate=0;
   }
   if(MSTemtemplate!=0){
     delete MSTemtemplate;
     MSTemtemplate=0;
   }
   if(MSTbmtemplate!=0){
     delete MSTbmtemplate;
     MSTbmtemplate=0;
   }

}

void NueModule::LoadBeamSummaryData(const char *bd)
{
    MSG("NueModule",Msg::kError)<<"Beam Summary ntuples are obsolete!"<<endl;
    return;


  DIR *dfd;
  if(!(dfd =  opendir(bd))){
    MSG("NueModule",Msg::kError)<<"Can not ls Beam Monitoring path "
                                <<bd<<" "<<dfd<<std::endl;
    return;
  }
}

BeamType::BeamType_t NueModule::DetermineBeamType(VldContext vc){

/*    BDSpillAccessor& sa = BDSpillAccessor::Get();
    VldContext evt_vldc = nr->GetHeader().GetVldContext();

    const BeamMonSpill* bms = sa.LoadSpill(vc.GetTimeStamp());
    VldTimeStamp bms_vts;
    if (!bms) {
       MSG("NueBeamMon",Msg::kError) << "No BeamMonSpill found for " << evt_vldc << endl;
       bms_vts = VldTimeStamp::GetEOT();
    }
    else {
        bms_vts = bms->SpillTime();
    }
    return bms->BeamType();
*/
    int time = vc.GetTimeStamp().GetSec();

    BeamType::BeamType_t beam = BeamType::kUnknown;

    if(time >= 1107216000 && time < 1109539850) beam = BeamType::kL100z200i;
    if(time >= 1109540615 && time < 1109899325) beam = BeamType::kL250z200i;
    if(time >= 1109899938 && time < 1110239564) beam = BeamType::kL100z200i;
    if(time >= 1110323763 && time < 1111622400) beam = BeamType::kL000z200i;
    if(time >= 1114892377 && time < 1115927583) beam = BeamType::kL100z200i;
    if(time >= 1115937438 && time < 1116604821) beam = BeamType::kL250z200i;
    if(time >= 1116618256 && time < 1122659668) beam = BeamType::kL010z185i;
    if(time >= 1122659886 && time < 1122922688) beam = BeamType::kL010z170i;
    if(time >= 1122922890 && time < 1123112674) beam = BeamType::kL010z200i;
    if(time >= 1123112803 && time < 1139605423) beam = BeamType::kL010z185i;
    if(time >= 1139605543 && time < 1140022084) beam = BeamType::kL010z000i;
    if(time >= 1140026702 && time < 1140908579) beam = BeamType::kL010z185i;
    // End of Run 1
    if(time >= 1149180600 && time < 1150047780) beam = BeamType::kL150z200i;
    if(time >= 1150047780 && time < 1151690460) beam = BeamType::kL250z200i;
    if(time >= 1153956600 && time < 1155510000) beam = BeamType::kL250z200i;
    if(time >= 1158004800 && time < 1158019870) beam = BeamType::kL010z200i;
    if(time >= 1158019870 && time < 1161892800) beam = BeamType::kL010z185i;
    if(time >= 1161892800 && time < 1184351737) beam = BeamType::kL010z185i;
    if(time >= 1184351737 && time < 1184708040) beam = BeamType::kL010z000i;
    //End of Run 2

    if(time >= 1184800000 ) beam = BeamType::kL010z185i;

    return beam;
}

BeamType::BeamType_t NueModule::DetermineBeamType(string file)
{
   BeamType::BeamType_t beam = BeamType::kUnknown;
   if(file.find("L010185")!=string::npos){ beam=BeamType::kL010z185i; }
   if(file.find("L100200")!=string::npos){ beam=BeamType::kL100z200i; }
   if(file.find("L250200")!=string::npos){ beam=BeamType::kL250z200i; }
   if(file.find("L150200")!=string::npos){ beam=BeamType::kL150z200i; }
   if(file.find("L010200")!=string::npos){ beam=BeamType::kL010z200i; }
   if(file.find("L010170")!=string::npos){ beam=BeamType::kL010z170i; }
   if(file.find("L010000")!=string::npos){ beam=BeamType::kL010z000i; }

   return beam;
}

bool NueModule::PassesBlindingCuts(NueRecord *nr)
{
  bool passes = true;

  if( kPidName != "None"){
     Selection::Selection_t pid = Selection::StringToEnum(kPidName.c_str());

     double pidVal = NueStandard::GetPIDValue(nr, pid);
     
     if(!ANtpDefVal::IsDefault(kPIDHighCut))
         passes = pidVal < kPIDHighCut;

     if(!ANtpDefVal::IsDefault(kPIDLowCut))
         passes = (passes) && pidVal > kPIDLowCut;

    if(pid != Selection::kANN6) nr->ann.pid_6inp = ANtpDefVal::kDouble;
    if(pid != Selection::kANN30) nr->ann.pid_30inp = ANtpDefVal::kDouble;
    if(pid != Selection::kSSPID) 
             nr->subshowervars.pid_daikon = ANtpDefVal::kDouble;
    if(pid != Selection::kMCNN) nr->mcnnv.fracCCy = ANtpDefVal::kDouble;
    if(pid != Selection::kCuts) nr->treepid.fCutPID = ANtpDefVal::kInt;

  }

  if(!ANtpDefVal::IsDefault(kCCPIDCut)){
     if(kCCPidName == "CC_DP")
         passes = passes &&  nr->mri.orig_cc_pid > kCCPIDCut;
     if(kCCPidName == "CC_AB")
         passes = passes &&  nr->mri.orig_abCCPID > kCCPIDCut;
  }

  NueConvention::NueEnergyCorrection(nr);

  if(!ANtpDefVal::IsDefault(kHighECut))
     passes = passes && nr->srevent.phNueGeV < kHighECut;

  if(!ANtpDefVal::IsDefault(kLowECut))
     passes =  passes && nr->srevent.phNueGeV > kLowECut;

  return passes;
}

---