From 8fae28b96a0a20935bfe095cd78e1c8419d0a25a Mon Sep 17 00:00:00 2001
From: fox <fox@foolishfox.cn>
Date: Mon, 25 Nov 2024 13:50:15 +0800
Subject: [PATCH] Fix: cutEnable, tof histogram draw, fill event hist only if
 there good tracks; Feat: draw option

---
 include/BluetAnalyzer.hh |  2 +-
 main.cc                  | 24 ++++++++++++----------
 src/BluetAnalyzer.cc     | 43 ++++++++++++++++++++++------------------
 3 files changed, 39 insertions(+), 30 deletions(-)

diff --git a/include/BluetAnalyzer.hh b/include/BluetAnalyzer.hh
index 9e818fd..6586f78 100644
--- a/include/BluetAnalyzer.hh
+++ b/include/BluetAnalyzer.hh
@@ -143,7 +143,7 @@ public:
   void readTreeData(vector<string> files);
   void FillHistograms();
   void postOperations();
-  void saveHistograms(TString);
+  void saveHistograms(TString filename);
   void drawHistograms();
 
 private:
diff --git a/main.cc b/main.cc
index a829b3a..4e23f0c 100644
--- a/main.cc
+++ b/main.cc
@@ -13,22 +13,25 @@ int main(int argc, char **argv) {
   string cfgfile;
   string datafile;
   string outputfile = "bluetAnalyzer.root";
+  bool draw = false;
   vector<string> ranges;
   vector<string> runfiles;
   Long_t rl, rr;
 
   auto cfg = value("config file", cfgfile);
-  auto data = (((required("-f", "--file").doc("data file path") &
-                 value("data file", datafile)) |
-                (required("-r", "--range").doc("data file number range") &
-                 values("number1 number2", ranges))),
-               option("-o", "--output")
-                       .doc("output file path, default: bluetAnalyzer.root") &
-                   value("output file", outputfile));
+  auto data = ((required("-f", "--file").doc("data file path") &
+                value("data file", datafile)) |
+               (required("-r", "--range").doc("data file number range") &
+                values("number1 number2", ranges)));
+  auto output = (option("-o", "--output")
+                         .doc("output file path, default: bluetAnalyzer.root") &
+                     value("output file", outputfile),
+                 option("-d", "--draw").doc("draw histograms").set(draw, true));
 
+  auto cmd = (cfg & data & output);
   auto fmt = doc_formatting{}.first_column(4);
-  if (!parse(argc, const_cast<char **>(argv), cfg & data)) {
-    cerr << make_man_page(cfg & data, "BluetAnalyzer", fmt)
+  if (!parse(argc, const_cast<char **>(argv), cmd)) {
+    cerr << make_man_page(cmd, "BluetAnalyzer", fmt)
                 .prepend_section("DESCRIPTION",
                                  "    Program for MTPC Data Analysis.")
          << endl;
@@ -60,7 +63,8 @@ int main(int argc, char **argv) {
   bluet->FillHistograms();
   bluet->postOperations();
   bluet->saveHistograms(outputfile);
-  bluet->drawHistograms();
+  if (draw)
+    bluet->drawHistograms();
 
   return 0;
 }
diff --git a/src/BluetAnalyzer.cc b/src/BluetAnalyzer.cc
index 8a2aac4..09405cc 100644
--- a/src/BluetAnalyzer.cc
+++ b/src/BluetAnalyzer.cc
@@ -144,7 +144,7 @@ void BluetAnalyzer::readCutParameters(string cfgfile) {
   dv0 = bluet::XMLToQuantity(gas.child("DriftVelocity")) * 1e4;
 
   pugi::xml_node cut = root.child("CutParameters");
-  cutEnable = bluet::XMLToQuantity(cut.child("CutEnable"));
+  cutEnable = bluet::XMLToQuantity(cut.child("cutEnable"));
   EcutMin = bluet::XMLToQuantity(cut.child("ECutMin"));
   EcutMax = bluet::XMLToQuantity(cut.child("ECutMax"));
   tofCutMin = bluet::XMLToQuantity(cut.child("TOFCutMin"));
@@ -165,7 +165,7 @@ void BluetAnalyzer::readCutParameters(string cfgfile) {
   lengthCutMax = bluet::XMLToQuantity(cut.child("LengthCutMax"));
   cosCutMin = bluet::XMLToQuantity(cut.child("cosCutMin"));
 
-  CenterPos[0] = -3.03631;
+  CenterPos[0] = 0;
   CenterPos[1] = 0.634384;
   sliceW = 2.5; // [mm]
   rCut = bluet::XMLToQuantity(cut.child("radiusCutEnable"));
@@ -474,26 +474,13 @@ void BluetAnalyzer::FillHistograms() {
     double ts = fEvent->TimeStamp;
     double tof = fEvent->TOF0;
     double tcath = fEvent->CathodeTime;
-    hts->Fill(ts, tof);
-
-    // convert time of flight to neutron energy using relativistic formula
-    double tof1 = fEvent->TOF0 - kickdT;
-    double vn = flighL / tof1 * 1e9;
-    double gamma = 1. / Sqrt(1. - vn * vn / (vc * vc));
-    double en = (gamma - 1) * mn * 1e6; // [MeV]->[eV]
-    hec->Fill(en);
-    hte->Fill(tof1, en);
-
-    for (size_t ihit = 0; ihit < fEvent->Hits.size(); ++ihit) {
-      BluetHit *hit = &(fEvent->Hits[ihit]);
-      hta->Fill(hit->HitTau, hit->HitAmp);
-      htt->Fill(hit->HitT * sampleT / 1e3, hit->HitTau);
-    }
 
+    bool goodTrack = false;
     for (size_t itrk = 0; itrk < fEvent->Tracks.size(); ++itrk) {
       if (!isGoodTrack(&(fEvent->Tracks[itrk])))
         continue;
 
+      goodTrack = true;
       double sumE = fEvent->Tracks[itrk].TrSumEdep;
       double length = fEvent->Tracks[itrk].TrRawLength;
       double traTmax = fEvent->Tracks[itrk].Tmax;
@@ -540,6 +527,24 @@ void BluetAnalyzer::FillHistograms() {
         htofE->Fill(tof, sumE);
       }
     }
+
+    if (goodTrack) {
+      hts->Fill(ts, tof);
+
+      // convert time of flight to neutron energy using relativistic formula
+      double tof1 = fEvent->TOF0 - kickdT;
+      double vn = flighL / tof1 * 1e9;
+      double gamma = 1. / Sqrt(1. - vn * vn / (vc * vc));
+      double en = (gamma - 1) * mn * 1e6; // [MeV]->[eV]
+      hec->Fill(en);
+      hte->Fill(tof1, en);
+
+      for (size_t ihit = 0; ihit < fEvent->Hits.size(); ++ihit) {
+        BluetHit *hit = &(fEvent->Hits[ihit]);
+        hta->Fill(hit->HitTau, hit->HitAmp);
+        htt->Fill(hit->HitT * sampleT / 1e3, hit->HitTau);
+      }
+    }
   }
   printf("\n===> Finish reading data ...\n");
 }
@@ -701,9 +706,9 @@ void BluetAnalyzer::drawTOF() {
   pad->cd();
   gPad->SetLogx();
   gPad->SetLogy();
-  htofr0->Draw("HIST");
+  htof->Draw("HIST");
+  htofr0->Draw("same HIST");
   htofr1->Draw("same HIST");
-  htof->Draw("same HIST");
   TLegend *leg = new TLegend(0.45, 0.75, 0.8, 0.9);
   leg->AddEntry(htofr0, Form("r < %.1f mm", rCutMin));
   leg->AddEntry(htofr1, Form("%.1f < r < %.1f mm", rCutMin, rCutMax));