SelectAntiNeutrino.cxx

Go to the documentation of this file.

// $Id: SelectAntiNeutrino.cxx,v 1.26 2009/02/28 21:46:16 gmieg Exp $

// C++
#include <cmath>
#include <cassert>
#include <climits>
#include <iostream>
#include <sstream>

// ROOT
#include "TDirectory.h"
#include "TMath.h"

// MINOS
#include "Registry/Registry.h"

// Local
#include "PhysicsNtuple/Default.h"
#include "PhysicsNtuple/Factory.h"

// Local
#include "SelectAntiNeutrino.h"
#include "SelectCount.h"

REGISTER_ANP_OBJECT(AlgEvent,SelectAntiNeutrino)
REGISTER_ANP_OBJECT(AlgSnarl,SelectAntiNeutrino)

using namespace std;

//---------------------------------------------------------------------------------------
Anp::SelectAntiNeutrino::SelectAntiNeutrino()
   :fCount(new SelectCount()),
    fAddData(false),
    fErase(false),
    fAddKey(5530),
    fAngKey(108),
    fPidKey(4280),
    fSelKey(5524),
    fSelect(1),
    fPidCut(0.5)
{
}

//------------------------------------------------------------------------------------------
Anp::SelectAntiNeutrino::~SelectAntiNeutrino()
{
}

//------------------------------------------------------------------------------------------
bool Anp::SelectAntiNeutrino::Run(Record &record)
{
   //
   // Run selection for individual events
   //

   EventIterator ievent = record.EventBegIterator();   
   while(ievent != record.EventEndIterator())
   {
      if(!SelectAntiNeutrino::Run(*ievent, record, true) && fErase)
      {
         ievent = record.Erase(ievent);
      }
      else
      {
         ++ievent;
      }
   }

   return true;
}

//------------------------------------------------------------------------------------------
bool Anp::SelectAntiNeutrino::Run(Event &event, const Record &record, const bool pass)
{
   //
   // Run \bar{\nu}_{\mu} selection algorithm
   //

   bool select = false;
   switch(fSelect)
   {
   case 1: 
      select = SelectAntiNeutrino::Run1(event, record, pass);
      break;
   case 2:
      select = SelectAntiNeutrino::Run2(event, record, pass);
      break;
   default:
      cerr << "SelectAntiNeutrino::Run - unknown selection case" << endl;
      break;
   }

   //
   // Event passed this selection - add corresponding key
   //
   if(select)
   {
      if(!event.KeyExists(fSelKey))
      {
         event.Add(fSelKey, +1.0);
      }
      else
      {
         cerr << "SelectAntiNeutrino::Run - key already exists: " << fSelKey << endl;
      }
   }

   return select;
}

//------------------------------------------------------------------------------------------
bool Anp::SelectAntiNeutrino::Run1(Event &event, const Record &record, const bool pass)
{
   //
   // Run standard anti-neutrino selection
   //
   if(pass)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }  

   //
   // select longest track, with largest number of scintillator planes
   //
   const TrackIter itrack = Anp::LongestTrack(event, record);   
   if(itrack != record.TrackEnd())
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }
   
   //
   // fill event data for variables used in a cut
   //
   if(fAddData) SelectAntiNeutrino::AddData(event, *itrack);

   //
   // check that track passed fitter algorithm
   //
   if(itrack -> PassFit())
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   //
   // select mu- tracks and tracks with small error on the fit momentum
   //
   if(itrack -> ErrorQP() > 0.0 && itrack -> QP()/itrack -> ErrorQP() > +1.0)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }
   
   //
   // select tracks that are at least 40 scintillator planes long
   //
   if(itrack -> GetBasic().NPlane() > 39)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   //
   // cut on directional charge id
   //
   if(itrack -> KeyExists(fAngKey) && (itrack -> DataAt(fAngKey) < 1.00 || 
                                      itrack -> DataAt(fAngKey) > 5.00))
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }
        
   //
   // cut on fit probability
   //
   const double chi2_prob = TMath::Prob(itrack -> Chi2(), int(itrack -> Ndof()));
   if(chi2_prob > 0.01)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   return true;
}

//------------------------------------------------------------------------------------------
bool Anp::SelectAntiNeutrino::Run2(Event &event, const Record &record, const bool pass)
{
   if(pass)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }  

   //
   // select longest track, with largest number of scintillator planes
   //
   const TrackIter itrack = Anp::LongestTrack(event, record);
   if(itrack != record.TrackEnd())
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   // fill event data for variables used in a cut
   if(fAddData) SelectAntiNeutrino::AddData(event, *itrack);
   
   // check that track passed fitter algorithm
   if(itrack -> PassFit())
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   // select mu- tracks and tracks with small error on the fit momentum
   if(itrack -> ErrorQP() > 0.0 && itrack -> QP()/itrack -> ErrorQP() > +1.0)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }
   
   // select tracks that pass kNN muon cut
   if(event.KeyExists(fPidKey) && event[fPidKey] > fPidCut)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   // cut on directional charge id
   if(itrack -> KeyExists(fAngKey) && (itrack -> DataAt(fAngKey) < 1.00 ||
                                      itrack -> DataAt(fAngKey) > 5.00))
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }       

   // cut on fit probability
   const double chi2_prob = TMath::Prob(itrack -> Chi2(), int(itrack -> Ndof()));
   if(chi2_prob > 0.01)
   {
      fCount -> Run(event, record, true);
   }
   else
   {
      fCount -> Run(event, record, false);
      return false;
   }

   return true;
}

//------------------------------------------------------------------------------------------
void Anp::SelectAntiNeutrino::Config(const Registry &reg)
{
   //
   // Configure self
   //

   Anp::Read(reg, "SelectAntiNeutrinoAddData", fAddData);
   Anp::Read(reg, "SelectAntiNeutrinoErase",   fErase);

   reg.Get("SelectAntiNeutrino",       fSelect);
   reg.Get("SelectAntiNeutrinoAddKey", fAddKey);
   reg.Get("SelectAntiNeutrinoAngKey", fAngKey);
   reg.Get("SelectAntiNeutrinoPidKey", fPidKey);
   reg.Get("SelectAntiNeutrinoSelKey", fSelKey);
   reg.Get("SelectAntiNeutrinoPidCut", fPidCut);

   assert(std::abs(fPidKey) < SHRT_MAX && "key is out of range");

   if(reg.KeyExists("PrintConfig"))
   {
      cout << "SelectAntiNeutrino::Config" << endl
           << "   AddData = " << fAddData << endl
           << "   Erase = " << fErase << endl
           << "   AddKey = " << fAddKey << endl
           << "   AngKey = " << fAngKey << endl
           << "   PidKey = " << fPidKey << endl
           << "   SelKey = " << fSelKey << endl
           << "   Select = " << fSelect << endl
           << "   PidCut = " << fPidCut << endl;
   }

   fCount -> Config(reg);
}

//------------------------------------------------------------------------------------------
void Anp::SelectAntiNeutrino::AddData(Event &event, const Track &track)
{   
   //
   // check that track passed fitter algorithm
   //
   if(track.PassFit())
   {
      event.Add(fAddKey + 0, 1.0);
   }

   //
   // select mu- tracks and tracks with small error on the fit momentum
   //
   if(track.ErrorQP() > 0.0)
   {
      event.Add(fAddKey + 1, track.QP()/track.ErrorQP());
   }
     
   //
   // select tracks that are at least 40 steel planes long
   //
   event.Add(fAddKey + 2, track.GetBasic().NPlane());

   //
   // cut on directional charge id
   //
   if(track.KeyExists(fAngKey))
   {
      event.Add(fAddKey + 3, track.DataAt(fAngKey));

      if(!event.KeyExists(fAngKey))
      {
         event.Add(fAngKey, track.DataAt(fAngKey));
      }
   }
     
   //
   // cut on fit probability
   //
   event.Add(fAddKey + 4, TMath::Prob(track.Chi2(), int(track.Ndof())));
}

//------------------------------------------------------------------------------------------
void Anp::SelectAntiNeutrino::End(const DataBlock &block)
{
   fCount -> End(block);
}

---