OpHitAna_module.cc

Go to the documentation of this file.

// Christie Chiu and Ben Jones, MIT, 2012
//
// This is an analyzer module which writes the raw optical
// detector pulses for each PMT to an output file
//

// LArSoft includes
#include "larana/OpticalDetector/OpDigiProperties.h"
#include "larcore/Geometry/WireReadout.h"
#include "lardataobj/RecoBase/OpHit.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art_root_io/TFileService.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// ROOT includes
#include "TH1D.h"
#include "TLorentzVector.h"
#include "TTree.h"

// C++ Includes
#include "math.h"

namespace opdet {

  class OpHitAna : public art::EDAnalyzer {
  public:
    // Standard constructor and destructor for an ART module.
    OpHitAna(const fhicl::ParameterSet&);

    // The analyzer routine, called once per event.
    void analyze(const art::Event&);

  private:
    // The stuff below is the part you'll most likely have to change to
    // go from this custom example to your own task.

    // The parameters we'll read from the .fcl file.
    std::string fInputModule; // Input tag for OpDet collection
    float fSampleFreq;        // in MHz
    float fTimeBegin;         // in us
    float fTimeEnd;           // in us
                              //  short fPEheight;                       // in ADC counts

    // Flags to enable or disable output of debugging TH1 / TH2s
    bool fMakeHistPerChannel;
    bool fMakeHistAllChannels;
    bool fMakeHitTree;

    // Output TTree and its branch variables
    TTree* fHitTree;
    Int_t fEventID;
    Int_t fOpChannel;
    Float_t fPeakTime;
    Float_t fNPe;
  };

}

// OpHitAna_module.cc

// This is a short program required by the ART framework.  It enables
// a program (OpHitAna, in this case) to be called as a module
// from a .fcl file. It is unlikely that you'll have to make any
// changes to this file.

namespace opdet {
  DEFINE_ART_MODULE(OpHitAna)
}

namespace opdet {

  //-----------------------------------------------------------------------
  // Constructor
  OpHitAna::OpHitAna(fhicl::ParameterSet const& pset) : EDAnalyzer(pset)
  {

    // Indicate that the Input Module comes from .fcl
    fInputModule = pset.get<std::string>("InputModule");

    art::ServiceHandle<OpDigiProperties> odp;
    fTimeBegin = odp->TimeBegin();
    fTimeEnd = odp->TimeEnd();
    fSampleFreq = odp->SampleFreq();

    fMakeHistPerChannel = pset.get<bool>("MakeHistPerChannel");
    fMakeHistAllChannels = pset.get<bool>("MakeHistAllChannels");
    fMakeHitTree = pset.get<bool>("MakeHitTree");

    // If required, make TTree for output

    if (fMakeHitTree) {
      art::ServiceHandle<art::TFileService const> tfs;
      fHitTree = tfs->make<TTree>("HitTree", "HitTree");
      fHitTree->Branch("EventID", &fEventID, "EventID/I");
      fHitTree->Branch("OpChannel", &fOpChannel, "OpChannel/I");
      fHitTree->Branch("PeakTime", &fPeakTime, "PeakTime/F");
      fHitTree->Branch("NPe", &fNPe, "NPe/F");
    }
  }

  //-----------------------------------------------------------------------
  void OpHitAna::analyze(const art::Event& evt)
  {
    // Create a handle for our vector of pulses
    art::Handle<std::vector<recob::OpHit>> HitHandle;

    // Create string for histogram name
    char HistName[50];

    // Read in HitHandle
    evt.getByLabel(fInputModule, HitHandle);

    // Access ART's TFileService, which will handle creating and writing
    // histograms for us.
    art::ServiceHandle<art::TFileService const> tfs;

    int NOpChannels = art::ServiceHandle<geo::WireReadout const>()->Get().NOpChannels();

    std::vector<TH1D*> HitHist;

    sprintf(HistName, "Event %d AllOpDets", evt.id().event());

    TH1D* AllHits = nullptr;
    if (fMakeHistAllChannels) {
      AllHits = tfs->make<TH1D>(HistName,
                                ";t (ns);",
                                int((fTimeEnd - fTimeBegin) * fSampleFreq),
                                fTimeBegin * 1000.,
                                fTimeEnd * 1000.);
    }

    for (int i = 0; i != NOpChannels; ++i) {

      sprintf(HistName, "Event %d OpDet %i", evt.id().event(), i);
      if (fMakeHistPerChannel)
        HitHist.push_back(tfs->make<TH1D>(HistName,
                                          ";t (ns);",
                                          int((fTimeEnd - fTimeBegin) * fSampleFreq),
                                          fTimeBegin * 1000.,
                                          fTimeEnd * 1000.));
    }

    fEventID = evt.id().event();

    // For every OpHit in the vector
    for (unsigned int i = 0; i < HitHandle->size(); ++i) {
      // Get OpHit
      art::Ptr<recob::OpHit> TheHitPtr(HitHandle, i);
      recob::OpHit TheHit = *TheHitPtr;

      fOpChannel = TheHit.OpChannel();
      fNPe = TheHit.PE();
      fPeakTime = TheHit.PeakTime();

      if (fMakeHitTree) fHitTree->Fill();

      if (fMakeHistPerChannel) {
        if (fOpChannel > int(HitHist.size())) {
          mf::LogError("OpHitAna") << "Error : Trying to fill channel out of range, " << fOpChannel;
        }
        HitHist[fOpChannel]->Fill(fPeakTime, fNPe);
      }

      if (fMakeHistAllChannels) AllHits->Fill(fPeakTime, fNPe);
    }
  }

} // namespace opdet