diff --git a/Armory/aarootscript.C b/Armory/aarootscript.C new file mode 100644 index 0000000..2c82588 --- /dev/null +++ b/Armory/aarootscript.C @@ -0,0 +1,34 @@ +#include "TFile.h" +#include "TTree.h" +#include "TGraph.h" +#include "TLegend.h" +#include "TCanvas.h" +#include "TH1D.h" +#include "TObjArray.h" +#include "TBranch.h" +#include + +void aarootscript() { + TFile *f = new TFile("SimAnasen1.root"); + TTree *tree1 = (TTree*)f->Get("tree1"); + TTree *tree2 = (TTree*)f->Get("tree2"); + + TTreeReader reader("tree"); + TTreeReaderValue branchValue(reader, "Tb"); + double totalSum = 0; + while (reader.Next()) { + totalSum += *branchValue; + } + std::cout << "Total Sum: " << totalSum << std::endl; + + TTreeReader reader2("tree2"); + TTreeReaderValue branchValue2(reader2, "Tb"); + double totalSum2 = 0; + while (reader2.Next()) { + totalSum2 += *branchValue2; + } + std::cout << "Total Sum: " << totalSum2 << std::endl; + + std::cout << "Difference: " << totalSum - totalSum2 << std::endl; + +} \ No newline at end of file diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index d2d1185..190de60 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -270,22 +270,6 @@ int main(int argc, char **argv){ visTrackHitPos.push_back(hitPos); visTrackWires.push_back({anodeID[0], cathodeID[0]}); } - /* - // apply energy loss from vertex to SX3 hit position (for light particle) - double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm) - double EkinLight = Pb.E() - Pb.M(); - double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx - double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling) - if (dE_light < EkinLight) { - Pb.SetE(Pb.E() - dE_light); - // update momentum to conserve direction - double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M()); - TVector3 dir_new = Pb.Vect().Unit() * p_new; - Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E()); - } - // update kinetic energy after loss - Tb = Pb.E() - Pb.M();*/ - // reconstruct track from PW readings + SX3 hit pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false); reTheta = pw->GetTrackTheta() * TMath::RadToDeg(); @@ -298,22 +282,26 @@ int main(int argc, char **argv){ z0 = pw->GetZ0(); tree->Fill(); + //printf("Event %d: Tb before energy loss = %.2f MeV\n", i, Tb); // apply energy loss from vertex to SX3 hit position (for light particle) - double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm) - double EkinLight = Pb.E() - Pb.M(); + double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm) + //printf("Event %d: dl = %f cm\n", i, dl); + double EkinLight = Pb.E() - Pb.M();// kinetic energy of light particle before loss double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx + //printf("Event %d: dE/dx = %f MeV/(mg/cm^2)\n", i, dedxLight); double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling) - if (dE_light < EkinLight) { - Pb.SetE(Pb.E() - dE_light); - // update momentum to conserve direction - double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M()); - TVector3 dir_new = Pb.Vect().Unit() * p_new; - Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E()); - } + //printf("Event %d: dE_light = %f MeV\n", i, dE_light); + //if (dE_light < EkinLight) { + Pb.SetE(Pb.E() - dE_light); // update energy after loss + double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M()); // update momentum to conserve direction + TVector3 dir_new = Pb.Vect().Unit() * p_new; // maintain original direction, adjust magnitude + Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E()); // update 4-vector with new energy and momentum + //} // update kinetic energy after loss - Tb = Pb.E() - Pb.M(); + Tb = Tb - dE_light; // fill tree2 with energy loss adjusted data tree2->Fill(); + //printf("Event %d: Tb after energy loss = %.2f MeV\n", i, Tb); }else{ // no valid SX3 hit: mark clearly invalid diff --git a/Armory/histcomp.C b/Armory/histcomp.C index bd3393b..66214fc 100644 --- a/Armory/histcomp.C +++ b/Armory/histcomp.C @@ -1,89 +1,107 @@ void histcomp() { + gROOT->SetBatch(kTRUE); -// Open file -TFile *f = new TFile("SimAnasen1.root"); + // 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"); + // 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 - tree1->Draw(name, "", "goff"); - double min = tree1->GetMinimum(name); - double max = tree1->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 - 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, 800, 600); - - h1->SetTitle(name + ";"+name+";Counts"); - h1->Draw("HIST"); - h2->Draw("HIST SAME"); - - // Legend - TLegend *leg = new TLegend(0.65,0.75,0.88,0.88); - leg->AddEntry(h1, "tree1", "l"); - leg->AddEntry(h2, "tree2", "l"); - leg->Draw(); - - // Save plot (overwrite each run) - TString filename = "plots/" + name + ".png"; - c->SaveAs(filename); - - // Clean up - delete c; - delete h1; - delete h2; - -} - -printf("Done! Plots saved in ./plots/\n"); + 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 + tree1->Draw(name, "", "goff"); + double min = tree1->GetMinimum(name); + double max = tree1->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 + 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); + + 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(); + + // Save plot (overwrite each run) + TString filename = "plots/" + name + ".png"; + c->SaveAs(filename); + + // Clean up + delete c; + delete h1; + delete h2; + + } + + printf("Done! Plots saved in ./plots/\n"); }