Energy loss loop

Armory/ClassAnasen.h

@@ -108,12 +108,12 @@ 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);
   //--- define some media
-  TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum); //name, number of materials, material
+  //TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum); //name, number of materials, material
   TGeoMedium *He = new TGeoMedium("He",2, matHe);
   TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
   TGeoMedium *Si = new TGeoMedium("Si",3, matSi);

Armory/anasenMS.cpp

@@ -33,21 +33,29 @@ int main(int argc, char **argv){
   if( argc >= 2 ) numEvent = atoi(argv[1]);
 
   // 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* elossLight = LoadELoss("../ELoss/Eloss_HeAlpha"); // for light particle (alpha)
   TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p");     // for heavy particle (proton)
-  double density = 0.0000876; // example for aluminum target, adjust as needed
+  double density = 0.0000861; // example for aluminum target, adjust as needed
   auto c1 = new TCanvas("c1", "Graph Example", 800, 600);
   auto g = elossLight;
   g->SetTitle("Energy Loss Table;Kinetic Energy (MeV);dE/dx (MeV/(mg/cm^{2}))");
   g->Draw("ALP");
   g->SetLineColor(kRed);
+  c1->SetLogy();
+  c1->SetLogx();
   c1->Print("eloss_light.png");
 
+  auto g2 = elossHeavy;
+  g2->SetTitle("Energy Loss Table;Kinetic Energy (MeV);dE/dx (MeV/(mg/cm^{2}))");
+  g2->Draw("ALP");
+  g2->SetLineColor(kBlue);
+  c1->Print("eloss_heavy.png");
+
   // 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.SetIncidentEnergyAngle(5.46, 0, 0); // 5.46 MeV beam, 0 polar and azimuthal angle
   transfer.Seta( 4, 2);              // identify reaction product a in internal indexing e.g., 4He (alpha)
   transfer.Setb( 1, 1);              // identify reaction product b e.g., 1H (proton)
 
@@ -289,26 +297,36 @@ int main(int argc, char **argv){
 
       z0 = pw->GetZ0();
       tree->Fill();
-      //printf("Event %d: Tb before energy loss = %.2f MeV\n", i, Tb);
+      //Energy loss
-      // apply energy loss from vertex to SX3 hit position (for light particle)
       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 EkinLight = Pb.E() - Pb.M();// kinetic energy of light particle before loss
+      //double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
-      double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
+      if (numEvent <= 10){
-      //printf("Event %d: dE/dx = %f MeV/(mg/cm^2)\n", i, dedxLight);
+        printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl);
-      double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
+      }
-      //printf("Event %d: dE_light = %f MeV\n", i, dE_light);
+      double tb_temp = Tb;
-      //if (dE_light < EkinLight) {
+      //double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
-      Pb.SetE(Pb.E() - dE_light); // update energy after loss
+      double dx = 0;
-      double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M()); // update momentum to conserve direction
+      double counter = 0;
-      TVector3 dir_new = Pb.Vect().Unit() * p_new; // maintain original direction, adjust magnitude
+      while(dx < dl){
-      Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E()); // update 4-vector with new energy and momentum
+        double step = 0.01; // cm, step size for energy loss calculation
-      //}
+        if(dx + step > dl) step = dl - dx; // adjust last step to end at hit position
-      // update kinetic energy after loss
+        double EkinStep = Tb; // kinetic energy at current step
-      Tb = Tb - dE_light;
+        double dedxStep = elossLight->Eval(EkinStep); // dE/dx at current energy
+        double dE_step = dedxStep * step * density; // energy loss for this step
+        if (numEvent <= 10 && fmod(counter, 10) == 0){ 
+          printf("Event %d: step = %f, Ekin = %f MeV, dE/dx = %f MeV/(g/cm^2), dE_step = %f MeV\n", i, counter, EkinStep, dedxStep, dE_step);
+        }
+        Pb.SetE(Pb.E() - dE_step); // update energy after loss for this step
+        dx += step;
+        counter++;
+        Tb = Tb - dE_step; // update kinetic energy for tree storage
+      }
       // fill tree2 with energy loss adjusted data
       tree2->Fill();
-      //printf("Event %d: Tb after energy loss = %.2f MeV\n", i, Tb);
+      if (numEvent <= 10){
+        printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb);
+      }
 
     }else{
       // no valid SX3 hit: mark clearly invalid
@@ -328,8 +346,8 @@ int main(int argc, char **argv){
       z0       = TMath::QuietNaN();
 
       // fill tree with original data (no energy loss for these events)
-      tree->Fill();
+      //tree->Fill();
-      tree2->Fill();
+      //tree2->Fill();
     }
 
     //#################################################################### Timer