ATLAS Offline Software
ReFitTrackWithTruth.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 // ReFitTrackWithTruth.cxx
7 // Implementation file for class ReFitTrackWithTruth
9 // version 1.0 03/10/19 Gabriel Facini (Rel21)
10 // version 2.0 22/03/23 Andrea Sciandra (Rel22/Rel23)
12 
13 //SiTBLineFitter includes
15 
16 // Gaudi includes
17 #include "GaudiKernel/TypeNameString.h"
18 
30 #include "TRandom3.h"
31 
32 #include <vector>
33 
34 // Constructor with parameters:
35 Trk::ReFitTrackWithTruth::ReFitTrackWithTruth(const std::string &name, ISvcLocator *pSvcLocator) :
36  AthAlgorithm(name,pSvcLocator)
37 {
38 
39 }
40 
41 // Initialize method:
43 {
44  ATH_MSG_INFO ("ReFitTrackWithTruth::initialize()");
45 
46  ATH_CHECK (m_siHitCollectionName.initialize());
47  ATH_CHECK (m_SDOContainerName.initialize());
48  ATH_CHECK (m_truthMapName.initialize());
49 
50  // get tracks
51  ATH_CHECK( m_inputTrackColName.initialize() );
52 
53  if (m_ITrackFitter.retrieve().isFailure()) {
54  ATH_MSG_FATAL ("Failed to retrieve tool "<<m_ITrackFitter.typeAndName());
55  return StatusCode::FAILURE;
56  } else {
57  ATH_MSG_INFO ("Retrieved general fitter " << m_ITrackFitter.typeAndName());
58  }
59 
60  if (m_trkSummaryTool.retrieve().isFailure()) {
61  ATH_MSG_FATAL ("Failed to retrieve tool " << m_trkSummaryTool);
62  return StatusCode::FAILURE;
63  } else {
64  ATH_MSG_INFO ("Retrieved tool " << m_trkSummaryTool.typeAndName());
65  }
66 
67  // needed if want to check on shared clusters
68  if ( m_assoTool.retrieve().isFailure() ) {
69  ATH_MSG_FATAL ("Failed to retrieve tool " << m_assoTool);
70  return StatusCode::FAILURE;
71  } else ATH_MSG_INFO("Retrieved tool " << m_assoTool);
72 
73 
74  // Configuration of the material effects
75  m_ParticleHypothesis = Trk::ParticleSwitcher::particle[m_matEffects];
76 
77  // Get ID Helper
78  ATH_CHECK(detStore()->retrieve(m_idHelper, "AtlasID"));
79 
80  // Get Pixel Helper
81  if (detStore()->retrieve(m_pixelID, "PixelID").isFailure()) {
82  ATH_MSG_FATAL ("Could not get Pixel ID helper");
83  return StatusCode::FAILURE;
84  }
85 
86  //set seed for random-number generator
87  m_random->SetSeed();
88 
89  // ensure the vector is the correct size in case the user does not want to smear
90  if( m_resolutionRPhi.size() < 4 ) { m_resolutionRPhi={0.,0.,0.,0.}; }
91  if( m_resolutionZ.size() < 4 ) { m_resolutionZ={0.,0.,0.,0.}; }
92 
93  if( m_errorRPhi.size() < 4 ) { m_errorRPhi={1.,1.,1.,1.}; }
94  if( m_errorZ.size() < 4 ) { m_errorZ={1.,1.,1.,1.}; }
95 
96  ATH_MSG_INFO (" Resolutions " << m_resolutionRPhi.size() << " " << m_resolutionZ.size());
97  for( unsigned int i=0; i<m_resolutionRPhi.size(); i++) {
98  ATH_MSG_INFO (" " << i << " " << m_resolutionRPhi[i] << "\t" << m_resolutionZ[i]);
99  }
100  ATH_MSG_INFO (" Error SR " << m_errorRPhi.size() << " " << m_errorZ.size());
101  for( unsigned int i=0; i<m_errorRPhi.size(); i++) {
102  ATH_MSG_INFO (" " << i << " " << m_errorRPhi[i] << "\t" << m_errorZ[i]);
103  }
104 
105  ATH_CHECK(m_outputTrackCollectionName.initialize());
106  ATH_MSG_DEBUG("m_outputTrackCollectionName: " << m_outputTrackCollectionName);
107 
108  return StatusCode::SUCCESS;
109 }
110 
111 // Execute method:
113 {
114  ATH_MSG_DEBUG ("ReFitTrackWithTruth::execute()");
115  std::unique_ptr<Trk::PRDtoTrackMap> prd_to_track_map(m_assoTool->createPRDtoTrackMap());
116  const EventContext& ctx = Gaudi::Hive::currentContext();
117 
118  SG::ReadHandle<TrackCollection> tracks(m_inputTrackColName, ctx);
119  if (!tracks.isValid()) {
120  ATH_MSG_ERROR(m_inputTrackColName.key() << " not found");
121  return StatusCode::FAILURE;
122  }
123 
124  //retrieve truth information needed (SiHitCollection)
125  SG::ReadHandle<SiHitCollection> siHits(m_siHitCollectionName, ctx);
126  if (!siHits.isValid()) {
127  ATH_MSG_WARNING( "Error retrieving SiHitCollection " << m_siHitCollectionName);
128  return StatusCode::FAILURE;
129  }
130 
131  //retrieve truth information needed (SDOCollection)
132  SG::ReadHandle<InDetSimDataCollection> sdoCollection(m_SDOContainerName, ctx);
133  if (!sdoCollection.isValid()) {
134  ATH_MSG_WARNING( "Error retrieving SDOCollection " << m_SDOContainerName );
135  return StatusCode::FAILURE;
136  }
137 
138  //retrieve truth map
139  SG::ReadHandle<TrackTruthCollection> truthMap(m_truthMapName, ctx);
140  if (!truthMap.isValid()) {
141  ATH_MSG_WARNING( "Error retrieving truth map " << m_truthMapName );
142  return StatusCode::FAILURE;
143  }
144 
145  // create new collection of tracks to write in storegate
146  std::vector<std::unique_ptr<Trk::Track> > newtracks;
147 
148  // loop over tracks
149  for (TrackCollection::const_iterator itr = tracks->begin(); itr != tracks->end(); ++itr) {
150  ATH_MSG_DEBUG("input track");
151 
152  // Get original parameters
153  const TrackParameters* origPerigee = (*itr)->perigeeParameters();
154  double od0(0);
155  double oz0(0);
156  double ophi0(0);
157  double otheta(0);
158  double oqOverP(0);
159  if (!origPerigee){
160  ATH_MSG_WARNING("Cannot get original parameters");
161  }
162  else if (msgLvl(MSG::DEBUG)) {
163  od0 = origPerigee->parameters()[Trk::d0];
164  oz0 = origPerigee->parameters()[Trk::z0];
165  ophi0 = origPerigee->parameters()[Trk::phi0];
166  otheta = origPerigee->parameters()[Trk::theta];
167  oqOverP = origPerigee->parameters()[Trk::qOverP];
168  ATH_MSG_DEBUG ("Original parameters " << od0 << " " << oz0 << " " << ophi0 << " " << otheta << " " << oqOverP);
169  }
170 
171  // get the unique ID of truth particle matched to the track
172  float minProb(0);
173 
174  // copy from DenseEnvironmentsAmbiguityProcessorTool.cxx
176  tracklink.setElement(const_cast<Trk::Track*>(*itr));
177  tracklink.setStorableObject(*tracks);
178  const ElementLink<TrackCollection> tracklink2=tracklink;
179 
180  TrackTruthCollection::const_iterator found = truthMap->find(tracklink2);
181  if ( found == truthMap->end() ) { continue; }
182  if ( !found->second.particleLink().isValid() ) { continue; }
183  if ( found->second.probability() < minProb ) { continue; }
184  int uniqueIdToMatch = HepMC::uniqueID(&(found->second.particleLink()));
185  ATH_MSG_DEBUG ("Unique ID to match " << uniqueIdToMatch);
186 
187  // now that have the track, loop through and build a new set of measurements for a track fit
188  std::vector<const Trk::MeasurementBase*> measurementSet;
189  std::vector<const Trk::MeasurementBase*> trash;
190 
191  ATH_MSG_DEBUG ("Loop over measurementsOnTrack " << (*itr)->measurementsOnTrack()->size() );
192  for (const auto *measurement : *((*itr)->measurementsOnTrack())) {
193  ATH_MSG_DEBUG ("Next Measurement: " << *measurement);
194 
195  // Get measurement as RIO_OnTrack
196  const Trk::RIO_OnTrack* rio = dynamic_cast <const Trk::RIO_OnTrack*>( measurement );
197  if (rio == nullptr) {
198  ATH_MSG_WARNING("Cannot get RIO_OnTrack from measurement. Hit will NOT be included in track fit.");
199  continue;
200  }
201 
202  // if not a pixel cluster, keep the measurement as is and press on
203  const Identifier& surfaceID = (rio->identify()) ;
204  if( !m_idHelper->is_pixel(surfaceID) ) {
205  measurementSet.push_back( measurement );
206  continue;
207  }
208 
209  // only PIXEL CLUSTERS from here on
210  const InDet::PixelCluster* pix = dynamic_cast<const InDet::PixelCluster*>(rio->prepRawData());
211  if (pix == nullptr) {
212  ATH_MSG_WARNING("Cannot get PixelCluster from RIO_OnTrack");
213  continue;
214  }
215 
216  const InDetDD::SiDetectorElement* element = pix->detectorElement();
217  const InDetDD::SiDetectorDesign* design = static_cast<const InDetDD::SiDetectorDesign*>(&element->design());
218  // To locate cluster module
219  Identifier clusterId = pix->identify();
220  int bec = m_pixelID->barrel_ec(clusterId);
221  int layer_disk = m_pixelID->layer_disk(clusterId);
222 
223  // Do any SDOs in reconstructed cluster match uniqueIdToMatch
224  // Should always return true for pseudotracks. For reco tracks could differ
225  bool hasSDOMatch = IsClusterFromTruth( pix, uniqueIdToMatch, *sdoCollection );
226 
227  // --- hit flags' logic
228  // m_saveWrongHits - only has impact on reco track, where hit might be from wrong particle
229  // m_fixWrongHits - fix will only work on hits that are not noise, noise hits behavior unchanged by this flag
230  // m_rejNoiseHits - if saving but don't want to save these
231 
232  // could be used already here to correct wrong hits
233  double maxEnergyDeposit(-1);
234  SiHit maxEDepSiHit;
235 
236  if ( !hasSDOMatch ) {
237  ATH_MSG_DEBUG ("No SDO matching cluster");
238  // save wrong hits
239  if (not m_saveWrongHits) continue;
240  const std::vector<SiHit> matchedSiHits = matchSiHitsToCluster( -999 , pix, siHits );
241  if( matchedSiHits.empty() ) { // then noise, and go to next hit
242  if( !m_rejNoiseHits ) { measurementSet.push_back( measurement ); }
243  continue;
244  } // is a real hit, just not from the particle in question -- will NOT go to next hit just yet
245  if( m_fixWrongHits ) { // if fix, use highest energy truth hit
246  for( const auto& siHit : matchedSiHits ) {
247  if (siHit.energyLoss() > maxEnergyDeposit) {
248  maxEnergyDeposit = siHit.energyLoss();
249  maxEDepSiHit = siHit;
250  }
251  }
252  } else { // save the wrong hit as is and go to next cluster
253  measurementSet.push_back( measurement );
254  continue;
255  }
256  } else { // hasSDOMatch
257  // Get All SiHits truth matched to the reconstruction cluster for the particle associated with the track
258  const std::vector<SiHit> matchedSiHits = matchSiHitsToCluster( uniqueIdToMatch, pix, siHits );
259  if( matchedSiHits.empty() ) {
260  ATH_MSG_WARNING ("No SiHit matching cluster");
261  continue; // should NOT HAPPEN for pseudotracks
262  }
263  ATH_MSG_DEBUG ("N SiHit matching cluster: " << matchedSiHits.size());
264 
265  // If multiple SiHits / cluster FOR THE SAME TRUTH PARTICLE,
266  // Take position of SiHit giving the most energy
267  for( const auto& siHit : matchedSiHits ) {
268  if (siHit.energyLoss() > maxEnergyDeposit) {
269  maxEnergyDeposit = siHit.energyLoss();
270  maxEDepSiHit = siHit;
271  }
272  }
273  } // if else hasSDOMatch
274 
275  // - Retrieve true position of cluster from entry/exit point of truth particle
276  //get average position of entry/exit point
277  HepGeom::Point3D<double> averagePosition = (maxEDepSiHit.localStartPosition() + maxEDepSiHit.localEndPosition()) * 0.5;
278  ATH_MSG_DEBUG (" Average position : " << averagePosition);
279 
280  // SiHit coordinate system i.e. localStartPosition
281  // z = eta -> this is y in ATLAS sensor local frame
282  // y = phi -> this is x in ATLAS sensor local frame
283  // USE hitLocalToLocal from SiDetectorElement
284 
285  HepGeom::Point3D<double> smearedPosition = smearTruthPosition( averagePosition, bec, layer_disk, design );
286  ATH_MSG_DEBUG (" Smeared position : " << smearedPosition );
287 
288  auto locparOrig = rio->localParameters();
289  ATH_MSG_DEBUG(" Original locpar " << locparOrig);
290 
291  Trk::LocalParameters locpar = element->hitLocalToLocal(smearedPosition.z(), smearedPosition.y()); // eta, phi
292  ATH_MSG_DEBUG(" locpar " << locpar);
293 
294  InDetDD::SiLocalPosition centroid(locpar.get(Trk::loc2), locpar.get(Trk::loc1), 0); // eta, phi, depth
295  const Amg::Vector3D& globPos = element->globalPosition(centroid);
296 
297  InDet::SiWidth pixWidth = pix->width();
298  Amg::MatrixX cov = pix->localCovariance();
299  // Original code took width / nrows (or columns)* 1/sqrt(12) to check if > 0...
300  // only need to check width to test against 0...
301  // disk one is layer 0
302  if(bec!=0) layer_disk++;
303  if(pixWidth.phiR()>0) {
304  float error(1.0);
305  error = getPhiPosResolution(layer_disk)*getPhiPosErrorFactor(layer_disk);
306  cov(0,0) = error*error;
307  } else {
308  ATH_MSG_WARNING("pixWidth.phiR not > 0");
309  }
310 
311  if(pixWidth.z()>0) {
312  float error(1.0);
313  error = getEtaPosResolution(layer_disk)*getEtaPosErrorFactor(layer_disk);
314  cov(1,1) = error*error;
315  } else {
316  ATH_MSG_WARNING("pixWidth.z not > 0");
317  }
318 
319  auto iH = element->identifyHash();
320 
321  InDet::PixelClusterOnTrack* pcot = new InDet::PixelClusterOnTrack(pix,std::move(locpar),std::move(cov),iH,globPos,
322  pix->gangedPixel(),
323  false);
324 
325  if(pcot) {
326  measurementSet.push_back( pcot);
327  trash.push_back(pcot);
328  } else {
329  ATH_MSG_WARNING("Could not make new PixelClusterOnTrack");
330  }
331 
332  } // loop over measurements on track
333 
334 
335 
336  ATH_MSG_DEBUG ("Fit new tracks with measurementSet : " << measurementSet.size());
337  std::unique_ptr<Trk::Track> newtrack;
338  try {
339  newtrack = m_ITrackFitter->fit(ctx,
340  measurementSet,
341  *origPerigee,
342  m_runOutlier,
343  m_ParticleHypothesis);
344  }
345  catch(const std::exception& e) {
346  ATH_MSG_ERROR ("Refit Logic Error. No new track. Message: " << e.what());
347  newtrack = nullptr;
348  }
349 
350  ATH_MSG_DEBUG ("Track fit is done!");
351 
352  if (msgLvl(MSG::DEBUG)) {
353  if (!newtrack) { ATH_MSG_DEBUG ("Refit Failed"); }
354  else {
355 
356  ATH_MSG_VERBOSE ("re-fitted track:" << *newtrack);
357  const Trk::Perigee* aMeasPer = newtrack->perigeeParameters();
358  if (aMeasPer==nullptr){
359  ATH_MSG_ERROR ("Could not get Trk::MeasuredPerigee");
360  } else {
361  double d0 = aMeasPer->parameters()[Trk::d0];
362  double z0 = aMeasPer->parameters()[Trk::z0];
363  double phi0 = aMeasPer->parameters()[Trk::phi0];
364  double theta = aMeasPer->parameters()[Trk::theta];
365  double qOverP = aMeasPer->parameters()[Trk::qOverP];
366  ATH_MSG_DEBUG ("Refitted parameters differences "
367  << (od0-d0)/od0 << " "
368  << (oz0-z0)/oz0 << " "
369  << (ophi0-phi0)/ophi0 << " "
370  << (otheta-theta)/otheta << " "
371  << (oqOverP-qOverP)/oqOverP );
372  } // aMeasPer exists
373  } // newtrack exists
374  } // if debug
375 
376  if (newtrack) { newtracks.push_back(std::move(newtrack)); }
377  else { ATH_MSG_WARNING ("Refit Failed"); }
378 
379  } // loop over tracks
380 
381  ATH_MSG_VERBOSE ("Add PRDs to assoc tool.");
382 
383  // recreate the summaries on the final track collection with correct PRD tool
384  for(const std::unique_ptr<Trk::Track> &new_track : newtracks ) {
385  if((m_assoTool->addPRDs(*prd_to_track_map, *new_track)).isFailure()) {ATH_MSG_WARNING("Failed to add PRDs to map");}
386  }
387 
388  ATH_MSG_VERBOSE ("Recalculate the summary");
389  // and copy tracks from vector of non-const tracks to collection of const tracks
390  std::unique_ptr<TrackCollection> new_track_collection = std::make_unique<TrackCollection>();
391  new_track_collection->reserve(newtracks.size());
392  for(std::unique_ptr<Trk::Track> &new_track : newtracks ) {
393  m_trkSummaryTool->computeAndReplaceTrackSummary(ctx, *new_track, false /* DO NOT suppress hole search*/);
394  new_track_collection->push_back(std::move(new_track));
395  }
396 
397  ATH_MSG_VERBOSE ("Save tracks");
398  ATH_CHECK(SG::WriteHandle<TrackCollection>(m_outputTrackCollectionName).record(std::move(new_track_collection)));
399 
400  ATH_MSG_INFO ("ReFitTrackWithTruth::execute() completed");
401  return StatusCode::SUCCESS;
402 }
403 
404 std::vector<SiHit> Trk::ReFitTrackWithTruth::matchSiHitsToCluster( const int uniqueIdToMatch,
405  const InDet::PixelCluster* pixClus,
406  SG::ReadHandle<AtlasHitsVector<SiHit>> &siHitCollection) const {
407 
408  // passing a negative unique ID value skip this requirement - can get multiple SiHits upon return from different particles
409 
410  ATH_MSG_VERBOSE( " Have " << (*siHitCollection).size() << " SiHits to look through" );
411  std::vector<SiHit> matchingHits;
412 
413  // Check if we have detector element -- needed to find the local position of the SiHits
414  const InDetDD::SiDetectorElement* de = pixClus->detectorElement();
415  if(!de) {
416  ATH_MSG_WARNING("Do not have detector element to find the local position of SiHits!");
417  return matchingHits;
418  }
419 
420  // To locate cluster module
421  Identifier clusterId = pixClus->identify();
422 
423  std::vector<const SiHit* > multiMatchingHits;
424 
425  // match SiHits to unique ID and make sure in same module as reco hit
426  for ( const auto& siHit : *siHitCollection) {
427 
428  if ( uniqueIdToMatch > 0 ) { // negative uniqueIdToMatch will keep all
429  if ( HepMC::uniqueID(&(siHit.particleLink())) != uniqueIdToMatch ) { continue; }
430  }
431 
432  // Check if it is a Pixel hit
433  if( !siHit.isPixel() ) { continue; }
434 
435  // Match to the cluster module
436  if( m_pixelID->barrel_ec(clusterId) != siHit.getBarrelEndcap() ) { continue; }
437  if( m_pixelID->layer_disk(clusterId)!= siHit.getLayerDisk() ) { continue; }
438  if( m_pixelID->phi_module(clusterId)!= siHit.getPhiModule() ) { continue; }
439  if( m_pixelID->eta_module(clusterId)!= siHit.getEtaModule() ) { continue; }
440 
441  // Have SiHits in the same module as the cluster at this point
442  ATH_MSG_DEBUG("Hit is on the same module");
443  multiMatchingHits.push_back(&siHit);
444 
445  } // loop over SiHitCollection
446 
447 
448  //Now we will now make 1 SiHit for each true particle if the SiHits "touch" other
449  std::vector<const SiHit* >::iterator siHitIter = multiMatchingHits.begin();
450  std::vector<const SiHit* >::iterator siHitIter2 = multiMatchingHits.begin();
451  ATH_MSG_DEBUG( "Found " << multiMatchingHits.size() << " SiHit " );
452 
453  // double loop - for each matching SiHit, consider all the SiHits _next_ in the collection
454  // to see if they overlap.
455  // if overlapping, combine and only consider new merged hits
456  for ( ; siHitIter != multiMatchingHits.end(); ++siHitIter) {
457  const SiHit* lowestXPos = *siHitIter;
458  const SiHit* highestXPos = *siHitIter;
459 
460 
461  // We will merge these hits
462  std::vector<const SiHit* > ajoiningHits;
463  ajoiningHits.push_back( *siHitIter );
464 
465  siHitIter2 = siHitIter+1;
466  while ( siHitIter2 != multiMatchingHits.end() ) {
467  // Need to come from the same truth particle
468 
469  // wasn't the unique ID match already done!?
470  if ( HepMC::uniqueID(&((*siHitIter)->particleLink())) != HepMC::uniqueID(&((*siHitIter2)->particleLink()))) {
471  ++siHitIter2;
472  continue;
473  }
474 
475  // Check to see if the SiHits are compatible with each other.
476  if (std::abs((highestXPos->localEndPosition().x()-(*siHitIter2)->localStartPosition().x()))<0.00005 &&
477  std::abs((highestXPos->localEndPosition().y()-(*siHitIter2)->localStartPosition().y()))<0.00005 &&
478  std::abs((highestXPos->localEndPosition().z()-(*siHitIter2)->localStartPosition().z()))<0.00005 )
479  {
480  highestXPos = *siHitIter2;
481  ajoiningHits.push_back( *siHitIter2 );
482  // Dont use hit more than once
483  siHitIter2 = multiMatchingHits.erase( siHitIter2 );
484  //--siHitIter2; // maybe
485  }else if (std::abs((lowestXPos->localStartPosition().x()-(*siHitIter2)->localEndPosition().x()))<0.00005 &&
486  std::abs((lowestXPos->localStartPosition().y()-(*siHitIter2)->localEndPosition().y()))<0.00005 &&
487  std::abs((lowestXPos->localStartPosition().z()-(*siHitIter2)->localEndPosition().z()))<0.00005)
488  {
489  lowestXPos = *siHitIter2;
490  ajoiningHits.push_back( *siHitIter2 );
491  // Dont use hit more than once
492  siHitIter2 = multiMatchingHits.erase( siHitIter2 );
493  // --siHitIter2; // maybe
494  } else {
495  ++siHitIter2;
496  }
497  } // loop over matching SiHits to see if any overlap
498 
499  if( ajoiningHits.empty()){
500  ATH_MSG_WARNING("This should really never happen");
501  continue;
502  }
503  if(ajoiningHits.size() == 1){
504  // Copy Si Hit ready to return
505  matchingHits.push_back( *ajoiningHits[0] );
506  continue;
507  }
508  // Build new SiHit and merge information together.
509  ATH_MSG_DEBUG("Merging " << ajoiningHits.size() << " SiHits together." );
510 
511 
512  float energyDep(0);
513  float time(0);
514  for( auto& siHit : ajoiningHits){
515  energyDep += siHit->energyLoss();
516  time += siHit->meanTime();
517  }
518  time /= (float)ajoiningHits.size();
519 
520  matchingHits.emplace_back(lowestXPos->localStartPosition(),
521  highestXPos->localEndPosition(),
522  energyDep,
523  time,
524  HepMC::uniqueID(&((*siHitIter)->particleLink())),
525  0, // 0 for pixel 1 for Pixel
526  (*siHitIter)->getBarrelEndcap(),
527  (*siHitIter)->getLayerDisk(),
528  (*siHitIter)->getEtaModule(),
529  (*siHitIter)->getPhiModule(),
530  (*siHitIter)->getSide() );
531  ATH_MSG_DEBUG("Finished Merging " << ajoiningHits.size() << " SiHits together." );
532  } // loop over all matching SiHits
533 
534  return matchingHits;
535 }
536 
538  const int uniqueIdToMatch,
539  const InDetSimDataCollection &sdoCollection) {
540 
541  // Should be true for all pseudotracks
542  // Can be false for reco tracks - misassigned hits
543 
544  bool match(false);
545 
546  // loop over reconstructed energy depoists in the cluster
547  for( const auto &hitIdentifier : pixClus->rdoList() ) {
548 
549  // find the rdo in the sdo collection
550  auto pos = sdoCollection.find(hitIdentifier);
551  if( pos == sdoCollection.end() ) { continue; }
552 
553  // get the unique ID from each deposit
554  for( const auto& deposit : pos->second.getdeposits() ){
555  if( !deposit.first ){ continue; } // if truthparticle(?) link doesn't exists? Energy deposit is still known
556  if( HepMC::uniqueID(&(deposit.first)) != uniqueIdToMatch ) { continue; }
557  match = true;
558  break;
559  }
560  if(match) { break; }
561  }
562 
563  return match;
564 }
565 
566 HepGeom::Point3D<double> Trk::ReFitTrackWithTruth::smearTruthPosition( const HepGeom::Point3D<double>& orig,
567  const int bec,
568  const int layer_disk,
569  const InDetDD::SiDetectorDesign* design) const {
570 
571  HepGeom::Point3D<double> smeared(0,0,0);
572 
573  smeared.setX(orig.x());
574 
575  if (bec == 0) {
576  double smearLocY = m_random->Gaus(0, getPhiPosResolution(layer_disk));
577  double smearLocZ = m_random->Gaus(0, getEtaPosResolution(layer_disk));
578  smeared.setY(orig.y() + smearLocY);
579  smeared.setZ(orig.z() + smearLocZ);
580 
581  } else {
582 
583  smeared.setY(orig.y());
584  smeared.setZ(orig.z());
585  }
586 
587  //check for module boundaries
588  if (smeared.y()>design->width()/2) {
589  smeared.setY(design->width()/2-1e-6);
590  } else if (smeared.y()<-design->width()/2) {
591  smeared.setY(-design->width()/2+1e-6);
592  }
593  if (smeared.z()>design->length()/2) {
594  smeared.setZ(design->length()/2-1e-6);
595  } else if (smeared.z()<-design->length()/2) {
596  smeared.setZ(-design->width()/2+1e-6);
597  }
598 
599 
600  return smeared;
601 }
602 
603 double Trk::ReFitTrackWithTruth::getPhiPosResolution(int layer) const { return m_resolutionRPhi[layer]; }
604 double Trk::ReFitTrackWithTruth::getEtaPosResolution(int layer) const { return m_resolutionZ[layer]; }
605 
606 double Trk::ReFitTrackWithTruth::getPhiPosErrorFactor(int layer) const { return m_errorRPhi[layer]; }
607 double Trk::ReFitTrackWithTruth::getEtaPosErrorFactor(int layer) const { return m_errorZ[layer]; }
DataVector::reserve
void reserve(size_type n)
Attempt to preallocate enough memory for a specified number of elements.
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::LocalParameters
Definition: LocalParameters.h:98
ITrackSummaryTool.h
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
Trk::ReFitTrackWithTruth::getPhiPosErrorFactor
double getPhiPosErrorFactor(int layer) const
Definition: ReFitTrackWithTruth.cxx:606
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:27
InDetDD::DetectorDesign::width
virtual double width() const =0
Method to calculate average width of a module.
SiHit::localEndPosition
HepGeom::Point3D< double > localEndPosition() const
Definition: SiHit.cxx:153
TrackParameters.h
Trk::ReFitTrackWithTruth::ReFitTrackWithTruth
ReFitTrackWithTruth(const std::string &name, ISvcLocator *pSvcLocator)
standard Algorithm constructor
Definition: ReFitTrackWithTruth.cxx:35
MeasurementBase.h
Trk::ReFitTrackWithTruth::execute
virtual StatusCode execute() override
Definition: ReFitTrackWithTruth.cxx:112
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
Trk::ReFitTrackWithTruth::IsClusterFromTruth
static bool IsClusterFromTruth(const InDet::PixelCluster *pixClus, const int uniqueIdToMatch, const InDetSimDataCollection &sdoCollection)
Definition: ReFitTrackWithTruth.cxx:537
Trk::Track
The ATLAS Track class.
Definition: Tracking/TrkEvent/TrkTrack/TrkTrack/Track.h:73
ParticleGun_SamplingFraction.bec
int bec
Definition: ParticleGun_SamplingFraction.py:89
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
Trk::ReFitTrackWithTruth::matchSiHitsToCluster
std::vector< SiHit > matchSiHitsToCluster(const int uniqueIdToMatch, const InDet::PixelCluster *pixClus, SG::ReadHandle< AtlasHitsVector< SiHit >> &siHitCollection) const
Definition: ReFitTrackWithTruth.cxx:404
Trk::ReFitTrackWithTruth::initialize
virtual StatusCode initialize() override
Definition: ReFitTrackWithTruth.cxx:42
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
AtlasHitsVector
Definition: AtlasHitsVector.h:33
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:201
InDetDD::DetectorDesign::length
virtual double length() const =0
Method to calculate length of a module.
Trk::PrepRawData::rdoList
const std::vector< Identifier > & rdoList() const
return the List of rdo identifiers (pointers)
Trk::z0
@ z0
Definition: ParamDefs.h:64
Trk::loc2
@ loc2
generic first and second local coordinate
Definition: ParamDefs.h:35
Trk::ReFitTrackWithTruth::getEtaPosErrorFactor
double getEtaPosErrorFactor(int layer) const
Definition: ReFitTrackWithTruth.cxx:607
InDetSimDataCollection
Definition: InDetSimDataCollection.h:25
Trk::RIO_OnTrack
Definition: RIO_OnTrack.h:70
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
ParamDefs.h
InDetDD::SolidStateDetectorElementBase::identifyHash
virtual IdentifierHash identifyHash() const override final
identifier hash (inline)
SiHit
Definition: SiHit.h:19
InDetDD::SiLocalPosition
Definition: SiLocalPosition.h:31
Trk::LocalParameters::get
double get(ParamDefs par) const
Retrieve specified parameter (const version).
AtlasDetectorID.h
This class provides an interface to generate or decode an identifier for the upper levels of the dete...
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TRT::Hit::layer
@ layer
Definition: HitInfo.h:79
calibdata.exception
exception
Definition: calibdata.py:496
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:116
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Trk::ReFitTrackWithTruth::getEtaPosResolution
double getEtaPosResolution(int layer) const
Definition: ReFitTrackWithTruth.cxx:604
TrackSummary.h
Trk::ParametersBase
Definition: ParametersBase.h:55
InDet::SiCluster::detectorElement
virtual const InDetDD::SiDetectorElement * detectorElement() const override final
return the detector element corresponding to this PRD The pointer will be zero if the det el is not d...
AthAlgorithm
Definition: AthAlgorithm.h:47
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
InDetDD::SolidStateDetectorElementBase::hitLocalToLocal
Amg::Vector2D hitLocalToLocal(double xEta, double xPhi) const
Simulation/Hit local frame to reconstruction local frame.
Definition: SolidStateDetectorElementBase.cxx:95
Trk::PrepRawData::identify
Identifier identify() const
return the identifier
python.PyKernel.detStore
detStore
Definition: PyKernel.py:41
Trk::Track::perigeeParameters
const Perigee * perigeeParameters() const
return Perigee.
Definition: Tracking/TrkEvent/TrkTrack/src/Track.cxx:163
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
Trk::d0
@ d0
Definition: ParamDefs.h:63
MagicNumbers.h
RIO_OnTrack.h
InDetDD::SiDetectorElement
Definition: SiDetectorElement.h:109
DataVector::push_back
value_type push_back(value_type pElem)
Add an element to the end of the collection.
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::RIO_OnTrack::prepRawData
virtual const Trk::PrepRawData * prepRawData() const =0
returns the PrepRawData (also known as RIO) object to which this RIO_OnTrack is associated.
LocalParameters.h
python.LumiBlobConversion.pos
pos
Definition: LumiBlobConversion.py:18
Trk::MeasurementBase::localParameters
const LocalParameters & localParameters() const
Interface method to get the LocalParameters.
Definition: MeasurementBase.h:132
InDet::PixelCluster
Definition: InnerDetector/InDetRecEvent/InDetPrepRawData/InDetPrepRawData/PixelCluster.h:49
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:76
ReFitTrackWithTruth.h
CondAlgsOpts.found
int found
Definition: CondAlgsOpts.py:101
InDet::PixelClusterOnTrack
Definition: PixelClusterOnTrack.h:51
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
InDetDD::SolidStateDetectorElementBase::globalPosition
HepGeom::Point3D< double > globalPosition(const HepGeom::Point3D< double > &localPos) const
transform a reconstruction local position into a global position (inline):
Trk::RIO_OnTrack::identify
Identifier identify() const
return the identifier -extends MeasurementBase
Definition: RIO_OnTrack.h:152
DEBUG
#define DEBUG
Definition: page_access.h:11
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
InDet::SiWidth
Definition: SiWidth.h:25
Trk::ReFitTrackWithTruth::smearTruthPosition
HepGeom::Point3D< double > smearTruthPosition(const HepGeom::Point3D< double > &orig, const int bec, const int layer_disk, const InDetDD::SiDetectorDesign *design) const
Definition: ReFitTrackWithTruth.cxx:566
InDet::SiWidth::phiR
double phiR() const
Definition: SiWidth.h:126
InDetDD::SiDetectorDesign
Definition: SiDetectorDesign.h:50
get_generator_info.error
error
Definition: get_generator_info.py:40
pix
Definition: PixelMapping.cxx:16
PixelClusterOnTrack.h
InDetDD::SiDetectorElement::design
virtual const SiDetectorDesign & design() const override final
access to the local description (inline):
Trk::loc1
@ loc1
Definition: ParamDefs.h:34
error
Definition: IImpactPoint3dEstimator.h:70
Trk::ReFitTrackWithTruth::getPhiPosResolution
double getPhiPosResolution(int layer) const
Definition: ReFitTrackWithTruth.cxx:603
ITrackFitter.h
readCCLHist.float
float
Definition: readCCLHist.py:83
InDet::SiWidth::z
double z() const
Definition: SiWidth.h:131
Trk::phi0
@ phi0
Definition: ParamDefs.h:65
NSWL1::centroid
Vertex centroid(const Polygon &p)
Definition: GeoUtils.cxx:59
SiHit::localStartPosition
HepGeom::Point3D< double > localStartPosition() const
Definition: SiHit.cxx:146
match
bool match(std::string s1, std::string s2)
match the individual directories of two strings
Definition: hcg.cxx:356
Identifier
Definition: IdentifierFieldParser.cxx:14