ATLAS Offline Software
MaterialMapper.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 // MaterialMapper.cxx, (c) ATLAS Detector software
8 
9 // Trk include
10 #include "TrkDetDescrTestTools/MaterialMapper.h"
11 #include "TrkDetDescrUtils/MaterialAssociationType.h"
12 #include "TrkGeometry/AssociatedMaterial.h"
13 #include "TrkGeometry/TrackingVolume.h"
14 #include "TrkGeometry/CylinderLayer.h"
15 #include "TrkGeometry/DiscLayer.h"
16 #include "TrkGeometry/MaterialLayer.h"
17 // Gaudi
18 #include "GaudiKernel/ITHistSvc.h"
19 
20 
21 // constructor
22 Trk::MaterialMapper::MaterialMapper(const std::string& t, const std::string& n, const IInterface* p)
23 : AthAlgTool(t,n,p),
24  m_incidentSvc("IncidentSvc", n),
25  m_materialAssociationType(1),
26  m_maxMappingEvents(100000),
27  m_processedEvents(0),
28  m_totalMaterialTree("TotalMaterial"),
29  m_validationTree(nullptr),
30  m_validationTreeName("MaterialMapper"),
31  m_validationTreeDescription("Material Effects Updator information"),
32  m_validationTreeFolder("/val/MaterialMapper"),
33  m_materialSteps(0),
34  m_averageEta{}, m_averagePhi{}, m_mappedPath{}, m_mappedPathInX0{},
35  m_mappedPathInL0{}, m_mappedPathRho{}, m_mappedPathZARho{},
36  m_unmappedPathInX0{},
37  m_volumeValidation(true),
38  m_volumeTreePrefix("VolumeMaterial"),
39  m_layerValidation(true),
40  m_layerTreePrefix("LayerMaterial_"),
41  m_surfaceValidation(true),
42  m_surfaceTreePrefix("SurfaceMaterial_")
43 {
44  declareInterface<IMaterialMapper>(this);
45 
46  // The maximum mapping events to control the output file size
47  declareProperty("MaximumMappingEvents", m_maxMappingEvents);
48  // specification for the Validation Tree
49  declareProperty("ValidationTreeName", m_validationTreeName);
50  declareProperty("ValidationTreeDescription", m_validationTreeDescription);
51  declareProperty("ValidationTreeFolder", m_validationTreeFolder);
52  // the total material statistics
53  declareProperty("TotalMaterialTreeName", m_totalMaterialTree);
54  // declare the association type
55  declareProperty("MaterialAssociationType", m_materialAssociationType);
56  // the volume validation
57  declareProperty("TrackingVolumeValidation", m_volumeValidation);
58  declareProperty("TrackingVolumeTreePrefix", m_volumeTreePrefix);
59  // the layer Validation
60  declareProperty("DetailedLayerValidation", m_layerValidation);
61  declareProperty("DetailedLayerTreePrefix", m_layerTreePrefix);
62  // the surface Validation
63  declareProperty("DetailedSurfaceValidation", m_surfaceValidation);
64  declareProperty("DetailedSurfaceTreePrefix", m_surfaceTreePrefix);
65 
66 }
67 
68 // destructor
69 Trk::MaterialMapper::~MaterialMapper()
70 {
71  // memory cleanup
72  if (m_volumeValidation) {
73  auto volTreeIter = m_volumeTrees.begin();
74  auto endVolTreeIter = m_volumeTrees.end();
75  for ( ; volTreeIter != endVolTreeIter; delete (volTreeIter->second), ++volTreeIter) {}
76  }
77 
78  if (m_layerValidation) {
79  // delete the layer hit TTree
80  auto layTreeIter = m_layerTrees.begin();
81  auto endLayTreeIter = m_layerTrees.end();
82  for ( ; layTreeIter != endLayTreeIter; delete (layTreeIter->second), ++layTreeIter) {}
83  // delete the layer per event TTree
84  layTreeIter = m_layerFullTrees.begin();
85  endLayTreeIter = m_layerFullTrees.end();
86  for ( ; layTreeIter != endLayTreeIter; delete (layTreeIter->second), ++layTreeIter) {}
87  }
88 
89  if (m_surfaceValidation) {
90  // finalize the Hit Tree
91  auto surfTreeIter = m_surfaceTrees.begin();
92  auto endSurfTreeIter = m_surfaceTrees.end();
93  for ( ; surfTreeIter != endSurfTreeIter; delete (surfTreeIter->second), ++surfTreeIter) {}
94  }
95 }
96 
97 // Athena standard methods
98 // initialize
99 StatusCode Trk::MaterialMapper::initialize()
100 {
101 
102  // Athena/Gaudi framework
103  if (m_incidentSvc.retrieve().isFailure()){
104  ATH_MSG_WARNING("Could not retrieve " << m_incidentSvc << ". Exiting.");
105  return StatusCode::FAILURE;
106  }
107  // register to the incident service: EndEvent for histogram filling and reset
108  m_incidentSvc->addListener( this, IncidentType::EndEvent);
109 
110  // book the Tree
111  bookValidationTree();
112  return StatusCode::SUCCESS;
113 }
114 
115 // finalize
116 StatusCode Trk::MaterialMapper::finalize()
117 {
118  ATH_MSG_INFO( "finalize() successful" );
119  return StatusCode::SUCCESS;
120 }
121 
122 void Trk::MaterialMapper::recordMaterialHit(const AssociatedMaterial& amhit, const Amg::Vector3D& projectedPosition) const
123 {
124 
125  if (m_processedEvents >= m_maxMappingEvents) return;
126 
127  // ----------------------------------- Recording Hits section --------------------------------
128  const Amg::Vector3D& position = amhit.materialPosition();
129 
130  // get the tracking volume
131  const Trk::TrackingVolume* tvol = amhit.associatedTrackingVolume();
132  const Trk::Layer* lay = amhit.associatedLayer();
133 
134  if (!tvol) {
135  ATH_MSG_WARNING("Could not associate TrackingVolume to material step!");
136  return;
137  }
138 
139  // ----------------------------------- The Material Validation code --------------------------
140  const Trk::TrackingVolume* itvol = tvol;
141  while (itvol) {
142  // get the volume tree object
143  ATH_MSG_VERBOSE( "getting volTreeObj for: " << itvol->volumeName() << " at: " << itvol );
144  Trk::VolumeTreeObject* volTreeObj = volumeTreeObject(lay,itvol);
145  if (volTreeObj) {
146  // add the quantities
147  (*volTreeObj).path += amhit.steplength();
148  (*volTreeObj).pathInX0 += amhit.steplengthInX0();
149  (*volTreeObj).pathInL0 += amhit.steplengthInL0();
150  if(amhit.A() != 0.0) {
151  (*volTreeObj).pathZARho += (amhit.Z()/amhit.A())*amhit.rho()*amhit.steplength();
152  }
153  } else
154  ATH_MSG_WARNING( "Cannot find/create VolumeTreeObject for volume '" << itvol->volumeName() << "'." );
155  itvol = itvol->getMotherVolume();
156  if (itvol) ATH_MSG_VERBOSE( "motherVolume found: " << itvol->volumeName() );
157  }
158 
159  // -------------------------------- main tree section -----------------------------------------
160  // eta / phi
161  m_averageEta += position.eta();
162  m_averagePhi += position.phi();
163  // effectively crossed Z/A*rho
164  double steplength = amhit.steplength();
165  // fill the variables for spatial information & event characteristics
166  if (lay) {
167  m_mappedPath += steplength;
168  m_mappedPathInX0 += amhit.steplengthInX0();
169  m_mappedPathInL0 += amhit.steplengthInL0();
170  m_mappedPathRho += amhit.rho()*steplength;
171  if(amhit.A() != 0.0) {
172  m_mappedPathZARho += (amhit.Z()/amhit.A())*amhit.rho()*steplength;
173  }
174  ATH_MSG_VERBOSE("[ MaterialMapper ] Accumulated Path in Rho (*Z/A) = " << m_mappedPathRho << " ( " << m_mappedPathZARho << " )");
175  } else
176  m_unmappedPathInX0 += amhit.steplengthInX0();
177 
178  // mapping information
179  m_mapped[m_materialSteps] = lay ? lay->layerIndex().value() : 0;
180  // the position & step information
181  m_materialAccumPathInX0[m_materialSteps] = m_mappedPathInX0;
182  m_materialAccumPathZARho[m_materialSteps] = m_mappedPathZARho;
183  m_materialStepPath[m_materialSteps] = steplength;
184  m_materialStepX0[m_materialSteps] = amhit.x0();
185  m_materialStepL0[m_materialSteps] = amhit.l0();
186  m_materialStepA[m_materialSteps] = amhit.A();
187  m_materialStepZ[m_materialSteps] = amhit.Z();
188  m_materialStepRho[m_materialSteps] = amhit.rho();
189  m_materialStepPositionX[m_materialSteps] = position.x();
190  m_materialStepPositionY[m_materialSteps] = position.y();
191  m_materialStepPositionZ[m_materialSteps] = position.z();
192  m_materialStepPositionR[m_materialSteps] = position.perp();
193  // record the projected position as well
194  m_materialProjPositionX[m_materialSteps] = projectedPosition.x();
195  m_materialProjPositionY[m_materialSteps] = projectedPosition.y();
196  m_materialProjPositionZ[m_materialSteps] = projectedPosition.z();
197  m_materialProjPositionR[m_materialSteps] = projectedPosition.perp();
198  m_materialProjDistance[m_materialSteps] = (position-projectedPosition).mag();
199  // and increase the number of steps
200  ++m_materialSteps;
201 
202 
203  // record the layer hit
204  if (m_layerValidation) recordLayerHit(amhit);
205 }
206 
207 
208 void Trk::MaterialMapper::finalizeVolumeHits(bool mapped) const
209 {
210  ATH_MSG_VERBOSE( "finalizeEvent() on Volume hits " );
211 
212  auto volTreeIter = mapped ? m_volumeTrees.begin() : m_volumeTreesUnmapped.begin();
213  auto endVolTreeIter = mapped ? m_volumeTrees.end() : m_volumeTreesUnmapped.end();
214 
215  // do the loop
216  for ( ; volTreeIter != endVolTreeIter; ++volTreeIter ) {
217 
218  Trk::VolumeTreeObject* volTreeObj = volTreeIter->second;
219 
220  if (!volTreeObj) continue;
221 
222  if ((*volTreeObj).path > 0.) {
223  // spatial mapping
224  (*volTreeObj).eta = m_averageEta;
225  (*volTreeObj).phi = m_averagePhi;
226  // fill the volume Tree Object
227  (*volTreeObj).tree->Fill();
228  }
229  // reset
230  (*volTreeObj).eta = 0.;
231  (*volTreeObj).phi = 0.;
232  (*volTreeObj).path = 0.;
233  (*volTreeObj).pathInX0 = 0.;
234  (*volTreeObj).pathInL0 = 0.;
235  (*volTreeObj).pathZARho = 0.;
236 
237  }
238 }
239 
240 
241 void Trk::MaterialMapper::recordLayerHit(const AssociatedMaterial& amhit, bool full) const
242 {
243 
244  if (m_processedEvents >= m_maxMappingEvents) return;
245 
246  const Trk::Layer* lay = amhit.associatedLayer();
247  if (m_layerValidation && lay) {
248  // try to find the histogram
249  Trk::LayerTreeObject* layTreeObj = layerTreeObject(*lay,full);
250  if (layTreeObj
251  && (*layTreeObj).layerHits < TRKDETDESCRTOOLS_MAXLAYERHITS
252  && amhit.steplengthInX0() > 0.)
253  {
254 
255  const Amg::Vector3D& pos = amhit.materialPosition();
256 
257  // hit positions
258  double posX = pos.x();
259  double posY = pos.y();
260  double posZ = pos.z();
261  double posR = pos.perp();
262  double posEta = pos.eta();
263 
264  // fill the hit position --------------------
265  (*layTreeObj).hitPositionX[(*layTreeObj).layerHits] = posX;
266  (*layTreeObj).hitPositionY[(*layTreeObj).layerHits] = posY;
267  (*layTreeObj).hitPositionZ[(*layTreeObj).layerHits] = posZ;
268  (*layTreeObj).hitPositionR[(*layTreeObj).layerHits] = posR;
269  (*layTreeObj).hitPositionEta[(*layTreeObj).layerHits] = posEta;
270 
271  // increase the layer Hits
272  ++(*layTreeObj).layerHits;
273  ++(*layTreeObj).densedHits;
274 
275  // fill the densed Hit position ---------------
276  (*layTreeObj).densedHitX += posX;
277  (*layTreeObj).densedHitY += posY;
278  (*layTreeObj).densedHitZ += posZ;
279  (*layTreeObj).densedHitR += posR;
280  (*layTreeObj).densedHitPhi += pos.phi();
281  (*layTreeObj).densedHitTheta += pos.theta();
282 
283  // fill the correction factor ------------------
284  (*layTreeObj).correctionFactor += amhit.correctionFactor();
285  (*layTreeObj).path += amhit.steplength();
286  (*layTreeObj).pathInX0 += amhit.steplengthInX0();
287  (*layTreeObj).pathInL0 += amhit.steplengthInL0();
288  (*layTreeObj).A += amhit.A() * amhit.rho() * amhit.steplength();
289  (*layTreeObj).Z += amhit.Z() * amhit.rho() * amhit.steplength();
290  (*layTreeObj).rho += amhit.rho() * amhit.steplength();
291 
292  }
293  }
294 }
295 
296 void Trk::MaterialMapper::finalizeLayerHits() const
297 {
298 
299  ATH_MSG_VERBOSE( "finalizeEvent() on Layer hits " );
300 
301  for (size_t iltm = 0; iltm < 2; ++iltm) {
302  // loop over all layer Tree objects and fill them
303  auto layObjIter = iltm ? m_layerFullTrees.begin() : m_layerTrees.begin();
304  auto endIter = iltm ? m_layerFullTrees.end() : m_layerTrees.end();
305  // do the loop
306  for ( ; layObjIter != endIter; ++layObjIter ) {
307 
308  Trk::LayerTreeObject* layTreeObj = layObjIter->second;
309 
310  if (!layTreeObj) continue;
311 
312  // get the number of single hits
313  int hits = (*layTreeObj).densedHits;
314 
315  if ((*layTreeObj).pathInX0 > 0.) {
316  if (hits>1) {
317  (*layTreeObj).densedHitX /= hits;
318  (*layTreeObj).densedHitY /= hits;
319  (*layTreeObj).densedHitZ /= hits;
320  (*layTreeObj).densedHitR /= hits;
321  (*layTreeObj).densedHitPhi /= hits;
322  (*layTreeObj).densedHitTheta /= hits;
323  (*layTreeObj).correctionFactor /= hits;
324  (*layTreeObj).A /= (*layTreeObj).rho;
325  (*layTreeObj).Z /= (*layTreeObj).rho;
326  (*layTreeObj).rho /= (*layTreeObj).path;
327  }
328  // fill the layer Tree Object
329  (*layTreeObj).tree->Fill();
330  }
331  // reset
332  (*layTreeObj).layerHits = 0;
333  (*layTreeObj).path = 0.;
334  (*layTreeObj).pathInX0 = 0.;
335  (*layTreeObj).pathInL0 = 0.;
336  (*layTreeObj).densedHits = 0;
337  (*layTreeObj).densedHitX = 0.;
338  (*layTreeObj).densedHitY = 0.;
339  (*layTreeObj).densedHitZ = 0.;
340  (*layTreeObj).densedHitR = 0.;
341  (*layTreeObj).densedHitPhi = 0.;
342  (*layTreeObj).densedHitTheta = 0.;
343  (*layTreeObj).correctionFactor = 0.;
344  (*layTreeObj).A = 0.;
345  (*layTreeObj).Z = 0.;
346  (*layTreeObj).rho = 0.;
347  }
348  }
349 }
350 
351 
352 void Trk::MaterialMapper::recordSurfaceHit(const Amg::Vector2D& locpos,
353  const AssociatedMaterial& amhit) const
354 {
355 
356  if (m_processedEvents >= m_maxMappingEvents) return;
357 
358  const Trk::Layer* lay = amhit.associatedLayer();
359  if (m_surfaceValidation && lay) {
360  // try to find the histogram
361  Trk::SurfaceTreeObject* surfTreeObj = surfaceTreeObject(*lay);
362  // now fill
363  if (surfTreeObj) {
364  // fill it
365  (*surfTreeObj).loc1 += locpos[0];
366  (*surfTreeObj).loc2 += locpos[1];
367  (*surfTreeObj).eta += amhit.materialPosition().eta();
368  (*surfTreeObj).correctionFactor += amhit.correctionFactor();
369  (*surfTreeObj).path += amhit.steplength();
370  (*surfTreeObj).pathInX0 += amhit.steplengthInX0();
371  (*surfTreeObj).pathInL0 += amhit.steplengthInL0();
372  (*surfTreeObj).A += amhit.A()*amhit.steplength()*amhit.rho();
373  (*surfTreeObj).Z += amhit.Z()*amhit.steplength()*amhit.rho();
374  (*surfTreeObj).rho += amhit.rho()*amhit.steplength();
375  // increase the surface hits
376  ++(*surfTreeObj).surfaceHits;
377 
378  }
379  }
380 }
381 
382 
383 void Trk::MaterialMapper::finalizeSurfaceHits() const
384 {
385 
386  ATH_MSG_VERBOSE( "finalizeEvent() on Surface hits " );
387 
388  // loop over all layer Tree objects and fill them
389  auto surfObjIter = m_surfaceTrees.begin();
390  auto endIter = m_surfaceTrees.end();
391  // do the loop
392  for ( ; surfObjIter != endIter; ++surfObjIter ) {
393  Trk::SurfaceTreeObject* surfTreeObj = surfObjIter->second;
394  // get the number of single hits
395  int hits = (*surfTreeObj).surfaceHits;
396  // fill only if hits have been there
397  if ((*surfTreeObj).pathInX0 > 0.) {
398  if (hits>1) {
399  (*surfTreeObj).loc1 /= hits;
400  (*surfTreeObj).loc2 /= hits;
401  (*surfTreeObj).eta /= hits;
402  (*surfTreeObj).correctionFactor /= hits;
403  (*surfTreeObj).A /= hits;
404  (*surfTreeObj).Z /= hits;
405  (*surfTreeObj).rho /= hits;
406  }
407  // fill the tree
408  (*surfTreeObj).tree->Fill();
409  }
410  // reset
411  (*surfTreeObj).loc1 = 0.;
412  (*surfTreeObj).loc2 = 0.;
413  (*surfTreeObj).eta = 0.;
414 
415  (*surfTreeObj).correctionFactor = 0.;
416  (*surfTreeObj).path = 0.;
417  (*surfTreeObj).pathInX0 = 0.;
418  (*surfTreeObj).pathInL0 = 0.;
419  (*surfTreeObj).rho = 0.;
420 
421  (*surfTreeObj).surfaceHits = 0;
422 
423  }
424 }
425 
426 
427 Trk::VolumeTreeObject* Trk::MaterialMapper::volumeTreeObject(const Trk::Layer* lay, const Trk::TrackingVolume* tvol) const
428 {
429  // try to find the histogram
430  auto endIter = lay ? m_volumeTrees.end() : m_volumeTreesUnmapped.end();
431  auto findIter = lay ? m_volumeTrees.end() : m_volumeTreesUnmapped.end();
432 
433  Trk::VolumeTreeObject* tvolTreeObj = nullptr;
434 
435  findIter = lay ? m_volumeTrees.find(tvol) : m_volumeTreesUnmapped.find(tvol);
436 
437  if (findIter == endIter) {
438 
439  TString tvolName = tvol->volumeName();
440  tvolName.ReplaceAll("::","_");
441 
442  TString treeName = tvolName;
443  treeName += "_";
444  treeName += m_volumeTreePrefix;
445  if (!lay) treeName += "_UNMAPPED";
446 
447  TString treeTitle = "TrackingVolume : ";
448  treeTitle += tvolName;
449  if (!lay) treeTitle += " - UNMAPPED";
450 
451  TString treeRegName = "/val/";
452  treeRegName += treeName;
453 
454  ATH_MSG_INFO( "No Tree found for " << treeTitle );
455  ATH_MSG_INFO( " -> Booking it now with register name : " << treeRegName.Data() );
456 
457  tvolTreeObj = new Trk::VolumeTreeObject(treeName, treeTitle);
458  if (lay)
459  m_volumeTrees.insert(std::make_pair(tvol,tvolTreeObj));
460  else
461  m_volumeTreesUnmapped.insert(std::make_pair(tvol,tvolTreeObj));
462 
463  // now register the Tree
464  ITHistSvc* tHistSvc = nullptr;
465  if (service("THistSvc",tHistSvc).isFailure()) {
466  ATH_MSG_ERROR( "initialize() Could not find Hist Service -> Switching Tree output for this volume off !" );
467  delete tvolTreeObj; tvolTreeObj = nullptr;
468  }
469  else if (tHistSvc && (tHistSvc->regTree(treeRegName.Data(), (*tvolTreeObj).tree)).isFailure()) {
470  ATH_MSG_ERROR( "initialize() Could not register the validation Tree -> Switching Tree output for this volume off !" );
471  delete tvolTreeObj; tvolTreeObj = nullptr;
472  }
473 
474  } else // a tree is found
475  tvolTreeObj = findIter->second;
476 
477  return tvolTreeObj;
478 }
479 
480 
481 Trk::LayerTreeObject* Trk::MaterialMapper::layerTreeObject(const Trk::Layer& lay, bool full) const
482 {
483  if (!lay.layerIndex().value()) return nullptr;
484 
485  // try to find the histogram
486 
487  Trk::LayerTreeObject* layTreeObj = nullptr;
488 
489  if (full) {
490  auto it = m_layerFullTrees.find(&lay);
491  if (it != m_layerFullTrees.end()) {
492  return it->second;
493  }
494  }
495  else {
496  auto it = m_layerTrees.find(&lay);
497  if (it != m_layerTrees.end()) {
498  return it->second;
499  }
500  }
501 
502  {
503  // check if it is a boundary MaterialLayer
504  const Trk::MaterialLayer* mLayer = dynamic_cast<const Trk::MaterialLayer*>(&lay);
505  if (mLayer)
506  ATH_MSG_INFO("MaterialLayer from BoundarySurface detected.");
507  // get the TrackingVolume
508  const Trk::TrackingVolume* enclosingVolume = lay.enclosingTrackingVolume();
509  TString tvolName = (enclosingVolume) ? enclosingVolume->volumeName() : "BoundaryLayers";
510  tvolName.ReplaceAll("::","_");
511 
512  TString treeName = tvolName;
513  treeName += "_";
514  treeName += m_layerTreePrefix;
515 
516  TString layerType = (lay.surfaceRepresentation().type() == Trk::SurfaceType::Cylinder) ?
517  "CylinderLayer_" : "DiscLayer_";
518  if (full) treeName += "full_";
519  treeName += layerType;
520  if (mLayer) treeName += "boundary_";
521  treeName += lay.layerIndex().value();
522 
523  TString treeTitle = "TrackingVolume :";
524  treeTitle += tvolName;
525 
526  TString treeRegName = "/val/";
527  treeRegName += treeName;
528 
529  ATH_MSG_INFO( "No Tree found for Layer " << lay.layerIndex().value() << " in Volume '" << tvolName << "'.");
530  ATH_MSG_INFO( " -> Booking it now with register name : " << treeRegName.Data() );
531 
532  layTreeObj = new Trk::LayerTreeObject(treeName, treeTitle);
533  if (full) m_layerFullTrees.insert(std::make_pair(&lay,layTreeObj));
534  else m_layerTrees.insert(std::make_pair(&lay,layTreeObj));
535 
536  // now register the Tree
537  ITHistSvc* tHistSvc = nullptr;
538  if (service("THistSvc",tHistSvc).isFailure()) {
539  ATH_MSG_ERROR( "initialize() Could not find Hist Service -> Switching Tree output for this layer off !" );
540  delete layTreeObj; layTreeObj = nullptr;
541  }
542  else if (tHistSvc && (tHistSvc->regTree(treeRegName.Data(), (*layTreeObj).tree)).isFailure()) {
543  ATH_MSG_ERROR( "initialize() Could not register the validation Tree -> Switching Tree output for this layer off !" );
544  delete layTreeObj; layTreeObj = nullptr;
545  }
546 
547  }
548 
549  return layTreeObj;
550 }
551 
552 
553 Trk::SurfaceTreeObject* Trk::MaterialMapper::surfaceTreeObject(const Trk::Layer& lay) const
554 {
555  // try to find the histogram
556  auto endIter = m_surfaceTrees.end();
557  auto findIter = m_surfaceTrees.end();
558 
559  Trk::SurfaceTreeObject* surfTreeObj = nullptr;
560 
561  findIter = m_surfaceTrees.find(&lay);
562  if (findIter == endIter) {
563  // get the TrackingVolume
564  const Trk::TrackingVolume* enclosingVolume = lay.enclosingTrackingVolume();
565  TString volumeName = (enclosingVolume) ? enclosingVolume->volumeName() : "Unknown";
566 
567  TString treeName = m_surfaceTreePrefix;
568  treeName += lay.layerIndex().value();
569 
570  TString treeTitle = "TrackingVolume :";
571  treeTitle += volumeName;
572 
573  TString treeRegName = "/val/";
574  treeRegName += treeName;
575 
576  ATH_MSG_INFO( "No Tree found for Layer " << lay.layerIndex().value()
577  << " in Volume '" << volumeName << "'." );
578  ATH_MSG_INFO( " -> Booking it now with register name : " << treeRegName.Data() );
579 
580  surfTreeObj = new Trk::SurfaceTreeObject(treeName, treeTitle);
581  m_surfaceTrees.insert(std::make_pair(&lay,surfTreeObj));
582 
583  // now register the Tree
584  ITHistSvc* tHistSvc = nullptr;
585  if (service("THistSvc",tHistSvc).isFailure()) {
586  ATH_MSG_INFO( "initialize() Could not find Hist Service -> Switching Tree output for this surface off !" );
587  delete surfTreeObj; surfTreeObj = nullptr;
588  }
589  else if (tHistSvc && (tHistSvc->regTree(treeRegName.Data(), (*surfTreeObj).tree)).isFailure()) {
590  ATH_MSG_INFO( "initialize() Could not register the validation Tree -> Switching Tree output for this surface off !" );
591  delete surfTreeObj; surfTreeObj = nullptr;
592  }
593 
594  } else // a tree is found
595  surfTreeObj = findIter->second;
596  return surfTreeObj;
597 }
598 
599 
600 void Trk::MaterialMapper::handle( const Incident& inc ) {
601  // check the incident type
602  if ( inc.type() == IncidentType::EndEvent && m_materialSteps){
603  // check if the hit collection already contains:
604  ATH_MSG_VERBOSE("EndEvent incident caught, finalize histogramns");
605  // average Eta / Phi
606  m_averageEta /= m_materialSteps;
607  m_averagePhi /= m_materialSteps;
608  // (a) the volume validation
609  finalizeVolumeHits(true); // trees with mapped hits
610  finalizeVolumeHits(false); // trees with unmapped hits
611  // (b) the layer validation
612  finalizeLayerHits();
613  // (c) the surface validation (reference material)
614  finalizeSurfaceHits();
615  // fill the main TTree & reset
616  m_validationTree->Fill();
617  // initialize the eta/phi
618  m_averageEta = 0.;
619  m_averagePhi = 0.;
620  // initialize the rest
621  m_materialSteps = 0;
622  m_mappedPath = 0.;
623  m_mappedPathInX0 = 0.;
624  m_mappedPathInL0 = 0.;
625  m_mappedPathRho = 0.;
626  m_mappedPathZARho = 0.;
627  // increate counter
628  ++m_processedEvents;
629  }
630  }
631 
632 
633 void Trk::MaterialMapper::bookValidationTree()
634 {
635 
636  ATH_MSG_INFO( "Booking the Validation Tree ... " );
637 
638  // now register the Tree
639  ITHistSvc* tHistSvc = nullptr;
640 
641  // (1) Main MaterialMapper TTree
642  // ------------- validation section ------------------------------------------
643  m_validationTree = new TTree(m_validationTreeName.c_str(), m_validationTreeDescription.c_str());
644 
645  // position coordinates of the update
646  m_validationTree->Branch("Eta", &m_averageEta, "averageEta/F");
647  m_validationTree->Branch("Phi", &m_averagePhi, "averagePhiF");
648  m_validationTree->Branch("Path", &m_mappedPath, "path/F");
649  m_validationTree->Branch("PathInX0", &m_mappedPathInX0, "pathInX0/F");
650  m_validationTree->Branch("PathInL0", &m_mappedPathInL0, "pathInL0/F");
651  m_validationTree->Branch("PathRho", &m_mappedPathRho, "pathRho/F");
652  m_validationTree->Branch("PathZARho", &m_mappedPathZARho, "pathZARho/F");
653  m_validationTree->Branch("UnmappedPathInX0", &m_unmappedPathInX0, "unmappedPathInX0/F");
654  m_validationTree->Branch("MaterialSteps", &m_materialSteps, "steps/I");
655  m_validationTree->Branch("Mapped", m_mapped, "mapped[steps]/I");
656  m_validationTree->Branch("MaterialAccumPathInX0", m_materialAccumPathInX0, "materialAccumPinX0[steps]/F");
657  m_validationTree->Branch("MaterialAccumPathZARho", m_materialAccumPathZARho, "materialAccumPZARho[steps]/F");
658  m_validationTree->Branch("MaterialStepPath", m_materialStepPath, "materialStepPath[steps]/F");
659  m_validationTree->Branch("MaterialStepX0", m_materialStepX0, "materialStepX0[steps]/F");
660  m_validationTree->Branch("MaterialStepL0", m_materialStepL0, "materialStepL0[steps]/F");
661  m_validationTree->Branch("MaterialStepZ", m_materialStepZ, "materialStepZ[steps]/F");
662  m_validationTree->Branch("MaterialStepA", m_materialStepA, "materialStepA[steps]/F");
663  m_validationTree->Branch("MaterialStepRho", m_materialStepRho, "materialStepRho[steps]/F");
664  m_validationTree->Branch("MaterialStepPositionX", m_materialStepPositionX , "materialStepX[steps]/F");
665  m_validationTree->Branch("MaterialStepPositionY", m_materialStepPositionY , "materialStepY[steps]/F");
666  m_validationTree->Branch("MaterialStepPositionZ", m_materialStepPositionZ , "materialStepZ[steps]/F");
667  m_validationTree->Branch("MaterialStepPositionR", m_materialStepPositionR , "materialStepR[steps]/F");
668  m_validationTree->Branch("MaterialProjPositionX", m_materialProjPositionX , "materialProjX[steps]/F");
669  m_validationTree->Branch("MaterialProjPositionY", m_materialProjPositionY , "materialProjY[steps]/F");
670  m_validationTree->Branch("MaterialProjPositionZ", m_materialProjPositionZ , "materialProjZ[steps]/F");
671  m_validationTree->Branch("MaterialProjPositionR", m_materialProjPositionR , "materialProjR[steps]/F");
672  m_validationTree->Branch("MaterialProjDistance", m_materialProjDistance , "materialProjD[steps]/F");
673 
674  // now register the Tree
675  if (service("THistSvc",tHistSvc).isFailure()) {
676  ATH_MSG_ERROR("initialize() Could not find Hist Service -> Switching ValidationMode Off !" );
677  delete m_validationTree; m_validationTree = nullptr;
678  return;
679  }
680  if ((tHistSvc->regTree(m_validationTreeFolder.c_str(), m_validationTree)).isFailure()) {
681  ATH_MSG_ERROR("initialize() Could not register the validation Tree -> Switching ValidationMode Off !" );
682  delete m_validationTree; m_validationTree = nullptr;
683  return;
684  }
685 
686  ATH_MSG_INFO( " ... successful." );
687 
688 }
689 
690 
Trk::MaterialMapper::finalizeSurfaceHits
void finalizeSurfaceHits() const
Finalize the SingleLayer Steps.
Definition: MaterialMapper.cxx:383
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Trk::MaterialMapper::finalizeVolumeHits
void finalizeVolumeHits(bool mapped=true) const
Finalize the Volume Steps.
Definition: MaterialMapper.cxx:208
Trk::AssociatedMaterial::steplengthInL0
double steplengthInL0() const
Access method : steplength.
Definition: AssociatedMaterial.h:182
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition: GeoPrimitives.h:48
Trk::AssociatedMaterial::steplength
double steplength() const
Access method : steplength.
Definition: AssociatedMaterial.h:167
Trk::TrackingVolume::getMotherVolume
const TrackingVolume * getMotherVolume() const
Return the MotherVolume - if it exists.
skel.it
it
Definition: skel.GENtoEVGEN.py:423
Trk::MaterialMapper::bookValidationTree
void bookValidationTree()
Validation : book the Tree.
Definition: MaterialMapper.cxx:633
DiscLayer.h
Trk::AssociatedMaterial::associatedLayer
const Trk::Layer * associatedLayer() const
Trivial Access methods.
Definition: AssociatedMaterial.h:155
beamspotman.posX
posX
Definition: beamspotman.py:1624
Trk::AssociatedMaterial::x0
double x0() const
Access method : material X0/A/Z/rho.
Definition: AssociatedMaterial.h:160
Trk::MaterialMapper::finalizeLayerHits
void finalizeLayerHits() const
Finalize the SingleLayer Steps.
Definition: MaterialMapper.cxx:296
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
Trk::SurfaceTreeObject
Definition: MaterialMapper.h:147
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::MaterialMapper::layerTreeObject
LayerTreeObject * layerTreeObject(const Layer &lay, bool full=false) const
find (&&,||) create the LayerTreeObject
Definition: MaterialMapper.cxx:481
Trk::LayerTreeObject::layerHits
int layerHits
Definition: MaterialMapper.h:90
Trk::MaterialMapper::initialize
StatusCode initialize()
AlgTool initialize method.
Definition: MaterialMapper.cxx:99
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
Trk::Layer::surfaceRepresentation
virtual const Surface & surfaceRepresentation() const =0
Transforms the layer into a Surface representation for extrapolation.
Trk::MaterialMapper::handle
void handle(const Incident &inc)
Handle the incident from the incident service.
Definition: MaterialMapper.cxx:600
TRKDETDESCRTOOLS_MAXLAYERHITS
#define TRKDETDESCRTOOLS_MAXLAYERHITS
Definition: MaterialMapper.h:25
beamspotman.n
n
Definition: beamspotman.py:731
beamspotman.posZ
posZ
Definition: beamspotman.py:1624
Trk::AssociatedMaterial::A
double A() const
Definition: AssociatedMaterial.h:195
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
LArG4ShowerLibProcessing.hits
hits
Definition: LArG4ShowerLibProcessing.py:136
MaterialMapper.h
dumpFileToPlots.treeName
string treeName
Definition: dumpFileToPlots.py:20
CylinderLayer.h
Trk::AssociatedMaterial::l0
double l0() const
Definition: AssociatedMaterial.h:189
Trk::LayerIndex::value
int value() const
layerIndex expressed in an integer
Definition: LayerIndex.h:71
Trk::MaterialMapper::surfaceTreeObject
SurfaceTreeObject * surfaceTreeObject(const Layer &lay) const
find (&&,||) create the LayerTreeObject
Definition: MaterialMapper.cxx:553
Trk::AssociatedMaterial::associatedTrackingVolume
const Trk::TrackingVolume * associatedTrackingVolume() const
Trivial Access methods.
Definition: AssociatedMaterial.h:149
Trk::AssociatedMaterial
Definition: AssociatedMaterial.h:33
Trk::TrackingVolume::volumeName
const std::string & volumeName() const
Returns the VolumeName - for debug reason, might be depreciated later.
Trk::MaterialMapper::volumeTreeObject
VolumeTreeObject * volumeTreeObject(const Layer *lay=0, const TrackingVolume *tvol=0) const
find (&&,||) create a VolumeTreObject
Definition: MaterialMapper.cxx:427
Trk::LayerTreeObject
Definition: MaterialMapper.h:69
Trk::AssociatedMaterial::Z
double Z() const
Definition: AssociatedMaterial.h:201
Trk::MaterialMapper::recordLayerHit
void recordLayerHit(const AssociatedMaterial &amhit, bool full=false) const
Record material hit - if various hits per uniform layer are recorded, or if you want to record one fu...
Definition: MaterialMapper.cxx:241
MaterialLayer.h
Trk::AssociatedMaterial::materialPosition
const Amg::Vector3D & materialPosition() const
Trivial Access methods.
Definition: AssociatedMaterial.h:139
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
find_data.full
full
Definition: find_data.py:27
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
Trk::MaterialMapper::recordSurfaceHit
void recordSurfaceHit(const Amg::Vector2D &locpos, const AssociatedMaterial &amhit) const
Record material hit on a surface.
Definition: MaterialMapper.cxx:352
beamspotman.posY
posY
Definition: beamspotman.py:1624
Trk::AssociatedMaterial::correctionFactor
double correctionFactor() const
Trivial Access methods.
Definition: AssociatedMaterial.h:144
TrackingVolume.h
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::VolumeTreeObject
Definition: MaterialMapper.h:39
Trk::Layer::enclosingTrackingVolume
const TrackingVolume * enclosingTrackingVolume() const
get the confining TrackingVolume
Trk::SurfaceType::Cylinder
@ Cylinder
Trk::AssociatedMaterial::steplengthInX0
double steplengthInX0() const
Access method : steplength.
Definition: AssociatedMaterial.h:174
mapped
std::vector< std::string > mapped
Definition: hcg.cxx:51
Trk::VolumeTreeObject::eta
float eta
Definition: MaterialMapper.h:43
AssociatedMaterial.h
Trk::AssociatedMaterial::rho
double rho() const
Definition: AssociatedMaterial.h:207
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
Trk::MaterialLayer
Definition: MaterialLayer.h:37
Trk::MaterialMapper::~MaterialMapper
virtual ~MaterialMapper()
Virtual destructor.
Definition: MaterialMapper.cxx:69
Trk::MaterialMapper::MaterialMapper
MaterialMapper(const std::string &, const std::string &, const IInterface *)
AlgTool like constructor.
Definition: MaterialMapper.cxx:22
Trk::MaterialMapper::finalize
StatusCode finalize()
AlgTool finalize method.
Definition: MaterialMapper.cxx:116
MaterialAssociationType.h
AthAlgTool
Definition: AthAlgTool.h:26
Trk::TrackingVolume
Definition: TrackingVolume.h:121
Trk::Surface::type
constexpr virtual SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
mag
Scalar mag() const
mag method
Definition: AmgMatrixBasePlugin.h:25
Trk::MaterialMapper::recordMaterialHit
void recordMaterialHit(const AssociatedMaterial &amhit, const Amg::Vector3D &projectedPosition) const
Record material hit along the recording.
Definition: MaterialMapper.cxx:122
Trk::Layer
Definition: Layer.h:73
Trk::Layer::layerIndex
const LayerIndex & layerIndex() const
get the layerIndex
Generated on Sun Jun 30 2024 21:21:10 for ATLAS Offline Software by   doxygen 1.8.18