22 #include "GaudiKernel/ITHistSvc.h"
23 #include "GaudiKernel/MsgStream.h"
24 #include "GaudiKernel/SystemOfUnits.h"
33 m_highestVolume(nullptr),
34 m_extrapolator(
"Trk::Extrapolator/AtlasExtrapolator"),
37 m_materialCollectionValidation(false),
38 m_validationTree(nullptr),
39 m_validationRunTree(nullptr),
40 m_validationTreeFolder(
"/val/EventTreeTG"),
41 m_validationTreeName(
"EventTreeTG"),
42 m_validationTreeDescription(
"Event output of the EnergyLossExtrapolationValidation Algorithm"),
43 m_validationRunTreeFolder(
"/val/RunTreeTG"),
44 m_validationRunTreeName(
"RunTreeTG"),
45 m_validationRunTreeDescription(
"Run stats of the EnergyLossExtrapolationValidation Algorithm"),
55 m_totalRecordedLayers(0),
56 m_avgRecordedLayers(0.),
77 m_collectedLayerForward(0),
78 m_collectedLayerBack(0),
81 m_theCylinders(
nullptr),
87 declareProperty(
"UseMaterialCollection" , m_materialCollectionValidation);
91 declareProperty(
"ValidationTreeDescription" , m_validationTreeDescription);
92 declareProperty(
"ValidationRunTreeDescription" , m_validationRunTreeDescription);
118 delete m_theCylinders;
131 ATH_MSG_DEBUG(
"initialize() m_materialCollectionValidation = " << m_materialCollectionValidation );
134 if (m_extrapolator.retrieve().isFailure()) {
135 ATH_MSG_FATAL(
"initialize() Could not retrieve Tool " << m_extrapolator <<
". Exiting." );
136 return StatusCode::FAILURE;
140 m_validationTree =
new TTree(m_validationTreeName.c_str(), m_validationTreeDescription.c_str());
141 m_validationRunTree =
new TTree(m_validationRunTreeName.c_str(), m_validationRunTreeDescription.c_str());
144 m_validationTree->Branch(
"Entries", &m_entries,
"entries/i");
145 m_validationTree->Branch(
"Energy", m_energy,
"energy[entries]/F");
146 m_validationTree->Branch(
"EnergyLoss", m_energyLoss,
"eLoss[entries]/F");
147 m_validationTree->Branch(
"tInX0", m_parameterX0,
"tinX0[entries]/F");
148 m_validationTree->Branch(
"Radius", m_radius,
"radius[entries]/F");
149 m_validationTree->Branch(
"PosX", m_positionX,
"posX[entries]/F");
150 m_validationTree->Branch(
"PosY", m_positionY,
"posY[entries]/F");
151 m_validationTree->Branch(
"PosZ", m_positionZ,
"posZ[entries]/F");
152 m_validationTree->Branch(
"Eta", m_parameterEta,
"eta[entries]/F");
153 m_validationTree->Branch(
"Phi", m_parameterPhi,
"phi[entries]/F");
154 m_validationTree->Branch(
"Layer", m_layer,
"layer[entries]/i");
156 m_validationRunTree->Branch(
"Layers", &m_cylinders,
"layers/i");
157 m_validationRunTree->Branch(
"CylR", m_cylinderR,
"cylR[layers]/F");
158 m_validationRunTree->Branch(
"CylZ", m_cylinderZ,
"cylZ[layers]/F");
159 m_validationRunTree->Branch(
"Momentum",&m_momentum,
"momentum/F");
160 m_validationRunTree->Branch(
"UsePt", &m_usePt,
"usePt/O");
161 m_validationRunTree->Branch(
"MinEta", &m_minEta,
"minEta/F");
162 m_validationRunTree->Branch(
"MaxEta", &m_maxEta,
"maxEta/F");
163 m_validationRunTree->Branch(
"PDG", &m_pdg,
"pdg/I");
164 m_validationRunTree->Branch(
"Events", &m_events,
"events/i");
165 m_validationRunTree->Branch(
"AvgRecordedLayers", &m_avgRecordedLayers,
"avgRecLayers/F");
168 ITHistSvc* tHistSvc =
nullptr;
169 if (service(
"THistSvc",tHistSvc).isFailure()){
170 ATH_MSG_ERROR(
"initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
171 delete m_validationTree; m_validationTree =
nullptr;
172 delete m_validationRunTree; m_validationRunTree =
nullptr;
174 if ((tHistSvc->regTree(m_validationTreeFolder, m_validationTree)).isFailure()
175 || (tHistSvc->regTree(m_validationRunTreeFolder, m_validationRunTree)).isFailure() ) {
176 ATH_MSG_ERROR(
"initialize() Could not register the validation Trees -> Switching ValidationMode Off !" );
177 delete m_validationTree; m_validationTree =
nullptr;
178 delete m_validationRunTree; m_validationRunTree =
nullptr;
183 m_gaussDist =
new Rndm::Numbers(randSvc(), Rndm::Gauss(0.,1.));
184 m_flatDist =
new Rndm::Numbers(randSvc(), Rndm::Flat(0.,1.));
191 ATH_MSG_INFO(
"initialize() cylinder dimensions vector from jobOptions :" );
192 for (
size_t lay=0; lay<m_cylinders+1; ++lay) {
193 ATH_MSG_INFO(
"initialize() m_cylinderVR[" << lay <<
"] = " << m_cylinderVR[lay] <<
"\t ... m_cylinderVZ[" << lay <<
"] = " << m_cylinderVZ[lay] );
196 ATH_MSG_INFO(
"initialize() cylinder dimensions array for algorithm and ROOT tree :" );
197 for (
size_t lay=0; lay<m_cylinders+1; ++lay) {
198 m_cylinderR[lay] = m_cylinderVR[lay] > 0 ? m_cylinderVR[lay] : s_cylInitR[lay];
199 m_cylinderZ[lay] = m_cylinderVZ[lay] > 0 ? m_cylinderVZ[lay] : s_cylInitZ[lay];
205 if (m_onion && lay>0 && (m_cylinderR[lay] < m_cylinderR[lay-1])) {
206 ATH_MSG_WARNING(
"initialize() layer " << lay <<
"dimensions are smaller than those of layer " << lay-1 <<
" - constraining m_cylinders to " << lay-1 );
208 ATH_MSG_INFO(
"initialize() m_cylinderR[" << lay <<
"] = " << m_cylinderR[lay] <<
"\t ... m_cylinderZ[" << lay <<
"] = " << m_cylinderZ[lay] );
212 ATH_MSG_INFO(
"initialize() m_cylinderR[" << lay <<
"] = " << m_cylinderR[lay] <<
"\t ... m_cylinderZ[" << lay <<
"] = " << m_cylinderZ[lay] );
219 for (
size_t lay=0; lay<m_cylinders+1; ++lay) {
221 ATH_MSG_INFO(
"initialize() Cylinder " << lay <<
": " << *m_theCylinders->at(lay) );
222 m_theDiscs1->push_back(
new Trk::DiscSurface(*createTransform(0.,0.,-m_cylinderZ[lay]), 0., m_cylinderR[lay]));
223 ATH_MSG_INFO(
"initialize() Disc1 " << lay <<
": " << *m_theDiscs1->at(lay) );
224 m_theDiscs2->push_back(
new Trk::DiscSurface(*createTransform(0.,0., m_cylinderZ[lay]), 0., m_cylinderR[lay]));
225 ATH_MSG_INFO(
"initialize() Disc2 " << lay <<
": " << *m_theDiscs2->at(lay) );
228 if (m_particleType==0) m_pdg = 999;
229 else if (m_particleType==1) m_pdg = 11;
230 else if (m_particleType==2) m_pdg = 13;
231 else if (m_particleType==3) m_pdg = 211;
232 else if (m_particleType==4) m_pdg = 321;
233 ATH_MSG_INFO(
"initialize() ParticleType = " << m_particleType <<
" ... PDG = " << m_pdg );
236 return StatusCode::SUCCESS;
244 ATH_MSG_INFO(
"finalize() ================== Output Statistics =========================" );
246 ATH_MSG_INFO(
"finalize() = - breaks fwd : " <<
double(m_breaksForward)/
double(m_triesForward)
247 <<
" (" << m_breaksForward <<
"/" << m_triesForward <<
")" );
248 ATH_MSG_INFO(
"finalize() = - breaks bwd : " <<
double(m_breaksBack)/
double(m_triesBack)
249 <<
" (" << m_breaksBack <<
"/" << m_triesBack <<
")" );
250 if (m_materialCollectionValidation){
252 ATH_MSG_INFO(
"finalize() = - layer collected fwd : " << m_collectedLayerForward );
253 ATH_MSG_INFO(
"finalize() = - layer collected bwd : " << m_collectedLayerBack );
255 ATH_MSG_INFO(
"finalize() ==============================================================" );
257 m_avgRecordedLayers = m_events ? (
float)m_totalRecordedLayers / (
float)m_events : 0;
259 if (m_validationRunTree)
260 m_validationRunTree->Fill();
263 return StatusCode::SUCCESS;
270 const EventContext& ctx = Gaudi::Hive::currentContext();
272 if (!m_highestVolume){
280 ATH_MSG_WARNING(
"execute() No highest TrackingVolume / no VolumeBounds ... pretty useless!" );
281 return StatusCode::SUCCESS;
293 m_parameterP[
par] = 0.;
295 m_energyLoss[
par] = 0.;
296 m_parameterX0[
par] = 0.;
298 m_positionX[
par] = 0.;
299 m_positionY[
par] = 0.;
300 m_positionZ[
par] = 0.;
301 m_parameterEta[
par] = 0.;
302 m_parameterPhi[
par] = 0.;
303 m_parameterTheta[
par] = 0.;
310 m_parameterPhi[0] =
M_PI * (2 * m_flatDist->shoot() - 1);
311 m_parameterEta[0] = m_minEta + m_flatDist->shoot()*(m_maxEta-m_minEta);
312 m_parameterTheta[0] = 2.*
atan(
exp(-m_parameterEta[0]));
315 m_parameterP[0] = m_momentum;
318 m_parameterP[0] /=
sin(m_parameterTheta[0]);
319 m_parameterQoverP[0] =
charge/m_parameterP[0];
323 double energy1 = sqrt(m_parameterP[0]*m_parameterP[0] +
mass*
mass);
333 m_parameterQoverP[0],
337 ATH_MSG_DEBUG(
"execute() Start Parameters : [phi,eta] = [ " << startParameters.momentum().phi() <<
", " << startParameters.eta() <<
" ]" );
345 if (!m_materialCollectionValidation) {
347 lastParameters = m_extrapolator->extrapolate(
350 *(m_theCylinders->at(0)),
358 const std::vector<const Trk::TrackStateOnSurface*>* collectedMaterial =
359 m_extrapolator->extrapolateM(
362 *(m_theCylinders->at(0)),
368 if (collectedMaterial && !collectedMaterial->empty()) {
372 m_collectedLayerForward += collectedMaterial->size();
374 std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter = collectedMaterial->begin();
375 std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIterEnd = collectedMaterial->end();
376 for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
377 newX0 += (*tsosIter)->materialEffectsOnTrack() ? (*tsosIter)->materialEffectsOnTrack()->thicknessInX0() : 0;
385 for (
size_t lay = 1; lay<m_cylinders+1; ++lay) {
387 if (!m_onion) newX0 = 0;
391 if (!lastParameters) {
392 ATH_MSG_WARNING(
"execute() Layer " << lay <<
": start parameters for cylinder NOT found - skip event !" );
395 ATH_MSG_VERBOSE(
"execute() Layer " << lay <<
": start parameters for cylinder found: " << *lastParameters );
399 newParameters =
nullptr;
400 if (!m_materialCollectionValidation) {
402 newParameters = m_extrapolator->extrapolate(
404 m_onion ? *lastParameters : startParameters,
405 *(m_theCylinders->at(lay)),
413 const std::vector<const Trk::TrackStateOnSurface*>* collectedMaterial =
414 m_extrapolator->extrapolateM(ctx,
415 m_onion ? *lastParameters : startParameters,
416 *(m_theCylinders->at(lay)),
422 if (collectedMaterial && !collectedMaterial->empty()){
427 m_collectedLayerForward += collectedMaterial->size();
429 m_collectedLayerForward = collectedMaterial->size();
431 std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter = collectedMaterial->begin();
432 std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIterEnd = collectedMaterial->end();
433 for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
434 newX0 += (*tsosIter)->materialEffectsOnTrack() ? (*tsosIter)->materialEffectsOnTrack()->thicknessInX0() : 0;
442 if (!newParameters) {
444 if (!m_materialCollectionValidation) {
446 newParameters = m_extrapolator->extrapolate(
448 m_onion ? *lastParameters : startParameters,
449 (m_parameterEta[0] < 0) ? *(m_theDiscs1->at(lay))
450 : *(m_theDiscs2->at(lay)),
458 const std::vector<const Trk::TrackStateOnSurface*>* collectedMaterial =
459 m_extrapolator->extrapolateM(ctx,
460 m_onion ? *lastParameters : startParameters,
461 (m_parameterEta[0] < 0) ? *(m_theDiscs1->at(lay)) : *(m_theDiscs2->at(lay)),
467 if (collectedMaterial && !collectedMaterial->empty()){
472 m_collectedLayerForward += collectedMaterial->size();
474 m_collectedLayerForward = collectedMaterial->size();
476 std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIter = collectedMaterial->begin();
477 std::vector<const Trk::TrackStateOnSurface*>::const_iterator tsosIterEnd = collectedMaterial->end();
478 for ( ; tsosIter != tsosIterEnd; ++tsosIter) {
479 newX0 += (*tsosIter)->materialEffectsOnTrack() ? (*tsosIter)->materialEffectsOnTrack()->thicknessInX0() : 0;
488 if (!newParameters) {
489 ATH_MSG_WARNING(
"execute() Layer " << lay <<
" intersection did not work !" );
492 else if (m_highestVolume && newParameters && !(m_highestVolume->inside(newParameters->
position()))) {
493 ATH_MSG_WARNING(
"execute() Layer " << lay <<
" intersection is outside the known world !" );
500 ATH_MSG_VERBOSE(
"execute() Track Parameters at layer " << lay <<
": " << *newParameters );
501 ATH_MSG_DEBUG(
"execute() Track Parameters at layer " << lay <<
": [r,z] = [ " << newPosition.perp() <<
", " << newPosition.z() );
506 m_parameterPhi[m_entries] = newParameters->parameters()[
Trk::phi];
507 m_parameterEta[m_entries] = newParameters->
momentum().eta();
508 m_parameterTheta[m_entries] = newParameters->parameters()[
Trk::theta];
509 m_parameterP[m_entries] = newParameters->
momentum().mag();
510 m_parameterX0[m_entries] = (
float)newX0;
511 ATH_MSG_DEBUG(
"execute() Layer " << lay <<
": cumulated X0 = " << m_parameterX0[m_entries] );
514 m_energy[m_entries] = sqrt(m_parameterP[m_entries]*m_parameterP[m_entries] +
mass*
mass);
515 m_energyLoss[m_entries] = energy1-m_energy[m_entries];
516 ATH_MSG_DEBUG(
"execute() Layer " << lay <<
": cumulated Energy Loss = " << m_energyLoss[m_entries] );
519 m_layer[m_entries] = lay;
521 m_radius[m_entries] = newPosition.perp();
522 m_positionX[m_entries] = newPosition.x();
523 m_positionY[m_entries] = newPosition.y();
524 m_positionZ[m_entries] = newPosition.z();
528 lastParameters = newParameters;
533 m_totalRecordedLayers += m_entries;
539 if (m_validationTree)
540 m_validationTree->Fill();
543 ATH_MSG_DEBUG(
"execute() deleting DataVector parameters ... " );
545 return StatusCode::SUCCESS;
550 double x,
double y,
double z,
double phi,
double theta,
double alphaZ)
565 Amg::Vector3D surfaceXdirection(surfaceYdirection.cross(surfaceZdirection));
568 surfaceRotation.col(0) = surfaceXdirection;
569 surfaceRotation.col(1) = surfaceYdirection;
570 surfaceRotation.col(2) = surfaceZdirection;
573 return std::make_unique<Amg::Transform3D>(surfaceRotation, surfacePosition);