LArG4Ana_module.cc

Go to the documentation of this file.

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"

// Framework includes
#include "art/Framework/Principal/Event.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 "larcore/Geometry/Geometry.h"
#include "larcore/Geometry/WireReadout.h"
#include "lardataobj/Simulation/SimChannel.h"
#include "lardataobj/Simulation/sim.h"
#include "larsim/MCCheater/ParticleInventoryService.h"
#include "nug4/ParticleNavigation/ParticleList.h"

// ROOT includes
#include "TH1.h"
#include "TLorentzVector.h"
#include "TProfile.h"
#include <TTree.h>

// C++ Includes
#include <cstring>

namespace larg4 {

  class LArG4Ana : public art::EDAnalyzer {
  public:
    explicit LArG4Ana(fhicl::ParameterSet const& pset);

  private:
    void analyze(const art::Event& evt) override;
    void beginJob() override;

    std::string fG4ModuleLabel;    
    std::string fTruthModuleLabel; 

    TH1D* fPDGCodes;
    TH1D* fPi0Momentum;
    TH1D* fnEnergy;
    TH1D* fnDist;
    TH1D* fnumChannels;       
    TProfile* fnumIDEs;       
    TH1D* fEventCharge;       
    TH1D* fEventEnergy;       
    TProfile* fChannelCharge; 
    TProfile* fChannelEnergy; 

    //    Int_t stringDim = 35;

    TTree* fTree;
    Int_t fTEvt;
    Int_t fTSub;
    Int_t fTRun;
    Int_t fTPdg;
    Int_t fTID;
    Int_t fTNdsOriginal;
    Int_t fTNds;
    Int_t fTNds4;
    Int_t* fTDID;
    Int_t* fTDPdg;
    Float_t* fTDWt;
    Char_t fTProcess[35];
    Char_t fTVolume[35];
    Char_t fTTVolume[35]; // Termination Volume
    Char_t fTMaterial[35];
    Char_t fTDProcess[200][35];
    Int_t fTParentID;
    Int_t fTStatus;
    Float_t fTWeight;
    Float_t* fT4Origin;
    Float_t* fT4DOrigin;
    Float_t* fT4Termination; // Termination Coordinates
    Float_t* fT4Momentum;
    Float_t* fT4DMomentum;
  };

} // namespace larg4

namespace larg4 {

  //-----------------------------------------------------------------------
  // Constructor
  LArG4Ana::LArG4Ana(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
    , fG4ModuleLabel{pset.get<std::string>("GeantModuleLabel")}
    , fTruthModuleLabel{pset.get<std::string>("TruthModuleLabel")}
    , fTNdsOriginal{pset.get<int>("Ndaughters")}
    , fTNds{fTNdsOriginal}
  {}

  //-----------------------------------------------------------------------
  void LArG4Ana::beginJob()
  {
    art::ServiceHandle<art::TFileService const> tfs;

    fPDGCodes = tfs->make<TH1D>("pdgcodes", ";PDG Code;", 5000, -2500, 2500);
    fPi0Momentum = tfs->make<TH1D>("pi0mom", ";#pi^{0} Momentum (GeV);", 1000, 0., 1000.);

    fTree = tfs->make<TTree>("MCTTree", "MCTTree");
    fnEnergy = tfs->make<TH1D>("nEnergy", ";n,#Lambda^{0},K^{0} Momentum (GeV);", 100, 0., 10.);
    fnDist =
      tfs->make<TH1D>("nDistance", ";n,#Lambda^{0},K^{0} Distance (m);", 200, -30000.0, +30000.);

    // Some histograms relating to drift electrons, active detector
    // channels and charge/energy on channels
    auto const nchannels = art::ServiceHandle<geo::WireReadout const>()->Get().Nchannels();
    fnumChannels = tfs->make<TH1D>(
      "fnumChannels", "Active channels;Active channels;# events", 256, 0, nchannels);
    fnumIDEs = tfs->make<TProfile>("fnumIDEs",
                                   "Drift Electrons per channel;Channel;Drift electrons",
                                   nchannels + 1,
                                   0,
                                   nchannels,
                                   0,
                                   1e4);
    fEventCharge = tfs->make<TH1D>(
      "fEventCharge", "Charge in event;Total charge per event;# events", 100, 0, 2.5e8);
    fEventEnergy = tfs->make<TH1D>(
      "fEventEnergy", "Energy in event;Total energy per event;# events", 100, 0, 1e4);
    fChannelCharge = tfs->make<TProfile>("fChannelCharge",
                                         "Charge on channel;Channel;Total charge per channel",
                                         nchannels + 1,
                                         0,
                                         nchannels,
                                         0,
                                         1e5);
    fChannelEnergy = tfs->make<TProfile>("fChannelEnergy",
                                         "Energy on channel;Channel;Total energy per channel",
                                         nchannels + 1,
                                         0,
                                         nchannels,
                                         0,
                                         1e3);

    fT4Origin = new Float_t[4];
    fT4DOrigin = new Float_t[fTNds * 4];
    fT4Termination = new Float_t[4];
    fT4Momentum = new Float_t[4];
    fT4DMomentum = new Float_t[fTNds * 4];
    fTDID = new Int_t[fTNds];
    fTDPdg = new Int_t[fTNds];
    fTDWt = new Float_t[fTNds];
    fTNds4 = fTNds * 4; //  TTree/Branch requirement to store this.

    fTree->Branch("MCEvt", &fTEvt, "MCEvt/I");
    fTree->Branch("MCSub", &fTSub, "MCSub/I");
    fTree->Branch("MCRun", &fTRun, "MCRun/I");
    fTree->Branch("MCWt", &fTWeight, "MCWt/F");
    fTree->Branch("MCPdg", &fTPdg, "MCPdg/I");
    fTree->Branch("MCID", &fTID, "MCID/I");
    fTree->Branch("MCParentID", &fTParentID, "MCParentID/I");
    fTree->Branch("MCNumDs", &fTNds, "MCNumDs/I");
    fTree->Branch("MCNumDs4", &fTNds4, "MCNumDs4/I");
    fTree->Branch("MCDID", fTDID, "MCDID[MCNumDs]/I");
    fTree->Branch("MCDPdg", fTDPdg, "MCDPdg[MCNumDs]/I");
    fTree->Branch("MCDWt", fTDWt, "MCDWt[MCNumDs]/I");
    fTree->Branch("MCProcess", fTProcess, "MCProcess/C");
    fTree->Branch("MCVolume", fTVolume, "MCVolume/C");
    fTree->Branch("MCTVolume", fTTVolume, "MCTVolume/C");
    fTree->Branch("MCMaterial", fTMaterial, "MCMaterial/C");
    fTree->Branch("MCDProcess", fTDProcess, "MCDProcess[MCNumDs]/C");
    fTree->Branch("MCStatus", &fTStatus, "MCStatus/I");
    fTree->Branch("MCOrigin", fT4Origin, "MCOrigin[4]/F");
    fTree->Branch("MCDOrigin", fT4DOrigin, "MCDOrigin[MCNumDs4]/F");
    fTree->Branch("MCTermination", fT4Termination, "MCTermination[4]/F");
    fTree->Branch("MCMomentum", fT4Momentum, "MCMomentum[4]/F");
    fTree->Branch("MCDMomentum", fT4DMomentum, "MCDMomentum[MCNumDs4]/F");
  }

  //-----------------------------------------------------------------------
  void LArG4Ana::analyze(const art::Event& evt)
  {
    //get the list of particles from this event
    art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
    const sim::ParticleList& plist = pi_serv->ParticleList();
    art::ServiceHandle<geo::Geometry const> geom;

    // loop over all sim::SimChannels in the event and make sure there are no
    // sim::IDEs with trackID values that are not in the sim::ParticleList
    std::vector<const sim::SimChannel*> sccol;
    evt.getView(fG4ModuleLabel, sccol);

    double totalCharge = 0.0;
    double totalEnergy = 0.0;
    fnumChannels->Fill(sccol.size());
    for (size_t sc = 0; sc < sccol.size(); ++sc) {
      double numIDEs = 0.0;
      double scCharge = 0.0;
      double scEnergy = 0.0;
      const auto& tdcidemap = sccol[sc]->TDCIDEMap();
      for (auto mapitr = tdcidemap.begin(); mapitr != tdcidemap.end(); mapitr++) {
        const std::vector<sim::IDE> idevec = (*mapitr).second;
        numIDEs += idevec.size();
        for (size_t iv = 0; iv < idevec.size(); ++iv) {
          if (plist.find(idevec[iv].trackID) == plist.end() &&
              idevec[iv].trackID != sim::NoParticleId)
            mf::LogWarning("LArG4Ana") << idevec[iv].trackID << " is not in particle list";
          totalCharge += idevec[iv].numElectrons;
          scCharge += idevec[iv].numElectrons;
          totalEnergy += idevec[iv].energy;
          scEnergy += idevec[iv].energy;
        }
      }
      fnumIDEs->Fill(sc, numIDEs);
      fChannelCharge->Fill(sc, scCharge);
      fChannelEnergy->Fill(sc, scEnergy);
    }
    fEventCharge->Fill(totalCharge);
    fEventEnergy->Fill(totalEnergy);

    // get the particles from the back tracker
    const sim::ParticleList& Particles = pi_serv->ParticleList();
    std::vector<const simb::MCParticle*> pvec;
    pvec.reserve(Particles.size());
    for (const auto& PartPair : Particles) {
      pvec.push_back(PartPair.second);
      fPDGCodes->Fill(PartPair.second->PdgCode());
    }

    // now look for pi0's that decay to 2 gammas
    int pi0loc = -1;
    int numpi0gamma = 0;
    for (unsigned int i = 0; i < pvec.size(); ++i) {
      if (pvec[i]->PdgCode() == 111) pi0loc = i;
      if (pvec[i]->Mother() == pi0loc + 1 && pi0loc > 0 && pvec[i]->PdgCode() == 22) {
        mf::LogInfo("LArG4Ana") << pvec[i]->E() << " gamma energy ";
        ++numpi0gamma;
      }

      // n,Lambda,K0s,K0L,K0
      if (pvec[i]->PdgCode() == 2112 || pvec[i]->PdgCode() == 3122 || pvec[i]->PdgCode() == 130 ||
          pvec[i]->PdgCode() == 310 || pvec[i]->PdgCode() == 311) {
        fnEnergy->Fill(pvec[i]->E(), pvec[i]->Weight());
        fnDist->Fill(pvec[i]->Vx(), pvec[i]->Weight());
      }

      fTPdg = pvec[i]->PdgCode();
      fTID = pvec[i]->TrackId();
      // 0 out strings, else there may be cruft in here from prev evt.
      for (unsigned int s = 0; s < 35; ++s) {
        *(fTProcess + s) = 0;
        *(fTProcess + s) = 0;
        *(fTMaterial + s) = 0;
        *(fTMaterial + s) = 0;
        *(fTVolume + s) = 0;
        *(fTVolume + s) = 0;
        *(fTTVolume + s) = 0;
        *(fTTVolume + s) = 0;
      }

      for (unsigned int s = 0; s < pvec[i]->Process().length(); ++s)
        *(fTProcess + s) = pvec[i]->Process()[s];

      TVector3 dum = pvec[i]->Position().Vect();

      using geo::vect::toPoint;
      auto const material_at_pos = geom->MaterialName(toPoint(pvec[i]->Position().Vect()));
      for (unsigned int s = 0; s < material_at_pos.length(); ++s)
        *(fTMaterial + s) = material_at_pos[s];

      auto const volume_at_pos = geom->VolumeName(toPoint(pvec[i]->Position().Vect()));
      for (unsigned int s = 0; s < volume_at_pos.length(); ++s)
        *(fTVolume + s) = volume_at_pos[s];

      auto const volume_at_end = geom->VolumeName(toPoint(pvec[i]->EndPosition().Vect()));
      for (unsigned int s = 0; s < volume_at_end.length(); ++s)
        *(fTTVolume + s) = volume_at_end[s];

      fTEvt = evt.id().event();
      fTSub = evt.subRun();
      fTRun = evt.run();
      fTParentID = pvec[i]->Mother();
      fTStatus = pvec[i]->StatusCode();
      int daughter = 9999;
      fTNds = TMath::Min(pvec[i]->NumberDaughters(), fTNdsOriginal);
      for (int d = 0; d < fTNds; d++) {
        daughter = pvec[i]->Daughter(d);
        fTDID[d] = daughter;
        // zero it out.
        for (unsigned int s = 0; s < 35; ++s)
          *(fTDProcess[d] + s) = 0;

        for (unsigned int jj = i; jj < pvec.size(); ++jj) { // Don't look below i.

          if (fTDID[d] == pvec[jj]->TrackId()) {
            fTDPdg[d] = pvec[jj]->PdgCode(); // get the pointer,
            fTDWt[d] = pvec[jj]->Weight();

            for (unsigned int s = 0; s < pvec[jj]->Process().length(); ++s)
              *(fTDProcess[d] + s) = pvec[jj]->Process()[s];

            for (unsigned int kk = 0; kk < 4; ++kk) {
              fT4DOrigin[d * 4 + kk] = pvec[jj]->Position()[kk];
              fT4DMomentum[d * 4 + kk] = pvec[jj]->Momentum()[kk];
            }
            break;
          }
        }
      } //end loop over d

      for (unsigned int ii = 0; ii < 4; ++ii) {
        fT4Termination[ii] = 1e9;
        fT4Origin[ii] = pvec[i]->Position()[ii];
        if (ii != 3) fT4Termination[ii] = pvec[i]->EndPosition()[ii];
        if (ii == 4) fT4Termination[ii] = pvec[i]->Momentum()[ii]; // yes, odd
        fT4Momentum[ii] = pvec[i]->Momentum()[ii];
      }

      fTWeight = pvec[i]->Weight();
      fTree->Fill();

    } // end loop on particles in list
    if (numpi0gamma == 2 && pi0loc > 0) {
      mf::LogInfo("LArG4Ana") << pvec[pi0loc]->E();
      fPi0Momentum->Fill(pvec[pi0loc]->E());
    }

    return;
  }

} // namespace larg4

namespace larg4 {

  DEFINE_ART_MODULE(LArG4Ana)

} // namespace LArG4