diff --git a/benchmarks/b0_tracker/forward_protons.sh b/benchmarks/b0_tracker/forward_protons.sh
index d78b81b7f0820cce03e96ab1e555209e6e68eb74..e5d3572bc15f076118ecfbbe671e079ba426370b 100755
--- a/benchmarks/b0_tracker/forward_protons.sh
+++ b/benchmarks/b0_tracker/forward_protons.sh
@@ -7,7 +7,7 @@ if [[ ! -n "${JUGGLER_DETECTOR}" ]] ; then
fi
if [[ ! -n "${JUGGLER_N_EVENTS}" ]] ; then
- export JUGGLER_N_EVENTS=100
+ export JUGGLER_N_EVENTS=1000
fi
export FILE_NAME_TAG="forward_protons"
diff --git a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
index b4984adaa3ef6c31788a38c6ba35f3d0f3e9e466..c6dd692cb973c2be830fa5f42b8b05ea1979d389 100644
--- a/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
+++ b/benchmarks/b0_tracker/scripts/b0_tracker_hits.cxx
@@ -36,25 +36,67 @@ void b0_tracker_hits(const char* fname = "./sim_output/sim_forward_protons.root"
ROOT::RDataFrame d0(*t);
- auto hits_eta = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
+ auto hits_theta = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
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 local_position = [&](const std::vector<dd4pod::TrackerHitData>& hits) {
+ std::vector<std::array<double, 2>> result;
+ for (const auto& h : hits) {
+
+ auto pos0 = (h.position);
+ result.push_back({pos0.x , pos0.y});
+ }
+ return result;
+ };
+
+
+ auto x_pos = [&](const std::vector<std::array<double, 2>>& xypos) {
+ std::vector<double> result;
+ for (const auto& h : xypos) {
+
+ result.push_back(h.at(0));
+ }
+ return result;
+ };
- auto h1 = d1.Histo1D({"h1", "hits_eta", 100, 0,20}, "hits_eta");
+ auto y_pos = [&](const std::vector<std::array<double, 2>>& xypos) {
+ std::vector<double> result;
+ for (const auto& h : xypos) {
+
+ result.push_back(h.at(1));
+ }
+ return result;
+ };
+
+
+ auto d1 = d0.Define("nhits", hits_theta, {"B0TrackerHits"})
+ .Define("xy_hit_pos", local_position, {"B0TrackerHits"})
+ .Define("x_pos", x_pos, {"xy_hit_pos"})
+ .Define("y_pos", y_pos, {"xy_hit_pos"});
+
+ auto h_local_pos = d1.Histo2D({"h_local_pos", ";x [mm]; y [mm] ", 100, -100.0, -200.0, 100, -100.0, 100.0}, "x_pos", "y_pos");
+
+ auto d2 = d0.Define("hits_theta", hits_theta, {"B0TrackerHits"});
+
+ auto h1 = d2.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");
+
+ h_local_pos->DrawCopy("colz");
+ c->SaveAs("results/b0_tracker_hits_occupancy_disk_1.png");
+ c->SaveAs("results/b0_tracker_hits_occupancy_disk_1.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 6cb9bd50e89e21bfc9e13a12f449c6d5fc8886bc..9e1a3405a700f8ce90c3f18cf92da1af5512cff0 100644
--- a/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
+++ b/benchmarks/b0_tracker/scripts/gen_forward_protons.cxx
@@ -20,12 +20,19 @@ using namespace HepMC3;
/** Generate electrons in the central region.
* This is for testing detectors in the "barrel" region.
*/
-void gen_forward_protons(int n_events = 100,
+void gen_forward_protons(int n_events = 1000,
const char* out_fname = "forward_protons.hepmc")
{
- double cos_theta_min = std::cos(0.5*(M_PI/180.0));
+
+ double crossingAngle = -0.025; //radiansc
+
+ // generate protons in B0 acceptance - roughly 5 - 20 mrad
+ double cos_theta_min = std::cos(0.02);
double cos_theta_max = std::cos(0.0*(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;
WriterAscii hepmc_output(out_fname);
@@ -42,12 +49,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(200.0, 275.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);
@@ -55,20 +62,12 @@ void gen_forward_protons(int n_events = 100,
Double_t py = p * std::sin(phi) * std::sin(th);
Double_t pz = p * std::cos(th);
-
ROOT::Math::XYZVector p0 = {px,py,pz};
//ROOT::Math::Rotation3D r = (-0.025);
ROOT::Math::RotationY r(-0.025);
auto p_rot = r*p0;
-
- // Generates random vectors, uniformly distributed over the surface of a
- // sphere of given radius, in this case momentum.
- // r1->Sphere(px, py, pz, p);
-
- //std::cout << std::sqrt(px*px + py*py + pz*pz) - p << " is zero? \n";
-
// type 1 is final state
// pdgid 11 - electron 0.510 MeV/c^2
GenParticlePtr p3 = std::make_shared<GenParticle>(