diff --git a/Armory/aarootscript.C b/Armory/aarootscript.C index 96e454f..258d030 100644 --- a/Armory/aarootscript.C +++ b/Armory/aarootscript.C @@ -46,10 +46,10 @@ void aarootscript(int argument = 0) { std::cout << "Zero Count: " << zeroCount << std::endl; - //std::cout << "Making histograms..." << std::endl; + std::cout << "Making histograms..." << std::endl; gErrorIgnoreLevel = 2001; - //gROOT->ProcessLine(".x histcomp.C"); + gROOT->ProcessLine(".x histcomp.C"); std::cout << "=========================================\n"; diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 9b322a1..c177818 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -12,6 +12,8 @@ #include "ClassTransfer.h" // Reaction kinematics and MC event generation #include "ClassAnasen.h" // ANASEN detector model classes (SX3, PW, etc.) #include "TCanvas.h" // ROOT canvas for drawing +#include +#include //======== Generate light particle based on reaction // calculate real and reconstructed tracks and Q-value uncertainty @@ -31,18 +33,22 @@ int main(int argc, char **argv){ // number of events can be overridden from command line int numEvent = 1000000; 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_HeAlpha"); // for light particle (alpha) + //double density = (2.1525e-7) * 1000; // example for aluminum target, adjust as needed (400 torr is 0.0000861) + //char command[256]; + //snprintf(command, sizeof(command), "python3 ../ELoss/EvXconverter.py %f", density); + //printf("Command: %s\n", command); + //system(command); // run the conversion script to generate energy loss tables with correct density + // load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm^2), density in mg/cm^3) + TGraph* elossLight = LoadELoss("../ELoss/E_vs_x_light"); // for light particle (alpha) TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p"); // for heavy particle (proton) - double density = (2.1525e-7) * 400; // example for aluminum target, adjust as needed (400 torr is 0.0000861) + 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->SetTitle("Energy Loss Table;cm;Kinetic Energy (MeV)"); g->Draw("ALP"); g->SetLineColor(kRed); - c1->SetLogy(); - c1->SetLogx(); + //c1->SetLogy(); + //c1->SetLogx(); c1->Print("eloss_light.png"); int ELossTotal = 0; @@ -129,12 +135,14 @@ int main(int argc, char **argv){ // outgoing particles in lab frame (light/heavy) double thetab, phib, Tb; double thetaB, phiB, TB; + double dEb; tree->Branch("thetab", &thetab, "thetab/D"); tree->Branch("phib", &phib, "phib/D"); tree->Branch("Tb", &Tb, "Tb/D"); tree->Branch("thetaB", &thetaB, "thetaB/D"); tree->Branch("phiB", &phiB, "phiB/D"); tree->Branch("TB", &TB, "TB/D"); + tree->Branch("dEb", &dEb, "dEb/D"); // excitation state identifiers int ExAID; @@ -297,6 +305,7 @@ int main(int argc, char **argv){ rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg(); z0 = pw->GetZ0(); + dEb = 0; tree->Fill(); //Energy loss double dl = (hitPos - vertex).Mag() / 10; // path length in cm (positions in mm) @@ -304,9 +313,10 @@ int main(int argc, char **argv){ printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl); } double tb_temp = Tb; - double dx = 0; - double counter = 0; - while(dx < dl){ + //double dx = 0; + //double counter = 0; + //energy loss loop + /*while(dx < dl){ double step = 0.1; // 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 @@ -324,6 +334,13 @@ int main(int argc, char **argv){ break; } } + dEb = tb_temp - Tb; // total energy loss*/ + TGraph *invg = new TGraph(elossLight->GetN(), elossLight->GetY(), elossLight->GetX()); + double x0 = invg->Eval(Tb); + //double x0 = elossLight->GetX(Tb); // range corresponding to final kinetic energy + x0 = x0 + dl; + Tb = elossLight->Eval(x0); // kinetic energy corresponding to range at hit position + dEb = tb_temp - Tb; // total energy loss // fill tree2 with energy loss adjusted data if (Tb != 0) { tree2->Fill(); @@ -349,6 +366,7 @@ int main(int argc, char **argv){ reTheta1 = TMath::QuietNaN(); rePhi1 = TMath::QuietNaN(); z0 = TMath::QuietNaN(); + dEb = TMath::QuietNaN(); //Tb = -12354567; // mark kinetic energy as invalid for no hit case // fill tree with original data (no energy loss for these events) //tree->Fill();