PhotonInf_module.cc

Go to the documentation of this file.

// Plugin Type: analyzer
// File:        PhotonInf_module.cc
// Description:
// takes the SimPhotonsCollection generated by LArG4's sensitive detectors
// and gives: x, y, z, num of det_ph on each OpDet
// Oct. 20, 2020 by Mu Wei wmu@fnal.gov

// FMWK 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 "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft includes
#include "larana/OpticalDetector/OpDetResponseInterface.h"
#include "larcore/Geometry/Geometry.h"
#include "larcorealg/Geometry/CryostatGeo.h"
#include "larcorealg/Geometry/OpDetGeo.h"
#include "lardataobj/Simulation/SimPhotons.h"
#include "lardataobj/Simulation/sim.h"
#include "larsim/MCCheater/ParticleInventoryService.h"

#include "nug4/ParticleNavigation/ParticleList.h"
#include "nusimdata/SimulationBase/MCParticle.h"

// ROOT includes
#include "RtypesCore.h"
#include "TLorentzVector.h"
#include "TTree.h"
#include "TVector3.h"

// C++ language includes
#include <cstring>
#include <iostream>

namespace opdet {

  class PhotonInf : public art::EDAnalyzer {
  public:
    PhotonInf(const fhicl::ParameterSet&);
    void analyze(art::Event const&);
    void beginJob();
    void endJob();

    std::ofstream logFile;

  private:
    bool FillTree;
    std::string LogFileName;               // Input tag for name of log file
    std::vector<std::string> InputModules; // Input tag for OpDet collection

    std::vector<int> num_detph; // number of directly detected photons
    std::vector<int> num_refph; // number of reflected photons
    std::vector<int> num_totph; // number of total detected photons
    double pos[3];

    TTree* PhInf;

    int OpChannel;
    int nOpChannels;
    int EventID;
    float Energy;
    //float Time;
  };

  PhotonInf::PhotonInf(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
    , FillTree{pset.get<bool>("FillTree")}
    , LogFileName{pset.get<std::string>("LogFileName")}
  {
    InputModules = pset.get<std::vector<std::string>>("InputModules", {"largeant"});

    art::ServiceHandle<geo::Geometry const> geo;
    nOpChannels = int(geo->Cryostat().NOpDet());

    num_detph.resize(nOpChannels, 0);
    num_refph.resize(nOpChannels, 0);
    num_totph.resize(nOpChannels, 0);

    if (FillTree) {
      art::ServiceHandle<art::TFileService const> tfs;
      PhInf = tfs->make<TTree>("PhInf", "PhInf");
      PhInf->Branch("X", &(pos[0]), "X/D");
      PhInf->Branch("Y", &(pos[1]), "Y/D");
      PhInf->Branch("Z", &(pos[2]), "Z/D");

      for (int channel = 0; channel < nOpChannels; ++channel) {
        PhInf->Branch(Form("Channel_%03d_det", channel),
                      &(num_detph[channel]),
                      Form("Channel_%03d_det/I", channel));
        PhInf->Branch(Form("Channel_%03d_ref", channel),
                      &(num_refph[channel]),
                      Form("Channel_%03d_ref/I", channel));
        PhInf->Branch(Form("Channel_%03d_tot", channel),
                      &(num_totph[channel]),
                      Form("Channel_%03d_tot/I", channel));
      }
    }
  }

  void PhotonInf::beginJob()
  {
    logFile.open(LogFileName);

    return;
  }

  void PhotonInf::endJob()
  {
    logFile.close();

    return;
  }

  void PhotonInf::analyze(art::Event const& evt)
  {
    art::EventNumber_t event = evt.id().event();
    EventID = int(event);

    num_detph.resize(nOpChannels, 0);
    num_refph.resize(nOpChannels, 0);
    num_totph.resize(nOpChannels, 0);

    pos[0] = 0.0;
    pos[1] = 0.0;
    pos[2] = 0.0;

    bool rec = false;
    float vuv = 6e-6;

    auto photon_handles = evt.getMany<std::vector<sim::SimPhotons>>();
    if (photon_handles.empty()) {
      throw art::Exception(art::errors::ProductNotFound)
        << "sim SimPhotons retrieved and you requested them.";
    }

    for (auto const& mod : InputModules) {
      for (auto const& ph_handle : photon_handles) {
        if (!ph_handle.isValid()) { continue; }
        if (ph_handle.provenance()->moduleLabel() != mod) { continue; }

        if (!ph_handle->empty()) {
          for (auto const& itOpDet : (*ph_handle)) {
            OpChannel = int(itOpDet.OpChannel());
            const sim::SimPhotons& TheHit = itOpDet;
            for (const sim::OnePhoton& Phot : TheHit) {
              Energy = Phot.Energy;

              num_totph[OpChannel] += 1;

              if (Energy < vuv) // reflected photons
              {
                num_refph[OpChannel] += 1;
              }
              else {
                num_detph[OpChannel] += 1;

                if (rec == false) {
                  pos[0] = Phot.InitialPosition.X() / 10.0; // cm
                  pos[1] = Phot.InitialPosition.Y() / 10.0;
                  pos[2] = Phot.InitialPosition.Z() / 10.0;

                  rec = true;
                }
              }
            }
          }
        }
      }
    }

    if (rec == true) {
      if (FillTree) { PhInf->Fill(); }

      logFile << pos[0] << ", " << pos[1] << ", " << pos[2];
      for (int channel = 0; channel < nOpChannels; channel++) {
        logFile << ", " << num_detph[channel] << ", " << num_refph[channel] << ", "
                << num_totph[channel];

        num_detph[channel] = 0;
        num_refph[channel] = 0;
        num_totph[channel] = 0;
      }
      logFile << std::endl;
      pos[0] = 0.0;
      pos[1] = 0.0;
      pos[2] = 0.0;
    }
    else {
      std::cout << "Warning: No record for event " << EventID << std::endl;
    }

    rec = false;
    return;
  }
}

DEFINE_ART_MODULE(opdet::PhotonInf)