TAnalysisLoop.cxx

Go to the documentation of this file.

#include <cmath>
#include <exception>
#include <iostream>

#include <TString.h>
#include <TND280Output.hxx>
#include <TOADatabase.hxx>
#include <nd280EventLoop.hxx>

#include "eventAnalysisEnums.hxx"

#include "TAnalysisModuleBase.hxx"
#include "TBasicDataQualityModule.hxx"
#include "TBasicHeaderModule.hxx"
#include "TBeamSummaryDataModule.hxx"
#include "TFgdOnlyModule.hxx"
#include "TGRooTrackerVtxModule.hxx"
#include "TGeometrySummaryModule.hxx"
#include "TGlobalReconModule.hxx"
#include "TLowLevelInfoModule.hxx"
#include "TND280Input.hxx"
#include "TNRooTrackerVtxModule.hxx"
#include "TP0DECALReconModule.hxx"
#include "TSFGReconModule.hxx"
#include "TTagsModule.hxx"
#include "TP0DReconModule.hxx"
#include "TReconPerformanceEvalModule.hxx"
#include "TSmrdReconModule.hxx"
#include "TTRExReconModule.hxx"
#include "TTrackerECALReconModule.hxx"
#include "TTrackerReconModule.hxx"
#include "TTruthTrajectoriesModule.hxx"
#include "TTruthVerticesModule.hxx"

#include "TAnalysisLoop.hxx"

namespace AnalysisLoopGlobals {
int EventNumber = -1;
std::string CurrentlyProcessingModule = "";
bool StopProcessing = false;
;
}

namespace ND {
//******************************************************************************
//                               Utils Funcs
//******************************************************************************

/// Used to sort Preselecting modules to the end of the vector.
///
///\note this does not obey strict weak ordering and so should only be used with
/// std::stable_sort
bool PreselectionSorter(
    std::pair<std::string, ND::TAnalysisModuleBase *> const &LHS,
    std::pair<std::string, ND::TAnalysisModuleBase *> const &RHS) {
  /// Only swap if the LHS is used for preselection and the RHS isn't
  return !(LHS.second->IsUsedForPreselection() &&
           (!RHS.second->IsUsedForPreselection()));
}

//******************************************************************************
//                         Analysis Loop Funcs
//******************************************************************************

TAnalysisLoop::TAnalysisLoop()
    : fEventCount(0),
      fSaveOriginalFullEvent(false),
      fSaveOutputOaEventTree(false),
      fDisableAll(false),
      fEnableAll(false),
      fProduction(false),
      fExit(false),
      fValidation(false),
      fGeometry(false) {
  fAnalysisModules.push_back(
      std::make_pair("TBasicHeaderModule", new ND::TBasicHeaderModule));
  /*fAnalysisModules.push_back(std::make_pair("TBasicDataQualityModule",
                                            new ND::TBasicDataQualityModule));*/
  fAnalysisModules.push_back(std::make_pair("TTruthTrajectoriesModule",
                                            new ND::TTruthTrajectoriesModule));
  fAnalysisModules.push_back(
      std::make_pair("TTruthVerticesModule", new ND::TTruthVerticesModule));
  fAnalysisModules.push_back(
      std::make_pair("TSFGReconModule", new ND::TSFGReconModule));
   fAnalysisModules.push_back(
       std::make_pair("TTagsModule", new ND::TTagsModule));
  fAnalysisModules.push_back(
      std::make_pair("TP0DReconModule", new ND::TP0DReconModule));
  fAnalysisModules.push_back(
      std::make_pair("TP0DECALReconModule", new ND::TP0DECALReconModule));
  fAnalysisModules.push_back(
      std::make_pair("TTrackerECALReconModule",
                     new ND::TTrackerECALReconModule::OutputManager));
  fAnalysisModules.push_back(
      std::make_pair("TGlobalReconModule", new ND::TGlobalReconModule));
  fAnalysisModules.push_back(
      std::make_pair("TTrackerReconModule", new ND::TTrackerReconModule));
  fAnalysisModules.push_back(
      std::make_pair("TSmrdReconModule", new ND::TSmrdReconModule));
  fAnalysisModules.push_back(
      std::make_pair("TLowLevelInfoModule", new ND::TLowLevelInfoModule));
  fAnalysisModules.push_back(
      std::make_pair("TFgdOnlyModule", new ND::TFgdOnlyModule));
  fAnalysisModules.push_back(
      std::make_pair("TBeamSummaryDataModule", new ND::TBeamSummaryDataModule));
  fAnalysisModules.push_back(std::make_pair(
      "TReconPerformanceEvalModule", new ND::TReconPerformanceEvalModule));
  fAnalysisModules.push_back(
      std::make_pair("TGeometrySummaryModule", new ND::TGeometrySummaryModule));
  fAnalysisModules.push_back(
      std::make_pair("TTRExReconModule", new ND::TTRExReconModule));
  fAnalysisModules.push_back(
      std::make_pair("TGRooTrackerVtxModule", new ND::TGRooTrackerVtxModule));
  fAnalysisModules.push_back(
      std::make_pair("TNRooTrackerVtxModule", new ND::TNRooTrackerVtxModule));

  // Set defaults for enabling or disabling modules.
  SetModuleDefaults();
}

TAnalysisLoop::~TAnalysisLoop() {}

void TAnalysisLoop::Usage(void) {
  std::cout
      << "\t-O inputdir=<path/to/dir> Set an input directory where the "
         "modules\n"
         "\t   can pick up non-nd280 files. For example use this if you need "
         "to\n"
         "\t   recover pass-through information from a given directory\n"
         "\t-O production Run eventAnalysis in production mode.  It will exit if "
         "no \n"
         "\t              RooTracker information is found for MC files or if "
         "any \n"
         "\t              modules enabled modules crash\n"
         "\t-O validation Run eventAnalysis in validation mode.  It will exit if "
         "any \n"
         "\t              enabled module crashes\n"
         "\t-O geometry Save the ND280 geometry to the eventAnalysis file\n"
         "\t            (not saved by default)\n"
         "\t-O save Save the original full TND280Event in analysis file\n\n"
         "=====================Module Level Options=====================\n"
         "\t-O enable=<modulename> Enable a particular module \n"
         "\t                        (overrides all other disable/enable)\n"
         "\t-O disable=<modulename> Disable a particular module \n"
         "\t                        (overrides disable/enable all)\n"
         "\t-O enable=all Enable all modules (overrides disable all)\n"
         "\t-O disable=all Disable all modules\n"
         "\t-O preselection=<modulename> Use module to output pre-selected\n"
         "\t   events in the ND280Events oaEvent tree\n"
         "\t-O modulename=<option> Configure module\n"
         "\t(for example. -O TTruthTrajectoriesModule=saveall)\n\n"
         "\tAvailable modules {default status} and options are:"
      << std::endl;

  for (ModDescripterator mod_it = fAnalysisModules.begin();
       mod_it != fAnalysisModules.end(); ++mod_it) {
    std::cout << "\t\t" << mod_it->first
              << (mod_it->second->IsEnabledByDefault() ? " {enabled}"
                                                       : " {disabled}")
              << std::endl;
    mod_it->second->SayAvailableOpts("\t\t  ");
  }
}

bool TAnalysisLoop::SetOption(std::string option, std::string value) {
  if (option == "production") {
    ND280NamedLog("TAnalysisLoop", "Running in production mode");
    fProduction = true;
    return true;
  }

  if (option == "validation") {
    ND280NamedLog("TAnalysisLoop", "Running in validation mode");
    fValidation = true;
    return true;
  }

  if (option == "geometry") {
    ND280NamedLog("TAnalysisLoop", "Saving the ND280 geometry");
    fGeometry = true;
    return true;
  }

  if (option == "disable" && value == "all") {
    ND280NamedLog("TAnalysisLoop",
                  "Disabling all modules due to user command line option");
    fDisableAll = true;
    return true;
  }
  if (option == "enable" && value == "all") {
    ND280NamedInfo("TAnalysisLoop", " " << value);
    ND280NamedLog("TAnalysisLoop",
                  "Enabling all modules due to user command line option "
                      << "(takes precedence over the disable=all option)");
    fEnableAll = true;
    return true;
  }

  if ((option == "disable") || (option == "enable") ||
      (option == "preselection")) {
    ModDescripterator mod_it = FindNamedModule(fAnalysisModules, value);
    if (mod_it == fAnalysisModules.end()) {
      ND280NamedError("TAnalysisLoop",
                      "Module name " << value << " does not match any module.");
      return false;
    }

    ND::TAnalysisModuleBase *module = mod_it->second;

    if (option == "disable") {
      ND280NamedLog("TAnalysisLoop", "Disabling module "
                                         << value
                                         << " due to user cmd line option.");
      // Enabling takes precedence over disabling so only set to
      // false if not already set to true.
      if (!fUserEnabled.count(value)) {
        fUserEnabled[value] = false;
      }
    } else if (option == "enable") {
      ND280NamedLog("TAnalysisLoop", "Enabling module "
                                         << value
                                         << " due to user cmd line option ");
      // Enabling takes precedence over disabling.
      fUserEnabled[value] = true;
    } else if (option == "preselection") {
      ND280NamedLog("TAnalysisLoop", "Enabling module "
                                         << value << ", and using it for");
      ND280NamedLog("TAnalysisLoop",
                    "\toaEvent output preselection, "
                    "due to user cmd line option.");
      module->SetEnabled();
      module->SetUsedForPreselection();
      fSaveOutputOaEventTree = true;
    }
    return true;

  } else if (option == "save") {
    fSaveOriginalFullEvent = true;
    fSaveOutputOaEventTree = true;
  }
  // Add the input directory to all the modules.
  else if (option == "inputdir") {
    void *dir = gSystem->OpenDirectory(value.c_str());
    if (dir != NULL) {
      for (ModDescripterator mod_it = fAnalysisModules.begin();
           mod_it != fAnalysisModules.end(); ++mod_it) {
        mod_it->second->SetInputDirectory(value);
      }
      ND280NamedLog(
        "TAnalysisLoop",
        "Set input directory as: " << value.c_str());
      return true;
    }
    ND280NamedError(
        "TAnalysisLoop",
        "Unable to open specified input directory: " << value.c_str());
    return false;
  } else {
    // for options of the type TTruthTrajectoriesModule=SaveAll
    // Test if the option name is that of a module
    ModDescripterator mod_it = FindNamedModule(fAnalysisModules, option);
    if (mod_it == fAnalysisModules.end()) {
      ND280NamedError("TAnalysisLoop", "Module name "
                                           << option
                                           << " does not match any module.");
      return false;
    }
    return mod_it->second->Configure(value);
  }
  return false;
}

bool TAnalysisLoop::operator()(ND::TND280Event &event) {
  // Enables reporting of current event number while exiting quickly.
  fEventCount++;
  AnalysisLoopGlobals::EventNumber = fEventCount;
  if (AnalysisLoopGlobals::StopProcessing) {
    return false;
  }
  if ((fEventCount % 100) == 1) {
    ND280NamedInfo("TAnalysisLoop", "Processed Events " << fEventCount);
  }

  // Need to make sure we can get geoemtry before going through analysis
  // modules.
  fGeom = ND::TOADatabase::Get().Geometry();

  bool saveThisFullOaEvent = false;

  for (ModDescripterator mod_it = fAnalysisModules.begin();
       mod_it != fAnalysisModules.end(); ++mod_it) {
    ND::TAnalysisModuleBase &module =
        static_cast<TAnalysisModuleBase &>(*mod_it->second);

    if (!module.IsEnabled()) {
      continue;
    }

    ND280NamedVerbose("TAnalysisLoop", "Processing with module "
                                           << module.GetName());

    AnalysisLoopGlobals::CurrentlyProcessingModule = mod_it->first;

    bool passThroughException = false;
    try {
      if (fEventCount == 1) {
        // As a special case, call ProcessFirstEvent() for the first event.
        // This is a good time to save persistent quantities in modules'
        // data members, which will be retrievable from the eventAnalysis output
        // file.
        // Not intended for filling in the tree with the first event, as
        // Process() will also be called.
        module.ProcessFirstEvent(event);
      }
      // Call the module's own Process, passing it the event. If returns true
      // and module is selected for preselection then save original nd280
      // event as well as the output analysis trees.
      saveThisFullOaEvent =
          module.Process(event) && module.IsUsedForPreselection();
    } catch (ND::EDataFile &ex) {
      ND280NamedVerbose("TAnalysisLoop",
                        "Analysing data file, will disable "
                        "Monte Carlo module: "
                            << module.GetName() << " module");
      module.SetDisabled();
      continue;
    } catch (ND::ENoKinematicsPassThroughTrees &ex) {
      ND280NamedWarn("TAnalysisLoop",
                     "Caught ND::ENoKinematicsPassThroughTr"
                     "ees, "
                         << ex.what());
      passThroughException = true;
    } catch (ND::EIncorrectInputGenerator &ex) {
      ND280NamedWarn("TAnalysisLoop", "Caught ND::EIncorrectInputGenerator, "
                                          << ex.what());
      passThroughException = true;
    } catch (ND::ENoBeamSummaryData &ex) {
      if (fProduction) {
        fExit = true;
        ND280NamedError(
            "TAnalysisLoop",
            "No BSD data available for this "
            "Spill event - running in production mode - program will exit");
      }
    } catch (std::exception &ex) {
      ND280NamedError("TAnalysisLoop", ex.what());
      ND280NamedError("TAnalysisLoop", "Disabling module " << module.GetName());
      if (fProduction || fValidation) {
        fExit = true;
        ND280NamedError("TAnalysisLoop",
                        "eventAnalysis is running in production"
                        " mode and has caught an exception so will exit");
      }
      module.SetDisabled();
      continue;
    } catch (...) {
      ND280NamedError("TAnalysisLoop", "Disabling module " << module.GetName());
      if (fProduction || fValidation) {
        fExit = true;
        ND280NamedError("TAnalysisLoop",
                        "eventAnalysis is running in production"
                        " mode and has caught an exception so will exit");
      }
      module.SetDisabled();
      continue;
    }

    if (passThroughException) {
      ModDescripterator groo_mod_it =
          FindNamedModule(fAnalysisModules, "TGRooTrackerVtxModule");
      ModDescripterator nroo_mod_it =
          FindNamedModule(fAnalysisModules, "TNRooTrackerVtxModule");

      bool NRooIsEnabled = (nroo_mod_it != fAnalysisModules.end()) &&
                           (nroo_mod_it->second->IsEnabled());
      bool GRooIsEnabled = (groo_mod_it != fAnalysisModules.end()) &&
                           (groo_mod_it->second->IsEnabled());

      // check that the module pointers match up
      if ((nroo_mod_it->second == mod_it->second) && NRooIsEnabled) {
        ND280NamedVerbose("TAnalysisLoop",
                          "No neutrino interaction pass-through information"
                          " found for the "
                              << module.GetName() << " module, so it will be "
                                                     "disabled");
        nroo_mod_it->second->SetDisabled();
        ND280NamedLog("TAnalysisLoop", "Deleting pass-through tree for module: "
                                           << nroo_mod_it->second->GetName());
        nroo_mod_it->second->KillOutputTree();
        NRooIsEnabled = false;
      }

      // check that the module pointers match up
      if ((groo_mod_it->second == mod_it->second) && GRooIsEnabled) {
        ND280NamedVerbose("TAnalysisLoop",
                          "No neutrino interaction pass-through information"
                          " found for the "
                              << module.GetName() << " module, so it will be "
                                                     "disabled");
        groo_mod_it->second->SetDisabled();
        ND280NamedLog("TAnalysisLoop", "Deleting pass-through tree for module: "
                                           << groo_mod_it->second->GetName());
        groo_mod_it->second->KillOutputTree();
        GRooIsEnabled = false;
      }

      if (!NRooIsEnabled && !GRooIsEnabled) {
        ND280NamedWarn(
            "TAnalysisLoop",
            "No pass through information could be found from the "
            "neutrino generators. Both RooTracker modules have therefore "
            "been disabled so no generator information will be available in "
            "the eventAnalysis file.");
        if (fProduction) {
          fExit = true;
          ND280NamedError(
              "TAnalysisLoop",
              "No pass through information has been found for this"
              " MC file and eventAnalysis is running in production mode, so "
              "will exit");
        }
      }
      continue;
    }
  }

  if (fExit) {
    throw EProductionException();
  }

  if (saveThisFullOaEvent || fSaveOriginalFullEvent) {
    return true;
  } else {
    return false;
  }
}

void TAnalysisLoop::Initialize(void) {
  // Check that the output file is writable first
  if (!gFile) {
    ND280NamedError(
        "TAnalysisLoop",
        "Output file not specified. Please use the -o filename option");
    exit(1);
  }

  // Override default enable/disable based on usr cmd line options.
  for (ModDescripterator mod_it = fAnalysisModules.begin();
       mod_it != fAnalysisModules.end(); ++mod_it) {
    ND::TAnalysisModuleBase &module =
        static_cast<TAnalysisModuleBase &>(*mod_it->second);

    std::string className(mod_it->first);

    // Ignore enable/disable options if selected for preselection.
    if (!module.IsUsedForPreselection()) {
      if (fDisableAll) module.SetDisabled();  // Overrides default values.
      if (fEnableAll) module.SetEnabled();    // Overrides disable all.
      if (fUserEnabled.count(className)) {
        // Override enable/disable all.
        if (fUserEnabled[className]) {
          module.SetEnabled();
        } else {
          module.SetDisabled();
        }
      }
    }
  }

  // Move modules with fIsUsedForPreselection set to the end of the
  // fAnalysisModules vector
  std::stable_sort(fAnalysisModules.begin(), fAnalysisModules.end(),
                   PreselectionSorter);

  for (ModDescripterator mod_it = fAnalysisModules.begin();
       mod_it != fAnalysisModules.end(); ++mod_it) {
    ND::TAnalysisModuleBase &module =
        static_cast<TAnalysisModuleBase &>(*mod_it->second);

    module.Print();

    std::string dirname = module.GetDirectoryName();
    gFile->cd("/");
    gDirectory->cd("/");
    if (gFile && !gFile->GetDirectory(dirname.c_str())) {
      ND280NamedVerbose("TAnalysisLoop", "Making " << dirname);
      gDirectory->mkdir(dirname.c_str());
      gDirectory->cd(dirname.c_str());
      gDirectory->cd("/");
    }
    if (module.IsEnabled()) {
      gDirectory->cd(dirname.c_str());
      try {
        // Call the module's own Initialize, passing it its own tree
        module.Initialize(new TTree(module.GetName(), module.GetTitle()));
      } catch (std::exception &ex) {
        ND280NamedError("TAnalysisLoop", ex.what());
        ND280NamedError("TAnalysisLoop",
                        "Disabling module (unable to initialise) "
                            << module.GetName());
        module.SetDisabled();
        continue;
      } catch (...) {
        ND280NamedError("TAnalysisLoop",
                        "Disabling module (unable to initialise) "
                            << module.GetName());
        module.SetDisabled();
        continue;
      }
    }
  }
}

void TAnalysisLoop::BeginFile(ND::TVInputFile *vinput) {
  ND280NamedLog("TAnalysisLoop", "Begin reading input file");
  ND::TND280Input *input = dynamic_cast<ND::TND280Input *>(vinput);
  if (input) {
    ND280NamedLog("TAnalysisLoop", "Successfully cast TVInputFile");
    TFile *file = input->GetFilePointer();
    if (file) {
      ND280NamedLog("TAnalysisLoop", "Setting RooTracker module input files.");
      ModDescripterator groo_mod_it =
          FindNamedModule(fAnalysisModules, "TGRooTrackerVtxModule");
      if (groo_mod_it != fAnalysisModules.end()) {
        static_cast<TRooTrackerVtxModuleBase *>(groo_mod_it->second)
            ->SetBeginFile(file);
      }
      ModDescripterator nroo_mod_it =
          FindNamedModule(fAnalysisModules, "TNRooTrackerVtxModule");
      if (nroo_mod_it != fAnalysisModules.end()) {
        static_cast<TRooTrackerVtxModuleBase *>(nroo_mod_it->second)
            ->SetBeginFile(file);
      }
    }
  }
}

void TAnalysisLoop::EndFile(ND::TVInputFile *) {}

void TAnalysisLoop::Finalize(ND::TND280Output *const output) {
  // Write geometry if needed
  if (fGeometry) {
    output->WriteGeometry(fGeom);
  }
  // Testing [lp208:2016-03-10]
  output->cd("/");
  std::string utilsdir("Config");
  ND280NamedLog("TAnalysisLoop", "Writing Utils directory: \"" << utilsdir
                                                               << "\"");

  output->mkdir(utilsdir.c_str());
  if (output->cd(utilsdir.c_str())) {
    TTree *SubDetectorTree =
        new TTree("ESubdetector",
                  "Tree containing eventAnalysis subdetector enum Definitions.");

    int sdet;
    TObjString *sdetName = 0;
    SubDetectorTree->Branch("ESubdetector_Value", &sdet);
    SubDetectorTree->Branch("ESubdetector_Name", &sdetName);

    for (eventAnalysisEnums::ESubdetector sdet_it = eventAnalysisEnums::kFGD1;
         sdet_it != eventAnalysisEnums::kNSubdetectors; sdet_it = Next(sdet_it)) {
      (*sdetName) = eventAnalysisEnums::DetectorName(sdet_it).c_str();
      sdet = static_cast<int>(sdet_it);
      SubDetectorTree->Fill();
    }

    TTree *ParticleCategoryTree = new TTree(
        "EParticlecategory",
        "Tree containing eventAnalysis particle category enum Definitions.");

    int part;
    TObjString *partName = 0;
    ParticleCategoryTree->Branch("EParticleCategory_Value", &part);
    ParticleCategoryTree->Branch("EParticleCategory_Name", &partName);

    for (eventAnalysisEnums::EParticleCategory part_it =
             eventAnalysisEnums::kChargedLepton;
         part_it != eventAnalysisEnums::kNParticleCategories;
         part_it = Next(part_it)) {
      (*partName) = eventAnalysisEnums::ParticleName(part_it).c_str();
      part = static_cast<int>(part_it);
      ParticleCategoryTree->Fill();
    }

    for (ModDescripterator mod_it = fAnalysisModules.begin();
         mod_it != fAnalysisModules.end(); ++mod_it) {
      ND::TAnalysisModuleBase &module =
          static_cast<TAnalysisModuleBase &>(*mod_it->second);

      TTree *modConfigTree =
          new TTree(module.GetName(),
                    (std::string(module.GetName()) + "-ConfigTree").c_str());
      module.FillConfigTree(modConfigTree);
    }

  } else {
    ND280NamedError("TAnalysisLoop", "Failed to make Config directory.");
  }

  // Call the modules destructor to clean up memory resident objects so they
  // are not saved in the output file. For some reason std::list::clear() does
  // not call the destructors so doing it explicitly.
  //[ed] It's a list/map of pointers, the list doesn't own the memory
  // and has no right calling the destructor, it will correctly destruct the
  // pointers that it does own.
  for (ModDescripterator mod_it = fAnalysisModules.begin();
       mod_it != fAnalysisModules.end(); ++mod_it) {
    delete mod_it->second;
  }
  fAnalysisModules.clear();

  ND280NamedLog("TAnalysisLoop",
                "Total number of events processed in Analysis "
                "Format: "
                    << fEventCount);
  output->cd("/");
  if (!fSaveOutputOaEventTree) {
    output->Delete("ND280Events;*");
  }
}

void TAnalysisLoop::SetModuleDefaults() {
  for (ModDescripterator mod_it = fAnalysisModules.begin();
       mod_it != fAnalysisModules.end(); ++mod_it) {
    mod_it->second->SetEnabled(mod_it->second->IsEnabledByDefault());
  }
}

TAnalysisLoop::CModDescripterator TAnalysisLoop::FindNamedModule(
    std::vector<TAnalysisLoop::ModDescriptor> const &ModMap,
    std::string ModuleName) {
  for (TAnalysisLoop::CModDescripterator mod_it = ModMap.begin();
       mod_it != ModMap.end(); ++mod_it) {
    if (mod_it->first == ModuleName) {
      return mod_it;
    }
  }
  return ModMap.end();
}

TAnalysisLoop::ModDescripterator TAnalysisLoop::FindNamedModule(
    std::vector<TAnalysisLoop::ModDescriptor> &ModMap, std::string ModuleName) {
  for (TAnalysisLoop::ModDescripterator mod_it = ModMap.begin();
       mod_it != ModMap.end(); ++mod_it) {
    if (mod_it->first == ModuleName) {
      return mod_it;
    }
  }
  return ModMap.end();
}
}