diff --git a/Armory/ClassPW.h b/Armory/ClassPW.h index bb33b80..fbedea3 100755 --- a/Armory/ClassPW.h +++ b/Armory/ClassPW.h @@ -509,10 +509,10 @@ inline void PW::CalTrack2(TVector3 sx3Pos, PWHitInfo hitInfo, double sigmaA, dou inline double PW::GetZ0() { - double x = trackPos.X(); - double y = trackPos.Y(); - double rho = TMath::Sqrt(x * x + y * y); - double theta = trackVec.Theta(); + [[maybe_unused]]double x = trackPos.X(); + [[maybe_unused]]double y = trackPos.Y(); + [[maybe_unused]]double rho = TMath::Sqrt(x * x + y * y); + [[maybe_unused]]double theta = trackVec.Theta(); return trackVec.Z(); } diff --git a/Armory/ClassTransfer.h b/Armory/ClassTransfer.h index c12c43f..ba89669 100644 --- a/Armory/ClassTransfer.h +++ b/Armory/ClassTransfer.h @@ -16,7 +16,6 @@ #include "Isotope.h" class ReactionConfig{ - public: ReactionConfig(){} @@ -47,9 +46,9 @@ public: std::vector beamEx; ///excitation_energy_of_A[MeV] - void SetReaction(int beamA, int beamZ, - int targetA, int targetZ, - int recoilA, int recoilZ, float beamEnergy_AMeV){ + void SetReaction(int beamA, int beamZ, // projectile + int targetA, int targetZ, // target + int recoilA, int recoilZ, float beamEnergy_AMeV){ // light recoil, e.g. alpha this->beamA = beamA; this->beamZ = beamZ; this->targetA = targetA; @@ -178,8 +177,8 @@ public: void Setb(int A, int Z); void SetB(int A, int Z); void SetIncidentEnergyAngle(double KEA, double theta, double phi); - void SetExA(double Ex); - void SetExB(double Ex); + void SetExA(double Ex); // excitation energy of A in MeV + void SetExB(double Ex); // excitation energy of B in MeV void SetReactionFromFile(string settingFile); TString GetReactionName(); @@ -247,7 +246,7 @@ TransferReaction::TransferReaction(){ Seta(4,2); Setb(1,1); SetB(27,13); - TA = 2.5; + TA = 2.5; // MeV/u T = TA * reaction.beamA; ExA = 0; @@ -311,7 +310,7 @@ void TransferReaction::SetB(int A, int Z){ isBSet = true; } -void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){ +void TransferReaction::SetIncidentEnergyAngle(double KEA, double theta, double phi){ // KEA in MeV/u, theta and phi in degree this->TA = KEA; this->T = TA * reaction.beamA; this->thetaIN = theta; diff --git a/Armory/aarootscript.C b/Armory/aarootscript.C index 3aa9a22..dd6c0ca 100644 --- a/Armory/aarootscript.C +++ b/Armory/aarootscript.C @@ -55,7 +55,7 @@ void aarootscript(int argument = 0) { //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/analyze.C b/Armory/analyze.C index 4f82470..fe0c2e1 100644 --- a/Armory/analyze.C +++ b/Armory/analyze.C @@ -51,7 +51,6 @@ void analyze(const char* filename = "SimAnasen1.root") double max = tree->GetMaximum("TB"); min = min - max*0.1; max = max * 1.1; - printf("Tb min: %f, TB max: %f\n", min, max); TH1D *hTb = new TH1D("hTb","Tb and TB;Energy (MeV);Counts",200,min,max); //arguments are name, title (with axis labels), number of bins, x-min, x-max TH1D *hTB = new TH1D("hTB","",200,min,max); @@ -89,9 +88,9 @@ void analyze(const char* filename = "SimAnasen1.root") // Tb vs TB correlation (with gate) - TCanvas *c3 = new TCanvas("c3","Tb vs TB",800,600); + TCanvas *c3 = new TCanvas("c3","dEb vs SX3z",800,600); - tree->Draw("TB:Tb>>h2(200,min,max,200,min,max)","sx3ID>=0","COLZ"); + tree->Draw("TB:Tb>>h2(200,min,max,200,min,max)","","COLZ"); //arguments are "y:x>>histogram(bins,xmin,xmax,bins,ymin,ymax)", "selection", "options" c3->SaveAs("Tb_vs_TB.png"); diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index 097b016..25b344d 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -41,10 +41,12 @@ bool IsDeadCathode(int id){ } bool IsDeadSX3(int id){ - static std::set dead = {}; // add dead SX3 IDs here, 0-23 + static std::set dead = {}; // add dead SX3 IDs here, 0-23 1,7,9,3 return dead.count(id); } +static std::set> ReactionProductb = { {1,1} }; // add reaction product b (light particle) A,Z pairs here, e.g. {1,1} for proton, {4,2} for alpha + int main(int argc, char **argv){ printf("=========================================\n"); @@ -56,43 +58,43 @@ 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^2), density in mg/cm^3) - TGraph* elossLight = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha) - TGraph* elossHeavy = LoadELoss("../ELoss/E_vs_x_proton.dat"); // for heavy particle (proton) - //TGraph* elossRecoil = LoadELoss("../ELoss/E_vs_x_recoil.dat"); // for recoil particle (if needed) - TGraph *invgLight = new TGraph(elossLight->GetN(), elossLight->GetY(), elossLight->GetX()); - TGraph *invgHeavy = new TGraph(elossHeavy->GetN(), elossHeavy->GetY(), elossHeavy->GetX()); + TGraph* elossAlpha = LoadELoss("../ELoss/E_vs_x_alpha.dat"); // for light particle (alpha) + TGraph* elossProton = LoadELoss("../ELoss/E_vs_x_proton.dat"); // for heavy particle (proton) + TGraph *invgAlpha = new TGraph(elossAlpha->GetN(), elossAlpha->GetY(), elossAlpha->GetX()); + TGraph *invgProton = new TGraph(elossProton->GetN(), elossProton->GetY(), elossProton->GetX()); - /* + //Plot energy loss tables (sanity check), vis will not work if this is ran without X11 display (e.g. on cluster), so comment out if running in headless mode auto c1 = new TCanvas("c1", "Graph Example", 800, 600); - auto g = elossLight; - g->SetTitle("Energy Loss Table (Light);cm;Kinetic Energy (MeV)"); + auto g = elossAlpha; + g->SetTitle("Energy Loss Table (Alpha);cm;Kinetic Energy (MeV)"); g->Draw("ALP"); g->SetLineColor(kRed); //c1->SetLogy(); //c1->SetLogx(); - c1->Print("eloss_light.png"); + c1->Print("eloss_alpha.png"); auto c2 = new TCanvas("c2", "Graph Example", 800, 600); - auto g2 = elossHeavy; - g2->SetTitle("Energy Loss Table (Heavy);cm;Kinetic Energy (MeV)"); + auto g2 = elossProton; + g2->SetTitle("Energy Loss Table (Proton);cm;Kinetic Energy (MeV)"); g2->Draw("ALP"); g2->SetLineColor(kBlue); - c2->Print("eloss_heavy.png");*/ + c2->Print("eloss_proton.png"); // Reaction setup: projectile + target configuration, energy, and product IDs TransferReaction transfer; - transfer.SetA(27,13, 0); // e.g., 24Mg (Z=12) with 0 excitation - transfer.SetIncidentEnergyAngle(72, 0, 0); // 5.46 MeV beam for alpha source, 0 polar and azimuthal angle, 72 for 27Al + transfer.SetA(14, 7, 0); // e.g., 24Mg (Z=12) with 0 excitation + transfer.SetIncidentEnergyAngle((42.82/14.0), 0, 0); // arguments are KEA in MeV/u, theta and phi in degree 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) + transfer.Setb(ReactionProductb.begin()->first, ReactionProductb.begin()->second); // identify reaction product b e.g., 1H (proton) + transfer.SetB(17, 8); // identify reaction product B e.g., 23Na (Z=11) // TODO add alpha source or alternative reaction channel selection // Excited state lists (target and projectile/excited products) - std::vector ExAList = {0}; // projectile excitation states in MeV, e.g., 0 for ground state, 1.37 for first excited state of 24Mg, etc. - std::vector ExList = {0, 1, 2}; // target excitation states in MeV, e.g., 0 for ground state, 1.37 for first excited state of 24Mg, etc. + std::vector ExAList = {0}; // projectile excitation states in MeV + std::vector ExList = {0}; // target excitation states in MeV // define vertex position uniform distribution ranges (mm) double vertexXRange[2] = { -5, 5}; // mm @@ -263,14 +265,14 @@ int main(int argc, char **argv){ thetab = Pb.Theta() * TMath::RadToDeg(); thetaB = PB.Theta() * TMath::RadToDeg(); - Tb = Pb.E() - Pb.M(); - TB = PB.E() - PB.M(); + Tb = (Pb.E() - Pb.M()); // kinetic energy of light particle at vertex (before energy loss) units of MeV + TB = (PB.E() - PB.M()); T[0] = Tb; T[1] = TB; - if (Tb < 1.5) { - //skip event if light particle energy after loss is below detection threshold of 1.5 MeV - continue; - } + //if (Tb < 1.5) { + // //skip event if light particle energy after loss is below detection threshold of 1.5 MeV + // continue; + //} phib = Pb.Phi() * TMath::RadToDeg(); phiB = PB.Phi() * TMath::RadToDeg(); @@ -350,39 +352,38 @@ int main(int argc, char **argv){ //Energy loss double dl = (hitPos - vertex).Mag(); // path length in units of cm - if (numEvent <= 100){ //printf("Event %d: Ekin before loss = %f MeV, distance = %f cm\n", i, Tb, dl); //printf("Total T before loss: %f MeV\n", T); } - double tb_temp = Tb; - double tB_temp = TB; dEb = tb_temp - Tb; // total energy loss - dEB = tB_temp - TB; // total energy loss for heavy particle, currently set equal to light particle loss for simplicity - - double x0light = invgLight->Eval(Tb); - double x0heavy = invgHeavy->Eval(TB); - - x0light = x0light + dl; - x0heavy = x0heavy + dl; - - Tb = elossLight->Eval(x0light); // kinetic energy corresponding to range at hit position - TB = elossHeavy->Eval(x0heavy); // kinetic energy for heavy particle, currently set equal to light particle loss for simplicity - - if (Tb < 1.5) { - //skip event if light particle energy after loss is below detection threshold of 1.5 MeV - continue; + if (ReactionProductb.count({4, 2})){ // if light particle is alpha, use alpha energy loss table + double x0b = invgAlpha->Eval(Tb); + x0b = x0b + dl; + Tb = elossAlpha->Eval(x0b); + } else if (ReactionProductb.count({1, 1})){ // if light particle is proton, use proton energy loss table + double x0b = invgProton->Eval(Tb); + x0b = x0b + dl; + Tb = elossProton->Eval(x0b); + } else { + // for other particle types, can add additional energy loss tables or use a generic approximation + // for now, we will just apply a simple linear energy loss as a placeholder + double dE_dx = 5; // MeV/cm, placeholder value for energy loss per unit length + Tb = Tb - dE_dx * dl; } + //if (Tb < 0) { + // Tb = TMath::QuietNaN(); + //} + dEb = tb_temp - Tb; // total energy loss - dEB = tB_temp - TB; // total energy loss for heavy particle, currently set equal to light particle loss for simplicity // fill tree2 with energy loss adjusted data //Fill T so it can make a histogram of both Tb and TB in root script T[0] = Tb; - T[1] = TB; + T[1] = 0; //to plot both as one histogram in root, can use tree2->Draw("T(0)"); for light particle and tree2->Draw("T(1)") for heavy particle tree2->Fill(); @@ -390,7 +391,7 @@ int main(int argc, char **argv){ if (numEvent <= 10){ //printf("Event %d: Tb after energy loss = %f MeV, energy loss = %f MeV\n", i, Tb, tb_temp - Tb); } //to give in scientific notation, use %e instead of %f in the printf format string. For example: printf("Event %d: Tb after energy loss = %e MeV, energy loss = %e MeV\n", i, Tb, tb_temp - Tb); - ELossTotal += (tb_temp - Tb) + (tB_temp - TB); + ELossTotal += (tb_temp - Tb); }else{ // no valid SX3 hit: mark clearly invalid @@ -446,8 +447,8 @@ int main(int argc, char **argv){ saveFile->Close(); printf("=============== done. saved as %s. tree entries: %d, tree2 entries: %d\n", saveFileName.Data(), count, count2); - printf("Total energy loss across all events: %e MeV\n", (double)ELossTotal); - printf("Average energy loss across events: %e MeV\n", (double)ELossTotal / count); + printf("Total energy loss across all events: %f MeV\n", (double)ELossTotal); + printf("Average energy loss across events: %f MeV\n", (double)ELossTotal / count); if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis" printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size()); @@ -508,8 +509,8 @@ int main(int argc, char **argv){ } delete anasen; - delete elossLight; - delete elossHeavy; + delete elossAlpha; + delete elossProton; return 0; diff --git a/Armory/histcomp.C b/Armory/histcomp.C index bca0e2f..a0b1d21 100644 --- a/Armory/histcomp.C +++ b/Armory/histcomp.C @@ -143,5 +143,18 @@ void histcomp() { } + // dEb on y, SX3z on x + TH2D *h2d = new TH2D("h2d", "dEb vs SX3z;SX3z (cm);dEb (MeV)", 500, -110, 110, 500, 0, 12); //arguments are (name, title, xbins, xlow, xup, ybins, ylow, yup) + tree2->Draw("dEb:sx3Z>>h2d", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" + TCanvas *c2d = new TCanvas("c2d", "dEb vs SX3z", 900, 600); + h2d->Draw("COLZ"); + c2d->SaveAs("plots/dEb_vs_SX3z.png"); + + TH2D *h2z = new TH2D("h2z", "dEb vs z0", 500, -1, 1, 500, 0, 12); + tree2->Draw("dEb:z0>>h2z", "", "goff"); // arguments are "y:x>>histogram", "selection", "options" + TCanvas *c2z = new TCanvas("c2z", "dEb vs z0", 900, 600); + h2z->Draw("COLZ"); + c2z->SaveAs("plots/dEb_vs_z0.png"); + printf("Done! Plots saved in ./plots/\n"); }