Add nested tube hit info to ntuple output

macros/vis_nestedTubeArrays.mac

@@ -44,7 +44,9 @@
 /rat/procset include_tracking 0
 /rat/procset include_mcparticles 1
 /rat/procset include_pmthits 1
+/rat/procset include_nestedtubehits 1
 /rat/procset include_untriggered_events 1
+#/rat/proc outroot
 /rat/physics/setOpWLS g4
 ##### GENERATORS #################
 /generator/add combo pbomb:point:poisson

ratdb/Validation/NestedTubeArray.geo

@@ -166,5 +166,6 @@ outer_material: "Fpolyethylene",
 inner_material: "PMMA",
 core_material: "WLSExample",
 #drawstyle: "solid",
-color: [0.0,0.8,0.0,0.2]
+color: [0.0,0.8,0.0,0.2],
+sensitive_detector: "/mydet/fibers"
 }

src/ds/include/RAT/DS/LinkDef.hh

@@ -68,6 +68,8 @@
 #pragma link C++ class vector < RAT::DS::MCTrackStep>;
 #pragma link C++ class vector < RAT::DS::MCPMT>;
 #pragma link C++ class vector < RAT::DS::MCPhoton>;
+#pragma link C++ class vector < RAT::DS::MCNestedTube>;
+#pragma link C++ class vector < RAT::DS::MCNestedTubeHit>;
 #pragma link C++ class vector < RAT::DS::Calib>;
 #pragma link C++ class vector < RAT::DS::EV>;
 #pragma link C++ class vector < RAT::DS::PMT>;
@@ -89,6 +91,8 @@
 #pragma link C++ class vector < RAT::DS::MCTrack *>;
 #pragma link C++ class vector < RAT::DS::MCPMT *>;
 #pragma link C++ class vector < RAT::DS::MCPhoton *>;
+#pragma link C++ class vector < RAT::DS::MCNestedTube *>;
+#pragma link C++ class vector < RAT::DS::MCNestedTubeHit *>;
 #pragma link C++ class vector < RAT::DS::MCTrackStep *>;
 
 #endif

src/io/include/RAT/OutNtupleProc.hh

@@ -51,6 +51,7 @@ class OutNtupleProc : public Processor {
     bool pmthits;
     bool untriggered;
     bool mchits;
+    bool nthits;
   };
   NtupleOptions options;
 
@@ -74,6 +75,13 @@ class OutNtupleProc : public Processor {
   std::vector<double> pmtU;
   std::vector<double> pmtV;
   std::vector<double> pmtW;
+  std::vector<int> ntId;
+  std::vector<double> ntX;
+  std::vector<double> ntY;
+  std::vector<double> ntZ;
+  std::vector<double> ntU;
+  std::vector<double> ntV;
+  std::vector<double> ntW;
   // Data Branches
   Int_t mcpdg;
   double mcx, mcy, mcz;
@@ -95,6 +103,14 @@ class OutNtupleProc : public Processor {
   int mcpecount;
   std::vector<int> mcpmtid;
   std::vector<int> mcpmtnpe;
+  // MCNestedTube
+  int mcnNTs;
+  int mcnNThits;
+  std::vector<int> mcNTid;
+  std::vector<double> mcNThittime;
+  std::vector<double> mcNThitx;
+  std::vector<double> mcNThity;
+  std::vector<double> mcNThitz;
   // MCPE
   std::vector<double> mcpehittime;
   std::vector<double> mcpefrontendtime;
@@ -125,11 +141,11 @@ class OutNtupleProc : public Processor {
   std::vector<double> hitPMTDigitizedTime;
   std::vector<double> hitPMTDigitizedCharge;
   std::vector<int> hitPMTNCrossings;
-  // Tracking
   std::map<std::string, int> processCodeMap;
   std::vector<int> processCodeIndex;
   std::vector<std::string> processName;
 
+  // Tracking
   std::vector<int> trackPDG;
   std::vector<std::vector<double>> trackPosX;
   std::vector<std::vector<double>> trackPosY;

src/io/src/OutNtupleProc.cc

@@ -48,12 +48,14 @@ OutNtupleProc::OutNtupleProc() : Processor("outntuple") {
     options.pmthits = table->GetZ("include_pmthits");
     options.untriggered = table->GetZ("include_untriggered_events");
     options.mchits = table->GetZ("include_mchits");
+    options.nthits = table->GetZ("include_nestedtubehits");
   } catch (DBNotFoundError &e) {
     options.tracking = false;
     options.mcparticles = false;
     options.pmthits = true;
     options.untriggered = false;
     options.mchits = true;
+    options.nthits = false;
   }
 }
 
@@ -129,6 +131,22 @@ bool OutNtupleProc::OpenFile(std::string filename) {
     outputTree->Branch("hitPMTDigitizedCharge", &hitPMTDigitizedCharge);
     outputTree->Branch("hitPMTNCrossings", &hitPMTNCrossings);
   }
+  if (options.nthits) {
+    outputTree->Branch("mcnNTs", &mcnNTs);
+    outputTree->Branch("mcnNThits", &mcnNThits);
+    outputTree->Branch("mcNTid", &mcNTid);
+    outputTree->Branch("mcNThittime", &mcNThittime);
+    outputTree->Branch("mcNThitx", &mcNThitx);
+    outputTree->Branch("mcNThity", &mcNThity);
+    outputTree->Branch("mcNThitz", &mcNThitz);
+    metaTree->Branch("ntId", &ntId);
+    metaTree->Branch("ntX", &ntX);
+    metaTree->Branch("ntY", &ntY);
+    metaTree->Branch("ntZ", &ntZ);
+    metaTree->Branch("ntU", &ntU);
+    metaTree->Branch("ntV", &ntV);
+    metaTree->Branch("ntW", &ntW);
+  }
   if (options.mchits) {
     // Save full MC PMT hit information
     outputTree->Branch("mcPMTID", &mcpmtid);
@@ -171,6 +189,7 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) {
   }
   runBranch = DS::RunStore::GetRun(ds);
   DS::PMTInfo *pmtinfo = runBranch->GetPMTInfo();
+  DS::NestedTubeInfo *ntinfo = runBranch->GetNestedTubeInfo();
   ULong64_t stonano = 1000000000;
   dsentries++;
   // Clear the previous vectors
@@ -329,6 +348,23 @@ Processor::Result OutNtupleProc::DSEvent(DS::Root *ds) {
       }
     }
   }
+  if (options.nthits) {
+    mcnNTs = mc->GetMCNestedTubeCount();
+    mcnNThits = 0;
+    for (int iNT = 0; iNT < mc->GetMCNestedTubeCount(); iNT++) {
+      DS::MCNestedTube *mcnt = mc->GetMCNestedTube(iNT);
+      mcnNThits += mcnt->GetMCNestedTubeHitCount();
+      mcNTid.push_back(mcnt->GetID());
+      G4ThreeVector position = ntinfo->GetPosition(mcnt->GetID());
+      for (int ih = 0; ih < mcnt->GetMCNestedTubeHitCount(); ih++) {
+        RAT::DS::MCNestedTubeHit *mcph = mcnt->GetMCNestedTubeHit(ih);
+        mcNThittime.push_back(mcph->GetHitTime());
+        mcNThitx.push_back(position.x());
+        mcNThity.push_back(position.y());
+        mcNThitz.push_back(position.z());
+      }
+    }
+  }
 
   // EV Branches
   for (subev = 0; subev < ds->GetEVCount(); subev++) {
@@ -463,6 +499,20 @@ OutNtupleProc::~OutNtupleProc() {
       pmtV.push_back(direction.Y());
       pmtW.push_back(direction.Z());
     }
+    if (options.nthits) {
+      DS::NestedTubeInfo *ntinfo = runBranch->GetNestedTubeInfo();
+      for (int id = 0; id < ntinfo->GetNestedTubeCount(); id++) {
+        G4ThreeVector position = ntinfo->GetPosition(id);
+        G4ThreeVector direction = ntinfo->GetDirection(id);
+        ntId.push_back(id);
+        ntX.push_back(position.x());
+        ntY.push_back(position.y());
+        ntZ.push_back(position.z());
+        ntU.push_back(direction.x());
+        ntV.push_back(direction.y());
+        ntW.push_back(direction.z());
+      }
+    }
     runId = runBranch->GetID();
     runType = runBranch->GetType();
     // Converting to unix time
@@ -526,6 +576,9 @@ void OutNtupleProc::SetI(std::string param, int value) {
   if (param == "include_pmthits") {
     options.pmthits = value ? true : false;
   }
+  if (param == "include_nestedtubehits") {
+    options.nthits = value ? true : false;
+  }
   if (param == "include_untriggered_events") {
     options.untriggered = value ? true : false;
   }