diff --git a/Armory/ClassAnasen.h b/Armory/ClassAnasen.h index ff7bd1e..1502d0d 100644 --- a/Armory/ClassAnasen.h +++ b/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); diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index c6176f9..d6a0b7c 100644 --- a/Armory/anasenMS.cpp +++ b/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); - // apply energy loss from vertex to SX3 hit position (for light particle) + //Energy loss 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) - //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 = Tb - dE_light; + //double EkinLight = Pb.E() - Pb.M();// kinetic energy of light particle before loss + //double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx + if (numEvent <= 10){ + printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl); + } + double tb_temp = Tb; + //double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling) + double dx = 0; + double counter = 0; + while(dx < dl){ + 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 + double EkinStep = Tb; // kinetic energy at current step + 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(); - tree2->Fill(); + //tree->Fill(); + //tree2->Fill(); } //#################################################################### Timer