fixed detached head ugh

Armory/anasenMS.cpp

@@ -34,15 +34,15 @@ int main(int argc, char **argv){
   // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm²), density in mg/cm³)
   TGraph* elossLight = LoadELoss("../ELoss/Eloss_alpha"); // for light particle (alpha)
   TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p");     // for heavy particle (proton)
-  double density = 2.7; // example for aluminum target, adjust as needed
+  double density = 0.0000876; // example for aluminum target, adjust as needed
 
   // Reaction setup: projectile + target configuration, energy, and product IDs
   TransferReaction transfer;
 
   transfer.SetA(24,12, 0);           // e.g., 24Mg (Z=12) with 0 excitation
   transfer.SetIncidentEnergyAngle(10, 0, 0); // 10 MeV beam, 0 polar and azimuthal angle
-  transfer.Seta( 4, 2);              // identify reaction product a in internal indexing
+  transfer.Seta( 4, 2);              // identify reaction product a in internal indexing e.g., 4He (alpha)
-  transfer.Setb( 1, 1);              // identify reaction product b
+  transfer.Setb( 1, 1);              // identify reaction product b e.g., 1H (proton)
 
   // TODO add alpha source or alternative reaction channel selection
 
@@ -100,7 +100,7 @@ int main(int argc, char **argv){
   printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
   TFile * saveFile = new TFile(saveFileName, "recreate");
   TTree * tree = new TTree("tree", "tree");
-  TTree* tree2 = tree->CloneTree(0); // for 2D histograms or alternative data structure if needed
+   
 
   // beam and CM variables saved in tree
   double KEA;
@@ -175,6 +175,8 @@ int main(int argc, char **argv){
   double z0;
   tree->Branch("z0",     &z0, "reconstucted_Z/D");
 
+  TTree* tree2 = tree->CloneTree(0);
+  tree2->SetName("tree2");
 
   //========timer
   TBenchmark clock;
@@ -268,10 +270,34 @@ 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();*/
 
-      // fill tree with original data before energy loss
+      // reconstruct track from PW readings + SX3 hit
+      pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
+      reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi   = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      // alternative track algorithm with uncertainty parameters
+      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
+      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
+      rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
+
+      z0 = pw->GetZ0();
       tree->Fill();
-
       // 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();
@@ -286,19 +312,6 @@ int main(int argc, char **argv){
       }
       // 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();
-      rePhi   = pw->GetTrackPhi() * TMath::RadToDeg();
-
-      // alternative track algorithm with uncertainty parameters
-      pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
-      reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
-      rePhi1   = pw->GetTrackPhi() * TMath::RadToDeg();
-
-      z0 = pw->GetZ0();
-
       // fill tree2 with energy loss adjusted data
       tree2->Fill();
 
@@ -321,6 +334,7 @@ int main(int argc, char **argv){
 
       // fill tree with original data (no energy loss for these events)
       tree->Fill();
+      tree2->Fill();
     }
 
     //#################################################################### Timer
@@ -343,8 +357,10 @@ int main(int argc, char **argv){
   }
 
   // write results to ROOT file and close
-  tree->Write();
+  //tree->Write();
-  tree2->Write();
+  //tree2->Write();
+  tree->Write("", TObject::kOverwrite);
+  tree2->Write("", TObject::kOverwrite);
   int count = tree->GetEntries();
   int count2 = tree2->GetEntries();
   saveFile->Close();

Armory/histcomp.C

@@ -0,0 +1,89 @@
+void histcomp() {
+
+// 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
+    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");
+}