From a2e0d2955964a9159217a1869b82a3ac094e7ee7 Mon Sep 17 00:00:00 2001
From: Alexander Jentsch <ajentsch@bnl.gov>
Date: Wed, 11 Aug 2021 09:06:41 -0400
Subject: [PATCH] Working on updating particle gun and hits reader for B0.
---
benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx | 14 +++++++-------
.../b0_tracker/scripts/gen_forward_protons.cxx | 13 +++++++++----
2 files changed, 16 insertions(+), 11 deletions(-)
diff --git a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
index b4984ada..880973d4 100644
--- a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
+++ b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
@@ -40,21 +40,21 @@ void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.root"
std::vector<double> result;
for (const auto& h : hits) {
ROOT::Math::XYZVector vec(h.position.x,h.position.y,h.position.z);
- result.push_back(vec.eta());
- std::cout << vec.eta() << "\n";
+ result.push_back(1000*vec.theta());
+ std::cout << 1000*vec.theta() << "\n";
}
return result;
};
- auto d1 = d0.Define("hits_eta", hits_eta, {"B0TrackerHits"});
+ auto d1 = d0.Define("hits_theta", hits_theta, {"B0TrackerHits"});
- auto h1 = d1.Histo1D({"h1", "hits_eta", 100, 0,20}, "hits_eta");
+ auto h1 = d1.Histo1D({"h1", "hits_theta", 100, 0,20}, "hits_theta");
TCanvas* c = new TCanvas();
h1->DrawCopy();
- c->SaveAs("results/b0_tracker_hits_eta.png");
- c->SaveAs("results/b0_tracker_hits_eta.pdf");
+ c->SaveAs("results/b0_tracker_hits_theta.png");
+ c->SaveAs("results/b0_tracker_hits_theta.pdf");
auto n1 = h1->GetMean();
- std::cout << "Pseudorapidity of hits: " << n1 << std::endl;
+ std::cout << "Polar angle of hits: " << n1 << std::endl;
//if (n1 < 5) {
// std::quick_exit(1);
diff --git a/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
index a8884ca1..1773423d 100644
--- a/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
+++ b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
@@ -20,8 +20,13 @@ using namespace HepMC3;
void gen_forward_protons(int n_events = 100,
const char* out_fname = "forward_protons.hepmc")
{
- double cos_theta_min = std::cos(1.0*(M_PI/180.0));
- double cos_theta_max = std::cos(0.0*(M_PI/180.0));
+
+ // generate protons in B0 acceptance - roughly 5 - 20 mrad
+ double cos_theta_min = std::cos(0.005*(M_PI/180.0));
+ double cos_theta_max = std::cos(0.020*(M_PI/180.0));
+
+ double partEnergyMin = 270.0; // xL 0.98
+ double partEnergyMax = 275.0; // top beam energy
const double M_p = common_bench::particleMap.at(2212).mass;
@@ -39,12 +44,12 @@ void gen_forward_protons(int n_events = 100,
// pdgid 111 - pi0
// pdgid 2212 - proton
GenParticlePtr p1 =
- std::make_shared<GenParticle>(FourVector(0.0, 0.0, 10.0, 10.0), 2212, 4);
+ std::make_shared<GenParticle>(FourVector(0.0, 0.0, partEnergyMax, partEnergyMax), 2212, 4);
GenParticlePtr p2 = std::make_shared<GenParticle>(
FourVector(0.0, 0.0, 0.0, M_p), 2212, 4);
// Define momentum
- Double_t p = r1->Uniform(1.0, 10.0);
+ Double_t p = r1->Uniform(partEnergyMin, partEnergyMax);
Double_t phi = r1->Uniform(0.0, 2.0 * M_PI);
Double_t costh = r1->Uniform(cos_theta_min, cos_theta_max);
Double_t th = std::acos(costh);
--
GitLab