From 2cde8e05c0b2bb5548beaceea2d418c65ca844f0 Mon Sep 17 00:00:00 2001
From: James Shen <jamesjieranshen@gmail.com>
Date: Thu, 30 May 2024 14:31:34 -0400
Subject: [PATCH] Add MC charge information to Ntuple

---
 src/io/include/RAT/OutNtupleProc.hh | 2 ++
 src/io/src/OutNtupleProc.cc         | 6 ++++++
 2 files changed, 8 insertions(+)

diff --git a/src/io/include/RAT/OutNtupleProc.hh b/src/io/include/RAT/OutNtupleProc.hh
index 33499982..72bb6da3 100644
--- a/src/io/include/RAT/OutNtupleProc.hh
+++ b/src/io/include/RAT/OutNtupleProc.hh
@@ -95,6 +95,7 @@ class OutNtupleProc : public Processor {
   int mcpecount;
   std::vector<int> mcpmtid;
   std::vector<int> mcpmtnpe;
+  std::vector<double> mcpmtcharge;
   // MCPE
   std::vector<double> mcpehittime;
   std::vector<double> mcpefrontendtime;
@@ -103,6 +104,7 @@ class OutNtupleProc : public Processor {
   std::vector<double> mcpex;
   std::vector<double> mcpey;
   std::vector<double> mcpez;
+  std::vector<double> mcpecharge;
   // MCParticles
   int mcpcount;
   std::vector<Int_t> pdgcodes;
diff --git a/src/io/src/OutNtupleProc.cc b/src/io/src/OutNtupleProc.cc
index 0710ee37..fd08785e 100644
--- a/src/io/src/OutNtupleProc.cc
+++ b/src/io/src/OutNtupleProc.cc
@@ -133,6 +133,7 @@ bool OutNtupleProc::OpenFile(std::string filename) {
     // Save full MC PMT hit information
     outputTree->Branch("mcPMTID", &mcpmtid);
     outputTree->Branch("mcPMTNPE", &mcpmtnpe);
+    outputTree->Branch("mcPMTCharge", &mcpmtcharge);
 
     outputTree->Branch("mcPEHitTime", &mcpehittime);
     outputTree->Branch("mcPEFrontEndTime", &mcpefrontendtime);
@@ -143,6 +144,7 @@ bool OutNtupleProc::OpenFile(std::string filename) {
     outputTree->Branch("mcPEx", &mcpex);
     outputTree->Branch("mcPEy", &mcpey);
     outputTree->Branch("mcPEz", &mcpez);
+    outputTree->Branch("mcPECharge", &mcpecharge);
   }
   if (options.tracking) {
     // Save particle tracking information
@@ -286,6 +288,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) {
   // MCPMT information
   mcpmtid.clear();
   mcpmtnpe.clear();
+  mcpmtcharge.clear();
 
   // MCPE information
   mcpehittime.clear();
@@ -295,6 +298,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) {
   mcpex.clear();
   mcpey.clear();
   mcpez.clear();
+  mcpecharge.clear();
 
   mcnhits = mc->GetMCPMTCount();
   mcpecount = mc->GetNumPE();
@@ -303,6 +307,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) {
       DS::MCPMT *mcpmt = mc->GetMCPMT(ipmt);
       mcpmtid.push_back(mcpmt->GetID());
       mcpmtnpe.push_back(mcpmt->GetMCPhotonCount());
+      mcpmtcharge.push_back(mcpmt->GetCharge());
       TVector3 position = pmtinfo->GetPosition(mcpmt->GetID());
       for (int ipe = 0; ipe < mcpmt->GetMCPhotonCount(); ipe++) {
         RAT::DS::MCPhoton *mcph = mcpmt->GetMCPhoton(ipe);
@@ -312,6 +317,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) {
         mcpex.push_back(position.X());
         mcpey.push_back(position.Y());
         mcpez.push_back(position.Z());
+        mcpecharge.push_back(mcph->GetCharge());
         if (mcph->IsDarkHit()) {
           mcpeprocess.push_back(noise);
           continue;