diff --git a/.DS_Store b/.DS_Store index 0e0ba8f..e8617b1 100644 Binary files a/.DS_Store and b/.DS_Store differ diff --git a/Armory/.DS_Store b/Armory/.DS_Store index 0920f30..43cac33 100644 Binary files a/Armory/.DS_Store and b/Armory/.DS_Store differ diff --git a/Armory/AutoHist2D.h b/Armory/AutoHist2D.h index e7405c4..5e928a2 100644 --- a/Armory/AutoHist2D.h +++ b/Armory/AutoHist2D.h @@ -1,6 +1,7 @@ #ifndef AUTOHIST2D_H #define AUTOHIST2D_H +#include #include #include #include @@ -8,6 +9,7 @@ #include #include #include +#include // Buffers (name -> x,y) pairs during the event loop and, on WriteAll(), // books a TH2F per name with a fixed number of bins per axis and a range @@ -49,6 +51,7 @@ public: TH2F h(name.c_str(), title.c_str(), nbins, xlo, xhi, nbins, ylo, yhi); for (size_t i = 0; i < buf.x.size(); ++i) h.Fill(buf.x[i], buf.y[i]); h.Write(); + SavePNG(h, name); std::cout << "AutoHist2D: wrote \"" << name << "\" (" << buf.x.size() << " entries, " << nbins << "x" << nbins << " bins, range [" @@ -58,11 +61,23 @@ public: } private: + static constexpr const char* kPlotDir = "Plots"; + struct Buffer { std::vector x, y; std::string xTitle, yTitle; }; + // Draws 'h' on a throwaway canvas and saves it under kPlotDir/.png, + // creating the directory first if it doesn't already exist; overwrites + // any existing file of that name. + static void SavePNG(TH2F& h, const std::string& name) { + std::filesystem::create_directories(kPlotDir); + TCanvas c(("c_" + name).c_str(), "", 800, 600); + h.Draw("COLZ"); + c.SaveAs((std::string(kPlotDir) + "/" + name + ".png").c_str()); + } + // Registered histograms, keyed by name. Function-local static avoids // needing a separate translation unit for a header-only class. static std::unordered_map& Registry() { diff --git a/Armory/ClassQQQ.h b/Armory/ClassQQQ.h index 8cfb3d3..ea29baf 100644 --- a/Armory/ClassQQQ.h +++ b/Armory/ClassQQQ.h @@ -274,7 +274,7 @@ inline void QQQ::CalQQQPos(unsigned short ID, id = ID; chBk = chBack; chDn = chDown; - chUp = chUp; + //chUp = chUp; } #endif \ No newline at end of file diff --git a/Armory/Isotope.h b/Armory/Isotope.h index d44cdc9..a81a4a7 100644 --- a/Armory/Isotope.h +++ b/Armory/Isotope.h @@ -109,8 +109,8 @@ private: lineMass200 = 2774; } - char * heliosPath; - bool isFindOnce; + //char * heliosPath; + //bool isFindOnce; }; diff --git a/Armory/Makefile b/Armory/Makefile index dd4cb87..690b172 100644 --- a/Armory/Makefile +++ b/Armory/Makefile @@ -69,6 +69,6 @@ EventBuilder: EventBuilder.cpp ClassData.h fsuReader.h Hit.h @echo "--------- making EventBuilder" $(CXX) $(CXXFLAGS) EventBuilder.cpp -o EventBuilder $(LDFLAGS) -AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h EnergyLoss.h +AnasenMS: anasenMS.cpp constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h EnergyLoss.h AutoHist2D.h @echo "--------- making ANASEN Monte Carlo" $(CXX) $(CXXFLAGS) anasenMS.cpp -o AnasenMS $(LDFLAGS) \ No newline at end of file diff --git a/Armory/Training_Data/.DS_Store b/Armory/Plots/.DS_Store similarity index 64% rename from Armory/Training_Data/.DS_Store rename to Armory/Plots/.DS_Store index 8c1cba7..a37e4d9 100644 Binary files a/Armory/Training_Data/.DS_Store and b/Armory/Plots/.DS_Store differ diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 3ac6dcc..c6644d8 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -53,7 +53,7 @@ bool IsDeadCathode(int id){ } bool IsDeadSX3(int id){ - static std::set dead = {};//{0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3 + static std::set dead = {0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23}; // add dead SX3 IDs here, 0-23 1,7,9,3 return dead.count(id); } @@ -73,7 +73,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){ bool IsDeadSX3BackChannel(int sx3ID, int chBk){ static std::set> dead = { - //{1, 10} + {1, 10} // {sx3ID, back-channel} }; return dead.count({sx3ID, chBk}); @@ -155,25 +155,26 @@ int main(int argc, char **argv){ transfer.SetB(30, 14); // 30Si* heavy product const ReactionConfig reactionConfig = transfer.GetRectionConfig(); const double beamA = reactionConfig.beamA; // mass number of 14N beam - const double beamE = 72 / beamA; // beam energy in MeV + //const double beamE = 72 / beamA; // beam energy in MeV // Excited state lists (projectile and heavy-product excitation states) std::vector ExAList = {0}; // Beam excited energy - std::vector ExList = {3.4}; // Heavy product excited energy + std::vector ExList = {0.0, 2.2, 3.4, 6.0}; // Heavy product excited energy const int kMBeam = reactionConfig.beamA; // mass number of beam const int kMTarget = reactionConfig.targetA; // mass number of target const int kMLight = reactionConfig.recoilLightA; // mass number of light ejectile const int kMHeavy = reactionConfig.recoilHeavyA; // mass number of heavy product - const int kZBeam = reactionConfig.beamZ; // atomic number of beam - const int kZTarget = reactionConfig.targetZ; // atomic number of target - const int kZLight = reactionConfig.recoilLightZ; // atomic number of light ejectile - const int kZHeavy = reactionConfig.recoilHeavyZ; // atomic number of heavy product - bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable - const int decayDaughterA = 20; - const int decayDaughterZ = 10; - const int decayEjectA = 1; - const int decayEjectZ = 1; + + //const int kZBeam = reactionConfig.beamZ; // atomic number of beam + //const int kZTarget = reactionConfig.targetZ; // atomic number of target + //const int kZLight = reactionConfig.recoilLightZ; // atomic number of light ejectile + //const int kZHeavy = reactionConfig.recoilHeavyZ; // atomic number of heavy product + //bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable + //const int decayDaughterA = 20; + //const int decayDaughterZ = 10; + //const int decayEjectA = 1; + //const int decayEjectZ = 1; std::string b; if (reactionConfig.recoilLightA == 1) { @@ -186,7 +187,7 @@ int main(int argc, char **argv){ b = "alpha"; } - TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat"); + //TGraph* elossLight = LoadELoss("../ELoss/HeLoss/E_vs_x_" + b + ".dat"); // define vertex position uniform distribution ranges (mm) double vertexXRange[2] = { -5, 5}; // mm - 5, 5 double vertexYRange[2] = { -5, 5}; // -5, 5 @@ -521,9 +522,9 @@ int main(int argc, char **argv){ trackDir.X() * trackDir.X() + trackDir.Y() * trackDir.Y(); if (transverseDirection2 > 0.0) { - const double pathToRhoMin = + /*const double pathToRhoMin = -(hitPos.X() * trackDir.X() + hitPos.Y() * trackDir.Y()) - / transverseDirection2; + / transverseDirection2;*/ //const TVector3 rhoMin = hitPos + pathToRhoMin * trackDir; const TVector3 rhoMin = vertex; @@ -566,14 +567,18 @@ int main(int argc, char **argv){ if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) { Esx3 = NAN; + beamEnergy = NAN; + Ex = NAN; + continue; } Edet = Esx3; + Eqqq = TMath::QuietNaN(); // mark QQQ energy as invalid for SX3 hit case AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)"); - AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(reTheta * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)"); + AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)"); tree1->Fill(); - }else if (false){//(qqqID >= 0){ + }else if (qqqID >= 0){ // handle QQQ hit case sx3Up = -1; sx3Dn = -1; @@ -639,10 +644,14 @@ int main(int argc, char **argv){ if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) { Eqqq = NAN; + Ex = NAN; + continue; } + Esx3 = TMath::QuietNaN(); // mark SX3 energy as invalid for QQQ hit case Edet = Eqqq; EPC = Eanode - Ecathode; AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)"); + AutoHist2D::Fill("EPC x sin(theta) vs Eqqq", Eqqq, EPC * sin(thetab * TMath::DegToRad()), "Eqqq (MeV)", "EPC x sin(theta) (MeV)"); beamEnergy = TMath::QuietNaN(); // mark beam energy as invalid for QQQ hit case tree1->Fill(); diff --git a/Armory/histcomp.C b/Armory/histcomp.C deleted file mode 100644 index 760f4b7..0000000 --- a/Armory/histcomp.C +++ /dev/null @@ -1,178 +0,0 @@ -void histcomp() { - gROOT->SetBatch(kTRUE); - - // Open file - TFile *f = new TFile("SimAnasen1.root"); - - // Get trees (MAKE SURE names are correct) - TTree *tree1 = (TTree*)f->Get("tree"); - TTree *tree2 = (TTree*)f->Get("tree2"); - - if (!tree1 || !tree2) { - printf("Error: could not find trees. Check names!\n"); - return; - } - - // Create output directory (overwrite-safe) - gSystem->Exec("mkdir -p plots"); - - // Get list of branches - TObjArray *branches = tree1->GetListOfBranches(); - int nBranches = branches->GetEntries(); - //int nBranches = 1; - - // Loop over branches - for (int i = 0; i < nBranches; i++) { - TBranch *br = (TBranch*)branches->At(i); - TString name = br->GetName(); - - //printf("Processing branch: %s\n", name.Data()); - - // Create histograms (auto-range using Draw first) - TString h1name = "h1_" + name; - TString h2name = "h2_" + name; - - // Temporary draw to get range - double min, max; - - 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); - min -= margin; - max += margin; - - TH1D *h1 = new TH1D(h1name, name, 100, min, max); - TH1D *h2 = new TH1D(h2name, name, 100, min, max); - - // Fill histograms - 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); - h1->SetLineWidth(2); - - h2->SetLineColor(kBlue); - h2->SetLineWidth(2); - - // Normalize (optional but useful) - //if (h1->GetEntries() > 0) h1->Scale(1.0 / h1->GetEntries()); - //if (h2->GetEntries() > 0) h2->Scale(1.0 / h2->GetEntries()); - - // Canvas - TCanvas *c = new TCanvas("c", name, 900, 600); //arguments are (name, title, width, height) - - c->SetRightMargin(0.18); - c->Modified(); - c->Update(); - - h1->SetTitle(name + ";"+name+";Counts"); - h1->Draw("HIST"); - h2->Draw("HIST SAME"); - - gPad->Update(); - - TPaveStats *st = (TPaveStats*)h1->FindObject("stats"); - - st->SetX1NDC(0.85); // New X start (left) - st->SetY1NDC(0.5); // New Y start (bottom) - st->SetX2NDC(0.98); // New X end (right) - st->SetY2NDC(0.8); // New Y end (top) - st->Draw(); - gPad->Modified(); - gPad->Update(); - - - // Legend - TLegend *leg = new TLegend(0.65 + .2,0.75 + .1,0.88 + .1,0.88 + .1); - 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); - // Optional: save log plots as well - - 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; - } - - - - } - - // dEb on y, SX3z on x - TH2D *h2d = new TH2D("h2d", "dEb vs SX3z;SX3z (cm);dEb (MeV)", 500, tree2->GetMinimum("sx3Z"), tree2->GetMaximum("sx3Z"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb")); //arguments are (name, title, xbins, xlow, xup, ybins, ylow, yup) - tree2->Draw("dEb:sx3Z>>h2d", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" - TCanvas *c2d = new TCanvas("c2d", "dEb vs SX3z", 900, 600); - h2d->Draw("COLZ"); - c2d->SaveAs("plots/dEb_vs_SX3z.png"); - - TH2D *h2z = new TH2D("h2z", "dEb vs z0", 500, tree2->GetMinimum("z0"), tree2->GetMaximum("z0"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb")); - tree2->Draw("dEb:z0>>h2z", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" - TCanvas *c2z = new TCanvas("c2z", "dEb vs z0", 900, 600); - h2z->Draw("COLZ"); - c2z->SaveAs("plots/dEb_vs_z0.png"); - - TH2D *h2theta = new TH2D("h2theta", "dEb vs reTheta", 500, tree2->GetMinimum("reTheta"), tree2->GetMaximum("reTheta"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb")); - tree2->Draw("dEb:reTheta>>h2theta", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" - TCanvas *c2theta = new TCanvas("c2theta", "dEb vs reTheta", 900, 600); - h2theta->Draw("COLZ"); - c2theta->SaveAs("plots/dEb_vs_reTheta.png"); - - TH2D *h2phi = new TH2D("h2phi", "dEb vs rePhi", 500, tree2->GetMinimum("rePhi"), tree2->GetMaximum("rePhi"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb")); - tree2->Draw("dEb:rePhi>>h2phi", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" - TCanvas *c2phi = new TCanvas("c2phi", "dEb vs rePhi", 900, 600); - h2phi->Draw("COLZ"); - c2phi->SaveAs("plots/dEb_vs_rePhi.png"); - - TH2D *h2dE = new TH2D("h2dE", "dEb vs Tb", 500, tree2->GetMinimum("Tb"), tree2->GetMaximum("Tb"), 500, tree2->GetMinimum("dEb"), tree2->GetMaximum("dEb")); - tree2->Draw("dEb:Tb>>h2dE", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" - TCanvas *c2dE = new TCanvas("c2dE", "dEb vs Tb", 900, 600); - h2dE->Draw("COLZ"); - c2dE->SaveAs("plots/dEb_vs_Tb.png"); - - printf("Done! Plots saved in ./plots/\n"); -} diff --git a/Armory/run_script.C b/Armory/run_script.C deleted file mode 100644 index e798df6..0000000 --- a/Armory/run_script.C +++ /dev/null @@ -1,4 +0,0 @@ -.L ANASEN_model.C -.L anasenMS_root.cpp+ -ANASEN_model(); -Run(10); diff --git a/ELoss/.DS_Store b/ELoss/.DS_Store index eb54b34..ab248a0 100644 Binary files a/ELoss/.DS_Store and b/ELoss/.DS_Store differ diff --git a/ELoss/AlSi.zip b/ELoss/AlSi.zip deleted file mode 100644 index 3f9e153..0000000 Binary files a/ELoss/AlSi.zip and /dev/null differ diff --git a/ELoss/EnergyLoss.py b/ELoss/EnergyLoss.py index e22a416..5312008 100644 --- a/ELoss/EnergyLoss.py +++ b/ELoss/EnergyLoss.py @@ -565,30 +565,6 @@ class MyInteractiveApp(cmd.Cmd): print("Please input particle, final energy from detector, and distance travelled") - def do_uproot_file(self, arg): - """Open a specific root file for inspection""" - - args = shlex.split(arg) - - if len(args) > 0: - filename = args[0] - else: - filename = self.rootFile - - try: - print(f"Opening {filename}") - - # Try Armory path first - try: - self.file = uproot.open(f"../Armory/{filename}") - except FileNotFoundError: - self.file = uproot.open(filename) - - print("File loaded successfully.") - print("Keys:", self.file.keys()) - - except Exception as e: - print("Error opening file:", e) if __name__ == "__main__":