NuCut.cxx

Go to the documentation of this file.

#include <cassert>
#include <iomanip>

#include "MessageService/MsgService.h"
#include "Conventions/SimFlag.h"
#include "Conventions/Detector.h"
#include "Conventions/Munits.h"
#include "DataUtil/infid.h"

#include "NtupleUtils/NuEvent.h"
#include "NtupleUtils/NuMCEvent.h"
#include "NtupleUtils/NuCut.h"

#include "TString.h"

CVSID("$Id: NuCut.cxx,v 1.14 2010/02/01 14:26:46 bckhouse Exp $");


NuCut::NuCut(TString name, const NuPlots *plots) : 
  fFidVol(), 
  fSelName(name), 
  fPlots(plots),
  fPassCuts(true)
{
  fEventCount = 0;
  fAnaVersion = AnaVersionHash(name);
}

NuCut::~NuCut()
{
  
}


// Does a full cut
Bool_t NuCut::MakeCuts(const NuEvent &nu)
{
  // Since we are asking for a full cut, reset the pass status
  ResetStatus();
  
  MakePreselectionCuts(nu);
  MakeSelectionCuts(nu);
  
  return fPassCuts;
}


// Apply just the preselection
Bool_t NuCut::MakePreselectionCuts(const NuEvent &nu)
{
  Preselection(nu);
  MAXMSG("NuCut",Msg::kDebug,50) << "Passed " << GetName() << " Preselection = " << (fPassCuts ? "Yes" : "No") << endl;
  return fPassCuts;
}


// apply just the main selection
Bool_t NuCut::MakeSelectionCuts(const NuEvent &nu)
{
  Selection(nu);
  return fPassCuts;
}


// Sets the cutter back to an 'uncut' state
void NuCut::ResetStatus(void)
{
  fPassCuts = true;
  
  // Push back the cut statii onto the central list. This shouldn't go here, probably
  // but has for now.
  Bool_t passed = true;
  for (UInt_t i = 0; i < fCutList.size(); ++i) {
    // Take each of these, and push them back onto the cut totals list.
    TString &cutname = fCutList[i].first;
    Bool_t passcut = fCutList[i].second;
    
    // Check if it already exists in the vector
    UInt_t sumindex;
    for (sumindex = 0; sumindex < fCutTotals.size(); ++sumindex)
    {
      if (fCutTotals[sumindex].first == cutname) break;
    }
    
    // Now check if we have to push_back a new entry
    if (sumindex == fCutTotals.size())
    {
      // We have to push back
      fCutTotals.push_back(pair<TString,Int_t>(cutname, 0));
    }
    
    // Now, shall we add to this?
    // Add to totals, if we have still passed after this cut
    passed = passed && passcut;
    if (passed)
    {
      fCutTotals[sumindex].second++;
    }
    
    
    // std::pair<TString,Int_t>
  }
  
  // If the cut list was not empty, add to the event count
  if (!fCutList.empty()) fEventCount++;
  
  // Now, clear the interim list
  fCutList.clear();
}


Bool_t NuCut::Cut_If(Bool_t condition, const TString &s)
{
  fPassCuts = fPassCuts && (!condition);

  // Store this cut, whatever for now
  fCutList.push_back(pair<TString,Int_t>(s,!condition));
  
  return fPassCuts;
}


Bool_t NuCut::Keep_If(Bool_t condition, const TString &s)
{
  fPassCuts = fPassCuts && condition;

  // Store this cut, whatever for now
  fCutList.push_back(pair<TString,Int_t>(s,condition));
  
  return fPassCuts;
}


// Only cuts if data
Bool_t NuCut::Keep_Data_If(Bool_t condition, const NuEvent &nu, const TString &name)
{
  // Only cuts if data
  if (nu.simFlag == SimFlag::kData)
  {
    return Keep_If(condition, name);
  }
  
  return fPassCuts;
}


// Only cuts if data
Bool_t NuCut::Cut_Data_If(Bool_t condition, const NuEvent &nu, const TString &name)
{
  // Only cuts if data
  if (nu.simFlag == SimFlag::kData)
  {
    return Cut_If(condition, name);
  }
  
  return fPassCuts;
}


// Set the fiducial volume
void NuCut::SetFidVol(const TString &fidvol)
{
  if (fFidVol != "" && fidvol == "") {
    MSG("NuCut", Msg::kWarning) << "Setting blank fiducial volume in " << GetName() << ", but one has already been set."
      << "Cannot reset - will continue using whatever is set (but not forcing it)" << endl;
  }
  fFidVol = fidvol;
  // Set the fiducial volume
  if (FidVol::getName() != fidvol) {
    choose_infid_set(fidvol.Data());
  }
}


Bool_t NuCut::InFidVol(const NuEvent &nu) const
{
  // Check we are still in the correct fiducial volume (and we have set one)
  if (FidVol::getName() != fFidVol) {
    MSG("NuCut",Msg::kInfo) << "Fiducial volume has changed in " << GetName() << ", resetting to \"" << fFidVol << "\"" << endl;
    // Reset the fiducial volume to outs
    if (FidVol::getName() != fFidVol) {
      choose_infid_set(fFidVol.Data());
    }
  }
  
  // assert(nu.ntrk > 0);  // No track means the track vertex is meaningless
  if (nu.ntrk == 0) {
    MAXMSG("NuCut",Msg::kWarning,3) << "Testing vertex against Fiducial volume, but no tracks!" << endl;
  }
  
  return infid(static_cast<Detector::Detector_t>(nu.detector),
    static_cast<SimFlag::SimFlag_t>(nu.simFlag),
    nu.xTrkVtx,nu.yTrkVtx,nu.zTrkVtx);
}


Bool_t NuCut::InFidVolEvt(const NuEvent &nu) const
{
  // Check we are still in the correct fiducial volume (and we have set one)
  if (FidVol::getName() != fFidVol) {
    MSG("NuCut",Msg::kInfo) << "Fiducial volume has changed in " << GetName() << ", resetting to \"" << fFidVol << "\"" << endl;
    // Reset the fiducial volume to outs
    if (FidVol::getName() != fFidVol) {
      choose_infid_set(fFidVol.Data());
    }
  }

  return infid(Detector::Detector_t(nu.detector),
               SimFlag::SimFlag_t(nu.simFlag),
               nu.xEvtVtx, nu.yEvtVtx, nu.zEvtVtx);
}
  

Bool_t NuCut::InFidVolTrueEvt(const NuMCEvent& mc)  const
{
  // Check we are still in the correct fiducial volume (and we have set one)
  if (FidVol::getName() != fFidVol) {
    MSG("NuCut",Msg::kInfo) << "Fiducial volume has changed in " << GetName() << ", resetting to \"" << fFidVol << "\"" << endl;
    // Reset the fiducial volume to outs
    if (FidVol::getName() != fFidVol) {
      choose_infid_set(fFidVol.Data());
    }
  }
  
  return infid(Detector::Detector_t(mc.detector),
               SimFlag::SimFlag_t(mc.simFlag),
               mc.vtxxMC,mc.vtxyMC,mc.vtxzMC);
}


void NuCut::Defer_Preselection(NuCut &cut, const NuEvent & nu)
{
  // Make sure the object is in a fresh state
  cut.ResetStatus();
  
  // Now, run the preselection cut
  cut.MakePreselectionCuts(nu);
  
  // // Now, combine out status - this is currently very simple, as we are
  // // not yet tracking named cuts. This will become a bit more involved 
  // // at that stage, but for now, this works!
  // //fPassCuts = fPassCuts && !cut.Failed();
  // 
  // 
  // // Not yet tracking individual cuts, but in the interm track
  // // what the whole deferred selection does
  // Keep_If(cut.Passed(), "DeferredPreselection");

  // Now copy out all the selections into our list!
  for (UInt_t i = 0; i < cut.fCutList.size(); ++i) {
    Keep_If(cut.fCutList[i].second, cut.fCutList[i].first);
  }
}


void NuCut::Defer_Selection(NuCut &cut, const NuEvent & nu)
{
  // Make sure the object is in a fresh state
  cut.ResetStatus();
  
  // Now, run the preselection cut
  cut.MakeSelectionCuts(nu);
  
  // Now copy out all the selections
  for (UInt_t i = 0; i < cut.fCutList.size(); ++i) {
    Keep_If(cut.fCutList[i].second, cut.fCutList[i].first);
  }
  
  // Now, combine out status - this is currently very simple, as we are
  // not yet tracking named cuts. This will become a bit more involved 
  // at that stage, but for now, this works!
  //fPassCuts = fPassCuts && !cut.Failed();
}


Bool_t NuCut::Passed(TString cutname) const
{
  for (UInt_t i = 0; i < i < fCutList.size(); ++i) {
    if (fCutList[i].first == cutname) {
      return fCutList[i].second;
    }
  }
  MSG("NuCut",Msg::kWarning) << "Cut " << cutname << " not found in selection " << GetName() << ".  Assuming failed." << endl;
  PrintCurrent();
  return false;
}


Bool_t NuCut::PassedExcept(TString cutname) const
{
  Bool_t found = false;
  Bool_t pass = true;
  for (UInt_t i = 0; i < fCutList.size(); ++i) {
    if (fCutList[i].first != cutname) {
      pass = pass && fCutList[i].second;
    } else {
      found = true;
    }
  }
  if (!found) {
    MSG("NuCut",Msg::kWarning) << "Cut " << cutname << " not found in selection " << GetName() << ".  Just returning selection, not N-1." << endl;
    PrintCurrent();
  }

  return pass;
}


void NuCut::PrintCurrent(void) const
{
  MSG("NuCut",Msg::kInfo) << "Selection " << GetName() << " has " << fCutList.size() << " cuts:" << endl;
  for (UInt_t i = 0; i < fCutList.size(); ++i) {
    MSG("NuCut",Msg::kInfo) << "  " << fCutList[i].first << " " << (fCutList[i].second?"Passed":"Failed") << endl;
  }
  if (Passed()) MSG("NuCut",Msg::kInfo) << "  -- Event Passed" << endl;
  else          MSG("NuCut",Msg::kInfo) << "  ** Event Failed" << endl;

}



void NuCut::PrintSummary(void)
{
  // Make sure we are reset and counted
  ResetStatus();
  
  // Find the longest name
  Int_t lsize = 12;
  for (UInt_t i = 0; i < fCutTotals.size(); ++i)
  {
    Int_t len = fCutTotals[i].first.Length();
    if (len > lsize) lsize = len;
  }
  // The 'total width'
  Int_t width = lsize + 20;
  // Print a summary of this events cuts
  MSG("NuCut",Msg::kInfo) << setw(width/2) << right << setfill('=') << (" "+GetName()) << setw(width/2+width%2) << left << " Summary " << '\n';
  MSG("NuCut",Msg::kInfo) << setw(lsize) << left << setfill(' ') << "Cut" << "  # Passed  % Passed\n";

  MSG("NuCut",Msg::kInfo) << setw(width) << left << setfill('=') << "" << "\n";
  MSG("NuCut",Msg::kInfo) << setw(lsize) << left << setfill(' ') << "Total Events" << "  " << setw(8) << right << fEventCount << "   100.00%\n";
  for (UInt_t i = 0; i < fCutTotals.size(); ++i)
  {
    // stat = stat && fCutList[i].second;
    MSG("NuCut",Msg::kInfo)
      << setw(lsize) << left << fCutTotals[i].first
      << "  " << setw(8) << right << fCutTotals[i].second
      << "  " << setw(7) << right << fixed << setprecision(2) << fCutTotals[i].second*100./fEventCount << "%\n";
  }
  MSG("NuCut",Msg::kInfo) << endl;
}


Int_t NuCut::AnaVersionHash(const TString name)
{
  TString upper(name);
  name.Hash();
  upper.ToUpper();
  Int_t hash = upper.Hash();
  // This will reduce the hash space by half, but will ensure that we never
  // get any conflict from real (old) anaVersion numbers
  if (hash > 0) hash *= -1;
  
  return hash;
}

---