Merge pull request #136 from JamesJieranShen/pr/ntuple-mc-charge

Add MC charge information to Ntuples
rat-pac · Jun 20, 2024 · 07d87fa · 07d87fa
2 parents 4ff882a + 2cde8e0
commit 07d87fa
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 0 deletions.
diff --git 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
@@ -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;