From 851d044f25d006031cb57815a6ecf91ec3c8be63 Mon Sep 17 00:00:00 2001 From: James Szalkie Date: Wed, 13 May 2026 13:48:41 -0400 Subject: [PATCH] Plot full kientic energy spectrum --- Armory/ClassAnasen.h | 4 +- Armory/ClassPW.h | 15 ++++- Armory/aarootscript.C | 21 +++++-- Armory/anasenMS.cpp | 127 +++++++++++++++++++++++++++--------------- Armory/histcomp.C | 68 ++++++++++++++++------ 5 files changed, 162 insertions(+), 73 deletions(-) diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h index 1502d0d..549392e 100644 --- a/Armory/ClassAnasen.h +++ b/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); diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index 0fd64bc..bb33b80 100755 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -7,6 +7,9 @@ #include #include +std::vector 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 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 nearestWire; // anode, cathode @@ -340,7 +343,8 @@ inline std::tuple 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 (phiS < phi && phi < phiL) + if(std::find(skipAnodes.begin(), skipAnodes.end(), i) != skipAnodes.end()) { // check if the current anode wire ID is in the skipAnodes vector + 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) diff --git a/Armory/aarootscript.C b/Armory/aarootscript.C index 5f92212..3aa9a22 100644 --- a/Armory/aarootscript.C +++ b/Armory/aarootscript.C @@ -24,20 +24,23 @@ void aarootscript(int argument = 0) { TTree *tree2 = (TTree*)f->Get("tree2"); TTreeReader reader("tree"); - TTreeReaderValue branchValue(reader, "Tb"); + TTreeReaderValue Tb1(reader, "Tb"); // this line will read the "Tb" branch from the tree and save it as + TTreeReaderValue 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 branchValue2(reader2, "Tb"); + TTreeReaderValue Tb2(reader2, "Tb"); double totalSum2 = 0; int zeroCount = 0; while (reader2.Next()) { - totalSum2 += *branchValue2; - if (*branchValue2 == 0) { + totalSum2 += *Tb2; + 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 Tb_val(r1, "Tb"); + TTreeReaderValue TB_val(r2, "TB"); TTreeReaderValue dEb_val(r2, "dEb"); std::vector x; // Tb (tree1) @@ -121,5 +127,8 @@ void aarootscript(int argument = 0) { delete c1; delete gr;*/ + + + } diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 000cd79..0623a8f 100644 --- a/Armory/anasenMS.cpp +++ b/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 #include +#include +#include "TLegend.h" +#include "TH1D.h" +#include "TObjArray.h" +#include "TBranch.h" +#include +#include //======== 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 dead = {}; // add dead anode IDs here, 0-23 + return dead.count(id); +} + +bool IsDeadCathode(int id){ + static std::set dead = {}; // add dead cathode IDs here, 0-23 + return dead.count(id); +} + +bool IsDeadSX3(int id){ + static std::set 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 - tree->Branch("thetab", &thetab, "thetab/D"); - tree->Branch("phib", &phib, "phib/D"); + double dEB; + std::array T; + 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("ExA", &ExA, "ExA/D"); + tree->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID + tree->Branch("ExA", &ExA, "ExA/D"); // projectile excitation energy in MeV int ExID; double Ex; - tree->Branch("ExID", &ExID, "ExID/I"); - tree->Branch("Ex", &Ex, "Ex/D"); + tree->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID + 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("vY", &vertexY, "VertexY/D"); - tree->Branch("vZ", &vertexZ, "VertexZ/D"); + tree->Branch("vX", &vertexX, "VertexX/D"); // true vertex X position in mm + tree->Branch("vY", &vertexY, "VertexY/D"); // true vertex Y position in mm + 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("sx3Y", &sx3Y, "sx3Y/D"); - tree->Branch("sx3Z", &sx3Z, "sx3Z/D"); + tree->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing) + tree->Branch("sx3Y", &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing) + 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("cID", cathodeID, "cathodeID/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"); // 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; - cathodeID[0] = hitInfo.nearestWire.second; - anodeID[1] = hitInfo.nextNearestWire.first; - cathodeID[1] = hitInfo.nextNearestWire.second; + anodeID[0] = hitInfo.nearestWire.first; // nearest anode wire ID + cathodeID[0] = hitInfo.nearestWire.second; // nearest cathode wire ID + anodeID[1] = hitInfo.nextNearestWire.first; // next nearest anode wire ID + cathodeID[1] = hitInfo.nextNearestWire.second; // next nearest cathode wire ID - anodeDist[0] = hitInfo.nearestDist.first; - cathodeDist[0] = hitInfo.nearestDist.second; - anodeDist[1] = hitInfo.nextNearestDist.first; - cathodeDist[1] = hitInfo.nextNearestDist.second; + anodeDist[1] = hitInfo.nextNearestDist.first; // distance to next nearest anode wire + cathodeDist[1] = hitInfo.nextNearestDist.second; // distance to next nearest cathode wire + + 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) { - tree2->Fill(); - } + T[0] = Tb; + 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(); } } diff --git a/Armory/histcomp.C b/Armory/histcomp.C index feab86f..bca0e2f 100644 --- a/Armory/histcomp.C +++ b/Armory/histcomp.C @@ -33,11 +33,24 @@ void histcomp() { TString h2name = "h2_" + name; // Temporary draw to get range - tree1->Draw(name, "", "goff"); - double min = fmin(tree1->GetMinimum(name), tree2->GetMinimum(name)); - double max = fmax(tree1->GetMaximum(name), tree2->GetMaximum(name)); + double min, max; - 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"); - tree2->Draw(name + ">>" + h2name, "", "goff"); + if(name == "T"){ + + // 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); - h1->SetTitle(name + " (log);"+name+";Counts"); - c->SaveAs("plots/" + name + "_log.png"); + + if (false) { // set to True to also save log plots + c->SetLogy(1); + 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; }