DemoMCFilter.cxx

Go to the documentation of this file.

/************************************************************************
  $Id: DemoMCFilter.cxx,v 1.6 2002/12/15 16:25:12 rhatcher Exp $

  A job module to select events based on their MC information

  messier@huhepl.harvard.edu
***********************************************************************/
#include "Demo/DemoMCFilter.h"
#include "MessageService/MsgService.h"
#include "JobControl/JobCommand.h"
#include "JobControl/JobCModuleRegistry.h" // JOBMODULE registration macro

// Required form unpacking reroot events
// #include "TClonesArray.h"
// #include "TObject.h"
#include "MINF_Classes/MINFast.h"
#include "REROOT_Classes/REROOT_Event.h"
#include "REROOT_Classes/REROOT_NeuKin.h"

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

CVSID("$Id: DemoMCFilter.cxx,v 1.6 2002/12/15 16:25:12 rhatcher Exp $");
JOBMODULE(DemoMCFilter,"MCFilter","Select MC events based on truth info.");

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

DemoMCFilter::DemoMCFilter() : fNuIntSelect(0), fModeSelect(0) { }

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

JobCResult DemoMCFilter::Ana(const MomNavigator* /* mom */) 
{
//======================================================================
// Analyze the event to see if the MC information meets the selection
// criteria.
//======================================================================
  // Get the pointer to the neutrino interaction data
  const REROOT_Event*  rrevt = gMINFast->GetREROOTEvent();
  const REROOT_NeuKin* nukin = 
    dynamic_cast<REROOT_NeuKin*>(rrevt->neukins()->First());
  
  // Energy cut
  double Enu = this->GetEnu(nukin);
  if (Enu < fEnuMin) return JobCResult::kFailed;
  if (Enu > fEnuMax) return JobCResult::kFailed;

  // y cut
  double y = this->GetY(nukin);
  if (y < fYmin) return JobCResult::kFailed;
  if (y > fYmax) return JobCResult::kFailed;
  
  // Interaction cut
  NuInt_t nuint = this->GetNuInt(nukin);
  if ((nuint & fNuIntSelect) == 0) return JobCResult::kFailed;

  // Interaction mode cut
  Mode_t numode = this->GetMode(nukin);
  if ((numode & fModeSelect) == 0) return JobCResult::kFailed;

  // Get here if all checks passed
  if (fDoPrint) this->Print(nukin);
  return JobCResult::kPassed;
}

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

const Registry& DemoMCFilter::DefaultConfig() const 
{
//======================================================================
// Supply a default configuration for this module
//======================================================================
  static Registry r;
  std::string name = this->GetName();
  name += ".config.default";
  
  r.UnLockValues();
  r.Set("EnuMin",      0.0);
  r.Set("EnuMax", 999999.0);
  r.Set("Ymin",        0.0);
  r.Set("Ymax",        1.0);

  r.Set("CCnue",     false);
  r.Set("CCnuebar",  false);
  r.Set("CCnumu",    false);
  r.Set("CCnumubar", false);
  r.Set("CCnutau",   false);
  r.Set("CCnutaubar",false);
  r.Set("NC",        false);
  r.Set("AllNuInt",  false);

  r.Set("Elastic",  false);
  r.Set("Resonant", false);
  r.Set("DIS",      false);
  r.Set("AllMode",  false);

  r.Set("DoPrint",   false);
  r.LockValues();

  return r;
}

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

void DemoMCFilter::Config(const Registry& r) 
{
//======================================================================
// Configure the module given the registry r
//======================================================================
  int    tmpb;
  double tmpd;

  if (r.Get("EnuMin", tmpd)) { fEnuMin = tmpd; }
  if (r.Get("EnuMax", tmpd)) { fEnuMax = tmpd; }
  if (r.Get("Ymin",   tmpd)) { fYmin   = tmpd; }
  if (r.Get("Ymax",   tmpd)) { fYmax   = tmpd; }
  
  if (r.Get("CCnue", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnue : fNuIntSelect &= ~kCCnue;
  }
  if (r.Get("CCnuebar", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnuebar : fNuIntSelect &= ~kCCnuebar;
  }
  if (r.Get("CCnumu", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnumu : fNuIntSelect &= ~kCCnumu;
  }
  if (r.Get("CCnumubar", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnumubar : fNuIntSelect &= ~kCCnumubar;
  }
  if (r.Get("CCnutau", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnutau : fNuIntSelect &= ~kCCnutau;
  }
  if (r.Get("CCnutaubar", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnutaubar : fNuIntSelect &= ~kCCnutaubar;
  }
  if (r.Get("NC", tmpb)) {
    tmpb ? fNuIntSelect |= kCCnumubar : fNuIntSelect &= ~kCCnumubar;
  }
  if (r.Get("AllNuInt", tmpb)) { 
    if (tmpb) fNuIntSelect = kAllNuInt; 
  }
  // if absolutely nothing is selected default to "all"
  if (!fNuIntSelect) fNuIntSelect = kAllNuInt;

  if (r.Get("Elastic", tmpb)) {
    tmpb ? fModeSelect |= kElastic : fModeSelect &= ~kElastic;
  }
  if (r.Get("Resonant", tmpb)) {
    tmpb ? fModeSelect |= kResonant : fModeSelect &= ~kResonant;
  }
  if (r.Get("DIS", tmpb)) {
    tmpb ? fModeSelect |= kDIS : fModeSelect &= ~kDIS;
  }
  if (r.Get("AllMode", tmpb)) {
    if (tmpb) fModeSelect = kAllMode;
  }
  // if absolutely nothing is selected default to "all"
  if (!fModeSelect) fModeSelect = kAllMode;

  if (r.Get("DoPrint", tmpb)) { fDoPrint = tmpb; }
}

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

void DemoMCFilter::HandleCommand(JobCommand* cmd) 
{
  const char* c = cmd->PopCmd();
  if (c) {
    string cs(c);
    if (cs == "Report") { 
      this->Report(); 
      return; 
    }
    if (cs == "Accept") {
      NuInt_t nuIntAccept = 0;
      Mode_t  modeAccept  = 0;
      const char* opt = 0;
      while ( (opt=cmd->PopOpt()) !=0 ) {
        string sopt(opt);

        // Which interactions to accept?
        if      (sopt == "CCnue")       nuIntAccept |=  kCCnue;
        else if (sopt == "CCnuebar")    nuIntAccept |=  kCCnuebar;
        else if (sopt == "CCnumu")      nuIntAccept |=  kCCnumu;
        else if (sopt == "CCnumubar")   nuIntAccept |=  kCCnumubar;
        else if (sopt == "CCnutau")     nuIntAccept |=  kCCnutau;
        else if (sopt == "CCnutaubar")  nuIntAccept |=  kCCnutaubar;
        else if (sopt == "NC")          nuIntAccept |=  kNC;
        else if (sopt == "AllNuInt")    nuIntAccept |=  kAllNuInt;

        if      (sopt == "Elastic")     modeAccept  |=  kElastic;
        else if (sopt == "Resonant")    modeAccept  |=  kResonant;
        else if (sopt == "DIS")         modeAccept  |=  kDIS;
        else if (sopt == "AllMode")     modeAccept  |=  kAllMode;

        // Which interaction modes to accept?   
        if (sopt == "Emin") {
          float e = cmd->PopFloatOpt();
          this->SetEnuMin(e);
        }

        // Which range of neutrino energy to accept?
        if (sopt == "Emax") {
          float e = cmd->PopFloatOpt();
          this->SetEnuMax(e);
        }

        // Which range of y to accept?
        if (sopt == "Ymin") {
          float y = cmd->PopFloatOpt();
          this->SetYmin(y);
        }
        if (sopt == "Ymax") {
          float y = cmd->PopFloatOpt();
          this->SetYmax(y);
        }
      }
      // Apply interaction selections
      if (nuIntAccept !=0 ) this->SetNuInt(nuIntAccept);
      if (modeAccept  !=0 ) this->SetNuMode(modeAccept);
      return;
    }
  }
}

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

void DemoMCFilter::Help() {}

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

void DemoMCFilter::Report() 
{
  MsgStream& msg = MSGSTREAM("Demo",Msg::kInfo);
  
  msg << 
    this->GetName() << ": Selecting events with:\n" << 
    "-> Enu in range [" << fEnuMin << "," << fEnuMax << "]\n" <<
    "-> y in range [" << fYmin << "," << fYmax << "]\n" <<
    "-> Interactions types: ( ";
  if (fNuIntSelect>=kAllNuInt) {
    msg << "All ";
  }
  else {
    if (fNuIntSelect&kCCnue)      msg << "CCnue ";
    if (fNuIntSelect&kCCnuebar)   msg << "CCnuebar ";
    if (fNuIntSelect&kCCnumu)     msg << "CCnumu ";
    if (fNuIntSelect&kCCnumubar)  msg << "CCnumubar ";
    if (fNuIntSelect&kCCnutau)    msg << "CCnutau ";
    if (fNuIntSelect&kCCnutaubar) msg << "CCnutaubar ";
    if (fNuIntSelect&kNC)         msg << "NC ";
  }
  msg << ")\n";
  msg << "-> Interaction modes: ( ";
  if (fModeSelect>=kAllMode) {
    msg << "All ";
  }
  else {
    if (fModeSelect&kElastic)  msg << "Elastic ";  
    if (fModeSelect&kResonant) msg << "Resonant ";
    if (fModeSelect&kDIS)      msg << "DIS ";
  }
  msg << ")\n";
}

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

void DemoMCFilter::SetNuInt(NuInt_t i)    { 
  if (i>=kAllNuInt) {
    fNuIntSelect = kAllNuInt;
  }
  else {
    fNuIntSelect = i;  
  }
}

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

void DemoMCFilter::SetNuMode(Mode_t i) 
{
  if (i>=kAllMode) {
    fModeSelect = kAllMode;
  }
  else {
    fModeSelect  = i;  
  }
}

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

void DemoMCFilter::SetEnuMin(double e)    { fEnuMin = e;        }
void DemoMCFilter::SetEnuMax(double e)    { fEnuMax = e;        }
void DemoMCFilter::SetYmin(double y)      { fYmin = y;          }
void DemoMCFilter::SetYmax(double y)      { fYmax = y;          }
void DemoMCFilter::SetPrint(bool doPrint) { fDoPrint = doPrint; }

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

DemoMCFilter::NuInt_t DemoMCFilter::GetNuInt(const REROOT_NeuKin* nukin)
{ 
  bool isNC = (nukin->IAction()==0);

  if (isNC) { return kNC; }

  int ntype = nukin->INu();
  if (ntype ==  12) { return kCCnue;      }
  if (ntype == -12) { return kCCnuebar;   }
  if (ntype ==  14) { return kCCnumu;     }
  if (ntype == -14) { return kCCnumubar;  }
  if (ntype ==  16) { return kCCnutau;    }
  if (ntype == -16) { return kCCnutaubar; }

  return 0;
}

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

DemoMCFilter::Mode_t DemoMCFilter::GetMode(const REROOT_NeuKin* /* nukin */)
{ 
  return kAllMode;
}

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

double DemoMCFilter::GetEnu(const REROOT_NeuKin* nukin) 
{ 
  return nukin->P4Neu()[3];
}

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

double DemoMCFilter::GetY(const REROOT_NeuKin* nukin)
{ 
  return nukin->Y();
}

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

void DemoMCFilter::Print(const REROOT_NeuKin* /* nukin */)
{ 
}

---