Plot full kientic energy spectrum

Armory/ClassAnasen.h

@@ -16,6 +16,8 @@
 #include "ClassSX3.h"
 #include "ClassPW.h"
 
+//to not include certain wires in the simulation, pass the anode and cathode IDs to the constructor, e.g. for anode wires 5-10, pass anodeID1 = 5, anodeID2 = 10, and for cathode wires 20-25, pass cathodeID1 = 20, cathodeID2 = 25. To include all wires, pass -1 for all IDs.
+
 class ANASEN{
 public:
   ANASEN();
@@ -108,7 +110,7 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1,
   geom = new TGeoManager("Detector", "ANASEN");
 
   //--- define some materials
-  TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0); //name, A, Z, density 
+  //TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0); //name, A, Z, density 
   TGeoMaterial *matHe = new TGeoMaterial("He", 4.0026, 2, 0.000861);
   TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
   TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085,14,2.33);

Armory/ClassPW.h

@@ -7,6 +7,9 @@
 #include <TVector3.h>
 #include <TRandom.h>
 
+std::vector<int> skipAnodes = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23}; 
+std::vector<int> skipCathodes = {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23};
+
 struct PWHitInfo
 {
   std::pair<short, short> nearestWire;   // anode, cathode
@@ -340,7 +343,8 @@ inline std::tuple<TVector3,double,double,double,double,double,double,double> PW:
 
 inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
 {
-
+//to skip certain wires, add their IDs to the skipAnodes and skipCathodes vectors above, then in this function, add a check to set the distance to a large number if the wire ID is in the skip list, e.g. if(std::find(skipAnodes.begin(), skipAnodes.end(), i) != skipAnodes.end()) { disA = 99999999; } to skip anode wire i, and if(std::find(skipCathodes.begin(), skipCathodes.end(), i) != skipCathodes.end()) { disC = 99999999; } to skip cathode wire i.
+//do this on the line before the line "if (phiS < phi && phi < phiL)" in the anode and cathode loops below, respectively, so that the wires will be marked as invalid and not used in track reconstruction
   hitInfo.Clear();
   double phi = direction.Phi();
 
@@ -355,8 +359,10 @@ inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
       phiL = phiL + TMath::TwoPi();
     if (phi < 0 && phiS > phiL)
       phiS = phiS - TMath::TwoPi();
-
+    if(std::find(skipAnodes.begin(), skipAnodes.end(), i) != skipAnodes.end()) { // check if the current anode wire ID is in the skipAnodes vector
-    if (phiS < phi && phi < phiL)
+      disA = 99999999;
+    }
+    if (phiS < phi && phi < phiL) // check if the track direction is within the angular range of the wire
     {
       disA = Distance(pos, pos + direction, An[i].first, An[i].second);
       if (disA < hitInfo.nearestDist.first)
@@ -370,6 +376,9 @@ inline void PW::FindWireID(TVector3 pos, TVector3 direction, bool verbose)
     phiS = Ca[i].second.Phi() - TMath::PiOver4();
     phiL = Ca[i].first.Phi() + TMath::PiOver4();
     // printf("C%2d: %f %f\n", i, phiS * TMath::RadToDeg(), phiL * TMath::RadToDeg());
+    if(std::find(skipCathodes.begin(), skipCathodes.end(), i) != skipCathodes.end()) { // check if the current cathode wire ID is in the skipCathodes vector
+      disC = 99999999;
+    }
     if (phi > 0 && phiS > phiL)
       phiL = phiL + TMath::TwoPi();
     if (phi < 0 && phiS > phiL)

Armory/aarootscript.C

@@ -24,20 +24,23 @@ void aarootscript(int argument = 0) {
     TTree *tree2 = (TTree*)f->Get("tree2");
 
     TTreeReader reader("tree");
-    TTreeReaderValue<double> branchValue(reader, "Tb");
+    TTreeReaderValue<double> Tb1(reader, "Tb"); // this line will read the "Tb" branch from the tree and save it as 
+    TTreeReaderValue<double> TB1(reader, "TB"); // this line will read the "TB" branch from the tree and save it as Tb1
+    //add Tb1 and TB1 together 
+
     double totalSum = 0;
     while (reader.Next()) {
-        totalSum += *branchValue;
+        totalSum += *Tb1;
     }
     std::cout << "Total Sum: " << totalSum << std::endl;
 
     TTreeReader reader2("tree2");
-    TTreeReaderValue<double> branchValue2(reader2, "Tb");
+    TTreeReaderValue<double> Tb2(reader2, "Tb");
     double totalSum2 = 0;
     int zeroCount = 0;
     while (reader2.Next()) {
-        totalSum2 += *branchValue2;
+        totalSum2 += *Tb2;
-        if (*branchValue2 == 0) {
+        if (*Tb2 == 0) {
             zeroCount++;
         }
     }
@@ -47,7 +50,9 @@ void aarootscript(int argument = 0) {
 
     std::cout << "Zero Count: " << zeroCount << std::endl;
 
-    std::cout << "Making histograms..." << std::endl;
+
+
+    //std::cout << "Making histograms..." << std::endl;
 
     gErrorIgnoreLevel = 2001;
     gROOT->ProcessLine(".x histcomp.C");
@@ -77,6 +82,7 @@ void aarootscript(int argument = 0) {
     TTreeReader r2(tree2);
 
     TTreeReaderValue<double> Tb_val(r1, "Tb");
+    TTreeReaderValue<double> TB_val(r2, "TB");
     TTreeReaderValue<double> dEb_val(r2, "dEb");
 
     std::vector<double> x; // Tb (tree1)
@@ -121,5 +127,8 @@ void aarootscript(int argument = 0) {
     delete c1;
     delete gr;*/
 
+
+
+
 }
 

Armory/anasenMS.cpp

@@ -11,9 +11,15 @@
 #include "TApplication.h" // ROOT app loop
 #include "ClassTransfer.h" // Reaction kinematics and MC event generation
 #include "ClassAnasen.h"   // ANASEN detector model classes (SX3, PW, etc.)
-#include "TCanvas.h"      // ROOT canvas for drawing
 #include <stdio.h>
 #include <stdlib.h>
+#include <set>
+#include "TLegend.h"
+#include "TH1D.h"
+#include "TObjArray.h"
+#include "TBranch.h"
+#include <iostream>
+#include <fstream>
 
 //======== Generate light particle based on reaction
 // calculate real and reconstructed tracks and Q-value uncertainty
@@ -24,6 +30,21 @@ TGraph* LoadELoss(const char* filename) {
   return g;
 }
 
+bool IsDeadAnode(int id){
+    static std::set<int> dead = {}; // add dead anode IDs here, 0-23
+    return dead.count(id);
+}
+
+bool IsDeadCathode(int id){
+    static std::set<int> dead = {}; // add dead cathode IDs here, 0-23
+    return dead.count(id);
+}
+
+bool IsDeadSX3(int id){
+    static std::set<int> dead = {}; // add dead SX3 IDs here, 0-8 for upstream, 9-17 for downstream, 18-23 for whole detector
+    return dead.count(id);
+}
+
 int main(int argc, char **argv){
 
   printf("=========================================\n");
@@ -33,33 +54,31 @@ int main(int argc, char **argv){
   // number of events can be overridden from command line
   int numEvent = 1000000;
   if( argc >= 2 ) numEvent = atoi(argv[1]);
-  //double density = (2.1525e-7) * 1000; // example for aluminum target, adjust as needed (400 torr is 0.0000861)
+
-  //char command[256];
-  //snprintf(command, sizeof(command), "python3 /home/jamesszalkie/.venv/EvXconverter.py %f", density);
-  //printf("Command: %s\n", command);
-  //system(command); // run the conversion script to generate energy loss tables with correct density
   // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm^2), density in mg/cm^3)
   TGraph* elossLight = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha)
   TGraph* elossHeavy = LoadELoss("../ELoss/E_vs_x_proton.dat");     // for heavy particle (proton)
-  TGraph* elossRecoil = LoadELoss("../ELoss/E_vs_x_recoil.dat");     // for recoil particle (if needed)
+  //TGraph* elossRecoil = LoadELoss("../ELoss/E_vs_x_recoil.dat");     // for recoil particle (if needed)
   TGraph *invgLight = new TGraph(elossLight->GetN(), elossLight->GetY(), elossLight->GetX());
   TGraph *invgHeavy = new TGraph(elossHeavy->GetN(), elossHeavy->GetY(), elossHeavy->GetX());
+
+  /*
+  //Plot energy loss tables (sanity check), vis will not work if this is ran without X11 display (e.g. on cluster), so comment out if running in headless mode
   auto c1 = new TCanvas("c1", "Graph Example", 800, 600);
   auto g = elossLight;
-  g->SetTitle("Energy Loss Table;cm;Kinetic Energy (MeV)");
+  g->SetTitle("Energy Loss Table (Light);cm;Kinetic Energy (MeV)");
   g->Draw("ALP");
   g->SetLineColor(kRed);
   //c1->SetLogy();
   //c1->SetLogx();
   c1->Print("eloss_light.png");
-  int ELossTotal = 0;
 
   auto c2 = new TCanvas("c2", "Graph Example", 800, 600);
   auto g2 = elossHeavy;
-  g2->SetTitle("Energy Loss Table;cm;Kinetic Energy (MeV)");
+  g2->SetTitle("Energy Loss Table (Heavy);cm;Kinetic Energy (MeV)");
   g2->Draw("ALP");
   g2->SetLineColor(kBlue);
-  c2->Print("eloss_heavy.png");
+  c2->Print("eloss_heavy.png");*/
 
   // Reaction setup: projectile + target configuration, energy, and product IDs
   TransferReaction transfer;
@@ -139,43 +158,45 @@ int main(int argc, char **argv){
   double thetab, phib, Tb;
   double thetaB, phiB, TB;
   double dEb;
-  double dEB; // placeholder for heavy particle energy loss, currently set equal to light particle loss for simplicity
+  double dEB;
-  tree->Branch("thetab", &thetab, "thetab/D");
+  std::array<double, 2> T;
-  tree->Branch("phib",     &phib, "phib/D");
+  tree->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
+  tree->Branch("phib",     &phib, "phib/D"); // azimuthal angle of light particle in lab frame
   tree->Branch("Tb",         &Tb, "Tb/D"); // kinetic energy of light particle at vertex (before energy loss)
   tree->Branch("thetaB", &thetaB, "thetaB/D");
   tree->Branch("phiB",     &phiB, "phiB/D");
   tree->Branch("TB",         &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss)
   tree->Branch("dEb",        &dEb, "dEb/D");
   tree->Branch("dEB",        &dEB, "dEB/D"); // placeholder for heavy particle energy loss, currently set equal to light particle loss for simplicity
+  tree->Branch("T",          &T, "T/D");   // placeholder for true Q-value, currently set to 0 for simplicity
 
   // excitation state identifiers
   int ExAID;
   double ExA;
-  tree->Branch("ExAID", &ExAID, "ExAID/I");
+  tree->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID
-  tree->Branch("ExA",     &ExA, "ExA/D");
+  tree->Branch("ExA",     &ExA, "ExA/D"); // projectile excitation energy in MeV
 
   int ExID;
   double Ex;
-  tree->Branch("ExID", &ExID, "ExID/I");
+  tree->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID
-  tree->Branch("Ex",     &Ex, "Ex/D");
+  tree->Branch("Ex",     &Ex, "Ex/D"); // target excitation energy in MeV
 
   // true vertex position in target volume
   double vertexX, vertexY, vertexZ;
-  tree->Branch("vX",     &vertexX, "VertexX/D");
+  tree->Branch("vX",     &vertexX, "VertexX/D"); // true vertex X position in mm
-  tree->Branch("vY",     &vertexY, "VertexY/D");
+  tree->Branch("vY",     &vertexY, "VertexY/D"); // true vertex Y position in mm
-  tree->Branch("vZ",     &vertexZ, "VertexZ/D");
+  tree->Branch("vZ",     &vertexZ, "VertexZ/D"); // true vertex Z position in mm
 
   // reconstructed SX3 hit position
   double sx3X, sx3Y, sx3Z;
-  tree->Branch("sx3X",     &sx3X, "sx3X/D");
+  tree->Branch("sx3X",     &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing)
-  tree->Branch("sx3Y",     &sx3Y, "sx3Y/D");
+  tree->Branch("sx3Y",     &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
-  tree->Branch("sx3Z",     &sx3Z, "sx3Z/D");
+  tree->Branch("sx3Z",     &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
 
   // PW nearest and next nearest wires
   int anodeID[2], cathodeID[2];
-  tree->Branch("aID", anodeID, "anodeID/I");
+  tree->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire
-  tree->Branch("cID", cathodeID, "cathodeID/I");
+  tree->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire
 
   // distances to nearest wires
   double anodeDist[2], cathodeDist[2];
@@ -213,6 +234,7 @@ int main(int argc, char **argv){
   clock.Reset();
   clock.Start("timer");
   shown = false;
+  int ELossTotal = 0;
 
   //================================= Calculate event loop
   for( int i = 0; i < numEvent ; i++){
@@ -243,6 +265,8 @@ int main(int argc, char **argv){
 
     Tb = Pb.E() - Pb.M();
     TB = PB.E() - PB.M();
+    T[0] = Tb;
+    T[1] = TB;
 
     phib = Pb.Phi() * TMath::RadToDeg();
     phiB = PB.Phi() * TMath::RadToDeg();
@@ -265,19 +289,25 @@ int main(int argc, char **argv){
 
     PWHitInfo hitInfo = pw->GetHitInfo();
 
-    // store nearest/next closest wire data
+    anodeID[0]   = hitInfo.nearestWire.first; // nearest anode wire ID
-    anodeID[0]   = hitInfo.nearestWire.first;
+    cathodeID[0] = hitInfo.nearestWire.second; // nearest cathode wire ID
-    cathodeID[0] = hitInfo.nearestWire.second;
+    anodeID[1]   = hitInfo.nextNearestWire.first; // next nearest anode wire ID
-    anodeID[1]   = hitInfo.nextNearestWire.first;
+    cathodeID[1] = hitInfo.nextNearestWire.second; // next nearest cathode wire ID
-    cathodeID[1] = hitInfo.nextNearestWire.second;
 
-    anodeDist[0]   = hitInfo.nearestDist.first;
+    anodeDist[1]   = hitInfo.nextNearestDist.first; // distance to next nearest anode wire
-    cathodeDist[0] = hitInfo.nearestDist.second;
+    cathodeDist[1] = hitInfo.nextNearestDist.second; // distance to next nearest cathode wire
-    anodeDist[1]   = hitInfo.nextNearestDist.first;
+
-    cathodeDist[1] = hitInfo.nextNearestDist.second;
+    if(IsDeadAnode(anodeID[0])) continue;
+    if(IsDeadCathode(cathodeID[0])) continue;
 
     // SX3 hit channel info and depth fraction
     sx3ID = sx3->GetID();
+
+    if(IsDeadSX3(sx3ID)) continue;
+
+    anodeDist[0]   = hitInfo.nearestDist.first; // distance to nearest anode wire
+    cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire
+
     if( sx3ID >= 0 ){
       sx3Up    = sx3->GetChUp();
       sx3Dn    = sx3->GetChDn();
@@ -316,34 +346,40 @@ int main(int argc, char **argv){
 
       //Energy loss
       double dl = (hitPos - vertex).Mag(); // path length in units of cm
-      if (numEvent <= 10){
+
-        printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
+      if (numEvent <= 100){
+        //printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
+        //printf("Total T before loss: %f MeV\n", T);
       }
+
       double tb_temp = Tb;
       double tB_temp = TB;
+
       dEb = tb_temp - Tb; // total energy loss
       dEB = tB_temp - TB; // total energy loss for heavy particle, currently set equal to light particle loss for simplicity
+
       double x0light = invgLight->Eval(Tb);
       double x0heavy = invgHeavy->Eval(TB);
+
       x0light = x0light + dl;
       x0heavy = x0heavy + dl;
+
       Tb = elossLight->Eval(x0light); // kinetic energy corresponding to range at hit position
       TB = elossHeavy->Eval(x0heavy); // kinetic energy for heavy particle, currently set equal to light particle loss for simplicity
+
       dEb = tb_temp - Tb; // total energy loss
       dEB = tB_temp - TB; // total energy loss for heavy particle, currently set equal to light particle loss for simplicity
       
       // fill tree2 with energy loss adjusted data
-      if (Tb != 0) {
+      T[0] = Tb;
-        tree2->Fill();
+      T[1] = TB;
-      }
+      tree2->Fill();
+      
       if (numEvent <= 10){
-        printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb);
+        //printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb);
       } //to give in scientific notation, use %e instead of %f in the printf format string. For example: printf("Event %d: Tb after energy loss = %e MeV, energy loss = %e MeV\n", i, Tb, tb_temp - Tb);
       ELossTotal += (tb_temp - Tb) + (tB_temp - TB);
 
-
-
-
     }else{
       // no valid SX3 hit: mark clearly invalid
       sx3Up = -1;
@@ -400,7 +436,7 @@ int main(int argc, char **argv){
   printf("Total energy loss across all events: %e MeV\n", (double)ELossTotal);
   printf("Average energy loss across events: %e MeV\n", (double)ELossTotal / count);
 
-  if(enableVis){
+  if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
     printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
     
     // Build full geometry with all wires
@@ -453,6 +489,7 @@ int main(int argc, char **argv){
     }
     
     if(app){
+      printf("Entering ROOT event loop\n");
       app->Run();
     }
   }

Armory/histcomp.C

@@ -33,11 +33,24 @@ void histcomp() {
         TString h2name = "h2_" + name;
 
         // Temporary draw to get range
-        tree1->Draw(name, "", "goff");
+        double min, max;
-        double min = fmin(tree1->GetMinimum(name), tree2->GetMinimum(name));
-        double max = fmax(tree1->GetMaximum(name), tree2->GetMaximum(name));
 
-        if (min == max) continue; // skip constant branches
+        if(name == "T"){
+            //Get minimum value of T[0] and use as min  
+            min = tree2->GetMinimum("Tb");
+            max = tree1->GetMaximum("TB");
+        }else{
+
+            tree1->Draw(name, "", "goff");
+
+            min = fmin(tree1->GetMinimum(name),
+                    tree2->GetMinimum(name));
+
+            max = fmax(tree1->GetMaximum(name),
+                    tree2->GetMaximum(name));
+        }
+
+        //if (min == max) continue; // skip constant branches
 
         // Expand range slightly
         double margin = 0.1 * (max - min);
@@ -48,8 +61,21 @@ void histcomp() {
         TH1D *h2 = new TH1D(h2name, name, 100, min, max);
 
         // Fill histograms
-        tree1->Draw(name + ">>" + h1name, "", "goff");
+        if(name == "T"){
-        tree2->Draw(name + ">>" + h2name, "", "goff");
+
+            // Fill both array elements into same histogram
+            tree1->Draw("Tb>>+" + h1name, "", "goff");
+            tree1->Draw("TB>>+" + h1name, "", "goff");
+
+            tree2->Draw("Tb>>+" + h2name, "", "goff");
+            tree2->Draw("TB>>+" + h2name, "", "goff");
+
+        }else{
+
+            tree1->Draw(name + ">>" + h1name, "", "goff");
+            tree2->Draw(name + ">>" + h2name, "", "goff");
+
+        }
 
         // Style
         h1->SetLineColor(kRed);
@@ -63,7 +89,7 @@ void histcomp() {
         //if (h2->GetEntries() > 0) h2->Scale(1.0 / h2->GetEntries());
 
         // Canvas
-        TCanvas *c = new TCanvas("c", name, 900, 600);
+        TCanvas *c = new TCanvas("c", name, 900, 600); //arguments are (name, title, width, height)
 
         c->SetRightMargin(0.18);
         c->Modified(); 
@@ -91,23 +117,29 @@ void histcomp() {
         leg->AddEntry(h1, "tree1", "l");
         leg->AddEntry(h2, "tree2", "l");
         leg->Draw();
-
+        //to plot both as one histogram in root, can use tree2->Draw("T(0)"); for light particle and tree2->Draw("T(1)") for heavy particle
         // Save plot (overwrite each run)
         TString filename = "plots/" + name + ".png";
         c->SaveAs(filename);
-        //c->SetLogx(0);
-        //c->SetLogy(0);
         // Optional: save log plots as well
-        c->SetLogx(1);
+
-        c->SetLogy(1);
+        if (false) { // set to True to also save log plots
-        h1->SetTitle(name + " (log);"+name+";Counts");
+            c->SetLogy(1);
-        c->SaveAs("plots/" + name + "_log.png");
+            h1->SetTitle(name + " (log);"+name+";Counts");
+            c->SaveAs("plots/" + name + "_logy.png");
+
+            c->SetLogy(0);
+            c->SetLogx(1);
+            h1->SetTitle(name + " (log);"+name+";Counts");
+            c->SaveAs("plots/" + name + "_logx.png");
+
+            // Clean up
+            delete c;
+            delete h1;
+            delete h2;
+        }
 
 
-        // Clean up
-        delete c;
-        delete h1;
-        delete h2;
 
     }