ATLAS Offline Software
Loading...
Searching...
No Matches
InDetTrackSplitterTool.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3*/
4
9
16//#include "TrkFitterUtils/FitterTypes.h"
17#include "GaudiKernel/MsgStream.h"
21//#include "TrkParameters/MeasuredPerigee.h" //use perigee
26
29#include "GaudiKernel/ThreadLocalContext.h"
30
71
72//namespace InDetTrackSplitterHelpers{
73// /** A method to determine which half track the hit belongs to
74// */
75// bool inline CompareYPosition(Trk::MeasurementBase const* mb1, Trk::MeasurementBase const* mb2){
76// return mb1->globalPosition().y() > mb2->globalPosition().y();
77// }
78//}
79
80
81InDet::InDetTrackSplitterTool::InDetTrackSplitterTool( std::string const & type, std::string const & name, IInterface const* parent ) :
82 AthAlgTool( type, name, parent ),
83 m_trkfitter("Trk::GlobalChi2Fitter/InDetTrackFitter")
84{
85
86 declareInterface<IInDetTrackSplitterTool>(this);
87 declareProperty("TrackFitter", m_trkfitter);
88 declareProperty("OutputUpperTracksName", m_outputUpperTracksName="splitUpperTracks");
89 declareProperty("OutputLowerTracksName", m_outputLowerTracksName="splitLowerTracks");
90 declareProperty("KeepMaterialStates",m_keepmaterial=true);
91}
92
93//=========================================================
95
96
97//=========================================================
99{
100
101 msg(MSG::DEBUG) << "In initialize()" << endmsg;
102
105 if ((detStore()->retrieve(m_trtid)).isFailure()) {
106 msg(MSG::WARNING) << "Problem retrieving TRTID helper" << endmsg;
107 return StatusCode::SUCCESS;
108 }
109
112 if ((detStore()->retrieve(m_sctid)).isFailure()) {
113 msg(MSG::WARNING) << "Problem retrieving SCT ID helper" << endmsg;
114 return StatusCode::SUCCESS;
115 }
116
119 if(m_trkfitter.retrieve().isFailure()) {
120 msg(MSG::WARNING) << "Could not find refit tool " << m_trkfitter << endmsg;
121 return StatusCode::SUCCESS;
122 }
123
124 ATH_CHECK( m_outputUpperTracksName.initialize() );
125 ATH_CHECK( m_outputLowerTracksName.initialize() );
126
127 msg(MSG::DEBUG) << "InDetTrackSplitterTool initialized" << endmsg;
128 return StatusCode::SUCCESS;
129
130}
131
132//=========================================================
134{
135 msg(MSG::DEBUG) << "In finalize() "<< endmsg;
136 return StatusCode::SUCCESS;
137}
138
142std::pair<Trk::Track*, Trk::Track*> InDet::InDetTrackSplitterTool::splitInUpperLowerTrack(Trk::Track const& input,
143 bool siliconHitsOnly) const {
144 //std::cout << "input: " << input << std::endl;
145 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " In splitInUpperLowerTrack" <<endmsg;
146
147 //**The returned tracks */
148 Trk::Track* upperTrack(nullptr);
149 Trk::Track* lowerTrack(nullptr);
150
152 Trk::Perigee const* originalPerigee = input.perigeeParameters();
153
155 Trk::ParticleHypothesis hypo = input.info().particleHypothesis();
156
158 auto uppertraj = std::make_unique<Trk::TrackStates>();
159 auto lowertraj = std::make_unique<Trk::TrackStates>();
160
161 unsigned int totalNumberHits = 0;
162
163 unsigned int totalNumberPixelHits = 0;
164 unsigned int totalNumberSCTHits = 0;
165 unsigned int totalNumberTRTHits = 0;
166
167 unsigned int numberUpperPixelHits = 0;
168 unsigned int numberUpperSCTHits = 0;
169 unsigned int numberUpperTRTHits = 0;
170 unsigned int numberUpperPseudoMeas = 0;
171
172 unsigned int numberLowerPixelHits = 0;
173 unsigned int numberLowerSCTHits = 0;
174 unsigned int numberLowerTRTHits = 0;
175 unsigned int numberLowerPseudoMeas = 0;
176
177 //DataVector<Trk::MeasurementBase const>::const_iterator meas = input.measurementsOnTrack()->begin();
178 //DataVector<Trk::MeasurementBase const>::const_iterator measEnd = input.measurementsOnTrack()->end();
179 //for(;meas != measEnd; ++meas){
180 //Trk::RIO_OnTrack const* rio = dynamic_cast<Trk::RIO_OnTrack const*>(*meas);
181 Trk::TrackStates::const_iterator tsosit = input.trackStateOnSurfaces()->begin();
182 Trk::TrackStates::const_iterator tsosEnd = input.trackStateOnSurfaces()->end();
183 bool perigeeseen=false;
184 for(;tsosit != tsosEnd; ++tsosit){
185 if ((**tsosit).type(Trk::TrackStateOnSurface::Outlier)) continue;
186 if (originalPerigee==(**tsosit).trackParameters()){
187 perigeeseen=true;
188 uppertraj->push_back((**tsosit).clone());
189 lowertraj->push_back((**tsosit).clone());
190 continue;
191 }
192
193 if (m_keepmaterial &&
194 ((**tsosit).type(Trk::TrackStateOnSurface::Scatterer) ||
195 (**tsosit).type(Trk::TrackStateOnSurface::BremPoint) ||
196 (**tsosit).type(Trk::TrackStateOnSurface::CaloDeposit))) {
197 if (!perigeeseen)
198 uppertraj->push_back((**tsosit).clone());
199 else
200 lowertraj->push_back((**tsosit).clone());
201 continue;
202 }
203 const Trk::RIO_OnTrack *rio = nullptr;
204 const Trk::MeasurementBase *measb=(**tsosit).measurementOnTrack();
205 if (measb) rio=dynamic_cast<const Trk::RIO_OnTrack*>(measb);
206
207 if(rio){
208 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "an rio..." <<endmsg;
209
210 totalNumberHits++;
211
212 Identifier const& surfaceid = (rio->identify());
213
214 if (m_trtid->is_pixel(surfaceid))
215 ++totalNumberPixelHits;
216 if (m_trtid->is_sct(surfaceid))
217 ++totalNumberSCTHits;
218 if (m_trtid->is_trt(surfaceid))
219 ++totalNumberTRTHits;
220
221 //if( (*meas)->globalPosition().y() > 0){
222 if( !perigeeseen){
223
224 if (m_trtid->is_pixel(surfaceid))
225 ++numberUpperPixelHits;
226 if (m_trtid->is_sct(surfaceid))
227 ++numberUpperSCTHits;
228 if (m_trtid->is_trt(surfaceid))
229 ++numberUpperTRTHits;
230
231 if (!siliconHitsOnly || m_trtid->is_sct(surfaceid) ||
232 m_trtid->is_pixel(surfaceid))
233 uppertraj->push_back((**tsosit).clone());
234 }
235
236 // if( (*meas)->globalPosition().y() < 0){
237 if (perigeeseen) {
238
239 if (m_trtid->is_pixel(surfaceid))
240 ++numberLowerPixelHits;
241 if (m_trtid->is_sct(surfaceid))
242 ++numberLowerSCTHits;
243 if (m_trtid->is_trt(surfaceid))
244 ++numberLowerTRTHits;
245
246 if (!siliconHitsOnly || m_trtid->is_sct(surfaceid) ||
247 m_trtid->is_pixel(surfaceid))
248 // m_lowerHits.push_back( *meas);
249 lowertraj->push_back((**tsosit).clone());
250 }
251 } else {
252 if (msgLvl(MSG::VERBOSE))
253 msg(MSG::VERBOSE) << "Not an rio" << endmsg;
254 // Trk::PseudoMeasurementOnTrack const* ps =
255 // dynamic_cast<Trk::PseudoMeasurementOnTrack const*>(*measb);
257 dynamic_cast<const Trk::PseudoMeasurementOnTrack*>(measb);
258
259 if (ps) {
260 if (!perigeeseen || totalNumberHits == totalNumberTRTHits) {
261 if (msgLvl(MSG::DEBUG))
262 msg(MSG::DEBUG) << "Adding an upper pseudoMeasurement" << endmsg;
263 ++numberUpperPseudoMeas;
264 uppertraj->push_back((**tsosit).clone());
265 }
266 if (perigeeseen || totalNumberHits == totalNumberTRTHits) {
267 if (msgLvl(MSG::DEBUG))
268 msg(MSG::DEBUG) << "Adding a lower pseudoMeasurement" << endmsg;
269 ++numberLowerPseudoMeas;
270 lowertraj->push_back((**tsosit).clone());
271 }
272 }
273
274 const Trk::CompetingRIOsOnTrack* crot =
275 dynamic_cast<const Trk::CompetingRIOsOnTrack*>(measb);
276 if (crot) {
277 if (!perigeeseen)
278 uppertraj->push_back((**tsosit).clone());
279 else
280 lowertraj->push_back((**tsosit).clone());
281 }
282 }
283 }
284
285 Trk::Track upperorigtrack(input.info(),std::move(uppertraj),nullptr);
286 Trk::Track lowerorigtrack(input.info(),std::move(lowertraj),nullptr);
287
289 if(isConstrained(numberUpperPixelHits,numberUpperSCTHits,numberUpperTRTHits,numberUpperPseudoMeas)){
290
291 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "before calling upper fit" << endmsg;
292 //m_upperTrack = m_trkfitter->fit(m_upperHits, *originalPerigee, true, hypo);
293 upperTrack = (m_trkfitter->fit(Gaudi::Hive::currentContext(),upperorigtrack, false, hypo)).release();
294 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "after calling upper fit" << endmsg;
295
296 if(!upperTrack){
297 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Upper Fit Failed!" << endmsg ;
298 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
299 << numberUpperPixelHits << " upper Pixel hits, "
300 << numberUpperSCTHits << " upper SCT hits, "
301 << numberUpperTRTHits << "upper TRT hits"
302 << numberUpperPseudoMeas << "upper Pseudomeasurements"
303 << endmsg;
304
305 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
306 << totalNumberPixelHits << " total Pixel hits, "
307 << totalNumberSCTHits << " total SCT hits, "
308 << totalNumberTRTHits << "total TRT hits"
309 << endmsg;
310
311 }
312
313 }else{
314 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Not enough measurements on upper track. Fit fails." << endmsg ;
315 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
316 << numberUpperPixelHits << " upper Pixel hits, "
317 << numberUpperSCTHits << " upper SCT hits, "
318 << numberUpperTRTHits << "upper TRT hits"
319 << numberUpperPseudoMeas << "upper Pseudomeasurements"
320 << endmsg;
321 }
322
324 if(isConstrained(numberLowerPixelHits,numberLowerSCTHits,numberLowerTRTHits,numberLowerPseudoMeas)){
325
326 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "before calling lower fit" << endmsg;
327 lowerTrack = (m_trkfitter->fit(Gaudi::Hive::currentContext(),lowerorigtrack, false, hypo)).release();
328 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "after calling lower fit" << endmsg;
329
330 if(!lowerTrack){
331 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Lower Fit Failed!" << endmsg ;
332 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
333 << numberLowerPixelHits << " upper Pixel hits, "
334 << numberLowerSCTHits << " upper SCT hits, "
335 << numberLowerTRTHits << "upper TRT hits"
336 << numberLowerPseudoMeas << "upper Pseudomeasurements"
337 << endmsg;
338
339 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
340 << totalNumberPixelHits << " total Pixel hits, "
341 << totalNumberSCTHits << " total SCT hits, "
342 << totalNumberTRTHits << "total TRT hits"
343 << endmsg;
344 }
345
346 }else{
347 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Not enough measurements on lower track. Fit fails." << endmsg ;
348 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
349 << numberLowerPixelHits << " upper Pixel hits, "
350 << numberLowerSCTHits << " upper SCT hits, "
351 << numberLowerTRTHits << "upper TRT hits"
352 << numberLowerPseudoMeas << "upper Pseudomeasurements"
353 << endmsg;
354 }
355
356 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Leaving splitInUpperLowerTrack" <<endmsg;
357 return std::make_pair(upperTrack, lowerTrack);
358}
359
363 bool removeSilicon, bool applyConstraint ) const {
364 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " In stripTrack" <<endmsg;
365
367 Trk::Track* outputTrack(nullptr);
368
370 if(removeSilicon){
371 outputTrack = stripSiFromTrack(input);
372 }
374 else{
375 outputTrack = stripTRTFromTrack(input, applyConstraint);
376 }
377
378 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Leaving stripTrack" <<endmsg;
379 return outputTrack;
380}
381
385 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "In stripSiFromTrack" <<endmsg;
386
388 Trk::Track* outputTrack(nullptr);
389
391 Trk::MeasurementSet TRTHits;
392
394 Trk::Perigee const* originalPerigee = input.perigeeParameters();
395
397 Trk::ParticleHypothesis hypo = input.info().particleHypothesis();
398
400 Trk::PseudoMeasurementOnTrack const* constraint = makeThetaZ0Constraint(originalPerigee);
401
403 TRTHits.push_back(constraint);
404
405 unsigned int totalNumberTRTHits = 0;
406
407 DataVector<Trk::MeasurementBase const>::const_iterator meas = input.measurementsOnTrack()->begin();
408 DataVector<Trk::MeasurementBase const>::const_iterator measEnd = input.measurementsOnTrack()->end();
409 for(;meas != measEnd; ++meas){
410 Trk::RIO_OnTrack const* rio = dynamic_cast<Trk::RIO_OnTrack const*>(*meas);
411
412 if(!rio){
413 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "No ROT skipping measurement." <<endmsg;
414 continue;
415 }
416
417 Identifier const& surfaceid = (rio->identify());
418
419 if(!m_trtid->is_trt(surfaceid)){
420 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Removing Non-TRT hit." <<endmsg;
421 continue;
422 }
423
424 ++totalNumberTRTHits;
425 TRTHits.push_back(*meas);
426
427 }
428
430 // Should be Constrained !!!
431 //if(isConstrained(numberUpperPixelHits,numberUpperSCTHits,numberUpperTRTHits,numberUpperPseudoMeas)){
432
433 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "before calling fit" << endmsg;
434 outputTrack = (m_trkfitter->fit(Gaudi::Hive::currentContext(),TRTHits, *originalPerigee, true, hypo)).release();
435 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "after calling fit" << endmsg;
436
437 if(!outputTrack){
438 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit Failed!" << endmsg ;
439 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There were: "
440 << totalNumberTRTHits << "TRT hits"
441 << endmsg;
442 }
443
444 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving stripSiFromTrack" <<endmsg;
445 return outputTrack;
446}
447
448
452 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "In stripTRTFromTrack" <<endmsg;
453
455 Trk::Track* outputTrack(nullptr);
456
458 Trk::MeasurementSet SiHits;
459
461 Trk::Perigee const* originalPerigee = input.perigeeParameters();
462
464 Trk::ParticleHypothesis hypo = input.info().particleHypothesis();
465
467 bool addedConstraint = false;
468
469 unsigned int totalNumberSiHits = 0;
470
471 DataVector<Trk::MeasurementBase const>::const_iterator meas = input.measurementsOnTrack()->begin();
472 DataVector<Trk::MeasurementBase const>::const_iterator measEnd = input.measurementsOnTrack()->end();
473 for(;meas != measEnd; ++meas){
474 Trk::RIO_OnTrack const* rio = dynamic_cast<Trk::RIO_OnTrack const*>(*meas);
475
476 if(!rio){
477 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "No ROT skipping measurement." <<endmsg;
478 continue;
479 }
480
481 Identifier const& surfaceid = (rio->identify());
482
483 if(applyConstraint && !addedConstraint && m_trtid->is_trt(surfaceid)){
484
489 Trk::StraightLineSurface const* trtSurf = dynamic_cast<Trk::StraightLineSurface const*>(&(rio->associatedSurface()));
490 if (not trtSurf){
491 ATH_MSG_DEBUG("Cast of rio associated surface to StraightLineSurface failed.");
492 continue;
493 }
495 Trk::PseudoMeasurementOnTrack const* constraint = makePConstraint(originalPerigee,trtSurf);
496
498 SiHits.push_back(constraint);
499
501 addedConstraint = true;
502 }
503
504 if(!m_trtid->is_sct(surfaceid) && !m_trtid->is_pixel(surfaceid)){
505 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Removing Non-Si hit." <<endmsg;
506 continue;
507 }
508
509 ++totalNumberSiHits;
510 SiHits.push_back(*meas);
511
512 }
513
515 // Should be Constrained !!!
516 //if(isConstrained(numberUpperPixelHits,numberUpperSCTHits,numberUpperTRTHits,numberUpperPseudoMeas)){
517
518 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "before calling fit" << endmsg;
519 outputTrack = (m_trkfitter->fit(Gaudi::Hive::currentContext(),SiHits, *originalPerigee, false, hypo)).release();
520 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "after calling fit" << endmsg;
521
522 if(!outputTrack){
523 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Fit Failed!" << endmsg ;
524 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There were: " << totalNumberSiHits << "TRT hits" << endmsg;
525 }
526
527 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving stripTRTFromTrack" <<endmsg;
528 return outputTrack;
529}
530
534
535 if( !perigee->covariance() ) return nullptr;
536
539 Trk::DefinedParameter inputPMeasurement(perigee->parameters()[Trk::qOverP],Trk::qOverP);
540 std::array<Trk::DefinedParameter,1> constraints = {inputPMeasurement};
541
543 Amg::MatrixX constraintErrMatrix(1,1);
544 constraintErrMatrix(0,0)=(*perigee->covariance())(Trk::qOverP,Trk::qOverP);
545
548 std::move(constraintErrMatrix) ,*trtSurf);
549
550 return psmom;
551}
552
556 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "In makeThetaZ0Constraint" <<endmsg;
557
558 if( !perigee->covariance() ) return nullptr;
559
562 Trk::DefinedParameter inputTheta0Measurement(perigee->parameters()[Trk::theta],Trk::theta);
563 Trk::DefinedParameter inputZ0Measurement(perigee->parameters()[Trk::z0],Trk::z0);
564 std::array<Trk::DefinedParameter,2> constraints = {inputZ0Measurement, inputTheta0Measurement};
565
568 Amg::MatrixX constraintErrMatrix(2,2);
569 constraintErrMatrix(0,0) = (*perigee->covariance())(Trk::z0,Trk::z0);
570 constraintErrMatrix(0,1) = (*perigee->covariance())(Trk::z0,Trk::theta);
571 constraintErrMatrix(1,0) = (*perigee->covariance())(Trk::theta, Trk::z0);
572 constraintErrMatrix(1,1) = (*perigee->covariance())(Trk::theta, Trk::theta);
573
577 const Trk::PerigeeSurface & surface = perigee->associatedSurface();
578
582 std::move(constraintErrMatrix) ,
583 surface);
584
587 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving makeThetaZ0Constraint" <<endmsg;
588 return psmom;
589}
590
595std::pair<Trk::Track*, Trk::Track*> InDet::InDetTrackSplitterTool::splitInOddEvenHitsTrack(Trk::Track const& input) const {
596 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " In splitInOddEvenHitsTrack" <<endmsg;
597
599 Trk::Track* oddTrack(nullptr);
600 Trk::Track* evenTrack(nullptr);
601
603 Trk::Perigee const* originalPerigee = input.perigeeParameters();
604
606 Trk::ParticleHypothesis hypo = input.info().particleHypothesis();
607
609 Trk::MeasurementSet oddHits;
610 Trk::MeasurementSet evenHits;
611
612 unsigned int totalNumberPixelHits = 0;
613 unsigned int totalNumberSCTHits = 0;
614 unsigned int totalNumberTRTHits = 0;
615 unsigned int totalNumberHits = 0;
616
617 unsigned int numberOddSCTHits = 0;
618 unsigned int numberOddPixelHits = 0;
619 unsigned int numberOddTRTHits = 0;
620 unsigned int numberOddPseudoMeas = 0;
621
622 unsigned int numberEvenPixelHits = 0;
623 unsigned int numberEvenSCTHits = 0;
624 unsigned int numberEvenTRTHits = 0;
625 unsigned int numberEvenPseudoMeas = 0;
626
628 std::vector<Trk::MeasurementBase const*> unusedSCTHits = getSCTHits(input);
629
630 DataVector<Trk::MeasurementBase const>::const_iterator meas = input.measurementsOnTrack()->begin();
631 DataVector<Trk::MeasurementBase const>::const_iterator measEnd = input.measurementsOnTrack()->end();
632 for(;meas != measEnd; ++meas){
633 Trk::RIO_OnTrack const* rio = dynamic_cast<Trk::RIO_OnTrack const*>(*meas);
634 if(rio){
635 Identifier const& surfaceid = (rio->identify());
636 if(m_trtid->is_pixel(surfaceid)) ++totalNumberPixelHits;
637 if(m_trtid->is_sct(surfaceid)) ++totalNumberSCTHits;
638 if(m_trtid->is_trt(surfaceid)) ++totalNumberTRTHits;
639
641 if( totalNumberHits % 2 == 1){
642
643 if(m_trtid->is_pixel(surfaceid)) ++numberOddPixelHits;
644 if(m_trtid->is_sct(surfaceid)) ++numberOddSCTHits;
645 if(m_trtid->is_trt(surfaceid)) ++numberOddTRTHits;
646
647 if(!m_trtid->is_sct(surfaceid)){
648 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "adding odd hit" <<endmsg;
649 oddHits.push_back( *meas);
650 ++totalNumberHits;
651 }else{
652
655 std::vector<Trk::MeasurementBase const* >::iterator result = find(unusedSCTHits.begin(), unusedSCTHits.end(), *meas);
656
658 if(result != unusedSCTHits.end()){
659
660 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "add the (odd) SCT hit" <<endmsg;
661 oddHits.push_back( *meas);
662 ++totalNumberHits;
663
665 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "removing the SCT hit from unused list" <<endmsg;
666 unusedSCTHits.erase(result);
667
669 result = findSCTHitsFromSameSpacePoint( *meas, unusedSCTHits);
670 if(result != unusedSCTHits.end()){
671 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "adding the other (odd) SCT hit in the spacepoint" <<endmsg;
672 oddHits.push_back( *result);
673 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "removing the SCT hit from unused list" <<endmsg;
674 unusedSCTHits.erase(result);
675 }
676 }
677 }
678 }else{
679 if(m_trtid->is_pixel(surfaceid))
680 ++numberEvenPixelHits;
681 if(m_trtid->is_sct(surfaceid))
682 ++numberEvenSCTHits;
683 if(m_trtid->is_trt(surfaceid))
684 ++numberEvenTRTHits;
685
686 if(!m_trtid->is_sct(surfaceid)){
687 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "add the (even) hit" <<endmsg;
688 evenHits.push_back( *meas);
689 ++totalNumberHits;
690 }else{
693 std::vector<Trk::MeasurementBase const* >::iterator result = find(unusedSCTHits.begin(), unusedSCTHits.end(), *meas);
694
696 if(result != unusedSCTHits.end()){
697
698 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "add the (even) SCT hit" <<endmsg;
699 evenHits.push_back( *meas);
700 ++totalNumberHits;
701
703 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "removing the SCT hit from unused list" <<endmsg;
704 unusedSCTHits.erase(result);
705
707 result = findSCTHitsFromSameSpacePoint( *meas, unusedSCTHits);
708 if(result != unusedSCTHits.end()){
709 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "adding the other (even) SCT hit in the spacepoint" <<endmsg;
710 evenHits.push_back( *result);
711 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "removing the SCT hit from unused list" <<endmsg;
712 unusedSCTHits.erase(result);
713 }
714 }
715 }
716 }
717 }else{
718 Trk::PseudoMeasurementOnTrack const* ps = dynamic_cast<Trk::PseudoMeasurementOnTrack const*>(*meas);
719
720 if(ps){
721 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Adding an odd pseudoMeasurement" <<endmsg;
722 ++numberOddPseudoMeas;
723 oddHits.push_back( ps);
724
725 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Adding an even pseudoMeasurement" <<endmsg;
726 ++numberEvenPseudoMeas;
727 evenHits.push_back( ps);
728 }
729 }
730 }
731
733 //std::sort(oddHits.begin(), oddHits.end(), InDetTrackSplitterHelpers::CompareYPosition );
734 //std::sort(evenHits.begin(), evenHits.end(), InDetTrackSplitterHelpers::CompareYPosition );
735
738 if(isConstrained(numberOddPixelHits,numberOddSCTHits,numberOddTRTHits,numberOddPseudoMeas)){
739
740 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "before calling odd fit" << endmsg;
741 oddTrack = (m_trkfitter->fit(Gaudi::Hive::currentContext(),oddHits, *originalPerigee, false, hypo)).release();
742 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "after calling odd fit" << endmsg;
743
744 if(!oddTrack){
745 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Odd Fit Failed!" << endmsg ;
746 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
747 << numberOddPixelHits << " odd Pixel hits, "
748 << numberOddSCTHits << " odd SCT hits, "
749 << numberOddTRTHits << " odd TRT hits"
750 << numberOddPseudoMeas << "odd Pseudomeasurements"
751 << endmsg;
752
753 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: "
754 << totalNumberPixelHits << " total Pixel hits, "
755 << totalNumberSCTHits << " total SCT hits, "
756 << totalNumberTRTHits << "total TRT hits"
757 << endmsg;
758 }
759
760 }else
761 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Not enough measurements on odd track. Fit fails." << endmsg ;
762
764 if(isConstrained(numberEvenPixelHits,numberEvenSCTHits,numberEvenTRTHits,numberEvenPseudoMeas)){
765
766 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "before calling even fit" << endmsg;
767 evenTrack = (m_trkfitter->fit(Gaudi::Hive::currentContext(),evenHits, *originalPerigee, false, hypo)).release();
768 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "after calling even fit" << endmsg;
769
770 if(!evenTrack){
771 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Even Fit Failed!" << endmsg ;
772 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There was: " << numberEvenSCTHits << " even si hits and "
773 << numberEvenTRTHits << "even TRT hits"<< endmsg;
774
775 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "There were: " << totalNumberSCTHits << " total Si hits and "
776 << totalNumberTRTHits << "total TRT hits"<< endmsg;
777 }
778
779 }else
780 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Not enough measurements on even track. Fit fails." << endmsg ;
781
782 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Leaving splitInOddEvenHitsTrack" <<endmsg;
783 return std::make_pair(oddTrack, evenTrack);
784}
785
795 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " In splitTracks" <<endmsg;
796
797 // const EventContext& ctx = Gaudi::Hive::currentContext();
799 if ( upperTracks.record( std::make_unique<TrackCollection>() ).isFailure() ) {
800 ATH_MSG_FATAL( "Failed to record upper tracks collection " << m_outputUpperTracksName.key() );
801 }
802
804 if ( lowerTracks.record( std::make_unique<TrackCollection>() ).isFailure() ) {
805 ATH_MSG_FATAL( "Failed to record upper tracks collection " << m_outputLowerTracksName.key() );
806 }
807
808 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "There are: " << inputTracks->size() << " input tracks."<< endmsg;
809
810 TrackCollection::const_iterator it = inputTracks->begin();
811 TrackCollection::const_iterator itE = inputTracks->end();
812 for (; it != itE; ++it){
813
815 if(trackIsCandidate( **it) ){
816 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Track is a candidate."<< endmsg;
817
818 std::pair<Trk::Track*, Trk::Track*> splitTracks = this->splitInUpperLowerTrack(**it);
819 //For debugging
820 //std::pair<Trk::Track*, Trk::Track*> splitTracks = this->splitInOddEvenHitsTrack(**it);
821
822 Trk::Track* upperTrack = splitTracks.first;
823 if(upperTrack){
824 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Upper track fit suceeded"<< endmsg;
825 upperTracks->push_back( upperTrack);
826 }else
827 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Upper track fit failed!"<< endmsg;
828
829 Trk::Track* lowerTrack = splitTracks.second;
830 if(lowerTrack){
831 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Lower track fit suceeded"<< endmsg;
832 lowerTracks->push_back( lowerTrack);
833 }else
834 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Lower track fit failed!"<< endmsg;
835
836 } else
837 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "Track is NOT a candidate."<< endmsg;
838 }
839
840 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Leaving splitTracks" <<endmsg;
841}
842
848 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "In trackIsCandidate."<< endmsg;
849 Trk::Perigee const* trackPerigee = inputTrack.perigeeParameters();
850
851 if (!trackPerigee)
852 {
853 msg(MSG::WARNING) << "Found track with invalid perigee parameters. Not splitting." << endmsg;
854 return false;
855 }
856
857 //Trk::Perigee const* measuredperigee = dynamic_cast<Trk::Perigee const*>( trackPerigee );
858 double const d0 = trackPerigee->parameters()[Trk::d0];
859
860 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "do is: "<< d0 << endmsg;
861
863
864 //width of the cavity in the TRT
865 if( fabs(d0) < 600){
866 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "is a candidate" << endmsg;
867 return true;
868 }
869
870 if(msgLvl(MSG::VERBOSE)) msg(MSG::VERBOSE) << "is NOT a candidate" << endmsg;
871 return false;
872}
873
877std::vector<Trk::MeasurementBase const*> InDet::InDetTrackSplitterTool::getSCTHits(Trk::Track const& input) const {
878 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "In getSCTHits " << endmsg;
879 std::vector<Trk::MeasurementBase const*> SCTHits;
880
881 DataVector<Trk::MeasurementBase const>::const_iterator meas = input.measurementsOnTrack()->begin();
882 DataVector<Trk::MeasurementBase const>::const_iterator measEnd = input.measurementsOnTrack()->end();
883 for(;meas != measEnd; ++meas){
884
885 Trk::RIO_OnTrack const* rio = dynamic_cast<Trk::RIO_OnTrack const*>(*meas);
886
887 if(rio){
888 Identifier const& surfaceid = (rio->identify());
889 if(m_trtid->is_sct(surfaceid)){
890 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "we've found an SCT hit " << endmsg;
891 SCTHits.push_back(*meas);
892 }
893 }
894 }
895
896 //if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "sorting the hits " << endmsg;
898 //std::sort(SCTHits.begin(), SCTHits.end(), InDetTrackSplitterHelpers::CompareYPosition );
899
900 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving getSCTHits " << endmsg;
901 return SCTHits;
902}
903
908std::vector<Trk::MeasurementBase const*>::iterator InDet::InDetTrackSplitterTool::findSCTHitsFromSameSpacePoint(Trk::MeasurementBase const* sctHit, std::vector<Trk::MeasurementBase const*>& listOfSCTHits) const {
909 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "in findSCTHitsFromSameSpacePoint " << endmsg;
910 std::vector<Trk::MeasurementBase const*>::iterator result = listOfSCTHits.end();
911
912 Trk::RIO_OnTrack const* rio = dynamic_cast<Trk::RIO_OnTrack const*>(sctHit);
913 if(rio){
914 Identifier const& surfaceid = (rio->identify());
915
916 int targetEta = m_sctid->eta_module(surfaceid);
917 int targetPhi = m_sctid->phi_module(surfaceid);
918
919 std::vector<Trk::MeasurementBase const*>::const_iterator meas = listOfSCTHits.begin();
920 std::vector<Trk::MeasurementBase const*>::const_iterator measEnd = listOfSCTHits.end();
921 for(;meas != measEnd; ++meas){
922 Trk::RIO_OnTrack const* candidateRio = dynamic_cast<Trk::RIO_OnTrack const*>(*meas);
923 if(candidateRio){
924 Identifier const& candidateSurfaceid = (candidateRio->identify());
925 if(m_sctid->eta_module(candidateSurfaceid) == targetEta && m_sctid->phi_module(candidateSurfaceid) == targetPhi){
926 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Found another hit in the SpacePoint " << endmsg;
927 result = find(listOfSCTHits.begin(), listOfSCTHits.end(), *meas);
928 }
929 }
930 }
931 }
932
933 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving findSCTHitsFromSameSpacePoint " << endmsg;
934 return result;
935}
936
940bool InDet::InDetTrackSplitterTool::isConstrained(unsigned int nPixelHits, unsigned int nSCTHits, unsigned int nTRTHits, unsigned int nPseudomeasurements) const{
941 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "In isConstrained " << endmsg;
942
945 if( (( nPixelHits + nSCTHits) > 1 && (2*nPixelHits + nSCTHits + nTRTHits) > 5) ||
946 (nPseudomeasurements >= 1 && nTRTHits > 10)){
947 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving isConstrained (true)" << endmsg;
948 return true;
949 }
950
951 if(msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Leaving isConstrained (false)" << endmsg;
952 return false;
953}
954
#define endmsg
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_FATAL(x)
#define ATH_MSG_DEBUG(x)
This is an Identifier helper class for the SCT subdetector.
This is an Identifier helper class for the TRT subdetector.
DataVector< Trk::Track > TrackCollection
This typedef represents a collection of Trk::Track objects.
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
const ServiceHandle< StoreGateSvc > & detStore() const
bool msgLvl(const MSG::Level lvl) const
MsgStream & msg() const
DataModel_detail::const_iterator< DataVector > const_iterator
Definition DataVector.h:838
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
size_type size() const noexcept
Returns the number of elements in the collection.
bool m_keepmaterial
Do we use the material on the input track.
Trk::PseudoMeasurementOnTrack const * makePConstraint(Trk::Perigee const *originialPerigee, Trk::StraightLineSurface const *trtSurf) const
Make the qOverP constraint.
SG::WriteHandleKey< TrackCollection > m_outputUpperTracksName
Output track collection names recorded in storgate.
virtual std::pair< Trk::Track *, Trk::Track * > splitInUpperLowerTrack(Trk::Track const &input, bool siliconHitsOnly=false) const
Splits a single input track into upper and lower parts (based on global y) returns a pair of track th...
bool trackIsCandidate(Trk::Track const &inputTrack) const
Track selection applied in the splitTracks method (for the moment this is just a d0 cut requiring the...
Trk::Track * stripSiFromTrack(Trk::Track const &input) const
Strip the Si hits, fit the remaining with a theta, z0 constraint.
std::vector< Trk::MeasurementBase const * > getSCTHits(Trk::Track const &input) const
Return a vector of the SCT hits on track.
virtual Trk::Track * stripTrack(Trk::Track const &input, bool removeSilicon=true, bool applyConstraint=true) const
Takes a combined ID track and either 1) if removeSilicon = true Strips off the Si hits.
virtual StatusCode initialize()
standard Athena-Algorithm methods
std::vector< Trk::MeasurementBaseconst * >::iterator findSCTHitsFromSameSpacePoint(Trk::MeasurementBase const *sctHit, std::vector< Trk::MeasurementBase const * > &listOfSCTHits) const
Logic to check if there is another SCT hit associated with the input hit, which forms a space point.
InDetTrackSplitterTool(std::string const &, std::string const &, IInterface const *)
default constructor
Trk::Track * stripTRTFromTrack(Trk::Track const &input, bool applyConstraint=true) const
Strip the TRT hits, fit the remaining with a qOverP constraint.
virtual void splitTracks(TrackCollection const *inputTracks) const
Takes a trackCollection, splits them according to upper and lower parts and fills two track collectio...
bool isConstrained(unsigned int nPixelHits, unsigned int nSCTHits, unsigned int nTRTHits, unsigned int nPseudomeasurements) const
Logic to check if the track is constrained given the number of various types of hits.
SG::WriteHandleKey< TrackCollection > m_outputLowerTracksName
Trk::PseudoMeasurementOnTrack const * makeThetaZ0Constraint(Trk::Perigee const *originialPerigee) const
Make the theta and z0 constraint.
ToolHandle< Trk::ITrackFitter > m_trkfitter
Helper Functions.
virtual ~InDetTrackSplitterTool()
default destructor
virtual std::pair< Trk::Track *, Trk::Track * > splitInOddEvenHitsTrack(Trk::Track const &input) const
Splits a single input track into odd and even parts (with logic to aviod splitting SCT space points)
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
Base class for all CompetingRIOsOnTack implementations, extends the common MeasurementBase.
This class is the pure abstract base class for all fittable tracking measurements.
virtual const S & associatedSurface() const override final
Access to the Surface method.
Class describing the Line to which the Perigee refers to.
Class to handle pseudo-measurements in fitters and on track objects.
Class to handle RIO On Tracks ROT) for InDet and Muons, it inherits from the common MeasurementBase.
Definition RIO_OnTrack.h:70
virtual const Surface & associatedSurface() const override=0
returns the surface for the local to global transformation
Identifier identify() const
return the identifier -extends MeasurementBase
Class for a StraightLineSurface in the ATLAS detector to describe dirft tube and straw like detectors...
@ BremPoint
This represents a brem point on the track, and so will contain TrackParameters and MaterialEffectsBas...
@ Outlier
This TSoS contains an outlier, that is, it contains a MeasurementBase/RIO_OnTrack which was not used ...
@ Scatterer
This represents a scattering point on the track, and so will contain TrackParameters and MaterialEffe...
@ CaloDeposit
This TSOS contains a CaloEnergy object.
const Perigee * perigeeParameters() const
return Perigee.
std::string find(const std::string &s)
return a remapped string
Definition hcg.cxx:138
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
std::vector< const MeasurementBase * > MeasurementSet
vector of fittable measurements
Definition FitterTypes.h:30
ParametersT< TrackParametersDim, Charged, PerigeeSurface > Perigee
@ theta
Definition ParamDefs.h:66
@ qOverP
perigee
Definition ParamDefs.h:67
@ d0
Definition ParamDefs.h:63
@ z0
Definition ParamDefs.h:64
std::pair< double, ParamDefs > DefinedParameter
Typedef to of a std::pair<double, ParamDefs> to identify a passed-through double as a specific type o...
ParticleHypothesis
Enumeration for Particle hypothesis respecting the interaction with material.