23 #include "GaudiKernel/ITHistSvc.h"
52 m_validationTree =
new TTree(m_validationTreeName.value().c_str(),
53 m_validationTreeDescription.value().c_str());
56 m_validationTree->Branch(
"Parameters", &m_parameters,
"params/I");
57 m_validationTree->Branch(
"ParametersLoc1", m_parameterLoc1,
"paramLoc1[params]/F");
58 m_validationTree->Branch(
"ParametersLoc2", m_parameterLoc2,
"paramLoc2[params]/F");
59 m_validationTree->Branch(
"ParametersPhi", m_parameterPhi,
"paramPhi[params]/F");
60 m_validationTree->Branch(
"ParametersTheta", m_parameterTheta,
"paramTheta[params]/F");
61 m_validationTree->Branch(
"ParametersEta", m_parameterEta,
"paramEta[params]/F");
62 m_validationTree->Branch(
"ParametersQoverP", m_parameterQoverP,
"paramQoverP[params]/F");
64 m_validationTree->Branch(
"CovarianceLoc1", m_covarianceLoc1,
"covLoc1[params]/F");
65 m_validationTree->Branch(
"CovarianceLoc2", m_covarianceLoc2,
"covLoc2[params]/F");
66 m_validationTree->Branch(
"CovariancePhi", m_covariancePhi,
"covPhi[params]/F");
67 m_validationTree->Branch(
"CovarianceTheta", m_covarianceTheta,
"covTheta[params]/F");
68 m_validationTree->Branch(
"CovarianceQoverP", m_covarianceQoverP,
"covQoverP[params]/F");
69 m_validationTree->Branch(
"CovarianceDeterminant", m_covarianceDeterminant,
"covDet[params]/F");
71 m_validationTree->Branch(
"StartMomentum", &m_startP,
"startP/F");
73 m_validationTree->Branch(
"StartSurfaceX", &m_startX,
"startX/F");
74 m_validationTree->Branch(
"StartSurfaceY", &m_startY,
"startY/F");
75 m_validationTree->Branch(
"StartSurfaceR", &m_startR,
"startR/F");
76 m_validationTree->Branch(
"StartSurfaceZ", &m_startZ,
"startZ/F");
78 m_validationTree->Branch(
"EstimationSurfaceX", &m_estimationX,
"estimateX/F");
79 m_validationTree->Branch(
"EstimationSurfaceY", &m_estimationY,
"estimateY/F");
80 m_validationTree->Branch(
"EstimationSurfaceR", &m_estimationR,
"estimateR/F");
81 m_validationTree->Branch(
"EstimationSurfaceZ", &m_estimationZ,
"estimateZ/F");
83 m_validationTree->Branch(
"DestinationSurfaceType", &m_destinationSurfaceType,
"surfaceType/I");
84 m_validationTree->Branch(
"DestinationSurfaceX", &m_destinationX,
"surfaceX/F");
85 m_validationTree->Branch(
"DestinationSurfaceY", &m_destinationY,
"surfaceY/F");
86 m_validationTree->Branch(
"DestinationSurfaceR", &m_destinationR,
"surfaceR/F");
87 m_validationTree->Branch(
"DestinationSurfaceZ", &m_destinationZ,
"surfaceZ/F");
90 SmartIF<ITHistSvc> tHistSvc{service(
"THistSvc")};
92 ATH_MSG_ERROR(
"initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
93 delete m_validationTree; m_validationTree =
nullptr;
95 if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()) {
96 ATH_MSG_ERROR(
"initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
97 delete m_validationTree; m_validationTree =
nullptr;
101 m_gaussDist =
new Rndm::Numbers(randSvc(), Rndm::Gauss(0.,1.));
102 m_flatDist =
new Rndm::Numbers(randSvc(), Rndm::Flat(0.,1.));
105 return StatusCode::SUCCESS;
113 ATH_MSG_INFO(
"================== Output Statistics =========================");
115 ATH_MSG_INFO(
"= - breaks fwd : " <<
static_cast<double>(m_breaksFront)/
static_cast<double>(m_triesFront)
116 <<
" (" << m_breaksFront <<
"/" << m_triesFront <<
")");
117 ATH_MSG_INFO(
"= - breaks bwd : " <<
static_cast<double>(m_breaksBack)/
static_cast<double>(m_triesBack)
118 <<
" (" << m_breaksBack <<
"/" << m_triesBack <<
")");
119 if (m_materialCollectionValidation){
121 ATH_MSG_INFO(
"= - layer collected fwd : " << m_collectedLayerFront );
122 ATH_MSG_INFO(
"= - layer collected bwd : " << m_collectedLayerBack );
125 ATH_MSG_INFO(
"==============================================================");
127 return StatusCode::SUCCESS;
134 const EventContext& ctx = Gaudi::Hive::currentContext();
136 if (!m_highestVolume){
144 ATH_MSG_WARNING(
"No highest TrackingVolume / no VolumeBounds ... pretty useless! ");
145 return StatusCode::SUCCESS;
154 m_destinationSurfaceType = 0;
172 m_parameterLoc1[m_parameters] = m_sigmaLoc * m_gaussDist->shoot();
173 m_parameterLoc2[m_parameters] = m_sigmaLoc * m_gaussDist->shoot();
175 m_parameterPhi[m_parameters] =
M_PI * m_flatDist->shoot();
176 m_parameterPhi[m_parameters] *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
177 m_parameterEta[m_parameters] = m_minEta + m_flatDist->shoot()*(m_maxEta-m_minEta);
178 m_parameterTheta[m_parameters] = 2.*
std::atan(
std::exp(-m_parameterEta[m_parameters]));
180 m_covarianceLoc1[m_parameters] = std::abs(m_parameterLoc1[m_parameters] * 0.1);
181 m_covarianceLoc2[m_parameters] = std::abs(m_parameterLoc2[m_parameters] * 0.1);
182 m_covariancePhi[m_parameters] = std::abs(m_parameterPhi[m_parameters] * 0.1);
183 m_covarianceTheta[m_parameters] = std::abs(m_parameterTheta[m_parameters] * 0.1);
186 double p = m_minP + m_flatDist->shoot()*(m_maxP-m_minP);
187 double charge = (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
188 m_parameterQoverP[m_parameters] =
charge/
p;
190 m_covarianceQoverP[m_parameters] = std::abs(m_parameterQoverP[m_parameters] * 0.1);
196 m_startR = std::abs(m_sigmaR * m_gaussDist->shoot());
197 double surfacePhi =
M_PI * m_flatDist->shoot();
198 surfacePhi *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
199 m_startX = m_startR*
cos(surfacePhi);
200 m_startY = m_startR*
sin(surfacePhi);
201 m_startZ = m_sigmaZ * m_gaussDist->shoot();
204 double alphaZ =
M_PI * m_flatDist->shoot();
205 alphaZ *= (m_flatDist->shoot() > 0.5 ) ? -1. : 1.;
211 m_parameterPhi[m_parameters],
212 m_parameterTheta[m_parameters],
219 covariance.setZero();
221 covariance(0,0) = m_covarianceLoc1[m_parameters];
223 covariance(1,1) = m_covarianceLoc2[m_parameters];
225 covariance(2,2) = m_covariancePhi[m_parameters];
227 covariance(3,3) = m_covarianceTheta[m_parameters];
229 covariance(4,4) = m_covarianceQoverP[m_parameters];
233 m_covarianceDeterminant[m_parameters] = covariance.determinant();
237 m_parameterLoc2[m_parameters],
238 m_parameterPhi[m_parameters],
239 m_parameterTheta[m_parameters],
240 m_parameterQoverP[m_parameters],
242 std::move(covariance));
245 if(startParameters.covariance())
ATH_MSG_VERBOSE(
"Start Covariance : \n" << *startParameters.covariance() );
249 m_estimationR = m_maximumR * m_flatDist->shoot();
253 CylTrf.setIdentity();
260 if (!estimationParameters) {
261 ATH_MSG_VERBOSE(
"Estimation of intersection did not work - skip event !" );
262 return StatusCode::SUCCESS;
264 else if (m_highestVolume && estimationParameters && !(m_highestVolume->inside(estimationParameters->
position()))){
265 ATH_MSG_VERBOSE(
"Estimation of intersection is outside the known world - skip event !" );
266 delete estimationParameters;
267 return StatusCode::SUCCESS;
275 m_parameterLoc1[m_parameters] = estimationParameters->parameters()[
Trk::loc1];
276 m_parameterLoc2[m_parameters] = estimationParameters->parameters()[
Trk::loc2];
277 m_parameterPhi[m_parameters] = estimationParameters->parameters()[
Trk::phi];
278 m_parameterEta[m_parameters] = estimationParameters->
eta();
279 m_parameterTheta[m_parameters] = estimationParameters->parameters()[
Trk::theta];
280 m_parameterQoverP[m_parameters] = estimationParameters->parameters()[
Trk::qOverP];
281 if(estimationParameters->covariance()){
282 m_covarianceLoc1[m_parameters] = (*estimationParameters->covariance())(0,0);
283 m_covarianceLoc2[m_parameters] = (*estimationParameters->covariance())(1,1);
284 m_covariancePhi[m_parameters] = (*estimationParameters->covariance())(2,2);
285 m_covarianceTheta[m_parameters] = (*estimationParameters->covariance())(3,3);
286 m_covarianceQoverP[m_parameters] = (*estimationParameters->covariance())(4,4);
287 m_covarianceDeterminant[m_parameters] = (estimationParameters->covariance())->determinant();
290 m_covarianceLoc1[m_parameters] = 0;
291 m_covarianceLoc2[m_parameters] = 0;
292 m_covariancePhi[m_parameters] = 0;
293 m_covarianceTheta[m_parameters] = 0;
294 m_covarianceQoverP[m_parameters] = 0;
295 m_covarianceDeterminant[m_parameters] = 0;
302 m_estimationX = estimatedPosition.x();
303 m_estimationY = estimatedPosition.y();
304 m_estimationZ = estimatedPosition.z();
307 delete estimationParameters; estimationParameters =
nullptr;
313 m_parameterPhi[m_parameters],
314 m_parameterTheta[m_parameters]), 10e5 , 10e5);
317 ATH_MSG_VERBOSE(
"Extrapolation to Destination Surface: " << destinationSurface );
322 if (!m_materialCollectionValidation && !m_direct)
323 destParameters = m_extrapolator->extrapolate(ctx,
332 const std::vector<const Trk::TrackStateOnSurface*>*
333 collectedMaterial = m_extrapolator->extrapolateM(ctx,
341 if (collectedMaterial && !collectedMaterial->empty()){
345 m_collectedLayerFront += collectedMaterial->size();
347 for (
const auto* tsos : *collectedMaterial) {
354 destParameters = m_extrapolator->extrapolateDirectly(ctx,
368 m_parameterLoc1[m_parameters] = destParameters->parameters()[
Trk::loc1];
369 m_parameterLoc2[m_parameters] = destParameters->parameters()[
Trk::loc2];
370 m_parameterPhi[m_parameters] = destParameters->parameters()[
Trk::phi];
371 m_parameterEta[m_parameters] = destParameters->
eta();
372 m_parameterTheta[m_parameters] = destParameters->parameters()[
Trk::theta];
373 m_parameterQoverP[m_parameters] = destParameters->parameters()[
Trk::qOverP];
374 if(destParameters->covariance()){
375 m_covarianceLoc1[m_parameters] = (*destParameters->covariance())(0,0);
376 m_covarianceLoc2[m_parameters] = (*destParameters->covariance())(1,1);
377 m_covariancePhi[m_parameters] = (*destParameters->covariance())(2,2);
378 m_covarianceTheta[m_parameters] = (*destParameters->covariance())(3,3);
379 m_covarianceQoverP[m_parameters] = (*destParameters->covariance())(4,4);
380 m_covarianceDeterminant[m_parameters] = (destParameters->covariance())->determinant();
383 m_covarianceLoc1[m_parameters] = 0;
384 m_covarianceLoc2[m_parameters] = 0;
385 m_covariancePhi[m_parameters] = 0;
386 m_covarianceTheta[m_parameters] = 0;
387 m_covarianceQoverP[m_parameters] = 0;
388 m_covarianceDeterminant[m_parameters] = 0;
392 m_destinationX = destinationPosition.x();
393 m_destinationY = destinationPosition.y();
394 m_destinationZ = destinationPosition.z();
395 m_destinationR = destinationPosition.perp();
400 if (!m_materialCollectionValidation && !m_direct)
401 backParameters = m_extrapolator->extrapolate(ctx,
410 const std::vector<const Trk::TrackStateOnSurface*>*
411 collectedBackMaterial = m_extrapolator->extrapolateM(ctx,
418 if (collectedBackMaterial && !collectedBackMaterial->empty()){
423 m_collectedLayerBack += collectedBackMaterial->size();
425 for (
const auto* tsos : *collectedBackMaterial) {
432 backParameters = m_extrapolator->extrapolateDirectly(ctx,
447 m_parameterLoc1[m_parameters] = backParameters->parameters()[
Trk::loc1];
448 m_parameterLoc2[m_parameters] = backParameters->parameters()[
Trk::loc2];
449 m_parameterPhi[m_parameters] = backParameters->parameters()[
Trk::phi];
450 m_parameterEta[m_parameters] = backParameters->
eta();
451 m_parameterTheta[m_parameters] = backParameters->parameters()[
Trk::theta];
452 m_parameterQoverP[m_parameters] = backParameters->parameters()[
Trk::qOverP];
453 if(backParameters->covariance()){
454 m_covarianceLoc1[m_parameters] = (*backParameters->covariance())(0,0);
455 m_covarianceLoc2[m_parameters] = (*backParameters->covariance())(1,1);
456 m_covariancePhi[m_parameters] = (*backParameters->covariance())(2,2);
457 m_covarianceTheta[m_parameters] = (*backParameters->covariance())(3,3);
458 m_covarianceQoverP[m_parameters] = (*backParameters->covariance())(4,4);
459 m_covarianceDeterminant[m_parameters] = (backParameters->covariance())->determinant();
462 m_covarianceLoc1[m_parameters] = 0;
463 m_covarianceLoc2[m_parameters] = 0;
464 m_covariancePhi[m_parameters] = 0;
465 m_covarianceTheta[m_parameters] = 0;
466 m_covarianceQoverP[m_parameters] = 0;
467 m_covarianceDeterminant[m_parameters] = 0;
470 delete backParameters;
474 delete destParameters;
481 if (m_validationTree)
482 m_validationTree->Fill();
488 return StatusCode::SUCCESS;
508 Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
510 double nx = 1./sqrt(surfaceXdirection[0]*surfaceXdirection[0]+surfaceXdirection[1]*surfaceXdirection[1]+surfaceXdirection[2]*surfaceXdirection[2]);
511 double ny = 1./sqrt(surfaceYdirection[0]*surfaceYdirection[0]+surfaceYdirection[1]*surfaceYdirection[1]+surfaceYdirection[2]*surfaceYdirection[2]);
512 surfaceXdirection[0]*=nx;
513 surfaceXdirection[1]*=nx;
514 surfaceXdirection[2]*=nx;
516 surfaceYdirection[0]*=ny;
517 surfaceYdirection[1]*=ny;
518 surfaceYdirection[2]*=ny;
521 surfaceRotation.col(0) = surfaceXdirection;
522 surfaceRotation.col(1) = surfaceYdirection;
523 surfaceRotation.col(2) = surfaceZdirection;