Project

General

Profile

Writing LArIAT Reconstruction and Analysis Code

The ART framework is based on a hierarchy of locations where data records are stored in a ROOT file. The lowest level is the art::Event, which in LArIAT corresponds to a single spill of the beam. The next level up is the art::SubRun which stores information for data taking where the conditions of the detector and DAQ are assumed to be static. The highest level is the art::Run, which corresponds to several SubRuns - minor changes can happen from sub-run to sub-run, but the general state of the detector and DAQ configuration is consistent through a run. The ART framework cannot randomly access the art::Event, art::SubRun or art::Run records. Instead, it has to move sequentially forward in time from one record to the next.

LArIAT data are different from the data of other experiments using LArSoft because our spills are much longer than other experiments, 4 seconds compared to a few microseconds for neutrino beam experiments, and because each spill can have many different trigger conditions. For this reason, we have a unique requirement for reconstructing and analyzing our data, namely the need to separate our data into different triggers in each spill. In this context, a trigger corresponds to all raw data fragments coming from the DAQ that come from the same time and satisfy a condition identified by the V1495.

The technical solution to separating out these triggers for reconstruction and analysis is to create a LArSoft raw::Trigger object for each trigger and then to associate each trigger to the LArSoft raw::RawDigit, raw::AuxDetDigit and raw::OpDetPulse objects that are "in-time" with each other. The raw::RawDigits come from the TPC, the raw::AuxDetDigits come from the various detectors in the beam line that are not in the cryostat, and the raw::OpDetPulses are for the photodetectors in the cryostat.

This solution allows LArIAT analyzers to have access to all the activity from a single spill, while at the same time dividing that activity up into relevant chunks based on timing and other trigger conditions. This kind of access is necessary for studying issues of pile up in the TPC where many different trigger conditions may leave energy deposition during the same readout window of the TPC. The ART framework cannot provide that kind of access if we split the triggers into different art::Event records, so instead we keep a single art::Event record for each spill and subdivide it ourselves in time using the raw::Trigger objects. Doing so we have random access of any raw::Trigger in the spill.

Taking DAQ Fragments to Triggers, Digits and Pulses

The DAQ writes objects known as Fragments into the record of each art::Event. Each file produced by the DAQ corresponds to a single spill or art::Event. The first step in the reconstruction and analysis is to convert these Fragments into the various LArSoft data products. The code that does this conversion is located in the LArIATSoft/RawDataUtilities package and the module is the source:RawDataUtilities/FragmentToDigit_module.cc. This module collects the Fragments into groupings determined by the timestamps of the fragments. It then makes a raw::Trigger object for the Fragments from a given time and converts the Fragments into the LArSoft data products. At the time the LArSoft data products are created, an ART Association is also created to link the raw::Trigger object to the other LArSoft data products.

After the FragmentToDigit module has been run for an art::Event, a single collection for each of the following data types, raw::Triggers, raw::RawDigits, raw::AuxDetDigits, and raw::OpDetPulses, is written into the art::Event record. Associations linking each raw::Trigger in the collection to the corresponding raw::RawDigits, raw::AuxDetDigits, and raw::OpDetPulses are also written into the art::Event. These Associations allow downstream modules to access the data products for each trigger individually.

Accessing Triggers, Digits and Pulses

As stated in the previous section, each art::Event contains a collection of raw::Triggers and associations between the raw::Triggers and other data products. A utility object, source:RawDataUtilities/TriggerDigitUtility.h, has been written to provide convenient access to the LArSoft raw data products in each art::Event. The TriggerDigitUtility gives a user access to all raw::Trigger, raw::RawDigit, raw:AuxDetDigit, and raw::OpDetPulse objects in an art::Event, as well as storing the mapping between each raw::Trigger and the other data products. The TriggerDigitUtility has interfaces for accessing both constant pointers to the various collections of objects and art::PtrVectors of the collections of objects.

Below is an example of how one can access the TriggerDigitUtility from within the art::EDProducer::produce, art::EDAnalyzer::analyze, or art::EDFilter::filter method.

// instantiate the TriggerDigitUtility object.  
// evt is the art::Event
// triggerModuleLabel is a std::string containing the label of the module creating the raw::Trigger objects
rdu::TriggerDigitUtility tdu(evt, triggerModuleLabel);

// To loop over triggers in an event and grab the raw::RawDigits for each, you would do
for(size_t t = 0; t < tdu.NTriggers(); ++t){
     const raw::Trigger* trig = tdu.EventTriggers()[t];

     // get the raw digits for this trigger
     std::vector<const raw::RawDigit*> rdvec = tdu.TriggerRawDigits(t);

     // do stuff with the RawDigits and Trigger
}

Practical Implications for LArIATSoft

The primary consideration for people wanting to reconstruct or analyze the data is that the LArSoft modules will not work correctly for LArIAT because those modules do not take different raw::Triggers into account. Instead, those modules assume that all the raw::RawDigits, raw::AuxDetDigits, and raw::OpDetPulses in a given art::Event record correspond to a single trigger condition. Fortunately, LArSoft modules are simply wrappers of other objects that do the reconstruction and analysis and those objects can be called for a single trigger. LArIAT analyzers will have to write new modules to wrap the LArSoft reconstruction objects and call them for individual triggers.

Example LArIAT EDProducer Module

Here is an example of how one could write an EDProducer module to

#ifndef MYPRODUCER_H
#define MYPRODUCER_H

// c++ includes
#include <iostream>

// Framework includes
#include "art/Framework/Core/ModuleMacros.h" 
#include "art/Framework/Core/EDProducer.h" 
#include "art/Framework/Principal/Event.h" 
#include "art/Framework/Principal/Handle.h" 
#include "art/Persistency/Common/Ptr.h" 
#include "art/Persistency/Common/PtrVector.h" 
#include "art/Framework/Services/Registry/ServiceHandle.h" 
#include "art/Framework/Services/Optional/TFileService.h" 
#include "art/Framework/Services/Optional/TFileDirectory.h" 
#include "fhiclcpp/ParameterSet.h" 
#include "messagefacility/MessageLogger/MessageLogger.h" 
#include "cetlib/exception.h" 
#include "cetlib/search_path.h" 

// LArSoft Includes
#include "MyProducts/MyProduct.h" 
#include "MyProducts/MyOtherProduct.h" 
#include "RawData/TriggerData.h" 
#include "Utilities/AssociationUtil.h" 

// ROOT Includes
#include "TH1.h" 

// LArIATSoft Includes
#include "RawDataUtilities/TriggerDigitUtility.h" 

namespace mp {

  ///class to produce a data product
  class MyProducer : public art::EDProducer {

  public:

    explicit MyProducer(fhicl::ParameterSet const& pset);  // the explicit tag tells the compiler that MyProducer is different from fhicl::ParameterSet
    virtual ~MyProducer();
    void produce(art::Event& evt);                         // makes the objects I want made
    void beginJob();                                      // does things like make histograms at the beginning of the job
    void beginRun(art::Run& run);                         // Happens before each run, so this is a good place to pick up run dependent variables like Lifetime etc. 
    void reconfigure(fhicl::ParameterSet const& pset);    // every module should have one of these to allow the event display to alter configurations

  private:

    double      fDouble;                      ///< some data member
    double      fRunProperty;             ///< a property that changes from run to run
    std::string fPrevModuleLabel;     ///< label of the module making objects that i need for this step
    std::string fTriggerModuleLabel; ///< label of the module making the trigger objects
    TH1D*       fHist;                         ///< a histogram data member

  }; // class MyProducer

  //-----------------------------------------------------
  MyProducer::MyProducer(fhicl::ParameterSet const& pset)
  {
    this->reconfigure(pset);

    // the next line tells the module what it is going to make, you must have one line for each type of collection to be made
    produces< std::vector<mp::MyOtherProduct> >(); 

    // make an association between something already in the event record and the new something you are making with this module 
    produces< art::Assns<mp::MyProduct, mp::MyOtherProdut> >(); 

     // make an association between the Trigger already in the event record and the new something you are making with this module
    produces< art::Assns<raw::Trigger, mp::MyOtherProdut> >();

  }

  //-----------------------------------------------------
  MyProducer::~MyProducer()
  {
  }

  //-----------------------------------------------------
  void MyProducer::reconfigure(fhicl::ParameterSet const& pset)
  {
    fDouble                     = pset.get< double >     ("TheDouble");
    fPrevModuleLabel     = pset.get< std::string >("PreviousModuleLabel");
    fTriggerModuleLabel = pset.get< std::string >("TriggerModuleLabel");

    return;
  }

  //-----------------------------------------------------
  void MyProducer::beginJob()
  {
    // get access to the TFile service
    art::ServiceHandle<art::TFileService> tfs;

    fHist    = tfs->make<TH1D>("channels",  ";channels;",  100, 0, 100);

    return;
  }

  //-----------------------------------------------------
  void MyProducer::beginRun(art::Run& run)
  {
  }

  //-----------------------------------------------------
  void MyProducer::produce(art::Event& evt)
  {

    // declare the collection you are eventually going to put into the art::Event
    std::unique_ptr<std::vector<mp::MyOtherProduct> > mopcol(new std::vector<mp::MyOtherProduct>);

   // declare the association you are eventually going to put into the art::Event
   std::unique_ptr< art::Assns<mp::MyProduct, mp::MyOtherProduct> > mpmopAssn(new art::Assns<mp::MyProduct, mp::MyOtherProduct>);
   std::unique_ptr< art::Assns<raw::Trigger, mp::MyOtherProduct> > trmopAssn(new art::Assns<mp::MyProduct, mp::MyOtherProduct>);

   // grab the triggers from the event
   rdu::TriggerDigitUtility tdu(evt, fTriggerModuleLabel);

   // we already made an association between raw::Triggers and mp::MyProducts in the module with the label 
   // fPrevModuleLabel
   // the calls to FindOne and FindMany should always happen outside of loops as they are time consuming
   art::FindManyP<mp::MyProduct> fmpmp(tdu.EventTriggersPtr(), evt, fPrevModuleLabel);

   // loop over the triggers
   auto trigPtrVec = tdu.EventTriggersPtr();
   for(size_t t = 0; t < trigPtrVec.size(); ++t ){

    // get the mp::MyProduct collection for this trigger
    // the type of mpCol is std::vector<art::Ptr<mp::MyProduct> >
    auto mpCol = fmpmp.at(t);

    // loop over them
    for(auto mpPtr : mpCol){

      // do something to turn out MyOtherProducts - not shown here

      // add the new MyOtherProduct to the collection
      mopcol->push_back(myOtherProd)

      // make an association between the trigger and the latest member of mopcol
      util::CreateAssn(*this, evt, *mopcol, trigPtrVec[t], *trmopAssn);

      // make an association between the latest member of mopcol and mpPtr
      util::CreateAssn(*this, evt, *mopcol, mpPtr, *mpmopAssn);
    } // end loop over mpCol
   } // end loop over triggers

    // put the newly created object collections and associations into the event record
    evt.put(std::move(mopcol));
    evt.put(std::move(trmopAssn));
    evt.put(std::move(mpmopAssn));

    // mopcol, trmopAssn, and mpmopAssn are no longer a valid pointer - don't use them any more!

    return;

  }

  // A macro required for a JobControl module.
  DEFINE_ART_MODULE(MyProducer);

} // namespace mp

#endif //MYPRODUCER_H