CRTEventProducer_module.cc

Go to the documentation of this file.

// Framework includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Event.h" 
#include "canvas/Persistency/Common/Ptr.h" 
#include "canvas/Persistency/Common/PtrVector.h" 
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art_root_io/TFileService.h" 
#include "art_root_io/TFileDirectory.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Persistency/Common/PtrMaker.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include <memory>
#include <iostream>
#include <map>
#include <iterator>
#include <algorithm>

// LArSoft
#include "lardata/Utilities/AssociationUtil.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardataobj/RawData/ExternalTrigger.h"
#include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"

// ROOT
#include "TTree.h"
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TVector3.h"

#include "sbnobj/ICARUS/CRT/CRTData.hh"
#include "icaruscode/CRT/CRTDecoder/CRTPreProcessTree.h"

using std::string;

namespace icarus {
namespace crt {

  class CRTEventProducer : public art::EDProducer {
  public:
 
    explicit CRTEventProducer(fhicl::ParameterSet const & p);

    // The destructor generated by the compiler is fine for classes
    // without bare pointers or other resource use.

    // Plugins should not be copied or assigned.
    CRTEventProducer(CRTEventProducer const &) = delete;
    CRTEventProducer(CRTEventProducer &&) = delete;
    CRTEventProducer & operator = (CRTEventProducer const &) = delete; 
    CRTEventProducer & operator = (CRTEventProducer &&) = delete;

    // Required functions.
    void produce(art::Event & e) override;

    // Selected optional functions.
    void beginJob() override;

    void endJob() override;

    void reconfigure(fhicl::ParameterSet const & p);

  private:

   std::string fInFile;
   float       fTimeWindow;
   icarus::crt::CRTPreProcessTree* fAnaTree;

   std::map<size_t,size_t> fEventEntryStart;
   std::map<size_t,size_t> fEventEntryEnd;
 
  }; // class CRTEventProducer

  CRTEventProducer::CRTEventProducer(fhicl::ParameterSet const & p)
  : EDProducer{p}
  // Initialize member data here, if know don't want to reconfigure on the fly
  {
 
   // Call appropriate produces<>() functions here.
    produces< std::vector<icarus::crt::CRTData> >();
    
    reconfigure(p);

    /*std::cout << "opening file " << fInFile << std::endl;
    TFile fin(fInFile.c_str());
    TTree* tree=(TTree*)fin.FindObjectAny("anaTree");
    if(tree==nullptr) std::cout << "anaTree not found!" << std::endl;
    fAnaTree = new icarus::crt::CRTPreProcessTree(tree);
    std::cout << "constructed CRTPreProcessTree with " << fAnaTree->GetNEntries() << " entries" << std::endl;
    std::cout << "earliest timestamp: " << fAnaTree->GetAbsTime(0) << std::endl;
*/
  } // CRTEventProducer()

  void CRTEventProducer::reconfigure(fhicl::ParameterSet const & p)
  {
    //fCrtModuleLabel = (p.get<art::InputTag> ("CrtModuleLabel"));
    fInFile = p.get<string>("infile");
    fTimeWindow = p.get<float>("time_window");

  } // CRTEventProducer::reconfigure()

  void CRTEventProducer::beginJob()
  {
        std::cout << "opening file " << fInFile << std::endl;
        TFile fin(fInFile.c_str());
        TTree* tree=(TTree*)fin.FindObjectAny("anaTree");
        if(tree==nullptr) std::cout << "anaTree not found!" << std::endl;
        fAnaTree = new icarus::crt::CRTPreProcessTree(tree);

        std::cout << "running 'beginJob'" << std::endl;
        size_t event=0;
        uint64_t ti = fAnaTree->GetAbsTime(0);
        fEventEntryStart[0] = 0;

        std::cout << "time slicing tree with " << fAnaTree->GetNEntries() << " entries" << std::endl;

        for(size_t ientry =0; ientry<fAnaTree->GetNEntries(); ientry++){
                if(fAnaTree->GetAbsTime(ientry)>=ti &&
                   fAnaTree->GetAbsTime(ientry)<ti+fTimeWindow*1.0e9)   
                        continue;

                else if(fAnaTree->GetAbsTime(ientry)>=ti+fTimeWindow*1.0e9 ||
                        ientry == fAnaTree->GetNEntries()-1){
                        ti += fTimeWindow*1.0e9;
                        fEventEntryEnd[event] = ientry-1;
                        event++;
                        if (ientry < fAnaTree->GetNEntries()-1)
                                fEventEntryStart[event] = ientry;
                        
                }
        }

        std::cout << "found " << event << " time slices" << std::endl;
        if(fEventEntryEnd.size()!=fEventEntryStart.size())
                std::cout << "event map start/end size mismatch!" << std::endl;

        fin.Close();
        delete fAnaTree;

  } // CRTEventProducer::beginJob()
 
  void CRTEventProducer::produce(art::Event & event)
  {
    TFile fin(fInFile.c_str());
    TTree* tree=(TTree*)fin.FindObjectAny("anaTree");
    if(tree==nullptr) std::cout << "anaTree not found!" << std::endl;
    fAnaTree = new icarus::crt::CRTPreProcessTree(tree);

    std::unique_ptr< std::vector<icarus::crt::CRTData> > crtdata( new std::vector<icarus::crt::CRTData>);

    size_t eveId = (size_t) event.event();
    std::cout << "processing event " << eveId << " with " << fEventEntryEnd[eveId]+1-fEventEntryStart[eveId] << " triggers" << std::endl;

    int entry=0;
    for(size_t ientry=fEventEntryStart[eveId]; ientry<=fEventEntryEnd[eveId]; ientry++){
        fAnaTree->Load(ientry);
        CRTData data;
        data.fMac5 = fAnaTree->Mac5();
        data.fTs0 = fAnaTree->GetAbsTime();
        data.fTs1 = data.fTs0;
        for(size_t ichan=0; ichan<32; ichan++){
                if(fAnaTree->Above(ichan)&&fAnaTree->Active(ichan)) {
                        data.fAdc[ichan]  = fAnaTree->PE(ichan);
                }//if channel active and above threshold
                else
                        data.fAdc[ichan] = 0;
        }//for channels

        crtdata->push_back(data);
        entry++;

    }//for entries in time slice

    delete fAnaTree;
    fin.Close();

    event.put(std::move(crtdata));

  } // CRTEventProducer::produce()

  void CRTEventProducer::endJob()
  {

  } // CRTEventProducer::endJob()

  DEFINE_ART_MODULE(icarus::crt::CRTEventProducer)

}
} //end namespace