13 #include "Acts/Definitions/Units.hpp"
14 #include "Acts/MagneticField/MagneticFieldContext.hpp"
56 auto beamSpotData = beamSpotHandle.cptr();
65 Acts::MagneticFieldContext magFieldContext(fieldCondObj);
72 Acts::MagneticFieldProvider::Cache magFieldCache = magneticField.
makeCache( magFieldContext );
73 Acts::Vector3 bField = *magneticField.
getField( Acts::Vector3(beamPos.x(), beamPos.y(), 0),
90 [] (
const auto* seed) ->
double
91 {
return seed->sp()[0]->
x(); });
94 [] (
const auto* seed) ->
double
95 {
return seed->sp()[0]->
y(); });
98 [] (
const auto* seed) ->
double
99 {
return seed->sp()[0]->
z(); });
102 [] (
const auto* seed) ->
double
104 const auto* sp = seed->sp()[0];
105 return std::sqrt(sp->x()*sp->x() + sp->y()*sp->y());
111 [] (
const auto* seed) ->
double
112 {
return seed->sp()[1]->
x(); });
115 [] (
const auto* seed) ->
double
116 {
return seed->sp()[1]->
y(); });
119 [] (
const auto* seed) ->
double
120 {
return seed->sp()[1]->
z(); });
123 [] (
const auto* seed) ->
double
125 const auto* sp = seed->sp()[1];
126 return std::sqrt(sp->x()*sp->x() + sp->y()*sp->y());
132 [] (
const auto* seed) ->
double
133 {
return seed->sp()[2]->
x(); });
136 [] (
const auto* seed) ->
double
137 {
return seed->sp()[2]->
y(); });
140 [] (
const auto* seed) ->
double
141 {
return seed->sp()[2]->
z(); });
144 [] (
const auto* seed) ->
double
146 const auto* sp = seed->sp()[2];
147 return std::sqrt(sp->x()*sp->x() + sp->y()*sp->y());
150 std::vector< std::array<float, 7> > parametersCollection;
151 parametersCollection.reserve(seed_collection->
size());
153 for (
const auto* seed : *seed_collection) {
154 parametersCollection.push_back(
estimateParameters(*seed, 300. * bField[2] / 1000.) );
158 [] (
const auto&
params) ->
float
161 [] (
const auto&
params) ->
float
164 [] (
const auto&
params) ->
float
167 [] (
const auto&
params) ->
float
171 [] (
const auto&
params) ->
float
174 [] (
const auto&
params) ->
float
179 [] (
const auto&
params) ->
float
186 auto monitor_event_number =
Monitored::Scalar<long>(
"event_number",
static_cast<long>(eventInfo->eventNumber()));
189 std::vector<int> vec_truthBarcode;
190 std::vector<double> vec_truthProb;
197 monitor_x1, monitor_y1, monitor_z1, monitor_r1,
198 monitor_x2, monitor_y2, monitor_z2, monitor_r2,
199 monitor_x3, monitor_y3, monitor_z3, monitor_r3,
200 monitor_param_pt, monitor_param_theta, monitor_param_eta, monitor_param_d0,
201 monitor_param_dzdr_b, monitor_param_dzdr_t,
202 monitor_param_penalty,
203 monitor_event_number, monitor_actual_mu,
204 monitor_truth_barcode, monitor_truth_prob);
206 return StatusCode::SUCCESS;
211 std::vector<int>& truthBarCodeVec,
212 std::vector<double>& truthProbVec)
const
231 Acts::MagneticFieldContext magFieldContext(fieldCondObj);
237 auto retrieveSurfaceFunction =
238 [
this, &detEle] (
const ActsTrk::Seed& seed) ->
const Acts::Surface&
240 const auto& els = seed.sp().front()->measurements();
241 const auto* cluster = els[0];
249 std::vector<bool> vec_pass;
250 vec_pass.reserve(seed_container.
size());
252 std::vector<double> estimated_pt;
253 std::vector<double> estimated_eta;
254 estimated_pt.reserve(seed_container.
size());
255 estimated_eta.reserve(seed_container.
size());
257 for (
const auto* seed : seed_container) {
258 std::optional<Acts::BoundTrackParameters> optTrackParams =
261 geo_context.context(),
263 retrieveSurfaceFunction);
265 if ( not optTrackParams.has_value() )
continue;
267 const auto param = optTrackParams.value();
268 estimated_pt.push_back( param.transverseMomentum() );
269 estimated_eta.push_back( -
std::log(
std::tan(0.5 * param.parameters()[Acts::eBoundTheta]) ) );
271 std::map<int, int> truthHits;
273 const auto& sps = seed->sp();
274 for (
const auto* sp : sps) {
276 for (
int cluster_number(0); cluster_number < number_of_clusters; cluster_number++) {
277 const auto& els = sp->measurements();
278 const auto* cluster = els[cluster_number];
288 truthBarCodeVec.push_back(
barcode);
289 truthProbVec.push_back(
prob);
299 monitor_estimated_pt, monitor_estimated_eta);
301 return StatusCode::SUCCESS;
314 return (*pixelLink)->identify();
327 return (*stripLink)->identify();
333 std::map<int, int>& countMap)
const {
334 auto n1 = prdTruth->count(
id);
337 auto nBC = countMap.count(bc);
344 using iprdt = PRD_MultiTruthCollection::const_iterator;
345 std::pair<iprdt, iprdt>
range = prdTruth->equal_range(
id);
346 for (iprdt itr =
range.first; itr !=
range.second; ++itr) {
347 auto bc = itr->second.barcode();
348 auto nBC = countMap.count(bc);
365 if (
count > bestCount) {
373 double prob = bestCount / nTotal;
375 return std::make_pair(bestBarcode,
prob);
380 float pTPerHelixRadius)
const
382 auto extractCoordinates =
385 std::array<float, 4> coordinates {
static_cast<float>(sp->x()),
386 static_cast<float>(sp->y()),
387 static_cast<float>(sp->z()),
388 static_cast<float>(std::sqrt(sp->x()*sp->x() + sp->y()*sp->y()))};
392 auto extractQuantities =
393 [] (
const std::array<float, 4>& sp,
394 const std::array<float, 4>& spM,
395 bool isBottom) -> std::array<float, 5>
397 auto& [xM, yM, zM, rM] = spM;
398 auto& [xO, yO, zO, rO] = sp;
400 float cosPhiM = xM / rM;
401 float sinPhiM = yM / rM;
408 float iDeltaR = std::sqrt(iDeltaR2);
410 float cot_theta =
deltaZ * iDeltaR * bottomFactor;
413 std::array<float, 5>
params =
425 const auto& sps = seed.sp();
426 const auto* bottom = sps[0];
427 const auto* medium = sps[1];
428 const auto*
top = sps[2];
430 auto coo_b = extractCoordinates(bottom);
431 auto coo_m = extractCoordinates(medium);
432 auto coo_t = extractCoordinates(
top);
435 auto [cotThetaB, Zob, iDeltaRB, Ub, Vb] = extractQuantities(coo_b, coo_m,
true);
436 auto [cotThetaT, Zot, iDeltaRT, Ut, Vt] = extractQuantities(coo_t, coo_m,
false);
438 float squarediDeltaR2B = iDeltaRB*iDeltaRB;
439 float squarediDeltaR2T = iDeltaRB*iDeltaRT;
440 float squarediDeltaR =
std::min(squarediDeltaR2B, squarediDeltaR2T);
442 auto& [xB, yB, zB, rB] = coo_b;
443 auto& [xM, yM, zM, rM] = coo_m;
444 auto& [xT, yT, zT, rT] = coo_t;
452 float xb = dxb *
ax + dyb *
ay;
453 float yb = dyb *
ax - dxb *
ay;
454 float dxyb =
xb *
xb + yb * yb;
458 float xt = dxt *
ax + dyt *
ay;
459 float yt = dyt *
ax - dxt *
ay;
462 float tzb = dzb * std::sqrt( 1./dxyb );
463 float tzt = dzt * std::sqrt( 1./dxyt );
465 float sTzb2 = std::sqrt(1 + tzb*tzb);
469 return {-1, -1, -1, -1, -1, -1, -1};
472 float A = (Vt - Vb) / dU;
473 float S2 = 1. +
A *
A;
474 float B = Vb -
A * Ub;
478 float dzdr_b = (zM - zB) / (rM - rB);
479 float dzdr_t = (zT - zM) / (rT - rM);
482 float cotThetaAvg2 = cotThetaB * cotThetaT;
483 if (cotThetaAvg2 <= 0) {
484 return {-1, -1, -1, -1, -1, -1, -1};
486 float theta =
std::atan(1. / std::sqrt(cotThetaAvg2));
490 float pt = pTPerHelixRadius * std::sqrt(
S2 / B2) / 2.;
493 float d0 = std::abs((
A - B * rM) * rM);
498 float penalty = std::abs((tzb - tzt) / (squarediDeltaR * sTzb2));
500 return {
pt, theta, eta,
d0, dzdr_b, dzdr_t, penalty};