diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointForSeedITK.h b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointForSeedITK.h
index be68fd5df794100406e5d534f8c0be1510cd6ad4..28d3ec7907d4a0e2f26ba7961516b28b18be237c 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointForSeedITK.h
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointForSeedITK.h
@@ -46,6 +46,11 @@ namespace InDet {
void set(const Trk::SpacePoint*const&,const float*,const float*);
void setQuality(float);
void setParam(const float&);
+ void setDR(const float&);
+ void setDZDR(const float&);
+ void setEta(const float&);
+ void setScorePenalty(const float&);
+ void setPt(const float&);
const Trk::SpacePoint* spacepoint{} ;
const float& x() const {return m_x;}
@@ -57,6 +62,11 @@ namespace InDet {
const float& covz() const {return m_covz;}
const float& param() const {return m_param;}
const float& quality() const {return m_q ;}
+ const float& dzdr() const {return m_dzdr;}
+ const float& eta() const {return m_eta;}
+ const float& pt() const {return m_pt;}
+ const float& scorePenalty() const {return m_scorePenalty;} /// penalty term in the seed score
+ const float& dR() const {return m_dR;} /// distance between top and central SP
const Trk::Surface* sur() const {return m_su;}
const Trk::Surface* sun() const {return m_sn;}
@@ -72,6 +82,12 @@ namespace InDet {
float m_covz{}; //
float m_param{};
float m_q{} ;
+ float m_scorePenalty{}; /// penalty term in the seed score
+ float m_dR{};
+ float m_eta{} ;
+ float m_pt{} ;
+ float m_dzdr{};
+
float m_b0[3]{};
float m_b1[3]{};
@@ -81,6 +97,7 @@ namespace InDet {
const Trk::Surface* m_su{};
const Trk::Surface* m_sn{};
};
+
} // end of name space
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointsSeedMakerEventData.h b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointsSeedMakerEventData.h
index b036fdd90f65ab2f87c93d3173b4c8fc327050ef..80531162ea40255bd1e00df46e9734b21e6ee0c9 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointsSeedMakerEventData.h
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/SiSPSeededTrackFinderData/SiSpacePointsSeedMakerEventData.h
@@ -23,6 +23,12 @@
namespace InDet {
+ class FloatInt {
+ public:
+ float Fl;
+ int In;
+ };
+
/**
@class InDet::SiSpacePointsSeedMakerEventData
@@ -95,7 +101,9 @@ namespace InDet {
float zmaxB{0.};
float ftrig{0.};
float ftrigW{0.};
- float maxScore{0.};
+ float maxScore{0.};
+ float RTmin{0.};
+ float RTmax{0.};
/**
* @name Beam geometry
@@ -136,6 +144,7 @@ namespace InDet {
std::vector<float> X;
std::vector<float> Y;
std::vector<float> Er; ///< error component on 1/tan(theta)==dz/dr from the position errors on the space-points
+ std::vector<FloatInt> Tn;
//@}
InDet::SiSpacePointsSeed seedOutput;
@@ -201,6 +210,7 @@ namespace InDet {
SP_ITK.resize(newSize, nullptr);
X.resize(newSize, 0.);
Y.resize(newSize, 0.);
+ Tn.resize(newSize);
} else {
SP.resize(newSize, nullptr);
}
diff --git a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiSpacePointForSeedITK.cxx b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiSpacePointForSeedITK.cxx
index b244b0bad3e8e8a99a86126eea70c55b58a4f8e0..f3501a90bada0861b1655553de74664a4bfaddd0 100644
--- a/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiSpacePointForSeedITK.cxx
+++ b/InnerDetector/InDetRecEvent/SiSPSeededTrackFinderData/src/SiSpacePointForSeedITK.cxx
@@ -28,6 +28,11 @@ namespace InDet {
m_q = 0.;
m_su = nullptr ;
m_sn = nullptr ;
+ m_scorePenalty = 0.;
+ m_dR = 0.;
+ m_eta = 0.;
+ m_pt = 0.;
+ m_dzdr = 0;
for(int i=0; i!=3; ++i) {m_b0[i]=0.; m_b1[i]=0.; m_dr[i]=0.; m_r0[i]=0.;}
}
@@ -35,16 +40,21 @@ namespace InDet {
(const SiSpacePointForSeedITK& sp)
{
if(&sp!=this) {
- spacepoint = sp.spacepoint;
- m_x = sp.m_x ;
- m_y = sp.m_y ;
- m_z = sp.m_z ;
- m_r = sp.m_r ;
- m_covr = sp.m_covr ;
- m_covz = sp.m_covz ;
- m_q = sp.m_q ;
- m_su = sp.m_su ;
- m_sn = sp.m_sn ;
+ spacepoint = sp.spacepoint ;
+ m_x = sp.m_x ;
+ m_y = sp.m_y ;
+ m_z = sp.m_z ;
+ m_r = sp.m_r ;
+ m_covr = sp.m_covr ;
+ m_covz = sp.m_covz ;
+ m_q = sp.m_q ;
+ m_su = sp.m_su ;
+ m_sn = sp.m_sn ;
+ m_scorePenalty = sp.m_scorePenalty ;
+ m_dR = sp.m_dR ;
+ m_eta = sp.m_eta ;
+ m_pt = sp.m_pt ;
+ m_dzdr = sp.m_dzdr ;
for(int i=0; i!=3; ++i) m_b0[i]=sp.m_b0[i];
for(int i=0; i!=3; ++i) m_b1[i]=sp.m_b1[i];
for(int i=0; i!=3; ++i) m_dr[i]=sp.m_dr[i];
@@ -123,6 +133,33 @@ namespace InDet {
m_su = &sp->clusterList().first->detectorElement()->surface();
}
+
+ void SiSpacePointForSeedITK::setDR(const float& dr)
+ {
+ m_dR = dr;
+ }
+
+ void SiSpacePointForSeedITK::setEta(const float& eta)
+ {
+ m_eta = eta;
+ }
+
+ void SiSpacePointForSeedITK::setDZDR(const float& dzdr)
+ {
+ m_dzdr = dzdr;
+ }
+
+ void SiSpacePointForSeedITK::setPt(const float& pt)
+ {
+ m_pt = pt;
+ }
+
+ void SiSpacePointForSeedITK::setScorePenalty(const float& score)
+ {
+ m_scorePenalty = score;
+ }
+
+
/////////////////////////////////////////////////////////////////////////////////
// Set with error correction
// sc[0] - barrel pixels error correction
diff --git a/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/CMakeLists.txt b/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/CMakeLists.txt
index e901b6a054ae2bcdc27249279ccfdcf6b16cae9f..e40f9bda22afe65df014bf52731179a40ea55192 100644
--- a/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/CMakeLists.txt
+++ b/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/CMakeLists.txt
@@ -1,5 +1,7 @@
# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+
# Declare the package name:
atlas_subdir( SiSpacePointsSeedTool_xk )
diff --git a/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/SiSpacePointsSeedTool_xk/SiSpacePointsSeedMaker_ITK.h b/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/SiSpacePointsSeedTool_xk/SiSpacePointsSeedMaker_ITK.h
index eaba5f80fa290fb4fac6ae2110d24904ce37f5b9..ddd97e1ed0bf4aeeaa2b5fe51db52ad53c7b989c 100644
--- a/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/SiSpacePointsSeedTool_xk/SiSpacePointsSeedMaker_ITK.h
+++ b/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/SiSpacePointsSeedTool_xk/SiSpacePointsSeedMaker_ITK.h
@@ -23,6 +23,14 @@
#include "TrkSpacePoint/SpacePointOverlapCollection.h"
#include "TrkEventUtils/PRDtoTrackMap.h"
+//for validation
+#include "GaudiKernel/ITHistSvc.h"
+#include "TFile.h"
+#include "TTree.h"
+
+
+
+
//////////////////////////////////////////////////////////////////////////////////////////////////////////////
// MagField cache
#include "MagFieldConditions/AtlasFieldCacheCondObj.h"
@@ -38,6 +46,19 @@ class MsgStream;
namespace InDet {
+ ///////////////////////////////////////////////////////////////////
+ // Object function for ordering space point in R coordinate order
+ ///////////////////////////////////////////////////////////////////
+
+ class SiSpacePointsITKComparison_R {
+
+ public:
+ bool operator () (InDet::SiSpacePointForSeedITK* s1,InDet::SiSpacePointForSeedITK* s2) {
+ return((*s1).radius() < (*s2).radius());
+ }
+ };
+
+
using EventData = SiSpacePointsSeedMakerEventData;
/**
@@ -112,7 +133,6 @@ namespace InDet {
//@}
virtual void writeNtuple(const SiSpacePointsSeed* seed, const Trk::Track* track, int seedType, long eventNumber) const override;
-
virtual bool getWriteNtupleBoolProperty() const override;
///////////////////////////////////////////////////////////////////
@@ -124,14 +144,18 @@ namespace InDet {
private:
/// enum for array sizes
- enum Size {SizeRF=53,
- SizeZ=11,
- SizeRFZ=SizeRF*SizeZ,
- SizeI=9,
- SizeRFV=100,
- SizeZV=3,
- SizeRFZV=SizeRFV*SizeZV,
- SizeIV=6};
+ /// Note that this stores the maximum capacities, the actual binnings
+ /// do not always use the full size. See data members below for the
+ /// actual binning paramaters, which are determined in buildFramework.
+ enum Size {arraySizePhi=200, ///< capacity of the 1D phi arrays
+ arraySizeZ=11, ///< capacity of the 1D z arrays
+ arraySizePhiZ=arraySizePhi*arraySizeZ, ///< capacity for the 2D phi-z arrays
+ arraySizeNeighbourBins=9, ///< array size to store neighbouring phi-z-regions in the seed finding
+ arraySizePhiV=100, ///< array size in phi for vertexing
+ arraySizeZV=3, ///< array size in z for vertexing
+ arraySizePhiZV=arraySizePhiV*arraySizeZV, ///< array size in phi-Z 2D for the vertexing
+ arraySizeNeighbourBinsVertex=6}; ///< array size to store neighbouring phi-z regions for the vertexing
+
///////////////////////////////////////////////////////////////////
// Private data and methods
@@ -166,13 +190,35 @@ namespace InDet {
FloatProperty m_zmin{this, "minZ", -250.};
FloatProperty m_zmax{this , "maxZ", +250.};
FloatProperty m_r_rmin{this, "radMin", 0.};
- FloatProperty m_r_rstep{this, "radStep", 2.};
+ FloatProperty m_binSizeR{this, "radStep", 2.};
FloatProperty m_dzver{this, "maxdZver", 5.};
FloatProperty m_dzdrver{this, "maxdZdRver", 0.02};
- FloatProperty m_diver{this, "maxdImpact", 10.};
- FloatProperty m_diversss{this, "maxdImpactSSS", 20.};
+ FloatProperty m_maxdImpact{this, "maxdImpact", 10.};
+ FloatProperty m_maxdImpactSSS{this, "maxdImpactSSS", 20.};
FloatProperty m_divermax{this, "maxdImpactForDecays", 20.};
FloatProperty m_dzmaxPPP{this, "dZmaxForPPPSeeds", 600.};
+
+ FloatProperty m_maxScore{this, "maximumAcceptedSeedScore", 100.};
+ IntegerProperty m_maxOneSizeSSS{this, "maxSeedsForSpacePointStrips", 5};
+ IntegerProperty m_maxOneSizePPP{this, "maxSeedsForSpacePointPixels", 5};
+ BooleanProperty m_alwaysKeepConfirmedPixelSeeds{this, "alwaysKeepConfirmedPixelSeeds", false};
+ BooleanProperty m_alwaysKeepConfirmedStripSeeds{this, "alwaysKeepConfirmedStripSeeds", false};
+ BooleanProperty m_fastTracking{this, "useFastTracking", false};
+ FloatProperty m_seedScoreBonusPPP{this, "seedScoreBonusPPP", -200.};
+ FloatProperty m_seedScoreBonusSSS{this, "seedScoreBonusSSS", -400.};
+ BooleanProperty m_optimisePhiBinning{this, "optimisePhiBinning", true};
+ FloatProperty m_rminSSS{this, "radMinSSS", 400.};
+ FloatProperty m_rmaxSSS{this, "radMaxSSS", 1000.};
+ BooleanProperty m_isLRT{this, "isLRT", false};
+ FloatProperty m_drminPPP{this, "mindRadiusPPP", 6.};
+ FloatProperty m_drmaxPPP{this, "maxdRadiusPPP", 140.};
+ FloatProperty m_drminSSS{this, "mindRadiusSSS", 20.};
+ FloatProperty m_drmaxSSS{this, "maxdRadiusSSS", 3000.};
+ FloatProperty m_dImpactCutSlopeUnconfirmedSSS{this, "dImpactCutSlopeUnconfirmedSSS", 1.0};
+ FloatProperty m_dImpactCutSlopeUnconfirmedPPP{this, "dImpactCutSlopeUnconfirmedPPP", 0.};
+ FloatProperty m_seedScoreBonusConfirmationSeed{this, "seedScoreBonusConfirmationSeed", -200.};
+ BooleanProperty m_useSeedConfirmation{this, "useSeedConfirmation", false};
+ FloatProperty m_rminPPPFast{this, "m_rminPPPFast", 70.};
//@}
/// @name Properties, which will be updated in initialize
@@ -211,15 +257,15 @@ namespace InDet {
//@{
bool m_initialized{false};
int m_outputlevel{0};
- int m_r_size{0};
+
int m_fNmax{0};
int m_fvNmax{0};
- int m_rfz_b[SizeRFZ]{};
- int m_rfz_t[SizeRFZ]{};
- int m_rfz_ib[SizeRFZ][SizeI]{};
- int m_rfz_it[SizeRFZ][SizeI]{};
- int m_rfzv_n[SizeRFZV]{};
- int m_rfzv_i[SizeRFZV][SizeIV]{};
+ int m_rfz_b[arraySizePhiZ]{};
+ int m_rfz_t[arraySizePhiZ]{};
+ int m_rfz_ib[arraySizePhiZ][arraySizeNeighbourBins]{};
+ int m_rfz_it[arraySizePhiZ][arraySizeNeighbourBins]{};
+ int m_rfzv_n[arraySizePhiZV]{};
+ int m_rfzv_i[arraySizePhiZV][arraySizeNeighbourBinsVertex]{};
float m_dzdrmin0{0.};
float m_dzdrmax0{0.};
float m_ipt{0.};
@@ -227,8 +273,83 @@ namespace InDet {
float m_COF{0.};
float m_sF{0.};
float m_sFv{0.};
+ float m_dzMaxFast {200.};
+ float m_R2MaxFast {2500.};
+ float m_zmaxPPP {2700.};
+ float m_rmaxPPP {140.};
+ float m_zmaxSSS {2700.};
+ float m_dzmaxSSS {900.};
+ float m_drminSeedConf{5.};
+ //@}
+
+ /// @name Properties of writeNTuple for validation purpose
+ //@{
+ Gaudi::Property<bool> m_writeNtuple {this, "WriteNtuple", false, "Flag to write Validation Ntuples"};
+ ///Flag to write validation ntuples. Turned off by default
+ ITHistSvc* m_thistSvc;
+ TTree* m_outputTree;
+ mutable std::mutex m_mutex;
+ mutable std::string m_treeName ATLAS_THREAD_SAFE;
+ mutable TString m_treeFolder ATLAS_THREAD_SAFE;
+ mutable float m_d0 ATLAS_THREAD_SAFE = 0;
+ mutable float m_z0 ATLAS_THREAD_SAFE = 0;
+ mutable float m_pt ATLAS_THREAD_SAFE = 0;
+ mutable float m_eta ATLAS_THREAD_SAFE = 0;
+ mutable double m_x1 ATLAS_THREAD_SAFE = 0;
+ mutable double m_x2 ATLAS_THREAD_SAFE = 0;
+ mutable double m_x3 ATLAS_THREAD_SAFE = 0;
+ mutable double m_y1 ATLAS_THREAD_SAFE = 0;
+ mutable double m_y2 ATLAS_THREAD_SAFE = 0;
+ mutable double m_y3 ATLAS_THREAD_SAFE = 0;
+ mutable double m_z1 ATLAS_THREAD_SAFE = 0;
+ mutable double m_z2 ATLAS_THREAD_SAFE = 0;
+ mutable double m_z3 ATLAS_THREAD_SAFE = 0;
+ mutable double m_r1 ATLAS_THREAD_SAFE = 0;
+ mutable double m_r2 ATLAS_THREAD_SAFE = 0;
+ mutable double m_r3 ATLAS_THREAD_SAFE = 0;
+ mutable float m_quality ATLAS_THREAD_SAFE = 0;
+ mutable int m_type ATLAS_THREAD_SAFE = 0;
+ mutable double m_dzdr_t ATLAS_THREAD_SAFE = 0;
+ mutable double m_dzdr_b ATLAS_THREAD_SAFE = 0;
+ mutable bool m_givesTrack ATLAS_THREAD_SAFE = 0;
+ mutable float m_trackPt ATLAS_THREAD_SAFE = 0;
+ mutable float m_trackEta ATLAS_THREAD_SAFE = 0;
+ mutable long m_eventNumber ATLAS_THREAD_SAFE = 0;
//@}
+
+ /// @name Binning parameters
+ int m_nBinsR{0}; ///< number of bins in the radial coordinate
+ int m_maxPhiBinPPP{0}; ///< number of bins in phi
+ float m_inverseBinSizePhiPPP{0}; ///< cache the inverse bin size in phi which we use - needed to evaluate phi bin locations
+ int m_maxPhiBinSSS{0};
+ float m_inverseBinSizePhiSSS{0};
+
+
+ /** Seed score thresholds defined based on the modifiers defined
+ * as configurables above.
+ * These allow to categorise the seeds based on their quality score.
+ * The modifiers above are much larger than the range of the raw
+ * unmodified scores would be, resulting in a grouping of the scores
+ * based on the modifiers applied. The thresholds are just below
+ * the value reachable without the additional modifier
+ **/
+ float m_seedScoreThresholdPPPConfirmationSeed{0.}; ///< max (score is assigned negative sign) score for PPP seeds with confirmation seed requirement.
+ float m_seedScoreThresholdSSSConfirmationSeed{0.}; ///< max (score is assigned negative sign) score for SSS seeds with confirmation seed requirement.
+
+
+ /// arrays associating bins to each other for SP formation
+ std::array<int,arraySizePhiZ> m_nNeighbourCellsBottomPPP; ///< number of neighbouring phi-z bins to consider when looking for "bottom SP" candidates for each phi-z bin
+ std::array<int,arraySizePhiZ> m_nNeighbourCellsTopPPP; ///< number of neighbouring phi-z bins to consider when looking for "top SP" candidates for each phi-z bin
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> m_neighbourCellsBottomPPP; ///< mapping of neighbour cells in the 2D phi-z binning to consider for the "bottom SP" search for central SPs in each phi-z bin. Number of valid entries stored in m_nNeighboursPhiZbottom
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> m_neighbourCellsTopPPP; ///< mapping of neighbour cells in the 2D phi-z binning to consider for the "top SP" search for central SPs in each phi-z bin. Number of valid entries stored in m_nNeighboursPhiZtop
+
+ std::array<int,arraySizePhiZ> m_nNeighbourCellsBottomSSS;
+ std::array<int,arraySizePhiZ> m_nNeighbourCellsTopSSS;
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> m_neighbourCellsBottomSSS;
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> m_neighbourCellsTopSSS;
+
+
///////////////////////////////////////////////////////////////////
// Private methods
///////////////////////////////////////////////////////////////////
@@ -243,10 +364,44 @@ namespace InDet {
MsgStream& dumpEvent(EventData& data, MsgStream& out) const;
void buildFrameWork();
+
+ void buildConnectionMaps(std::array<int, arraySizePhiZ>& nNeighbourCellsBottom,
+ std::array<int, arraySizePhiZ>& nNeighbourCellsTop,
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ>& neighbourCellsBottom,
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ>& neighbourCellsTop,
+ int maxPhiBin, bool isSSS);
+
void buildBeamFrameWork(EventData& data) const;
- SiSpacePointForSeedITK* newSpacePoint
- (EventData& data, const Trk::SpacePoint*const&) const;
+ /** Determine the expected azimuthal trajectory displacement
+ * in phi in presence of the magnetic field for a particle
+ * with momentum pTmin and impact parameter maxd0,
+ * moving from a radial coordinate Rmin outward to Rmax.
+ * This method is used to determine the optimal binning
+ * of the phi-z regions we consider in the seed making,
+ * to ensure we contain the hits from our softest tracks
+ * in a set of consecutive bins.
+ * @param[in] pTmin: minimum pt cut applied in MeV
+ * @param[in] maxD0: maximum d0 allowed
+ * @param[in] Rmin: starting radius for trajectory displacement
+ * @param[in] Rmax: end radius for trajectory displacement
+ **/
+ float azimuthalStep(const float pTmin,const float maxd0,const float Rmin,const float Rmax) const;
+
+
+ /** Create a SiSpacePointForSeed from the space point.
+ * This will also add the point to the data object's
+ * l_spforseed list and update its i_spforseed iterator
+ * to point to the entry after the new SP
+ * for further additions.
+ * Returns a nullptr if the SP fails the eta cut,
+ * should we apply one
+ * @param[in,out] data: Provides beam spot location, receives updates to the l_spforseed and i_spforseed members
+ * @param[in] sp: Input space point.
+ **/
+ SiSpacePointForSeedITK* newSpacePoint(EventData& data, const Trk::SpacePoint*const& sp) const;
+ SiSpacePointForSeedITK* newSpacePoint(EventData& data, const Trk::SpacePoint*const& sp, float* r, bool usePixSctInform=false) const;
+
void newSeed
(EventData& data,
SiSpacePointForSeedITK*&,SiSpacePointForSeedITK*&,float) const;
@@ -262,31 +417,95 @@ namespace InDet {
void fillSeeds(EventData& data) const;
void fillLists(EventData& data) const;
+ void fillListsFast(const EventContext& ctx, EventData& data) const;
+ void pixInform(const Trk::SpacePoint* sp, float* r) const;
+ void sctInform(EventData& data,const Trk::SpacePoint* sp, float* r) const;
void erase(EventData& data) const;
void production2Sp(EventData& data) const;
void production3Sp(EventData& data) const;
- void production3SpSSS
- (EventData& data,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- int,int,int&) const;
- void production3SpPPP
- (EventData& data,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- int,int,int&) const;
- void production3SpTrigger
- (EventData& data,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator*,
- int,int,int&) const;
+
+ /** \brief: Seed production from space points.
+ *
+ * This method will try to find 3-SP combinations within a
+ * local phi-z region in the detector.
+ *
+ * The central SP of the seed will be taken from this region
+ * (technically via the first entry of the bottom candidate array,
+ * which always points to the phi-z bin of interest itself).
+ *
+ * The top SP is allowed to come from the same or one of several close-by
+ * phi-Z bins, as is the bottom SP.
+ *
+ * All SP collections are expected to be internally sorted in the radial coordinate.
+ *
+ * @param[in,out] data: Event data
+ * @param[in,out] iter_bottomCands: collection of iterators over SP collections for up to 9 phi-z cells to consider for the bottom space-point search
+ * @param[in,out] iter_endBottomCands: collection of end-iterators over the
+ * SP collections for up to 9 phi-z cells to consider for the bottom space-point search
+ * @param[in,out] iter_topCands: collection of iterators over SP collections for up to 9 phi-z cells to consider for the top space-point search
+ * @param[in,out] iter_endTopCands: collection of end-iterators over the
+ * SP collections for up to 9 phi-z cells to consider for the top space-point search
+ * @param[in] numberBottomCells: Number of bottom cells to consider. Determines how many entries in iter_(end)bottomCands are expected to be valid.
+ * @param[in] numberTopCells: Number of top cells to consider.Determines how many entries in iter_(end)topCands are expected to be valid.
+ * @param[out] nseed: Number of seeds found
+ **/
+ void production3SpSSS
+ (EventData& data,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_bottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_endBottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_topCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_endTopCands,
+ const int numberBottomCells, const int numberTopCells, int& nseed) const;
+
+ void production3SpPPP
+ (EventData& data,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_bottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_endBottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_topCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & iter_endTopCands,
+ const int numberBottomCells, const int numberTopCells, int& nseed) const;
+
+ /// as above, but for the trigger
+ void production3SpTrigger
+ (EventData& /*data*/,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & /*rb*/,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & /*rbe*/,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & /*rt*/,
+ std::array<std::list<InDet::SiSpacePointForSeedITK*>::iterator, arraySizeNeighbourBins> & /*rte*/,
+ const int /*numberBottomCells*/, const int /*numberTopCells*/, int& /*nseed*/) const;
+
+ /** This creates all possible seeds with the passed central and bottom SP, using all top SP
+ * candidates which are stored in the data.CmSp member. Seeds are scored by a quality score
+ * seeded by abs(d0), and modified if there is a second-seed confirmation or in case of PPP/SSS
+ * topologies. Then, they are written out via the newOneSeed method.
+ * @param[in] SPb Bottom Space point for the seed creation
+ * @param[in] SP0 Central Space point for the seed creation
+ * @param[in] Zob z0 estimate
+ **/
+ void newOneSeedWithCurvaturesComparisonSSS
+ (EventData& data, SiSpacePointForSeedITK*& SPb, SiSpacePointForSeedITK*& SP0, float Zob) const;
+ void newOneSeedWithCurvaturesComparisonPPP
+ (EventData& data, SiSpacePointForSeedITK*& SPb, SiSpacePointForSeedITK*& SP0, float Zob) const;
+ void newOneSeedWithCurvaturesComparisonSeedConfirmation
+ (EventData& data, SiSpacePointForSeedITK*& SPb, SiSpacePointForSeedITK*& SP0, float Zob) const;
+
+ /** Helper method to determine if a seed
+ * is 'confirmed' - this means that a second
+ * seed exists with compatible curvature,
+ * the same bottom and central SP, but a
+ * different third SP. This information
+ * is stored in a modification of the
+ * seed quality, which we check here.
+ * @param[in] bottomSP: bottom space point
+ * @param[in] topSP: top space point
+ * @param[in] quality: seed quality
+ * @return true if the seed is confirmed, false otherwise
+ **/
+ bool isConfirmedSeed(const InDet::SiSpacePointForSeedITK* bottomSP, const InDet::SiSpacePointForSeedITK* topSP, float quality) const;
+
+
+ void sort(std::vector<FloatInt>& s, int start, int size) const;
bool newVertices(EventData& data, const std::list<Trk::Vertex>&) const;
void findNext(EventData& data) const;
bool isZCompatible(EventData& data, float&,float&,float&) const;
@@ -329,7 +548,20 @@ namespace InDet {
return false;
}
-}
+///////////////////////////////////////////////////////////////////
+ // The procedure sorts the elements into ascending order.
+ ///////////////////////////////////////////////////////////////////
+
+ inline
+ void SiSpacePointsSeedMaker_ITK::sort(std::vector<FloatInt>& s, int start, int size) const
+ {
+ //QuickSort for fast tracking currently buggy
+ //TBC if really faster than std::sort
+ //Using std::sort in all cases for now
+ std::sort(s.begin()+start,s.begin()+start+size,[](const FloatInt a,const FloatInt b)->bool {return a.Fl < b.Fl;});
+ }
+
+}
#endif // SiSpacePointsSeedMaker_ITK_H
diff --git a/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/src/SiSpacePointsSeedMaker_ITK.cxx b/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/src/SiSpacePointsSeedMaker_ITK.cxx
index 5bb089e57196c1681ba098d1601f32a4906c0602..8e4f02d77b81e0d629d0074125c7a2c12ebebf67 100644
--- a/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/src/SiSpacePointsSeedMaker_ITK.cxx
+++ b/InnerDetector/InDetRecTools/SiSpacePointsSeedTool_xk/src/SiSpacePointsSeedMaker_ITK.cxx
@@ -1,6 +1,6 @@
/*
- Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
-*/
+ Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+ */
///////////////////////////////////////////////////////////////////
// Implementation file for class SiSpacePointsSeedMaker_ITK
@@ -16,19 +16,33 @@
#include "InDetPrepRawData/SiCluster.h"
+
+//for validation
+#include "TrkTrack/Track.h"
+#include "TrkParameters/TrackParameters.h"
+
#include <cmath>
#include <iomanip>
#include <ostream>
+
///////////////////////////////////////////////////////////////////
// Constructor
///////////////////////////////////////////////////////////////////
-InDet::SiSpacePointsSeedMaker_ITK::SiSpacePointsSeedMaker_ITK
-(const std::string& t, const std::string& n, const IInterface* p)
- : base_class(t, n, p)
+InDet::SiSpacePointsSeedMaker_ITK::SiSpacePointsSeedMaker_ITK(const std::string &t, const std::string &n, const IInterface *p)
+ : base_class(t, n, p),
+ m_thistSvc(nullptr),
+ m_outputTree(nullptr),
+ m_treeName(""),
+ m_treeFolder("/valNtuples/")
{
+ if (m_maxOneSize > 0)
+ {
+ m_maxOneSizeSSS = m_maxOneSize;
+ m_maxOneSizePPP = m_maxOneSize;
+ }
}
///////////////////////////////////////////////////////////////////
@@ -47,29 +61,68 @@ StatusCode InDet::SiSpacePointsSeedMaker_ITK::initialize()
//
ATH_CHECK(m_beamSpotKey.initialize());
- ATH_CHECK( m_fieldCondObjInputKey.initialize());
+ ATH_CHECK(m_fieldCondObjInputKey.initialize());
// PRD-to-track association (optional)
- ATH_CHECK( m_prdToTrackMap.initialize( !m_prdToTrackMap.key().empty()));
+ ATH_CHECK(m_prdToTrackMap.initialize(!m_prdToTrackMap.key().empty()));
- if (m_r_rmax < 1100.) m_r_rmax = 1100.;
-
// Build framework
//
buildFrameWork();
// Get output print level
//
- m_outputlevel = msg().level()-MSG::DEBUG;
- if (m_outputlevel<=0) {
+ m_outputlevel = msg().level() - MSG::DEBUG;
+
+ if (msgLvl(MSG::DEBUG)) {
EventData data;
initializeEventData(data);
data.nprint=0;
- dump(data, msg(MSG::DEBUG));
+ dump(data, msg(MSG::DEBUG));
+ }
+
+ m_umax = 100. - std::abs(m_umax) * 300.;
+
+ if (m_writeNtuple) {
+
+ ATH_CHECK( service("THistSvc",m_thistSvc) );
+
+ m_treeName = (std::string("SeedTree_")+name());
+ std::replace( m_treeName.begin(), m_treeName.end(), '.', '_' );
+
+ m_outputTree = new TTree( m_treeName.c_str() , "SeedMakerValTool");
+
+ m_outputTree->Branch("eventNumber", &m_eventNumber);
+ m_outputTree->Branch("d0", &m_d0);
+ m_outputTree->Branch("z0", &m_z0);
+ m_outputTree->Branch("pt", &m_pt);
+ m_outputTree->Branch("eta", &m_eta);
+ m_outputTree->Branch("x1", &m_x1);
+ m_outputTree->Branch("x2", &m_x2);
+ m_outputTree->Branch("x3", &m_x3);
+ m_outputTree->Branch("y1", &m_y1);
+ m_outputTree->Branch("y2", &m_y2);
+ m_outputTree->Branch("y3", &m_y3);
+ m_outputTree->Branch("z1", &m_z1);
+ m_outputTree->Branch("z2", &m_z2);
+ m_outputTree->Branch("z3", &m_z3);
+ m_outputTree->Branch("r1", &m_r1);
+ m_outputTree->Branch("r2", &m_r2);
+ m_outputTree->Branch("r3", &m_r3);
+ m_outputTree->Branch("quality", &m_quality);
+ m_outputTree->Branch("seedType", &m_type);
+ m_outputTree->Branch("givesTrack", &m_givesTrack);
+ m_outputTree->Branch("dzdr_b", &m_dzdr_b);
+ m_outputTree->Branch("dzdr_t", &m_dzdr_t);
+ m_outputTree->Branch("track_pt", &m_trackPt);
+ m_outputTree->Branch("track_eta", &m_trackEta);
+
+ TString fullTreeName = m_treeFolder + m_treeName;
+
+ ATH_CHECK( m_thistSvc->regTree( fullTreeName.Data(), m_outputTree ) );
+
}
- m_umax = 100.-std::abs(m_umax)*300.;
- m_initialized = true;
return sc;
}
@@ -84,157 +137,287 @@ StatusCode InDet::SiSpacePointsSeedMaker_ITK::finalize()
}
///////////////////////////////////////////////////////////////////
-// Initialize tool for new event
+// Initialize tool for new event
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::newEvent(const EventContext& ctx, EventData& data, int iteration) const
+void InDet::SiSpacePointsSeedMaker_ITK::newEvent(const EventContext &ctx, EventData &data, int iteration) const
{
- if (not data.initialized) initializeEventData(data);
- data.iteration0 = iteration;
+ /// if not done so, book the arrays etc inside the event data object
+ if (not data.initialized)
+ initializeEventData(data);
+
data.trigger = false;
- if (!m_pixel && !m_sct) return;
+ if (!m_pixel && !m_sct)
+ return;
- iteration <=0 ? data.iteration = 0 : data.iteration = iteration;
+ /// pass the iteration info into our data object
+ data.iteration = iteration;
+ if (iteration <= 0)
+ data.iteration = 0;
+ /// Erase any existing entries in the data object
erase(data);
+
data.dzdrmin = m_dzdrmin0;
data.dzdrmax = m_dzdrmax0;
+ data.maxScore = m_maxScore; ///< max score, where low scores are "better".
+
+ /// in the first iteration, initialise the beam framework - beam spot position and direction
+ if (data.iteration == 0)
+ {
- if (!data.iteration) {
buildBeamFrameWork(data);
- double f[3], gP[3] ={10.,10.,0.};
+ /// Read the field information
+ double magField[3]{0, 0, 0};
+ double globalPos[3] = {10., 10., 0.};
- MagField::AtlasFieldCache fieldCache;
+ MagField::AtlasFieldCache fieldCache;
- // Get field cache object
SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
- const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
- if (fieldCondObj == nullptr) {
+ const AtlasFieldCacheCondObj *fieldCondObj{*readHandle};
+ if (fieldCondObj == nullptr)
+ {
ATH_MSG_ERROR("SiSpacePointsSeedMaker_ITK: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
return;
}
- fieldCondObj->getInitializedCache (fieldCache);
-
- if (fieldCache.solenoidOn()) {
- fieldCache.getFieldZR(gP, f);
-
- data.K = 2./(300.*f[2]);
- } else {
- data.K = 2./(300.* 5. );
+ fieldCondObj->getInitializedCache(fieldCache);
+
+ if (fieldCache.solenoidOn())
+ {
+ /// retrieve field
+ fieldCache.getFieldZR(globalPos, magField);
+ /**
+ * Knowing the field (note that the field cache returns the field in units of kiloTesla!)
+ * allows to set the circle-radius to pT conversion factor.
+ *
+ * See for example ATLAS-CONF-2010-072
+ * R[mm] =pT[GeV] / (3·10−4×B[T]) = pT[MeV] / (300 *Bz[kT])
+ *
+ * We actually estimate the circle diameter, 2R, in the seeding.
+ * So what we want is: 2R = pT[MeV] x 2 / (300 x Bz) = K x pT[MeV].
+ **/
+ data.K = 2. / (300. * magField[2]);
}
-
- data.ipt2K = m_ipt2/(data.K*data.K);
- data.ipt2C = m_ipt2*m_COF;
- data.COFK = m_COF*(data.K*data.K);
+ else
+ {
+ data.K = 2. / (300. * 5.);
+ }
+ /** helper variables allowing us to directly apply our pt cut on the variables
+ * available at seed level.
+ * ipt2K is 1 / (K * 0.9 * pt cut)²
+ **/
+ data.ipt2K = m_ipt2 / (data.K * data.K);
+ /// related to the mysterious magic number, m_COF{134*.05*9}
+ data.ipt2C = m_ipt2 * m_COF;
+ data.COFK = m_COF * (data.K * data.K);
+
+ /// set the spacepoint iterator to the beginning of the space-point list
data.i_spforseed_ITK = data.l_spforseed_ITK.begin();
- } else {
- data.r_first = 0;
+ // Set the seed multiplicity strategy of the event data to the one configured
+ // by the user for strip seeds
+ data.maxSeedsPerSP = m_maxOneSizeSSS;
+ data.keepAllConfirmedSeeds = m_alwaysKeepConfirmedStripSeeds;
+ } ///< end if-statement for iteration 0
+ else
+ { /// for the second iteration (PPP pass), don't redo the full init required the first time
+ data.r_first = 0; ///< reset the first radial bin
+ // Set the seed multiplicity strategy of the event data to the one configured
+ // by the user for pixel seeds
+ data.maxSeedsPerSP = m_maxOneSizePPP;
+ data.keepAllConfirmedSeeds = m_alwaysKeepConfirmedPixelSeeds;
+
+ /// call fillLists to repopulate the candidate space points and exit
fillLists(data);
return;
}
+ if (m_fastTracking)
+ {
+ fillListsFast(ctx, data);
+ return;
+ }
+
+ /// the following will only happen in the first iteration
data.checketa = data.dzdrmin > 1.;
- float irstep = 1.f/m_r_rstep;
- int irmax = m_r_size-1;
- for (int i=0; i<data.nr; ++i) {
+ /// build the r-binning.
+ float oneOverBinSizeR = 1. / m_binSizeR;
+ int maxBinR = m_nBinsR - 1;
+
+ /** This cleans up remaining entries in the data object.
+ * In standard execution, we only run this in the first
+ * iterations on a newly created data object,
+ * in which case this loop will not ever be entered.
+ * Leaving it in place for nonstandard use cases.
+ **/
+ for (int i = 0; i < data.nr; ++i)
+ {
int n = data.r_index[i];
data.r_map[n] = 0;
data.r_Sorted_ITK[n].clear();
}
data.ns = data.nr = 0;
- SG::ReadHandle<Trk::PRDtoTrackMap> prd_to_track_map;
+ SG::ReadHandle<Trk::PRDtoTrackMap> prd_to_track_map;
const Trk::PRDtoTrackMap *prd_to_track_map_cptr = nullptr;
- if (!m_prdToTrackMap.key().empty()) {
- prd_to_track_map=SG::ReadHandle<Trk::PRDtoTrackMap>(m_prdToTrackMap, ctx);
- if (!prd_to_track_map.isValid()) {
+ if (!m_prdToTrackMap.key().empty())
+ {
+ prd_to_track_map = SG::ReadHandle<Trk::PRDtoTrackMap>(m_prdToTrackMap, ctx);
+ if (!prd_to_track_map.isValid())
+ {
ATH_MSG_ERROR("Failed to read PRD to track association map: " << m_prdToTrackMap.key());
}
prd_to_track_map_cptr = prd_to_track_map.cptr();
}
- // Get pixels space points containers from store gate
- //
+ /////////////////////////////
+ /// Now, we will populate the space point list in the event data object once.
+ /////////////////////////////
+
+ /// reset the first r index
data.r_first = 0;
- if (m_pixel) {
+ /// Get pixels space points containers from store gate
+ if (m_pixel)
+ {
SG::ReadHandle<SpacePointContainer> spacepointsPixel{m_spacepointsPixel, ctx};
- if (spacepointsPixel.isValid()) {
-
- for (const SpacePointCollection* spc: *spacepointsPixel) {
- for (const Trk::SpacePoint* sp: *spc) {
-
- if ((prd_to_track_map_cptr && isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax || sp->r() < m_r_rmin ) continue;
-
- InDet::SiSpacePointForSeedITK* sps = newSpacePoint(data, sp);
- if (!sps) continue;
-
- int ir = static_cast<int>(sps->radius()*irstep);
- if (ir>irmax) ir = irmax;
- data.r_Sorted_ITK[ir].push_back(sps);
- ++data.r_map[ir];
- if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
- if (ir > data.r_first) data.r_first = ir;
- ++data.ns;
- }
- }
- }
- ++data.r_first;
- }
-
- // Get sct space points containers from store gate
- //
- if (m_sct) {
+ if (spacepointsPixel.isValid())
+ {
+ /// loop over the pixel space points
+ for (const SpacePointCollection *spc : *spacepointsPixel)
+ {
+ for (const Trk::SpacePoint *sp : *spc)
+ {
+
+ /// if we use the PRD to track map and this SP has already been used in a track, bail out
+ /// also skip any SP outside the r binning
+ if ((prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr)) || sp->r() > m_r_rmax || sp->r() < m_r_rmin)
+ continue;
+
+ /// Remove DBM space points
+ ///
+ const InDetDD::SiDetectorElement *de =
+ static_cast<const InDetDD::SiDetectorElement *>(sp->clusterList().first->detectorElement());
+ if (!de || de->isDBM())
+ continue;
+
+ /** create a SiSpacePointForSeed from the space point.
+ * This will also add the point to the l_spforseed list and update
+ * the i_spforseed iterator
+ **/
+ InDet::SiSpacePointForSeedITK *sps = newSpacePoint(data, sp);
+ /// this can occur if we fail the eta cut
+ if (!sps)
+ continue;
+
+ /// determine the r-bin of this SP.
+ /// done by dividing the radius by the bin size.
+ int radiusBin = static_cast<int>(sps->radius() * oneOverBinSizeR);
+ /// catch outliers
+ if (radiusBin > maxBinR)
+ radiusBin = maxBinR;
+
+ /// now add the SP to the r-binned vector
+ data.r_Sorted_ITK[radiusBin].push_back(sps);
+ /// increment the counter for this bin
+ ++data.r_map[radiusBin];
+ /// if this is the first time we see this bin in use, we update the index map for this bin
+ /// to the radius bin index
+ if (data.r_map[radiusBin] == 1)
+ data.r_index[data.nr++] = radiusBin;
+ /// if this is the highest bin we saw so far, update the r_first member of the data object to this bin
+ if (radiusBin > data.r_first)
+ data.r_first = radiusBin;
+ /// update the space point counter
+ ++data.ns;
+ } ///< end loop over space points in collection
+ } ///< end loop over pixel SP collections
+ } ///< end if-statement on valid pixel SP container
+ ++data.r_first; ///< increment r_first past the last occupied bin we saw
+ } ///< end pixel case
+
+ /// Get sct space points containers from store gate
+ if (m_sct)
+ {
SG::ReadHandle<SpacePointContainer> spacepointsSCT{m_spacepointsSCT, ctx};
- if (spacepointsSCT.isValid()) {
-
- for (const SpacePointCollection* spc: *spacepointsSCT) {
- for (const Trk::SpacePoint* sp: *spc) {
-
- if ((prd_to_track_map_cptr && isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax || sp->r() < m_r_rmin ) continue;
-
- InDet::SiSpacePointForSeedITK* sps = newSpacePoint(data, sp);
- if (!sps) continue;
-
- int ir = static_cast<int>(sps->radius()*irstep);
- if (ir>irmax) ir = irmax;
- data.r_Sorted_ITK[ir].push_back(sps);
- ++data.r_map[ir];
- if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
- ++data.ns;
- }
+ if (spacepointsSCT.isValid())
+ {
+
+ for (const SpacePointCollection *spc : *spacepointsSCT)
+ {
+ for (const Trk::SpacePoint *sp : *spc)
+ {
+ /// as for the pixel, veto already used SP if we are using the PRD to track map in later passes of track finding.
+ /// Also, veto SP outside the maximum and minimum radii
+ if ((prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr)) || sp->r() > m_r_rmax || sp->r() < m_r_rmin)
+ continue;
+ /// create a space point and write it into the data object's list of points
+ InDet::SiSpacePointForSeedITK *sps = newSpacePoint(data, sp);
+ if (!sps)
+ continue;
+
+ /// as for PIX, determine the radial bin.
+ /// Note that for the SCT we do not update data.r_first.
+ int radiusBin = static_cast<int>(sps->radius() * oneOverBinSizeR);
+ if (radiusBin > maxBinR)
+ radiusBin = maxBinR;
+ /// again store the SP in the r-binned vectors
+ data.r_Sorted_ITK[radiusBin].push_back(sps);
+ /// update the count of SP in the given bin
+ ++data.r_map[radiusBin];
+ /// update the r_index map and data.nr if needed
+ if (data.r_map[radiusBin] == 1)
+ data.r_index[data.nr++] = radiusBin;
+ /// and increment the SP count too.
+ ++data.ns;
+ }
}
}
- // Get sct overlap space points containers from store gate
- //
- if (m_useOverlap) {
+ /// Get sct overlap space points containers from store gate
+ if (m_useOverlap && !data.checketa)
+ {
SG::ReadHandle<SpacePointOverlapCollection> spacepointsOverlap{m_spacepointsOverlap, ctx};
- if (spacepointsOverlap.isValid()) {
-
- for (const Trk::SpacePoint* sp: *spacepointsOverlap) {
-
- if ((prd_to_track_map_cptr && isUsed(sp,*prd_to_track_map_cptr)) || sp->r() > m_r_rmax || sp->r() < m_r_rmin) continue;
-
- InDet::SiSpacePointForSeedITK* sps = newSpacePoint(data, sp);
- if (!sps) continue;
-
- int ir = static_cast<int>(sps->radius()*irstep);
- if (ir>irmax) ir = irmax;
- data.r_Sorted_ITK[ir].push_back(sps);
- ++data.r_map[ir];
- if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
- ++data.ns;
- }
+ if (spacepointsOverlap.isValid())
+ {
+ for (const Trk::SpacePoint *sp : *spacepointsOverlap)
+ {
+ /// usual rejection of SP used in previous track finding passes if we run with the PRT to track map + check of the max and min radii
+ if ((prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr)) || sp->r() > m_r_rmax || sp->r() < m_r_rmin)
+ continue;
+
+ /// SP creation, entry into list of the data object
+ InDet::SiSpacePointForSeedITK *sps = newSpacePoint(data, sp);
+ if (!sps)
+ continue;
+
+ /// radial bin determination
+ int radiusBin = static_cast<int>(sps->radius() * oneOverBinSizeR);
+ if (radiusBin > maxBinR)
+ radiusBin = maxBinR;
+ /// insert into the "histogram" vector
+ data.r_Sorted_ITK[radiusBin].push_back(sps);
+ /// update the counter for each bin content
+ ++data.r_map[radiusBin];
+ /// update the bin index list and occupied bin counter
+ if (data.r_map[radiusBin] == 1)
+ data.r_index[data.nr++] = radiusBin;
+ /// and the total SP count too.
+ ++data.ns;
+ }
}
}
}
- if (iteration < 0) data.r_first = 0;
+ /// negative iterations are not used in the current ITk reco
+ if (iteration < 0)
+ data.r_first = 0;
+
+ /// populates the phi-z sorted histograms using the spacepoint lists and r-binning.
+ /// after this call, we have a 3D binning
fillLists(data);
}
@@ -242,112 +425,151 @@ void InDet::SiSpacePointsSeedMaker_ITK::newEvent(const EventContext& ctx, EventD
// Initialize tool for new region
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::newRegion
-(const EventContext& ctx, EventData& data,
- const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT) const
+void InDet::SiSpacePointsSeedMaker_ITK::newRegion(const EventContext &ctx, EventData &data,
+ const std::vector<IdentifierHash> &vPixel, const std::vector<IdentifierHash> &vSCT) const
{
- if (not data.initialized) initializeEventData(data);
+ if (not data.initialized)
+ initializeEventData(data);
data.iteration = 0;
data.trigger = false;
erase(data);
- if (!m_pixel && !m_sct) return;
+ if (!m_pixel && !m_sct)
+ return;
data.dzdrmin = m_dzdrmin0;
data.dzdrmax = m_dzdrmax0;
+ data.maxScore = m_maxScore;
buildBeamFrameWork(data);
- double f[3], gP[3] ={10.,10.,0.};
+ double magField[3]{0, 0, 0};
+ double globalPos[3] = {10., 10., 0.};
- MagField::AtlasFieldCache fieldCache;
+ MagField::AtlasFieldCache fieldCache;
// Get field cache object
SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey, ctx};
- const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
- if (fieldCondObj == nullptr) {
+ const AtlasFieldCacheCondObj *fieldCondObj{*readHandle};
+ if (fieldCondObj == nullptr)
+ {
ATH_MSG_ERROR("SiSpacePointsSeedMaker_ITK: Failed to retrieve AtlasFieldCacheCondObj with key " << m_fieldCondObjInputKey.key());
return;
}
- fieldCondObj->getInitializedCache (fieldCache);
+ fieldCondObj->getInitializedCache(fieldCache);
- if (fieldCache.solenoidOn()) {
- fieldCache.getFieldZR(gP, f);
+ if (fieldCache.solenoidOn())
+ {
+ fieldCache.getFieldZR(globalPos, magField);
- data.K = 2./(300.*f[2]);
- } else {
- data.K = 2./(300.* 5. );
+ data.K = 2. / (300. * magField[2]);
+ }
+ else
+ {
+ data.K = 2. / (300. * 5.);
}
- data.ipt2K = m_ipt2/(data.K*data.K);
- data.ipt2C = m_ipt2*m_COF;
- data.COFK = m_COF*(data.K*data.K);
+ data.ipt2K = m_ipt2 / (data.K * data.K);
+ data.ipt2C = m_ipt2 * m_COF;
+ data.COFK = m_COF * (data.K * data.K);
data.i_spforseed_ITK = data.l_spforseed_ITK.begin();
- float irstep = 1.f/m_r_rstep;
- int irmax = m_r_size-1;
+ float oneOverBinSizeR = 1. / m_binSizeR; //was float irstep = 1.f/m_binSizeR;
+ int maxBinR = m_nBinsR - 1; //was int irmax = m_nBinsR-1;
data.r_first = 0;
data.checketa = false;
- for (int i=0; i<data.nr; ++i) {
+ for (int i = 0; i < data.nr; ++i)
+ {
int n = data.r_index[i];
data.r_map[n] = 0;
data.r_Sorted_ITK[n].clear();
}
data.ns = data.nr = 0;
- // Get pixels space points containers from store gate
+ SG::ReadHandle<Trk::PRDtoTrackMap> prd_to_track_map;
+ const Trk::PRDtoTrackMap *prd_to_track_map_cptr = nullptr;
+ if (!m_prdToTrackMap.key().empty())
+ {
+ prd_to_track_map = SG::ReadHandle<Trk::PRDtoTrackMap>(m_prdToTrackMap, ctx);
+ if (!prd_to_track_map.isValid())
+ {
+ ATH_MSG_ERROR("Failed to read PRD to track association map: " << m_prdToTrackMap.key());
+ }
+ prd_to_track_map_cptr = prd_to_track_map.cptr();
+ }
+
+ // Get pixels space points containers from store gate
//
- if (m_pixel && !vPixel.empty()) {
+ if (m_pixel && !vPixel.empty())
+ {
SG::ReadHandle<SpacePointContainer> spacepointsPixel{m_spacepointsPixel, ctx};
- if (spacepointsPixel.isValid()) {
+ if (spacepointsPixel.isValid())
+ {
+ SpacePointContainer::const_iterator spce = spacepointsPixel->end();
// Loop through all trigger collections
//
- for (const IdentifierHash& l: vPixel) {
- const auto *w = spacepointsPixel->indexFindPtr(l);
- if (w==nullptr) continue;
- for (const Trk::SpacePoint* sp: *w) {
- float r = sp->r();
- if (r > m_r_rmax || r < m_r_rmin) continue;
- InDet::SiSpacePointForSeedITK* sps = newSpacePoint(data, sp);
- int ir = static_cast<int>(sps->radius()*irstep);
- if (ir>irmax) ir = irmax;
- data.r_Sorted_ITK[ir].push_back(sps);
+ for (const IdentifierHash &l : vPixel)
+ {
+ auto w = spacepointsPixel->indexFind(l);
+ if (w == spce)
+ continue;
+ for (const Trk::SpacePoint *sp : **w)
+ {
+ float r = sp->r();
+ if ((prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr)) || r > m_r_rmax || r < m_r_rmin)
+ continue;
+ InDet::SiSpacePointForSeedITK *sps = newSpacePoint(data, sp);
+ int ir = static_cast<int>(sps->radius() * oneOverBinSizeR);
+ if (ir > maxBinR)
+ ir = maxBinR;
+ data.r_Sorted_ITK[ir].push_back(sps);
++data.r_map[ir];
- if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
- ++data.ns;
- }
+ if (data.r_map[ir] == 1)
+ data.r_index[data.nr++] = ir;
+ ++data.ns;
+ }
}
}
}
- // Get sct space points containers from store gate
+ // Get sct space points containers from store gate
//
- if (m_sct && !vSCT.empty()) {
+ if (m_sct && !vSCT.empty())
+ {
SG::ReadHandle<SpacePointContainer> spacepointsSCT{m_spacepointsSCT, ctx};
- if (spacepointsSCT.isValid()) {
+ if (spacepointsSCT.isValid())
+ {
+
+ SpacePointContainer::const_iterator spce = spacepointsSCT->end();
// Loop through all trigger collections
//
- for (const IdentifierHash& l: vSCT) {
- const auto *w = spacepointsSCT->indexFindPtr(l);
- if (w==nullptr) continue;
- for (const Trk::SpacePoint* sp: *w) {
- float r = sp->r();
- if (r > m_r_rmax || r < m_r_rmin) continue;
- InDet::SiSpacePointForSeedITK* sps = newSpacePoint(data, sp);
- int ir = static_cast<int>(sps->radius()*irstep);
- if (ir>irmax) ir = irmax;
- data.r_Sorted_ITK[ir].push_back(sps);
+ for (const IdentifierHash &l : vSCT)
+ {
+ auto w = spacepointsSCT->indexFind(l);
+ if (w == spce)
+ continue;
+ for (const Trk::SpacePoint *sp : **w)
+ {
+ float r = sp->r();
+ if ((prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr)) || r > m_r_rmax || r < m_r_rmin)
+ continue;
+ InDet::SiSpacePointForSeedITK *sps = newSpacePoint(data, sp);
+ int ir = static_cast<int>(sps->radius() * oneOverBinSizeR);
+ if (ir > maxBinR)
+ ir = maxBinR;
+ data.r_Sorted_ITK[ir].push_back(sps);
++data.r_map[ir];
- if (data.r_map[ir]==1) data.r_index[data.nr++] = ir;
- ++data.ns;
- }
+ if (data.r_map[ir] == 1)
+ data.r_index[data.nr++] = ir;
+ ++data.ns;
+ }
}
}
}
@@ -358,62 +580,70 @@ void InDet::SiSpacePointsSeedMaker_ITK::newRegion
// Initialize tool for new region
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::newRegion
-(const EventContext& ctx, EventData& data,
- const std::vector<IdentifierHash>& vPixel, const std::vector<IdentifierHash>& vSCT, const IRoiDescriptor& IRD) const
+void InDet::SiSpacePointsSeedMaker_ITK::newRegion(const EventContext &ctx, EventData &data,
+ const std::vector<IdentifierHash> &vPixel, const std::vector<IdentifierHash> &vSCT, const IRoiDescriptor &IRD) const
{
- if (not data.initialized) initializeEventData(data);
+ constexpr float twoPi = 2. * M_PI;
+
+ if (not data.initialized)
+ initializeEventData(data);
newRegion(ctx, data, vPixel, vSCT);
data.trigger = true;
- double dzdrmin = 1.f/std::tan(2.f*std::atan(std::exp(-IRD.etaMinus())));
- double dzdrmax = 1.f/std::tan(2.f*std::atan(std::exp(-IRD.etaPlus ())));
-
- data.zminB = IRD.zedMinus()-data.zbeam[0]; // min bottom Z
- data.zmaxB = IRD.zedPlus ()-data.zbeam[0]; // max bottom Z
- data.zminU = data.zminB+550.*dzdrmin;
- data.zmaxU = data.zmaxB+550.*dzdrmax;
- double fmax = IRD.phiPlus ();
- double fmin = IRD.phiMinus();
- if (fmin > fmax) fmin-=(2.*M_PI);
- data.ftrig = (fmin+fmax)*.5;
- data.ftrigW = (fmax-fmin)*.5;
+ double dzdrmin = 1.f / std::tan(2.f * std::atan(std::exp(-IRD.etaMinus())));
+ double dzdrmax = 1.f / std::tan(2.f * std::atan(std::exp(-IRD.etaPlus())));
+
+ data.zminB = IRD.zedMinus() - data.zbeam[0]; // min bottom Z
+ data.zmaxB = IRD.zedPlus() - data.zbeam[0]; // max bottom Z
+ data.zminU = data.zminB + 550. * dzdrmin;
+ data.zmaxU = data.zmaxB + 550. * dzdrmax;
+ double fmax = IRD.phiPlus();
+ double fmin = IRD.phiMinus();
+ if (fmin > fmax)
+ fmin -= twoPi;
+ data.ftrig = (fmin + fmax) * .5;
+ data.ftrigW = (fmax - fmin) * .5;
}
+
///////////////////////////////////////////////////////////////////
// Methods to initilize different strategies of seeds production
// with two space points with or without vertex constraint
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::find2Sp(EventData& data, const std::list<Trk::Vertex>& lv) const
+void InDet::SiSpacePointsSeedMaker_ITK::find2Sp(EventData &data, const std::list<Trk::Vertex> &lv) const
{
- if (not data.initialized) initializeEventData(data);
+ if (not data.initialized)
+ initializeEventData(data);
data.zminU = m_zmin;
data.zmaxU = m_zmax;
int mode = 0;
- if (lv.begin()!=lv.end()) mode = 1;
+ if (lv.begin() != lv.end())
+ mode = 1;
bool newv = newVertices(data, lv);
-
- if (newv || !data.state || data.nspoint!=2 || data.mode!=mode || data.nlist) {
+
+ if (newv || !data.state || data.nspoint != 2 || data.mode != mode || data.nlist)
+ {
data.i_seede_ITK = data.l_seeds_ITK.begin();
- data.state = 1;
+ data.state = 1;
data.nspoint = 2;
- data.nlist = 0;
- data.mode = mode;
+ data.nlist = 0;
+ data.mode = mode;
data.endlist = true;
- data.fvNmin = 0;
- data.fNmin = 0;
- data.zMin = 0;
+ data.fvNmin = 0;
+ data.fNmin = 0;
+ data.zMin = 0;
production2Sp(data);
}
data.i_seed_ITK = data.l_seeds_ITK.begin();
-
- if (m_outputlevel<=0) {
- data.nprint=1;
+
+ if (m_outputlevel <= 0)
+ {
+ data.nprint = 1;
dump(data, msg(MSG::DEBUG));
}
}
@@ -423,34 +653,45 @@ void InDet::SiSpacePointsSeedMaker_ITK::find2Sp(EventData& data, const std::list
// with three space points with or without vertex constraint
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
+void InDet::SiSpacePointsSeedMaker_ITK::find3Sp(const EventContext &, EventData &data, const std::list<Trk::Vertex> &lv) const
{
- if (not data.initialized) initializeEventData(data);
+ if (not data.initialized)
+ initializeEventData(data);
+ /// reset the Z interval stored in the data object
data.zminU = m_zmin;
data.zmaxU = m_zmax;
-
+ /// mode 2 if we have no vertices in the list, otherwise mode 3
int mode = 2;
- if (lv.begin()!=lv.end()) mode = 3;
+ if (lv.begin() != lv.end())
+ mode = 3;
+ /** copy the vertices into the data object, if we have any.
+ * Note that by construction, newv will ALWAYS be false, also if we
+ * pass vertices.
+ **/
bool newv = newVertices(data, lv);
-
- if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
+ /// update the data object's config
+ if (newv || !data.state || data.nspoint != 3 || data.mode != mode || data.nlist)
+ {
data.i_seede_ITK = data.l_seeds_ITK.begin();
- data.state = 1;
+ data.state = 1;
data.nspoint = 3;
- data.nlist = 0;
- data.mode = mode;
+ data.nlist = 0;
+ data.mode = mode;
data.endlist = true;
- data.fvNmin = 0;
- data.fNmin = 0;
- data.zMin = 0;
- production3Sp(data);
+ data.fvNmin = 0;
+ data.fNmin = 0;
+ data.zMin = 0;
+ production3Sp(data); ///< This performs the actual seed finding
}
+
+ /// reset the i_seed_ITK iterator - this is used to return the seeds to the
+ /// consumer when they call next()
data.i_seed_ITK = data.l_seeds_ITK.begin();
- data.seed_ITK = data.seeds_ITK.begin();
- if (m_outputlevel<=0) {
- data.nprint=1;
+ if (msgLvl(MSG::DEBUG))
+ {
+ data.nprint = 1;
dump(data, msg(MSG::DEBUG));
}
}
@@ -460,36 +701,51 @@ void InDet::SiSpacePointsSeedMaker_ITK::find3Sp(const EventContext&, EventData&
// with three space points with or without vertex constraint
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::find3Sp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv, const double* ZVertex) const
+void InDet::SiSpacePointsSeedMaker_ITK::find3Sp(const EventContext &, EventData &data, const std::list<Trk::Vertex> &lv, const double *ZVertex) const
{
- if (not data.initialized) initializeEventData(data);
+ if (not data.initialized)
+ initializeEventData(data);
+
+ /// Update the data object's Z interval based on the interval passed as arg to this
+ /// function.
data.zminU = ZVertex[0];
- if (data.zminU < m_zmin) data.zminU = m_zmin;
+ if (data.zminU < m_zmin)
+ data.zminU = m_zmin; ///< don't go beyond user-specified intervals
data.zmaxU = ZVertex[1];
- if (data.zmaxU > m_zmax) data.zmaxU = m_zmax;
+ if (data.zmaxU > m_zmax)
+ data.zmaxU = m_zmax; ///< don't go beyond user-specified intervals
+ /// mode 2 when working with the Z constraint
int mode = 2;
- if (lv.begin()!=lv.end()) mode = 3;
+ if (lv.begin() != lv.end())
+ mode = 3;
+ /** copy the vertices into the data object, if we have any.
+ * Note that by construction, newv will ALWAYS be false, also if we
+ * pass vertices.
+ **/
bool newv = newVertices(data, lv);
-
- if (newv || !data.state || data.nspoint!=3 || data.mode!=mode || data.nlist) {
+ /// update the data object's config
+ if (newv || !data.state || data.nspoint != 3 || data.mode != mode || data.nlist)
+ {
data.i_seede_ITK = data.l_seeds_ITK.begin();
- data.state = 1;
+ data.state = 1;
data.nspoint = 3;
- data.nlist = 0;
- data.mode = mode;
+ data.nlist = 0;
+ data.mode = mode;
data.endlist = true;
- data.fvNmin = 0;
- data.fNmin = 0;
- data.zMin = 0;
- production3Sp(data);
+ data.fvNmin = 0;
+ data.fNmin = 0;
+ data.zMin = 0;
+ production3Sp(data); ///< This performs the actual seed finding
}
+ /// reset the i_seed_ITK iterator - this is used to return the seeds to the
+ /// consumer when they call next()
data.i_seed_ITK = data.l_seeds_ITK.begin();
- data.seed_ITK = data.seeds_ITK.begin();
- if (m_outputlevel<=0) {
- data.nprint=1;
+ if (msgLvl(MSG::DEBUG))
+ {
+ data.nprint = 1;
dump(data, msg(MSG::DEBUG));
}
}
@@ -500,35 +756,38 @@ void InDet::SiSpacePointsSeedMaker_ITK::find3Sp(const EventContext&, EventData&
// Variable means (2,3,4,....) any number space points
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::findVSp(const EventContext&, EventData& data, const std::list<Trk::Vertex>& lv) const
+void InDet::SiSpacePointsSeedMaker_ITK::findVSp(const EventContext &, EventData &data, const std::list<Trk::Vertex> &lv) const
{
- if (not data.initialized) initializeEventData(data);
+
+ if (not data.initialized)
+ initializeEventData(data);
data.zminU = m_zmin;
data.zmaxU = m_zmax;
int mode = 5;
- if (lv.begin()!=lv.end()) mode = 6;
+ if (lv.begin() != lv.end())
+ mode = 6;
bool newv = newVertices(data, lv);
-
- if (newv || !data.state || data.nspoint!=4 || data.mode!=mode || data.nlist) {
+ if (newv || !data.state || data.nspoint != 4 || data.mode != mode || data.nlist)
+ {
data.i_seede_ITK = data.l_seeds_ITK.begin();
- data.state = 1;
+ data.state = 1;
data.nspoint = 4;
- data.nlist = 0;
- data.mode = mode;
+ data.nlist = 0;
+ data.mode = mode;
data.endlist = true;
- data.fvNmin = 0;
- data.fNmin = 0;
- data.zMin = 0;
+ data.fvNmin = 0;
+ data.fNmin = 0;
+ data.zMin = 0;
production3Sp(data);
}
data.i_seed_ITK = data.l_seeds_ITK.begin();
- data.seed_ITK = data.seeds_ITK.begin();
- if (m_outputlevel<=0) {
- data.nprint=1;
+ if (msgLvl(MSG::DEBUG))
+ {
+ data.nprint = 1;
dump(data, msg(MSG::DEBUG));
}
}
@@ -537,11 +796,13 @@ void InDet::SiSpacePointsSeedMaker_ITK::findVSp(const EventContext&, EventData&
// Dumps relevant information into the MsgStream
///////////////////////////////////////////////////////////////////
-MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dump(EventData& data, MsgStream& out) const
+MsgStream &InDet::SiSpacePointsSeedMaker_ITK::dump(EventData &data, MsgStream &out) const
{
- if (not data.initialized) initializeEventData(data);
+ if (not data.initialized)
+ initializeEventData(data);
- if (data.nprint) return dumpEvent(data, out);
+ if (data.nprint)
+ return dumpEvent(data, out);
return dumpConditions(data, out);
}
@@ -549,17 +810,24 @@ MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dump(EventData& data, MsgStream& o
// Dumps conditions information into the MsgStream
///////////////////////////////////////////////////////////////////
-MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpConditions(EventData& data, MsgStream& out) const
+MsgStream &InDet::SiSpacePointsSeedMaker_ITK::dumpConditions(EventData &data, MsgStream &out) const
{
int n = 42-m_spacepointsPixel.key().size();
- std::string s2; for (int i=0; i<n; ++i) s2.append(" "); s2.append("|");
+ std::string s2;
+ for (int i=0; i<n; ++i) s2.append(" ");
+ s2.append("|");
n = 42-m_spacepointsSCT.key().size();
- std::string s3; for (int i=0; i<n; ++i) s3.append(" "); s3.append("|");
+ std::string s3;
+ for (int i=0; i<n; ++i) s3.append(" ");
+ s3.append("|");
n = 42-m_spacepointsOverlap.key().size();
- std::string s4; for (int i=0; i<n; ++i) s4.append(" "); s4.append("|");
+ std::string s4;
+ for (int i=0; i<n; ++i) s4.append(" ");
+ s4.append("|");
n = 42-m_beamSpotKey.key().size();
- std::string s5; for (int i=0; i<n; ++i) s5.append(" "); s5.append("|");
-
+ std::string s5;
+ for (int i=0; i<n; ++i) s5.append(" ");
+ s5.append("|");
out<<"|---------------------------------------------------------------------|"
<<endmsg;
@@ -586,17 +854,11 @@ MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpConditions(EventData& data, Ms
out<<"| pTmin (mev) | "
<<std::setw(12)<<std::setprecision(5)<<m_ptmin
<<" |"<<endmsg;
- out<<"| |rapidity| <= | "
- <<std::setw(12)<<std::setprecision(5)<<m_rapcut
- <<" |"<<endmsg;
out<<"| max radius SP | "
<<std::setw(12)<<std::setprecision(5)<<m_r_rmax
<<" |"<<endmsg;
- out<<"| min radius SP | "
- <<std::setw(12)<<std::setprecision(5)<<m_r_rmin
- <<" |"<<endmsg;
out<<"| radius step | "
- <<std::setw(12)<<std::setprecision(5)<<m_r_rstep
+ <<std::setw(12)<<std::setprecision(5)<<m_binSizeR
<<" |"<<endmsg;
out<<"| min Z-vertex position | "
<<std::setw(12)<<std::setprecision(5)<<m_zmin
@@ -604,42 +866,18 @@ MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpConditions(EventData& data, Ms
out<<"| max Z-vertex position | "
<<std::setw(12)<<std::setprecision(5)<<m_zmax
<<" |"<<endmsg;
- out<<"| min radius first SP(3) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r1min
- <<" |"<<endmsg;
- out<<"| min radius second SP(3) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r2min
- <<" |"<<endmsg;
- out<<"| min radius last SP(3) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r3min
- <<" |"<<endmsg;
- out<<"| max radius first SP(3) | "
- <<std::setw(12)<<std::setprecision(4)<<m_r1max
- <<" |"<<endmsg;
- out<<"| max radius second SP(3) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r2max
- <<" |"<<endmsg;
- out<<"| max radius last SP(3) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r3max
- <<" |"<<endmsg;
- out<<"| min radius first SP(2) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r1minv
- <<" |"<<endmsg;
- out<<"| min radius second SP(2) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r2minv
- <<" |"<<endmsg;
- out<<"| max radius first SP(2) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r1maxv
- <<" |"<<endmsg;
- out<<"| max radius second SP(2) | "
- <<std::setw(12)<<std::setprecision(5)<<m_r2maxv
- <<" |"<<endmsg;
- out<<"| min space points dR | "
- <<std::setw(12)<<std::setprecision(5)<<m_drmin
- <<" |"<<endmsg;
- out<<"| max space points dR | "
- <<std::setw(12)<<std::setprecision(5)<<m_drmax
- <<" |"<<endmsg;
+ out<<"| min space points dR SSS | "
+ <<std::setw(12)<<std::setprecision(5)<<m_drminSSS
+ <<" |"<<std::endl;
+ out<<"| max space points dR SSS | "
+ <<std::setw(12)<<std::setprecision(5)<<m_drmaxSSS
+ <<" |"<<std::endl;
+ out<<"| min space points dR PPP | "
+ <<std::setw(12)<<std::setprecision(5)<<m_drminPPP
+ <<" |"<<std::endl;
+ out<<"| max space points dR PPP | "
+ <<std::setw(12)<<std::setprecision(5)<<m_drmaxPPP
+ <<" |"<<std::endl;
out<<"| max dZ impact | "
<<std::setw(12)<<std::setprecision(5)<<m_dzver
<<" |"<<endmsg;
@@ -647,13 +885,10 @@ MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpConditions(EventData& data, Ms
<<std::setw(12)<<std::setprecision(5)<<m_dzdrver
<<" |"<<endmsg;
out<<"| max impact | "
- <<std::setw(12)<<std::setprecision(5)<<m_diver
- <<" |"<<endmsg;
- out<<"| max impact pps | "
- <<std::setw(12)<<std::setprecision(5)<<m_diverpps
+ <<std::setw(12)<<std::setprecision(5)<<m_maxdImpact
<<" |"<<endmsg;
out<<"| max impact sss | "
- <<std::setw(12)<<std::setprecision(5)<<m_diversss
+ <<std::setw(12)<<std::setprecision(5)<<m_maxdImpactSSS
<<" |"<<endmsg;
out<<"|---------------------------------------------------------------------|"
<<endmsg;
@@ -684,23 +919,25 @@ MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpConditions(EventData& data, Ms
out<<"|---------------------------------------------------------------------|"
<<endmsg;
return out;
+
+
}
///////////////////////////////////////////////////////////////////
// Dumps event information into the MsgStream
///////////////////////////////////////////////////////////////////
-MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpEvent(EventData& data, MsgStream& out) const
+MsgStream &InDet::SiSpacePointsSeedMaker_ITK::dumpEvent(EventData &data, MsgStream &out) const
{
out<<"|---------------------------------------------------------------------|"
<<endmsg;
- out<<"| data.ns | "
+ out<<"| ns | "
<<std::setw(12)<<data.ns
<<" |"<<endmsg;
- out<<"| data.nsaz | "
+ out<<"| nsaz | "
<<std::setw(12)<<data.nsaz
<<" |"<<endmsg;
- out<<"| data.nsazv | "
+ out<<"| nsazv | "
<<std::setw(12)<<data.nsazv
<<" |"<<endmsg;
out<<"| seeds | "
@@ -709,51 +946,65 @@ MsgStream& InDet::SiSpacePointsSeedMaker_ITK::dumpEvent(EventData& data, MsgStre
out<<"|---------------------------------------------------------------------|"
<<endmsg;
return out;
+
}
///////////////////////////////////////////////////////////////////
// Find next set space points
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::findNext(EventData& data) const
+void InDet::SiSpacePointsSeedMaker_ITK::findNext(EventData &data) const
{
- if (data.endlist) return;
+ if (data.endlist)
+ return;
data.i_seede_ITK = data.l_seeds_ITK.begin();
- if (data.mode==0 || data.mode==1) production2Sp(data);
- else if (data.mode==2 || data.mode==3) production3Sp(data);
- else if (data.mode==5 || data.mode==6) production3Sp(data);
+ if (data.mode == 2 || data.mode == 3 || data.mode == 5 || data.mode == 6)
+ production3Sp(data);
data.i_seed_ITK = data.l_seeds_ITK.begin();
- data.seed_ITK = data.seeds_ITK.begin();
++data.nlist;
-}
+}
///////////////////////////////////////////////////////////////////
// New and old list vertices comparison
///////////////////////////////////////////////////////////////////
-bool InDet::SiSpacePointsSeedMaker_ITK::newVertices(EventData& data, const std::list<Trk::Vertex>& lV) const
+bool InDet::SiSpacePointsSeedMaker_ITK::newVertices(EventData &data, const std::list<Trk::Vertex> &lV) const
{
+
unsigned int s1 = data.l_vertex.size();
unsigned int s2 = lV.size();
+ /// reset the isvertex flag
data.isvertex = false;
- if (s1==0 && s2==0) return false;
+ /// if we had no vertices before and have none now,
+ /// we can exit right away
+ if (s1 == 0 && s2 == 0)
+ return false;
+ /// clean up the vertex list
data.l_vertex.clear();
- if (s2 == 0) return false;
+ /// if we have no vertices now, we can exit
+ if (s2 == 0)
+ return false;
+ /// otherwise, update the data with the new vertices
data.isvertex = true;
- for (const Trk::Vertex& v: lV) {
+ for (const Trk::Vertex &v : lV)
+ {
data.l_vertex.insert(static_cast<float>(v.position().z()));
}
- data.zminU = (*data.l_vertex. begin())-20.;
- if ( data.zminU < m_zmin) data.zminU = m_zmin;
- data.zmaxU = (*data.l_vertex.rbegin())+20.;
- if ( data.zmaxU > m_zmax) data.zmaxU = m_zmax;
+ /// and also update the z interval, adding 20mm before/after the first/last vertex in z
+ /// make sure not to extend the interval beyond the user-configured z interval.
+ data.zminU = (*data.l_vertex.begin()) - 20.;
+ if (data.zminU < m_zmin)
+ data.zminU = m_zmin;
+ data.zmaxU = (*data.l_vertex.rbegin()) + 20.;
+ if (data.zmaxU > m_zmax)
+ data.zmaxU = m_zmax;
return false;
}
@@ -762,330 +1013,734 @@ bool InDet::SiSpacePointsSeedMaker_ITK::newVertices(EventData& data, const std::
// Initiate frame work for seed generator
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::buildFrameWork()
+void InDet::SiSpacePointsSeedMaker_ITK::buildFrameWork()
{
m_ptmin = std::abs(m_ptmin);
-
- if (m_ptmin < 100.) m_ptmin = 100.;
-
- if (m_diversss < m_diver ) m_diversss = m_diver;
- if (m_divermax < m_diversss) m_divermax = m_diversss;
-
- if (std::abs(m_etamin) < .1) m_etamin = -m_etamax;
- m_dzdrmax0 = 1.f/std::tan(2.f*std::atan(exp(-m_etamax)));
- m_dzdrmin0 = 1.f/std::tan(2.f*std::atan(exp(-m_etamin)));
-
- m_COF = 134*.05*9.;
- m_ipt = 1.f/std::abs(.9f*m_ptmin);
- m_ipt2 = m_ipt*m_ipt;
-
- // Build radius sorted containers
- //
- m_r_size = static_cast<int>((m_r_rmax+.1)/m_r_rstep);
-
- // Build radius-azimuthal sorted containers
- //
- constexpr float pi2 = 2.*M_PI;
- const int NFmax = SizeRF;
- const float sFmax = static_cast<float>(NFmax)/pi2;
- const float sFmin = 100./60.;
- float ptm = 400.;
- if (m_ptmin < ptm) ptm = m_ptmin;
+ if (m_ptmin < 100.)
+ m_ptmin = 100.;
+
+ /// ensure consistency in the transverse IP cuts
+ if (m_maxdImpactSSS < m_maxdImpact)
+ m_maxdImpactSSS = m_maxdImpact;
+
+ /// symmetrise eta cut if eta min not explicitly set
+ if (std::abs(m_etamin) < .1)
+ m_etamin = -m_etamax;
+ /// set dz/dr cut values based on eta cuts
+ m_dzdrmax0 = 1. / std::tan(2. * std::atan(std::exp(-m_etamax)));
+ m_dzdrmin0 = 1. / std::tan(2. * std::atan(std::exp(-m_etamin)));
+
+ /// cache inverse pt cuts
+ m_ipt = 1. / std::abs(m_ptmin);
+ m_ipt2 = m_ipt * m_ipt;
+
+ /// set up the score thresholds based on the user-supplied properties
+ /// The score of any seeds will always be >= the bonus applied,
+ /// since the starting value is |d0|.
+ /// Hence, by subtracting one, we get all seeds which have
+ /// received an additional bonus in addition to the PPP/SSS one,
+ /// which is the confirmation one by construction.
+ m_seedScoreThresholdPPPConfirmationSeed = m_seedScoreBonusPPP - 1.;
+ m_seedScoreThresholdSSSConfirmationSeed = m_seedScoreBonusSSS - 1.;
+
+ /// Build radius sorted containers
+ m_nBinsR = static_cast<int>((m_r_rmax + .1) / m_binSizeR);
+
+ /**
+ * Now we construct the radius-azimuthal sorted containers.
+ * Binning is determined semi-dynamically based on the tool config
+ **/
+
+ /// determine the phi binning
+ constexpr float twoPi = 2. * M_PI;
+
+ /// The max Nb. of bins possible is given by the binning enum
+ const int nPhiBinsMax = arraySizePhi;
+ const float inverseSizePhiMax = static_cast<float>(nPhiBinsMax) / twoPi; ///< for run-3: 200 : 6.28 ~ 31.8, bin size 0.0314
+ constexpr float inverseSizePhiMin = 10. / twoPi;
+
+ /// derive the optimum bin size in phi for our needs.
+ /// We do this by checking the size needed for each of the passes we run,
+ /// and then pick the larger one (more conservative).
+ /// This is slightly less optimal than using separate settings per pass, but avoids
+ /// having to book and store in memory two distinct R-phi-Z maps.
+
+ if (m_optimisePhiBinning)
+ {
+ /// case 1: PPP seeds, if we use them
+ const float radiusPixelStart = m_fastTracking ? 50. : 40.; /// approximate lowest R location of pixel hits (driven by barrel)
+ const float radiusPixelEnd = m_fastTracking ? 250. : 320.; /// approximate largest R location of pixel hits (driven by endcap)
+ const float binSizePhi_PPP = m_pixel ? azimuthalStep(m_ptmin, m_maxdImpact, radiusPixelStart, radiusPixelEnd) / 3. : 0.;
+ /// case 2: SSS seeds, if we use them
+ const float binSizePhi_SSS = m_sct ? azimuthalStep(m_ptmin, m_maxdImpactSSS, m_rminSSS, m_rmaxSSS) / 3. : 0.;
+ m_inverseBinSizePhiPPP = 1. / binSizePhi_PPP;
+ m_inverseBinSizePhiSSS = 1. / binSizePhi_SSS;
+ }
+ else
+ {
+ /// this is the default phi binning as operated in release 21 - optimised for
+ /// a trajectory with 400 MeV, from the origin, and Rmin = 0 / Rmax = 600mm float ptm = 400.;
+ float ptm = 400.;
+ /// if we cut below 400 MeV, adapt the ptm
+ if (m_ptmin < ptm)
+ ptm = m_ptmin;
+ m_inverseBinSizePhiPPP = m_inverseBinSizePhiSSS = ptm / 60.;
+ }
- m_sF = ptm /60.;
- if (m_sF > sFmax) m_sF = sFmax;
- else if (m_sF < sFmin) m_sF = sFmin;
- m_fNmax = static_cast<int>(pi2*m_sF);
- if (m_fNmax >=NFmax) m_fNmax = NFmax-1;
+ /// truncate the bin size to fall within our thresholds
+ if (m_inverseBinSizePhiPPP > inverseSizePhiMax)
+ m_inverseBinSizePhiPPP = inverseSizePhiMax;
+ else if (m_inverseBinSizePhiPPP < inverseSizePhiMin)
+ m_inverseBinSizePhiPPP = inverseSizePhiMin;
+ if (m_inverseBinSizePhiSSS > inverseSizePhiMax)
+ m_inverseBinSizePhiSSS = inverseSizePhiMax;
+ else if (m_inverseBinSizePhiSSS < inverseSizePhiMin)
+ m_inverseBinSizePhiSSS = inverseSizePhiMin;
+
+ /// now we can determine the number of bins by dividing the interval by the bin size
+ m_maxPhiBinPPP = static_cast<int>(twoPi * m_inverseBinSizePhiPPP);
+ /// additional protection against too many bins. Should not happen given constraints above
+ if (m_maxPhiBinPPP >= nPhiBinsMax)
+ m_maxPhiBinPPP = nPhiBinsMax - 1;
+ m_maxPhiBinSSS = static_cast<int>(twoPi * m_inverseBinSizePhiSSS);
+ if (m_maxPhiBinSSS >= nPhiBinsMax)
+ m_maxPhiBinSSS = nPhiBinsMax - 1;
+ /// recompute inverse bin size, taking into account rounding + truncation
+ m_inverseBinSizePhiPPP = m_maxPhiBinPPP / twoPi;
+ m_inverseBinSizePhiSSS = m_maxPhiBinSSS / twoPi;
+
+ buildConnectionMaps(m_nNeighbourCellsBottomPPP, m_nNeighbourCellsTopPPP,
+ m_neighbourCellsBottomPPP, m_neighbourCellsTopPPP,
+ m_maxPhiBinPPP, false);
+
+ buildConnectionMaps(m_nNeighbourCellsBottomSSS, m_nNeighbourCellsTopSSS,
+ m_neighbourCellsBottomSSS, m_neighbourCellsTopSSS,
+ m_maxPhiBinSSS, true);
+}
- // Build radius-azimuthal-Z sorted containers for Z-vertices
- //
- const int NFtmax = SizeRFV;
- const float sFvmax = static_cast<float>(NFtmax)/pi2;
- m_sFv = m_ptmin/120.;
- if (m_sFv>sFvmax) m_sFv = sFvmax;
- m_fvNmax = static_cast<int>(pi2*m_sFv);
- if (m_fvNmax>=NFtmax) m_fvNmax = NFtmax-1;
-
- // Build maps for radius-azimuthal-Z sorted collections
- //
- for (int f=0; f<=m_fNmax; ++f) {
+void InDet::SiSpacePointsSeedMaker_ITK::buildConnectionMaps(std::array<int, arraySizePhiZ> &nNeighbourCellsBottom,
+ std::array<int, arraySizePhiZ> &nNeighbourCellsTop,
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> &neighbourCellsBottom,
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> &neighbourCellsTop,
+ int maxPhiBin, bool isSSS)
+{
- int fb = f-1; if (fb<0 ) fb=m_fNmax;
- int ft = f+1; if (ft>m_fNmax) ft=0;
-
- // For each azimuthal region loop through all Z regions
- //
- for (int z=0; z<SizeZ; ++z) {
-
- int a = f *SizeZ+z;
- int b = fb*SizeZ+z;
- int c = ft*SizeZ+z;
- m_rfz_b [a] = 3; m_rfz_t [a] = 3;
- m_rfz_ib[a][0] = a; m_rfz_it[a][0] = a;
- m_rfz_ib[a][1] = b; m_rfz_it[a][1] = b;
- m_rfz_ib[a][2] = c; m_rfz_it[a][2] = c;
- if (z==5) {
-
- m_rfz_t [a] = 9;
- m_rfz_it[a][3] = a+1;
- m_rfz_it[a][4] = b+1;
- m_rfz_it[a][5] = c+1;
- m_rfz_it[a][6] = a-1;
- m_rfz_it[a][7] = b-1;
- m_rfz_it[a][8] = c-1;
- } else if (z> 5) {
-
- m_rfz_b [a] = 6;
- m_rfz_ib[a][3] = a-1;
- m_rfz_ib[a][4] = b-1;
- m_rfz_ib[a][5] = c-1;
-
- if (z<10) {
-
- m_rfz_t [a] = 6;
- m_rfz_it[a][3] = a+1;
- m_rfz_it[a][4] = b+1;
- m_rfz_it[a][5] = c+1;
- }
- } else {
-
- m_rfz_b [a] = 6;
- m_rfz_ib[a][3] = a+1;
- m_rfz_ib[a][4] = b+1;
- m_rfz_ib[a][5] = c+1;
-
- if (z>0) {
-
- m_rfz_t [a] = 6;
- m_rfz_it[a][3] = a-1;
- m_rfz_it[a][4] = b-1;
- m_rfz_it[a][5] = c-1;
- }
+ /// Build maps for radius-azimuthal-Z sorted collections.
+ /// Here, we associate which bins are 'connected' to a given phi-Z bin
+ /// for the seeding
+
+ for (int phiBin = 0; phiBin <= maxPhiBin; ++phiBin)
+ {
+
+ int phiBelow = phiBin - 1;
+ if (phiBelow < 0)
+ phiBelow = maxPhiBin; ///< loop around at edges of range
+
+ int phiAbove = phiBin + 1;
+ if (phiAbove > maxPhiBin)
+ phiAbove = 0; ///< loop around at edges of range
+
+ /// For each azimuthal region loop through all Z regions
+ for (int z = 0; z < arraySizeZ; ++z)
+ {
+
+ /// we always include the two neighbouring phi bins for the top / bottom
+ /// space point search
+
+ int twoDbinSamePhi = phiBin * arraySizeZ + z;
+ int twoDbinLowerPhi = phiBelow * arraySizeZ + z;
+ int twoDbinHigherPhi = phiAbove * arraySizeZ + z;
+
+ nNeighbourCellsBottom[twoDbinSamePhi] = 3;
+ nNeighbourCellsTop[twoDbinSamePhi] = 3;
+
+ neighbourCellsBottom[twoDbinSamePhi][0] = twoDbinSamePhi;
+ neighbourCellsTop[twoDbinSamePhi][0] = twoDbinSamePhi;
+
+ neighbourCellsBottom[twoDbinSamePhi][1] = twoDbinLowerPhi;
+ neighbourCellsTop[twoDbinSamePhi][1] = twoDbinLowerPhi;
+
+ neighbourCellsBottom[twoDbinSamePhi][2] = twoDbinHigherPhi;
+ neighbourCellsTop[twoDbinSamePhi][2] = twoDbinHigherPhi;
+
+ /** in addition, we usually add at least one neighbouring slice
+ * in Z. This depends on where we are in the detector.
+ * Guide for the following:
+ * z == 5: central z region, |z|<250mm
+ * 0 1 2 3 4 5 6 7 8 9 10 z bin index
+ * --------------------------------------> Z[mm]
+ * Z=-2500 IP,Z=0 Z=+2500
+ **/
+ if (z == 5)
+ {
+ nNeighbourCellsTop[twoDbinSamePhi] = 9;
+ // in the central z region, we include the two neighbouring
+ // z slices for the top neighbour search
+
+ neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi + 1;
+ neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi + 1;
+ neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi + 1;
+ neighbourCellsTop[twoDbinSamePhi][6] = twoDbinSamePhi - 1;
+ neighbourCellsTop[twoDbinSamePhi][7] = twoDbinLowerPhi - 1;
+ neighbourCellsTop[twoDbinSamePhi][8] = twoDbinHigherPhi - 1;
}
-
- if (z==3) {
- m_rfz_b[a] = 9;
- m_rfz_ib[a][6] = a+2;
- m_rfz_ib[a][7] = b+2;
- m_rfz_ib[a][8] = c+2;
- } else if (z==7) {
- m_rfz_b[a] = 9;
- m_rfz_ib[a][6] = a-2;
- m_rfz_ib[a][7] = b-2;
- m_rfz_ib[a][8] = c-2;
+ // z > 5: positive z values, |z| > 250mm
+ else if (z > 5)
+ {
+ // for the bottom SP search in positive non-central z, we include the
+ // neighbouring Z region on the left (towards the IP) in the bottom
+ // neighbour search
+ nNeighbourCellsBottom[twoDbinSamePhi] = 6;
+ neighbourCellsBottom[twoDbinSamePhi][3] = twoDbinSamePhi - 1;
+ neighbourCellsBottom[twoDbinSamePhi][4] = twoDbinLowerPhi - 1;
+ neighbourCellsBottom[twoDbinSamePhi][5] = twoDbinHigherPhi - 1;
+
+ if (z < 10)
+ {
+ /** for the top SP search in positive z,
+ * if we are not in the outermost z region,
+ * we include the neighbouring Z region on the right (away from the IP).
+ * In z = 10, the most forward, we do not have such a 'right side' neighbour we can add
+ **/
+ nNeighbourCellsTop[twoDbinSamePhi] = 6;
+ neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi + 1;
+ neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi + 1;
+ neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi + 1;
+ }
+ }
+ // z < 5: negative z values, |z| > 250mm
+ else
+ {
+ /** for the bottom SP in negative non-central z, we include
+ * the neighbouring z region on the right (towards the IP) in the
+ * bottom neighbour search
+ **/
+ nNeighbourCellsBottom[twoDbinSamePhi] = 6;
+ neighbourCellsBottom[twoDbinSamePhi][3] = twoDbinSamePhi + 1;
+ neighbourCellsBottom[twoDbinSamePhi][4] = twoDbinLowerPhi + 1;
+ neighbourCellsBottom[twoDbinSamePhi][5] = twoDbinHigherPhi + 1;
+
+ if (z > 0)
+ {
+ // if there is a z region on the left (away from the IP), we include it in the top
+ // neighbour search
+ nNeighbourCellsTop[twoDbinSamePhi] = 6;
+ neighbourCellsTop[twoDbinSamePhi][3] = twoDbinSamePhi - 1;
+ neighbourCellsTop[twoDbinSamePhi][4] = twoDbinLowerPhi - 1;
+ neighbourCellsTop[twoDbinSamePhi][5] = twoDbinHigherPhi - 1;
+ }
}
- }
- }
-
- // Build maps for radius-azimuthal-Z sorted collections for Z
- //
- for (int f=0; f<=m_fvNmax; ++f) {
- int fb = f-1;
- if (fb<0) fb=m_fvNmax;
- int ft = f+1;
- if (ft>m_fvNmax) ft=0;
-
- // For each azimuthal region loop through central Z regions
- //
- for (int z=0; z<SizeZV; ++z) {
-
- int a = f *SizeZV+z;
- int b = fb*SizeZV+z;
- int c = ft*SizeZV+z;
- m_rfzv_n[a] = 3;
- m_rfzv_i[a][0] = a;
- m_rfzv_i[a][1] = b;
- m_rfzv_i[a][2] = c;
- if (z>1) {
- m_rfzv_n[a] = 6;
- m_rfzv_i[a][3] = a-1;
- m_rfzv_i[a][4] = b-1;
- m_rfzv_i[a][5] = c-1;
- } else if (z<1) {
- m_rfzv_n[a] = 6;
- m_rfzv_i[a][3] = a+1;
- m_rfzv_i[a][4] = b+1;
- m_rfzv_i[a][5] = c+1;
+ /**
+ * z bins 3 / 7: 450mm < |z| < 925mm.:
+ * also include the central z region in the bottom SP search.
+ * likely for PPP seeds with hits in pixel barrel + endcaps
+ **/
+ /// Only used for strip seeds
+ if (isSSS)
+ {
+ if (z == 3)
+ {
+ nNeighbourCellsBottom[twoDbinSamePhi] = 9;
+ neighbourCellsBottom[twoDbinSamePhi][6] = twoDbinSamePhi + 2;
+ neighbourCellsBottom[twoDbinSamePhi][7] = twoDbinLowerPhi + 2;
+ neighbourCellsBottom[twoDbinSamePhi][8] = twoDbinHigherPhi + 2;
+ }
+ else if (z == 7)
+ {
+ nNeighbourCellsBottom[twoDbinSamePhi] = 9;
+ neighbourCellsBottom[twoDbinSamePhi][6] = twoDbinSamePhi - 2;
+ neighbourCellsBottom[twoDbinSamePhi][7] = twoDbinLowerPhi - 2;
+ neighbourCellsBottom[twoDbinSamePhi][8] = twoDbinHigherPhi - 2;
+ }
}
}
}
}
+float InDet::SiSpacePointsSeedMaker_ITK::azimuthalStep(const float pTmin, const float maxd0, const float Rmin, const float Rmax) const
+{
+ /// here we approximate the largest curvature
+ /// that can be expected for the seeds we build
+ /// using R[mm] ~ pT[MeV] / (0.3 * B[T]), with B == 2T
+ float Rm = pTmin / .6;
+
+ /// for LRT, the maximum allowed d0 may be larger than
+ /// the radius at which the innermost hit is possible.
+ /// In that case, our "worst case" trajectory
+ /// generating the largest phi displacement
+ /// will be the one with a hit at Rmin and
+ /// d0 == Rmin.
+ float worstCaseD0 = maxd0;
+ if (maxd0 > Rmin)
+ worstCaseD0 = Rmin;
+
+ float sI = std::abs(std::asin(worstCaseD0 / Rmin) - std::asin(worstCaseD0 / Rmax));
+ float sF = std::abs(std::asin(std::min(1., Rmax / (2. * Rm))) -
+ std::asin(std::min(1., Rmin / (2. * Rm))));
+ return sI + sF;
+}
+
///////////////////////////////////////////////////////////////////
// Initiate beam frame work for seed generator
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::buildBeamFrameWork(EventData& data) const
-{
- SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle { m_beamSpotKey };
+void InDet::SiSpacePointsSeedMaker_ITK::buildBeamFrameWork(EventData &data) const
+{
+ SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle{m_beamSpotKey};
const Amg::Vector3D &cb = beamSpotHandle->beamPos();
double tx = std::tan(beamSpotHandle->beamTilt(0));
double ty = std::tan(beamSpotHandle->beamTilt(1));
- double ph = atan2(ty,tx);
- double th = acos(1./sqrt(1.+tx*tx+ty*ty));
- double sint = sin(th);
- double cost = cos(th);
- double sinp = sin(ph);
- double cosp = cos(ph);
-
- data.xbeam[0] = static_cast<float>(cb.x());
- data.xbeam[1] = static_cast<float>(cost*cosp*cosp+sinp*sinp);
- data.xbeam[2] = static_cast<float>(cost*sinp*cosp-sinp*cosp);
- data.xbeam[3] = -static_cast<float>(sint*cosp );
-
- data.ybeam[0] = static_cast<float>(cb.y());
- data.ybeam[1] = static_cast<float>(cost*cosp*sinp-sinp*cosp);
- data.ybeam[2] = static_cast<float>(cost*sinp*sinp+cosp*cosp);
- data.ybeam[3] = -static_cast<float>(sint*sinp );
-
- data.zbeam[0] = static_cast<float>(cb.z());
- data.zbeam[1] = static_cast<float>(sint*cosp);
- data.zbeam[2] = static_cast<float>(sint*sinp);
- data.zbeam[3] = static_cast<float>(cost);
+ double phi = std::atan2(ty, tx);
+ double theta = std::acos(1. / std::sqrt(1. + tx * tx + ty * ty));
+ double sinTheta = std::sin(theta);
+ double cosTheta = std::cos(theta);
+ double sinPhi = std::sin(phi);
+ double cosPhi = std::cos(phi);
+
+ data.xbeam[0] = static_cast<float>(cb.x());
+ data.xbeam[1] = static_cast<float>(cosTheta * cosPhi * cosPhi + sinPhi * sinPhi);
+ data.xbeam[2] = static_cast<float>(cosTheta * sinPhi * cosPhi - sinPhi * cosPhi);
+ data.xbeam[3] = -static_cast<float>(sinTheta * cosPhi);
+
+ data.ybeam[0] = static_cast<float>(cb.y());
+ data.ybeam[1] = static_cast<float>(cosTheta * cosPhi * sinPhi - sinPhi * cosPhi);
+ data.ybeam[2] = static_cast<float>(cosTheta * sinPhi * sinPhi + cosPhi * cosPhi);
+ data.ybeam[3] = -static_cast<float>(sinTheta * sinPhi);
+
+ data.zbeam[0] = static_cast<float>(cb.z());
+ data.zbeam[1] = static_cast<float>(sinTheta * cosPhi);
+ data.zbeam[2] = static_cast<float>(sinTheta * sinPhi);
+ data.zbeam[3] = static_cast<float>(cosTheta);
}
///////////////////////////////////////////////////////////////////
// Initiate beam frame work for seed generator
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::convertToBeamFrameWork
-(EventData& data, const Trk::SpacePoint*const& sp,float* r) const
+void InDet::SiSpacePointsSeedMaker_ITK::convertToBeamFrameWork(EventData &data, const Trk::SpacePoint *const &sp, float *r) const
{
- r[0] = static_cast<float>(sp->globalPosition().x())-data.xbeam[0];
- r[1] = static_cast<float>(sp->globalPosition().y())-data.ybeam[0];
- r[2] = static_cast<float>(sp->globalPosition().z())-data.zbeam[0];
+ r[0] = static_cast<float>(sp->globalPosition().x()) - data.xbeam[0];
+ r[1] = static_cast<float>(sp->globalPosition().y()) - data.ybeam[0];
+ r[2] = static_cast<float>(sp->globalPosition().z()) - data.zbeam[0];
- if (!sp->clusterList().second) return;
+ //not in 21.9 ITk New implementation
+ /*
+ if (!sp->clusterList().second) return;
- // Only for SCT space points
- //
- const InDet::SiCluster* c0 = static_cast<const InDet::SiCluster*>(sp->clusterList().first );
- const InDet::SiCluster* c1 = static_cast<const InDet::SiCluster*>(sp->clusterList().second);
-
- Amg::Vector2D lc0 = c0->localPosition();
- Amg::Vector2D lc1 = c1->localPosition();
-
- std::pair<Amg::Vector3D, Amg::Vector3D > e0 =
- (c0->detectorElement()->endsOfStrip(InDetDD::SiLocalPosition(lc0.y(),lc0.x(),0.)));
- std::pair<Amg::Vector3D, Amg::Vector3D > e1 =
- (c1->detectorElement()->endsOfStrip(InDetDD::SiLocalPosition(lc1.y(),lc1.x(),0.)));
+ // Only for SCT space points
+ //
+ const InDet::SiCluster* c0 = static_cast<const InDet::SiCluster*>(sp->clusterList().first );
+ const InDet::SiCluster* c1 = static_cast<const InDet::SiCluster*>(sp->clusterList().second);
+
+ Amg::Vector2D lc0 = c0->localPosition();
+ Amg::Vector2D lc1 = c1->localPosition();
+
+ std::pair<Amg::Vector3D, Amg::Vector3D > e0 =
+ (c0->detectorElement()->endsOfStrip(InDetDD::SiLocalPosition(lc0.y(),lc0.x(),0.)));
+ std::pair<Amg::Vector3D, Amg::Vector3D > e1 =
+ (c1->detectorElement()->endsOfStrip(InDetDD::SiLocalPosition(lc1.y(),lc1.x(),0.)));
+
+ Amg::Vector3D b0 (e0.second-e0.first);
+ Amg::Vector3D b1 (e1.second-e1.first);
+ Amg::Vector3D d02(e0.first -e1.first);
+
+ // b0
+ r[ 3] = static_cast<float>(b0[0]);
+ r[ 4] = static_cast<float>(b0[1]);
+ r[ 5] = static_cast<float>(b0[2]);
+
+ // b1
+ r[ 6] = static_cast<float>(b1[0]);
+ r[ 7] = static_cast<float>(b1[1]);
+ r[ 8] = static_cast<float>(b1[2]);
+
+ // r0-r2
+ r[ 9] = static_cast<float>(d02[0]);
+ r[10] = static_cast<float>(d02[1]);
+ r[11] = static_cast<float>(d02[2]);
+
+ // r0
+ r[12] = static_cast<float>(e0.first[0])-data.xbeam[0];
+ r[13] = static_cast<float>(e0.first[1])-data.ybeam[0];
+ r[14] = static_cast<float>(e0.first[2])-data.zbeam[0];
+ */
+}
- Amg::Vector3D b0 (e0.second-e0.first);
- Amg::Vector3D b1 (e1.second-e1.first);
- Amg::Vector3D d02(e0.first -e1.first);
+///////////////////////////////////////////////////////////////////
+// Initiate space points seed maker
+///////////////////////////////////////////////////////////////////
- // b0
- r[ 3] = static_cast<float>(b0[0]);
- r[ 4] = static_cast<float>(b0[1]);
- r[ 5] = static_cast<float>(b0[2]);
-
- // b1
- r[ 6] = static_cast<float>(b1[0]);
- r[ 7] = static_cast<float>(b1[1]);
- r[ 8] = static_cast<float>(b1[2]);
+void InDet::SiSpacePointsSeedMaker_ITK::fillLists(EventData &data) const
+{
+ constexpr float twoPi = 2. * M_PI;
+
+ int firstRadialBin = 0;
+ int lastRadialBin = 0;
+ bool endcap = false;
+
+ /**
+ * The following is done separately for each iteration.
+ * We sort the hits in our radially sorted lists into the
+ * z-phi binning, keeping only those we want to consider to reduce
+ * combinatorics
+ *
+ * Note that the use of r_first to start the loop is what ensures that
+ * the first iteration is a pure SSS pass.
+ * In newEvent, r_first is set to the bin after the last
+ * radial bin containing Pixel space points.
+ * For the second iteration, we reset it to zero and thus start in the pixels.
+ **/
+
+ const std::map<float, int> ztoBin{
+ {-2500., 0},
+ {-1400., 1},
+ {-925., 2},
+ {-450., 3},
+ {-250, 4},
+ {250, 5},
+ {450, 6},
+ {925, 7},
+ {1400, 8},
+ {2500, 9},
+ {100000, 10}, ///< if we encounter Si hits at z > +100m, we are probably not in ATLAS anymore...
+ };
+
+ int nPhiBins = data.iteration ? m_maxPhiBinPPP : m_maxPhiBinSSS;
+ float inverseBinSizePhi = data.iteration ? m_inverseBinSizePhiPPP : m_inverseBinSizePhiSSS;
+
+ for (int radialBin = data.r_first; radialBin < m_nBinsR; ++radialBin)
+ {
+
+ /// skip empty radial bins
+ if (!data.r_map[radialBin])
+ continue;
+
+ // Stop when we reach strip SP in PPP iteration #1
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator SP_first = data.r_Sorted_ITK[radialBin].begin();
+ if (data.iteration && (*SP_first)->spacepoint->clusterList().second)
+ break;
+
+ /// remember the first non-empty bin we encounter
+ if (firstRadialBin == 0)
+ firstRadialBin = radialBin;
+ lastRadialBin = radialBin;
+
+ // loop over the space points in the r-bin and sort them into the 2d phi-z binning
+ for (InDet::SiSpacePointForSeedITK *SP : data.r_Sorted_ITK[radialBin])
+ {
+
+ /// Azimuthal angle sort
+ /// find the bin by dividing phi in 0...2pi by the bin size
+ float Phi = SP->phi();
+ if (Phi < 0.)
+ Phi += twoPi; // phi is defined in [0..2pi] for the binning
+ int phiBin = static_cast<int>(Phi * inverseBinSizePhi);
+ /// handle overflows
+ if (phiBin < 0)
+ {
+ phiBin = nPhiBins;
+ }
+ else if (phiBin > nPhiBins)
+ {
+ phiBin = 0;
+ }
- // r0-r2
- r[ 9] = static_cast<float>(d02[0]);
- r[10] = static_cast<float>(d02[1]);
- r[11] = static_cast<float>(d02[2]);
+ float Z = SP->z();
+ endcap = (std::abs(Z) > 1490);
+ /** z-coordinate sort.
+ * Here, we have a variable bin size.
+ * Use a map to replace 5 levels of nested ternaries
+ * for a somewhat more readable notation while retaining
+ * a logN lookup speed
+ **/
+ int zBin{0};
+ auto bound = ztoBin.lower_bound(Z);
+ /// some protection in case things go REALLY wrong
+ if (bound == ztoBin.end())
+ {
+ --bound; ///< z beyond the max: return last bin
+ }
+ zBin = bound->second;
- // r0
- r[12] = static_cast<float>(e0.first[0])-data.xbeam[0];
- r[13] = static_cast<float>(e0.first[1])-data.ybeam[0];
- r[14] = static_cast<float>(e0.first[2])-data.zbeam[0];
+ /// 2D bin index - computed from the 1D using standard 2D array bin arithmetics
+ int twoDbin = phiBin * arraySizeZ + zBin;
+ /// increment total counter of space points.
+ /// This is not reset between iterations.
+ ++data.nsaz;
+ // push our space point into the 2D binned array
+ data.rfz_Sorted_ITK[twoDbin].push_back(SP);
+ /// the conditional seems to always be true. The rfz_index vector stores
+ /// the 2D bin for each SP in the radius-sorted map. This way,
+ /// we obtain effectively a *3D binning* in r(via the r-sorted vector), phi and z (via the 2D index)
+ if (!data.rfz_map[twoDbin]++)
+ data.rfz_index[data.nrfz++] = twoDbin;
+ }
+ }
+
+ data.state = 0;
+ if(data.iteration){ // PPP
+ data.RTmin = m_binSizeR*firstRadialBin+10. ;
+ data.RTmax = m_binSizeR*lastRadialBin-10.;
+ }else{ //SSS
+ if ( endcap and m_isLRT) {
+ data.RTmin = m_binSizeR*firstRadialBin+10. ;
+ data.RTmax = m_binSizeR*lastRadialBin-10.;
+ }
+ else{
+ data.RTmin = m_binSizeR*firstRadialBin+30. ;
+ data.RTmax = m_binSizeR*lastRadialBin-150.;
+ }
+ }
+
}
-
+
///////////////////////////////////////////////////////////////////
-// Initiate space points seed maker
+/// Initiate space points seed maker for fast tracking
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::fillLists(EventData& data) const
+void InDet::SiSpacePointsSeedMaker_ITK::fillListsFast(const EventContext &ctx, EventData &data) const
{
- constexpr float pi2 = 2.*M_PI;
- std::list<InDet::SiSpacePointForSeedITK*>::iterator r, re;
+ constexpr float twoPi = 2. * M_PI;
+
+ /// Erase any existing entries in the data object
+ erase(data);
+
+ SG::ReadHandle<Trk::PRDtoTrackMap> prd_to_track_map;
+ const Trk::PRDtoTrackMap *prd_to_track_map_cptr = nullptr;
+ if (!m_prdToTrackMap.key().empty())
+ {
+ prd_to_track_map = SG::ReadHandle<Trk::PRDtoTrackMap>(m_prdToTrackMap, ctx);
+ if (!prd_to_track_map.isValid())
+ {
+ ATH_MSG_ERROR("Failed to read PRD to track association map: " << m_prdToTrackMap.key());
+ }
+ prd_to_track_map_cptr = prd_to_track_map.cptr();
+ }
- int ir0 =0;
-
- for (int i=data.r_first; i!=m_r_size; ++i) {
+ /**
+ * The following is done separately for each iteration.
+ * We sort the hits in our radially sorted lists into the
+ * z-phi binning, keeping only those we want to consider to reduce
+ * combinatorics
+ *
+ * Note that the use of r_first to start the loop is what ensures that
+ * the first iteration is a pure SSS pass.
+ * In newEvent, r_first is set to the bin after the last
+ * radial bin containing Pixel space points.
+ * For the second iteration, we reset it to zero and thus start in the pixels.
+ **/
+
+ const std::map<float, int> ztoBin{
+ {-2500., 0},
+ {-1400., 1},
+ {-925., 2},
+ {-450., 3},
+ {-250, 4},
+ {250, 5},
+ {450, 6},
+ {925, 7},
+ {1400, 8},
+ {2500, 9},
+ {100000, 10}, ///< if we encounter Si hits at z > +100m, we are probably not in ATLAS anymore...
+ };
+
+ SG::ReadHandle<SpacePointContainer> spacepoints;
+ if (m_sct)
+ {
+ SG::ReadHandle<SpacePointContainer> spacepointsSCT{m_spacepointsSCT, ctx};
+ spacepoints = spacepointsSCT;
+ }
+ else if (m_pixel)
+ {
+ SG::ReadHandle<SpacePointContainer> spacepointsPixel{m_spacepointsPixel, ctx};
+ spacepoints = spacepointsPixel;
+ }
- if (!data.r_map[i]) continue;
- r = data.r_Sorted_ITK[i].begin();
- re = data.r_Sorted_ITK[i].end();
- if (!ir0) ir0 = i;
+ int nPhiBins = m_pixel ? m_maxPhiBinPPP : m_maxPhiBinSSS;
+ float inverseBinSizePhi = m_pixel ? m_inverseBinSizePhiPPP : m_inverseBinSizePhiSSS;
- if (data.iteration && (*r)->spacepoint->clusterList().second) break;
+ if (spacepoints.isValid())
+ {
+ /// loop over the space points
+ for (const SpacePointCollection *spc : *spacepoints)
+ {
- for (; r!=re; ++r) {
-
- // Azimuthal angle sort
- //
- float F = (*r)->phi();
- if (F<0.) F+=pi2;
+ SpacePointCollection::const_iterator spFirst = spc->begin();
+ float r[15];
+ if (m_pixel)
+ pixInform(*spFirst, r);
+ else if (m_sct)
+ sctInform(data, *spFirst, r);
- int f = static_cast<int>(F*m_sF);
- if (f < 0) f = m_fNmax;
- else if (f > m_fNmax) f = 0;
+ for (const Trk::SpacePoint *sp : *spc)
+ {
- int z;
- float Z = (*r)->z();
+ if (prd_to_track_map_cptr && isUsed(sp, *prd_to_track_map_cptr))
+ continue;
- // Azimuthal angle and Z-coordinate sort
- //
- if (Z>0.) {
- Z< 250.?z=5:Z< 450.?z=6:Z< 925.?z=7:Z< 1400.?z=8:Z< 2500.?z=9:z=10;
- } else {
- Z>-250.?z=5:Z>-450.?z=4:Z>-925.?z=3:Z>-1400.?z=2:Z>-2500.?z=1:z= 0;
- }
+ /** create a SiSpacePointForSeedITK from the space point.
+ * This will also add the point to the l_spforseed list and update
+ * the i_spforseed iterator
+ **/
+ InDet::SiSpacePointForSeedITK *sps = newSpacePoint(data, sp, r);
+ if (!sps)
+ continue;
- int n = f*SizeZ+z;
- ++data.nsaz;
- data.rfz_Sorted_ITK[n].push_back(*r);
- if (!data.rfz_map[n]++) data.rfz_index[data.nrfz++] = n;
-
- if (!data.iteration && (*r)->spacepoint->clusterList().second == 0 && z>=3 && z<=7) {
- z<=4 ? z=0 : z>=6 ? z=2 : z=1;
-
- // Azimuthal angle and Z-coordinate sort for fast vertex search
- //
- f = static_cast<int>(F*m_sFv);
- if (f < 0) f += m_fvNmax;
- else if (f> m_fvNmax) f -= m_fvNmax;
-
- n = f*3+z;
- ++data.nsazv;
- data.rfzv_Sorted_ITK[n].push_back(*r);
- if (!data.rfzv_map[n]++) data.rfzv_index[data.nrfzv++] = n;
+ /// Azimuthal angle sort
+ /// find the bin by dividing phi in 0...2pi by the bin size
+ float Phi = sps->phi();
+ if (Phi < 0.)
+ Phi += twoPi; // phi is defined in [0..2pi] for the binning
+ int phiBin = static_cast<int>(Phi * inverseBinSizePhi);
+ /// handle overflows
+ if (phiBin < 0)
+ {
+ phiBin = nPhiBins;
+ }
+ else if (phiBin > nPhiBins)
+ {
+ phiBin = 0;
+ }
+
+ float Z = sps->z();
+ /** z-coordinate sort.
+ * Here, we have a variable bin size.
+ * Use a map to replace 5 levels of nested ternaries
+ * for a somewhat more readable notation while retaining
+ * a logN lookup speed
+ **/
+ int zBin{0};
+ auto bound = ztoBin.lower_bound(Z);
+ /// some protection in case things go REALLY wrong
+ if (bound == ztoBin.end())
+ {
+ --bound; ///< z beyond the max: return last bin
+ }
+ zBin = bound->second;
+
+ /// 2D bin index - computed from the 1D using standard 2D array bin arithmetics
+ int twoDbin = phiBin * arraySizeZ + zBin;
+ /// increment total counter of space points.
+ /// This is not reset between iterations.
+ ++data.nsaz;
+ // push our space point into the 2D binned array
+ data.rfz_Sorted_ITK[twoDbin].push_back(sps);
+ /// the conditional seems to always be true. The rfz_index vector stores
+ /// the 2D bin for each SP in the radius-sorted map. This way,
+ /// we obtain effectively a *3D binning* in r(via the r-sorted vector), phi and z (via the 2D index)
+ if (!data.rfz_map[twoDbin]++)
+ data.rfz_index[data.nrfz++] = twoDbin;
}
}
}
- data.state = 0;
+
+ // Loop through all RZ collections and sort them in radius order
+ //
+ for (int twoDbin(0); twoDbin != arraySizePhiZ; ++twoDbin)
+ {
+ if (data.rfz_Sorted_ITK[twoDbin].size() > 1)
+ {
+ data.rfz_Sorted_ITK[twoDbin].sort(InDet::SiSpacePointsITKComparison_R());
+ }
+ }
+
+ if(m_pixel){
+ data.RTmin = m_rminPPPFast ;
+ //m_RTmax set per zBin in production3Sp
+ }
+ else if(m_sct){
+ data.RTmin = m_rminSSS ;
+ data.RTmax = m_rmaxSSS ;
+ }
}
-
+
+///////////////////////////////////////////////////////////////////
+// Pixels information
+///////////////////////////////////////////////////////////////////
+
+void InDet::SiSpacePointsSeedMaker_ITK::pixInform(const Trk::SpacePoint *sp, float *r) const
+{
+ const InDet::SiCluster *cl = static_cast<const InDet::SiCluster *>(sp->clusterList().first);
+ const InDetDD::SiDetectorElement *de = cl->detectorElement();
+ const Amg::Transform3D &Tp = de->surface().transform();
+ r[3] = float(Tp(0, 2));
+ r[4] = float(Tp(1, 2));
+ r[5] = float(Tp(2, 2));
+}
+
+///////////////////////////////////////////////////////////////////
+// SCT information
+///////////////////////////////////////////////////////////////////
+
+void InDet::SiSpacePointsSeedMaker_ITK::sctInform(EventData &data, const Trk::SpacePoint *sp, float *r) const
+{
+ const InDet::SiCluster *c0 = static_cast<const InDet::SiCluster *>(sp->clusterList().first);
+ const InDet::SiCluster *c1 = static_cast<const InDet::SiCluster *>(sp->clusterList().second);
+ const InDetDD::SiDetectorElement *d0 = c0->detectorElement();
+ const InDetDD::SiDetectorElement *d1 = c1->detectorElement();
+
+ Amg::Vector2D lc0 = c0->localPosition();
+ Amg::Vector2D lc1 = c1->localPosition();
+
+ std::pair<Amg::Vector3D, Amg::Vector3D> e0 =
+ (d0->endsOfStrip(InDetDD::SiLocalPosition(lc0.y(), lc0.x(), 0.)));
+ std::pair<Amg::Vector3D, Amg::Vector3D> e1 =
+ (d1->endsOfStrip(InDetDD::SiLocalPosition(lc1.y(), lc1.x(), 0.)));
+
+ Amg::Vector3D s0(.5 * (e0.first + e0.second));
+ Amg::Vector3D s1(.5 * (e1.first + e1.second));
+
+ Amg::Vector3D b0(.5 * (e0.second - e0.first));
+ Amg::Vector3D b1(.5 * (e1.second - e1.first));
+ Amg::Vector3D d02(s0 - s1);
+
+ // b0
+ r[3] = float(b0[0]);
+ r[4] = float(b0[1]);
+ r[5] = float(b0[2]);
+
+ // b1
+ r[6] = float(b1[0]);
+ r[7] = float(b1[1]);
+ r[8] = float(b1[2]);
+
+ // r0-r2
+ r[9] = float(d02[0]);
+ r[10] = float(d02[1]);
+ r[11] = float(d02[2]);
+
+ // r0
+ r[12] = float(s0[0]) - data.xbeam[0];
+ r[13] = float(s0[1]) - data.ybeam[0];
+ r[14] = float(s0[2]) - data.zbeam[0];
+}
+
///////////////////////////////////////////////////////////////////
// Erase space point information
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::erase(EventData& data) const
+void InDet::SiSpacePointsSeedMaker_ITK::erase(EventData &data) const
{
- for (int i=0; i<data.nrfz; ++i) {
+ for (int i = 0; i < data.nrfz; ++i)
+ {
int n = data.rfz_index[i];
data.rfz_map[n] = 0;
data.rfz_Sorted_ITK[n].clear();
}
-
- for (int i=0; i<data.nrfzv; ++i) {
+
+ for (int i = 0; i < data.nrfzv; ++i)
+ {
int n = data.rfzv_index[i];
data.rfzv_map[n] = 0;
data.rfzv_Sorted_ITK[n].clear();
}
data.state = 0;
- data.nsaz = 0;
+ data.nsaz = 0;
data.nsazv = 0;
- data.nrfz = 0;
+ data.nrfz = 0;
data.nrfzv = 0;
}
@@ -1093,95 +1748,119 @@ void InDet::SiSpacePointsSeedMaker_ITK::erase(EventData& data) const
// 2 space points seeds production
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::production2Sp(EventData& data) const
+void InDet::SiSpacePointsSeedMaker_ITK::production2Sp(EventData &data) const
{
- if (data.nsazv<2) return;
+ if (data.nsazv < 2)
+ return;
- std::list<InDet::SiSpacePointForSeedITK*>::iterator r0,r0e,r,re;
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator r0, r0e, r, re;
int nseed = 0;
// Loop thorugh all azimuthal regions
//
- for (int f=data.fvNmin; f<=m_fvNmax; ++f) {
+ for (int f = data.fvNmin; f <= m_fvNmax; ++f)
+ {
// For each azimuthal region loop through Z regions
//
int z = 0;
- if (!data.endlist) z = data.zMin;
- for (; z<SizeZV; ++z) {
-
- int a = f*SizeZV+z;
- if (!data.rfzv_map[a]) continue;
- r0 = data.rfzv_Sorted_ITK[a].begin();
- r0e = data.rfzv_Sorted_ITK[a].end ();
-
- if (!data.endlist) {
+ if (!data.endlist)
+ z = data.zMin;
+ for (; z < arraySizeZV; ++z)
+ {
+
+ int a = f * arraySizeZV + z;
+ if (!data.rfzv_map[a])
+ continue;
+ r0 = data.rfzv_Sorted_ITK[a].begin();
+ r0e = data.rfzv_Sorted_ITK[a].end();
+
+ if (!data.endlist)
+ {
r0 = data.rMin_ITK;
data.endlist = true;
}
// Loop through trigger space points
//
- for (; r0!=r0e; ++r0) {
- float X = (*r0)->x();
- float Y = (*r0)->y();
- float R = (*r0)->radius();
- if (R<m_r2minv) continue;
- if (R>m_r2maxv) break;
- float Z = (*r0)->z();
- float ax = X/R;
- float ay = Y/R;
-
- // Bottom links production
- //
- int NB = m_rfzv_n[a];
- for (int i=0; i<NB; ++i) {
- int an = m_rfzv_i[a][i];
- if (!data.rfzv_map[an]) continue;
-
- r = data.rfzv_Sorted_ITK[an].begin();
- re = data.rfzv_Sorted_ITK[an].end ();
-
- for (; r!=re; ++r) {
- float Rb =(*r)->radius();
- if (Rb<m_r1minv) continue;
- if (Rb>m_r1maxv) break;
- float dR = R-Rb;
- if (dR<m_drminv) break;
- if (dR>m_drmax) continue;
- float dZ = Z-(*r)->z();
- float Tz = dZ/dR;
- if (Tz<data.dzdrmin || Tz>data.dzdrmax) continue;
- float Zo = Z-R*Tz;
-
- // Comparison with vertices Z coordinates
- //
- if (!isZCompatible(data, Zo, Rb, Tz)) continue;
-
- // Momentum cut
- //
- float dx =(*r)->x()-X;
- float dy =(*r)->y()-Y;
- float x = dx*ax+dy*ay;
- float y =-dx*ay+dy*ax;
- float xy = x*x+y*y; if (xy == 0.) continue;
- float r2 = 1.f/xy;
- float Ut = x*r2;
- float Vt = y*r2;
- float UR = Ut*R+1.f; if (UR == 0.) continue;
- float A = Vt*R/UR;
- float B = Vt-A*Ut;
- if (std::abs(B*data.K) > m_ipt*sqrt(1.f+A*A)) continue;
+ for (; r0 != r0e; ++r0)
+ {
+ float X = (*r0)->x();
+ float Y = (*r0)->y();
+ float R = (*r0)->radius();
+ if (R < m_r2minv)
+ continue;
+ if (R > m_r2maxv)
+ break;
+ float Z = (*r0)->z();
+ float ax = X / R;
+ float ay = Y / R;
+
+ // Bottom links production
+ //
+ int NB = m_rfzv_n[a];
+ for (int i = 0; i < NB; ++i)
+ {
+ int an = m_rfzv_i[a][i];
+ if (!data.rfzv_map[an])
+ continue;
+
+ r = data.rfzv_Sorted_ITK[an].begin();
+ re = data.rfzv_Sorted_ITK[an].end();
+
+ for (; r != re; ++r)
+ {
+ float Rb = (*r)->radius();
+ if (Rb < m_r1minv)
+ continue;
+ if (Rb > m_r1maxv)
+ break;
+ float dR = R - Rb;
+ if (dR < m_drminv)
+ break;
+ if (dR > m_drmax)
+ continue;
+ float dZ = Z - (*r)->z();
+ float Tz = dZ / dR;
+ if (Tz < data.dzdrmin || Tz > data.dzdrmax)
+ continue;
+ float Zo = Z - R * Tz;
+
+ // Comparison with vertices Z coordinates
+ //
+ if (!isZCompatible(data, Zo, Rb, Tz))
+ continue;
+
+ // Momentum cut
+ //
+ float dx = (*r)->x() - X;
+ float dy = (*r)->y() - Y;
+ float x = dx * ax + dy * ay;
+ float y = -dx * ay + dy * ax;
+ float xy = x * x + y * y;
+ if (xy == 0.)
+ continue;
+ float r2 = 1.f / xy;
+ float Ut = x * r2;
+ float Vt = y * r2;
+ float UR = Ut * R + 1.f;
+ if (UR == 0.)
+ continue;
+ float A = Vt * R / UR;
+ float B = Vt - A * Ut;
+ if (std::abs(B * data.K) > m_ipt * sqrt(1.f + A * A))
+ continue;
++nseed;
- newSeed(data, (*r), (*r0), Zo);
- }
- }
- if (nseed < m_maxsize) continue;
- data.endlist=false;
+ newSeed(data, (*r), (*r0), Zo);
+ }
+ }
+ if (nseed < m_maxsize)
+ continue;
+ data.endlist = false;
data.rMin_ITK = (++r0);
- data.fvNmin=f;
- data.zMin=z;
- return;
+ data.fvNmin = f;
+ data.zMin = z;
+ return;
}
}
}
@@ -1189,686 +1868,1321 @@ void InDet::SiSpacePointsSeedMaker_ITK::production2Sp(EventData& data) const
}
///////////////////////////////////////////////////////////////////
-// Production 3 space points seeds
+// Production 3 space points seeds
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::production3Sp(EventData& data) const
-{
- if (data.nsaz<3) return;
- data.seeds_ITK.clear();
+void InDet::SiSpacePointsSeedMaker_ITK::production3Sp(EventData &data) const
+{
+ /// we need at least 3 SP in our phi-z binning, otherwise can't build 3-SP seeds
+ if (data.nsaz < 3)
+ return;
+
+ /**
+ * This method will run a separate seed formation round
+ * for each phi-Z region, taking the central SP from there
+ * and allowing the top/bottom SP to come from either the same
+ * or certain neighbouring bins.
+ *
+ * The search in each region is performed in the overload of this
+ * method with the extended signature below.
+ * Here, we implement the loop over the 2D regions
+ **/
+
+ /**
+ * Order how we walk across z.
+ * 0-4 are negative z, 5 is central z, 6-10 are positive z.
+ * 0 1 2 3 4 5 6 7 8 9 10 z bin index
+ * --------------------------------------> Z[mm]
+ * Z=-2500 IP,Z=0 Z=+2500
+ **/
+ const std::array<int, arraySizeZ> zBinIndex_SSS{5, 6, 4, 7, 3, 8, 2, 9, 1, 10, 0};
+ const std::array<int, arraySizeZ> zBinIndex_PPP_fast{0, 10, 1, 9, 2, 8, 5, 3, 7, 4, 6};
+ const std::array<int, arraySizeZ> zBinIndex_PPP_long{0, 1, 2, 3, 10, 9, 8, 7, 5, 4, 6};
+ const auto zBinIndex_PPP = m_fastTracking ? zBinIndex_PPP_fast : zBinIndex_PPP_long;
+ // Fast tracking runs a single iteration, either pixel or strip
+ // Default tracking runs a 0-th iteration for strip then a 1-st for pixel
+ bool isPixel = (m_fastTracking && m_pixel) || data.iteration == 1;
+ const auto zBinIndex = isPixel ? zBinIndex_PPP : zBinIndex_SSS;
+
+ const float RTmax[11] = {80., 100., 150., 200., 250., 250., 250., 200., 150., 100., 80.};
+
+ /// prepare arrays to store the iterators over the SP containers for all
+ /// neighbouring cells we wish to consider in the seed formation
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> iter_topCands;
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> iter_endTopCands;
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> iter_bottomCands;
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> iter_endBottomCands;
+
+ int nPhiBins;
+ std::array<int, arraySizePhiZ> nNeighbourCellsBottom;
+ std::array<int, arraySizePhiZ> nNeighbourCellsTop;
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> neighbourCellsBottom;
+ std::array<std::array<int, arraySizeNeighbourBins>, arraySizePhiZ> neighbourCellsTop;
+
+ if (isPixel)
+ {
+ nPhiBins = m_maxPhiBinPPP;
+ nNeighbourCellsBottom = m_nNeighbourCellsBottomPPP;
+ nNeighbourCellsTop = m_nNeighbourCellsTopPPP;
+ neighbourCellsBottom = m_neighbourCellsBottomPPP;
+ neighbourCellsTop = m_neighbourCellsTopPPP;
+ }
+ else
+ {
+ nPhiBins = m_maxPhiBinSSS;
+ nNeighbourCellsBottom = m_nNeighbourCellsBottomSSS;
+ nNeighbourCellsTop = m_nNeighbourCellsTopSSS;
+ neighbourCellsBottom = m_neighbourCellsBottomSSS;
+ neighbourCellsTop = m_neighbourCellsTopSSS;
+ }
- const int ZI[SizeZ]= {5,6,7,8,9,10,4,3,2,1,0};
- std::list<InDet::SiSpacePointForSeedITK*>::iterator rt[9],rte[9],rb[9],rbe[9];
+ /// counter for the found
int nseed = 0;
+ /// prevent another pass from being run when we run out of Seeds
+ data.endlist = true;
- // Loop thorugh all azimuthal regions
- //
- for (int f=data.fNmin; f<=m_fNmax; ++f) {
-
- // For each azimuthal region loop through all Z regions
- //
- int z = 0;
- if (!data.endlist) z = data.zMin;
-
- for (; z<SizeZ; ++z) {
-
- int a = f *SizeZ+ZI[z];
- if (!data.rfz_map[a]) continue;
- int NB = 0, NT = 0;
- for (int i=0; i<m_rfz_b[a]; ++i) {
- int an = m_rfz_ib[a][i];
- if (!data.rfz_map[an]) continue;
- rb [NB] = data.rfz_Sorted_ITK[an].begin();
- rbe[NB++] = data.rfz_Sorted_ITK[an].end();
- }
- for (int i=0; i<m_rfz_t[a]; ++i) {
- int an = m_rfz_it[a][i];
- if (!data.rfz_map[an]) continue;
- rt [NT] = data.rfz_Sorted_ITK[an].begin();
- rte[NT++] = data.rfz_Sorted_ITK[an].end();
- }
-
- if (data.iteration == 0 && data.iteration0 ==0) production3SpSSS(data, rb, rbe, rt, rte, NB, NT, nseed);
- else production3SpPPP(data, rb, rbe, rt, rte, NB, NT, nseed);
-
- if (!data.endlist) {
- data.fNmin = f;
- data.zMin = z;
- return;
+ /// Loop through all azimuthal regions
+ for (int phiBin = data.fNmin; phiBin <= nPhiBins; ++phiBin)
+ {
+
+ /// For each azimuthal region loop through all Z regions
+ int z = (m_fastTracking && m_pixel) ? 2 : 0;
+ /// If we had to abort a previous run, continue where we left off
+ if (!data.endlist)
+ z = data.zMin;
+
+ /// note that this loop follows the order within 'zBinIndex',
+ /// not the ascending order of z regions. We start in the centre,
+ /// not at -2500 mm, and then move outward.
+ for (; z < arraySizeZ; ++z)
+ {
+
+ if (m_fastTracking && m_pixel)
+ {
+ data.RTmax = RTmax[ zBinIndex[z] ];
}
- }
- }
- data.endlist = true;
-}
-///////////////////////////////////////////////////////////////////
-// Production 3 pixel space points seeds for full scan
-///////////////////////////////////////////////////////////////////
+ int phiZbin = phiBin * arraySizeZ + zBinIndex[z];
-void InDet::SiSpacePointsSeedMaker_ITK::production3SpPPP
-(EventData& data,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rb ,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rbe,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rt ,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rte,
- int NB, int NT, int& nseed) const
-{
- std::list<InDet::SiSpacePointForSeedITK*>::iterator r0=rb[0], r;
- if (!data.endlist) {
- r0 = data.rMin_ITK;
- data.endlist = true;
+ /// can skip the rest if this particular 2D bin is empty
+ if (!data.rfz_map[phiZbin])
+ continue;
+
+ /// count how many non-emtpy cells should be searched for the
+ /// top and bottom neighbour
+ int numberBottomCells = 0;
+ int numberTopCells = 0;
+
+ /// walk through the cells in phi-z we wish to consider for the bottom SP search.
+ /// Typically, this will be 3 adjacent phi bins (including the one of the central SP)
+ /// and possibly neighbours in z on side towards the IP or on both sides,
+ /// depdending on the z region we are in
+ for (int neighbourCellNumber = 0; neighbourCellNumber < nNeighbourCellsBottom[phiZbin]; ++neighbourCellNumber)
+ {
+
+ int theNeighbourCell = neighbourCellsBottom[phiZbin][neighbourCellNumber];
+ /// only do something if this cell is populated
+ if (!data.rfz_map[theNeighbourCell])
+ continue;
+ /// plug the begin and end iterators to the SP in the cell into our array
+ iter_bottomCands[numberBottomCells] = data.rfz_Sorted_ITK[theNeighbourCell].begin();
+ iter_endBottomCands[numberBottomCells++] = data.rfz_Sorted_ITK[theNeighbourCell].end();
+ }
+
+ /// walk through the cells in phi-z we wish to consider for the top SP search.
+ /// Typically, this will be 3 adjacent phi bins (including the one of the central SP)
+ /// and possibly neighbours in z on the side opposed to the IP or on both sides,
+ /// depdending on the z region we are in
+ for (int neighbourCellNumber = 0; neighbourCellNumber < nNeighbourCellsTop[phiZbin]; ++neighbourCellNumber)
+ {
+
+ int theNeighbourCell = neighbourCellsTop[phiZbin][neighbourCellNumber];
+ /// only do something if this cell is populated
+ if (!data.rfz_map[theNeighbourCell])
+ continue;
+ /// plug the begin and end iterators to the SP in the cell into our array
+ iter_topCands[numberTopCells] = data.rfz_Sorted_ITK[theNeighbourCell].begin();
+ iter_endTopCands[numberTopCells++] = data.rfz_Sorted_ITK[theNeighbourCell].end();
+ }
+
+ /// now run the seed search for the current phi-z bin.
+ if (!data.trigger)
+ {
+ if (isPixel)
+ production3SpPPP(data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed);
+ else
+ production3SpSSS(data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed);
+ }
+ else
+ production3SpTrigger(data, iter_bottomCands, iter_endBottomCands, iter_topCands, iter_endTopCands, numberBottomCells, numberTopCells, nseed);
+ }
+
+ /** If we exceed the seed capacity, we stop here.
+ * Save where we were in z and phi, and set endlist to false.
+ * This will trigger another run of production3Sp when
+ * The client calls next() after processing all vertices seen
+ * so far (freeing up capacity).
+ **/
+ if (nseed >= m_maxsize)
+ {
+ data.endlist = false;
+ data.fNmin = phiBin + 1;
+ return;
+ }
}
+ /// Processed all seeds there are without aborting - no re-run needed!
+ data.endlist = true;
+}
- float ipt2K = data.ipt2K;
- float ipt2C = data.ipt2C;
- float COFK = data.COFK;
- float imaxp = m_diver;
- float imaxs = m_divermax;
+///////////////////////////////////////////////////////////////////
+// Production 3 pixel space points seeds for full scan
+///////////////////////////////////////////////////////////////////
+
+void InDet::SiSpacePointsSeedMaker_ITK::production3SpPPP(EventData &data,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_bottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_endBottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_topCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_endTopCands,
+ const int numberBottomCells, const int numberTopCells, int &nseed) const
+{
+ /**
+ * This method implements the seed search for a single phi-Z region of the detector.
+ * The central SP is taken from the region, while the top and bottom SP are allowed
+ * to come from either the same or a range of neighbouring cells.
+ **/
+
+ /// iterator across the candidates for the central space point.
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator iter_centralSP = iter_bottomCands[0];
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator iter_otherSP; ///< will be used for iterating over top/bottom SP
+
+ /**
+ * Next, we work out where we are within the ATLAS geometry.
+ * This will help us identify which space-points we need to
+ * consider as potential central points of a seed.
+ **/
+
+ /// find the first central SP candidate above the minimum radius.
+ for (; iter_centralSP != iter_endBottomCands[0]; ++iter_centralSP)
+ {
+ if((*iter_centralSP)->radius() > data.RTmin) break;
+ }
+ /// for the top candidates in the central phi-Z bin, we do not need to start at a smaller
+ /// radius than the lowest-r valid central SP candidate
+ iter_topCands[0] = iter_centralSP;
+ ++iter_topCands[0];
+
+ /// prepare cut values
+ const float &ipt2K = data.ipt2K;
+ const float &ipt2C = data.ipt2C;
+ const float &COFK = data.COFK;
+ const float &maxd0cut = m_maxdImpact;
+ const float &zmax = data.zmaxU;
+ const float &dzdrmax = data.dzdrmax;
data.CmSp_ITK.clear();
- // Loop through all trigger space points
- //
- for (; r0!=rbe[0]; ++r0) {
+ /// keep track of the SP storace capacity.
+ /// Extend it needed (should rarely be the case)
+ size_t SPcapacity = data.SP_ITK.size();
- data.nOneSeeds = 0;
- data.mapOneSeeds_ITK.clear();
+ /// Loop through all central space point candidates
+ for (; iter_centralSP != iter_endBottomCands[0]; ++iter_centralSP)
+ {
- float R = (*r0)->radius();
+ const float &R = (*iter_centralSP)->radius();
+
+ if(R > data.RTmax) break; ///< stop if we have moved outside our radial region of interest.
+
+ /// global coordinates of the central SP
+ const float &X = (*iter_centralSP)->x();
+ const float &Y = (*iter_centralSP)->y();
+ const float &Z = (*iter_centralSP)->z();
+
+ /// for the central SP, we veto locations on the last disk -
+ /// there would be no "outer" hits to complete a seed.
+ double absZ = std::abs(Z);
+ if (!m_fastTracking && absZ > m_zmaxPPP)
+ continue;
+
+ float covr0 = (*iter_centralSP)->covr();
+ float covz0 = (*iter_centralSP)->covz();
+ float ax = X / R;
+ float ay = Y / R;
+ float VR = maxd0cut / (R * R);
+ size_t Ntm = 2;
+ if (R > m_rmaxPPP)
+ Ntm = 1;
+
+ /// initialise a counter for found bottom links
+ /// This also serves as an index in the data.SP vector
+ size_t Nt = 0;
+
+ /// Top links production
+ /// Loop over all the cells where we expect to find such SP
+ for (int cell = 0; cell < numberTopCells; ++cell)
+ {
+
+ if (iter_otherSP == iter_endTopCands[cell])
+ continue;
+
+ /// loop over each SP in each cell
+ for (iter_otherSP = iter_topCands[cell]; iter_otherSP != iter_endTopCands[cell]; ++iter_otherSP)
+ {
+
+ /// evaluate the radial distance,
+ float Rt = (*iter_otherSP)->radius();
+ float dR = Rt - R;
+ /// and continue if we are too close
+ if (dR < m_drminPPP)
+ {
+ iter_topCands[cell] = iter_otherSP;
+ continue;
+ }
- const Trk::Surface* sur0 = (*r0)->sur();
- const Trk::Surface* surn = (*r0)->sun();
- float X = (*r0)->x();
- float Y = (*r0)->y();
- float Z = (*r0)->z();
- int Nb = 0;
+ /// if we are too far, the next ones will be even farther, so abort
+ if (!m_fastTracking && dR > m_drmaxPPP)
+ break;
- // Bottom links production
- //
- for (int i=0; i<NB; ++i) {
- for (r=rb[i]; r!=rbe[i]; ++r) {
- float Rb =(*r)->radius();
- float dR = R-Rb;
+ const float dz = (*iter_otherSP)->z() - Z;
+ const float dZdR = dz / dR;
+ /// Comparison with vertices Z coordinates
+ /// straight line extrapolation to r=0
+ const float z0 = Z - R * dZdR;
+ if (std::abs(z0) > zmax)
+ continue;
+
+ float dx = (*iter_otherSP)->x() - X;
+ float dy = (*iter_otherSP)->y() - Y;
+ float x = dx * ax + dy * ay;
+ float y = dy * ax - dx * ay;
+ float dxy = x * x + y * y;
+ float r2 = 1. / dxy;
+ float u = x * r2;
+ float v = y * r2;
+
+ if (std::abs(R * y) > maxd0cut * x)
+ {
+ float V0;
+ y < 0. ? V0 = VR : V0 = -VR;
+ float A = (v - V0) / (u + 1. / R);
+ float B = V0 + A / R;
+ if ((B * B) > (ipt2K * (1. + A * A)))
+ continue;
+ }
- if (dR > m_drmax) {
- rb[i]=r;
+ const float dr = std::sqrt(r2);
+ const float tz = dz * dr;
+ /// this is effectively a segment-level eta cut - exclude too shallow seed segments
+ if (std::abs(tz) > dzdrmax)
continue;
+
+ /// add SP to the list
+ data.SP_ITK[Nt] = (*iter_otherSP);
+ data.R[Nt] = dr; ///< inverse distance to central SP
+ data.U[Nt] = u; ///< transformed U coordinate
+ data.V[Nt] = v; ///< transformed V coordinate
+ data.Er[Nt] = ((covz0 + (*iter_otherSP)->covz()) + (tz * tz) * (covr0 + (*iter_otherSP)->covr())) * r2; ///<squared Error on 1/tan theta coming from the space-point position errors
+ data.SP_ITK[Nt]->setDR(std::sqrt(dxy + dz * dz));
+ data.SP_ITK[Nt]->setDZDR(dZdR);
+ data.Tn[Nt].Fl = tz;
+ data.Tn[Nt].In = Nt;
+
+ /// if we are exceeding the SP capacity of our data object,
+ /// make it resize its vectors. Will add 50 slots by default,
+ /// so rarely should happen more than once per event.
+ if (++Nt == SPcapacity)
+ {
+ data.resizeSPCont();
+ SPcapacity = data.SP_ITK.size();
}
- if (dR < m_drmin) break;
- if ((*r)->sur()==sur0 || (surn && surn==(*r)->sun())) continue;
-
- float Tz = (Z-(*r)->z())/dR, aTz =std::abs(Tz);
-
- if (aTz < data.dzdrmin || aTz > data.dzdrmax) continue;
-
- // Comparison with vertices Z coordinates
- //
- float Zo = Z-R*Tz;
- if (!isZCompatible(data, Zo, Rb, Tz)) continue;
- data.SP_ITK[Nb] = (*r);
- if (++Nb==m_maxsizeSP) goto breakb;
- }
- }
- breakb:
- if (!Nb || Nb==m_maxsizeSP) continue;
- int Nt = Nb;
-
- // Top links production
- //
- for (int i=0; i<NT; ++i) {
- for (r=rt[i]; r!=rte[i]; ++r) {
- float Rt =(*r)->radius();
- float dR = Rt-R;
-
- if (dR<m_drmin) {
- rt[i]=r;
+ } ///< end of loop over SP within top candidate cell
+ } ///< end of loop over top candidate cells
+
+ if (Nt < Ntm)
+ continue;
+
+ /// now continue with the bottom SP search.
+ /// Make the counter start at Nt, as this serves as a running
+ /// index in the SP list for this seed.
+ size_t Nb = Nt;
+
+ /// Bottom links production
+ /// Loop over all the cells where we expect to find such SP
+ for (int cell = 0; cell < numberBottomCells; ++cell)
+ {
+ /// in each cell, loop over the space points
+ for (iter_otherSP = iter_bottomCands[cell]; iter_otherSP != iter_endBottomCands[cell]; ++iter_otherSP)
+ {
+
+ /// evaluate the radial distance between the central and bottom SP
+ const float &Rb = (*iter_otherSP)->radius();
+ float dR = R - Rb;
+
+ /// if the bottom SP is too far, remember this for future iterations and
+ /// don't bother starting from the beginning again
+ if (dR > m_drmaxPPP)
+ {
+ iter_bottomCands[cell] = iter_otherSP;
continue;
}
- if (dR>m_drmax) break;
+ /// if the points are too close in r, abort (future ones will be even closer).
+ if (dR < m_drminPPP)
+ break;
- if ( (*r)->sur()==sur0 || (surn && surn==(*r)->sun())) continue;
+ const float dz = Z - (*iter_otherSP)->z();
+ const float dZdR = dz / dR;
+ /// Comparison with vertices Z coordinates
+ /// straight line extrapolation to r=0
+ const float z0 = Z - R * dZdR;
+ if (std::abs(z0) > zmax)
+ continue;
- float Tz = ((*r)->z()-Z)/dR, aTz =std::abs(Tz);
+ float dx = (*iter_otherSP)->x() - X;
+ float dy = (*iter_otherSP)->y() - Y;
+ float x = dx * ax + dy * ay;
+ float y = dy * ax - dx * ay;
+ float dxy = x * x + y * y;
+ float r2 = 1. / dxy;
+ float u = x * r2;
+ float v = y * r2;
+
+ if (std::abs(R * y) > -maxd0cut * x)
+ {
+ float V0;
+ y > 0. ? V0 = VR : V0 = -VR;
+ float A = (v - V0) / (u + 1. / R);
+ float B = V0 + A / R;
+ if ((B * B) > (ipt2K * (1. + A * A)))
+ continue;
+ }
- if (aTz < data.dzdrmin || aTz > data.dzdrmax) continue;
+ const float dr = std::sqrt(r2);
+ const float tz = dz * dr;
+ /// this is effectively a segment-level eta cut - exclude too shallow seed segments
+ if (std::abs(tz) > dzdrmax)
+ continue;
+ if (m_fastTracking && (*iter_otherSP)->radius() < 50. && std::abs(tz) > 1.5)
+ continue;
- // Comparison with vertices Z coordinates
- //
- float Zo = Z-R*Tz;
- if (!isZCompatible(data, Zo, R, Tz)) continue;
- data.SP_ITK[Nt] = (*r);
- if (++Nt==m_maxsizeSP) goto breakt;
- }
- }
-
- breakt:
- if (!(Nt-Nb)) continue;
- float covr0 = (*r0)->covr ();
- float covz0 = (*r0)->covz ();
- float ax = X/R;
- float ay = Y/R;
-
- for (int i=0; i<Nt; ++i) {
- InDet::SiSpacePointForSeedITK* sp = data.SP_ITK[i];
-
- float dx = sp->x()-X;
- float dy = sp->y()-Y;
- float dz = sp->z()-Z;
- float x = dx*ax+dy*ay;
- float y = dy*ax-dx*ay;
- float r2 = 1.f/(x*x+y*y);
- float dr = std::sqrt(r2);
- float tz = dz*dr;
- if (i < Nb) tz = -tz;
-
- data.Tz[i] = tz;
- data.Zo[i] = Z-R*tz;
- data.R [i] = dr;
- data.U [i] = x*r2;
- data.V [i] = y*r2;
- data.Er[i] = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;
- }
- covr0 *= .5;
- covz0 *= 2.;
-
- // Three space points comparison
- //
- for (int b=0; b<Nb; ++b) {
- float Zob = data.Zo[b];
- float Tzb = data.Tz[b];
- float Rb2r = data.R [b]*covr0;
- float Rb2z = data.R [b]*covz0;
- float Erb = data.Er[b];
- float Vb = data.V [b];
- float Ub = data.U [b];
- float Tzb2 = (1.f+Tzb*Tzb);
- float sTzb2 = std::sqrt(Tzb2);
- float CSA = Tzb2*COFK;
- float ICSA = Tzb2*ipt2C;
- float imax = imaxp;
- if (data.SP_ITK[b]->spacepoint->clusterList().second) imax = imaxs;
-
- for (int t=Nb; t<Nt; ++t) {
- float dT = ((Tzb-data.Tz[t])*(Tzb-data.Tz[t])-data.R[t]*Rb2z-(Erb+data.Er[t]))-(data.R[t]*Rb2r)*((Tzb+data.Tz[t])*(Tzb+data.Tz[t]));
- if (dT > ICSA) continue;
-
- float dU = data.U[t]-Ub;
- if (dU == 0.) continue;
- float A = (data.V[t]-Vb)/dU;
- float S2 = 1.f+A*A;
- float B = Vb-A*Ub;
- float B2 = B*B;
- if (B2 > ipt2K*S2 || dT*S2 > B2*CSA) continue;
-
- float Im = std::abs((A-B*R)*R);
-
- if (Im <= imax) {
- float dr = data.R[b];
- if (data.R[t] < data.R[b]) dr = data.R[t];
- Im+=std::abs((Tzb-data.Tz[t])/(dr*sTzb2));
- data.CmSp_ITK.emplace_back(std::make_pair(B/std::sqrt(S2), data.SP_ITK[t]));
- data.SP_ITK[t]->setParam(Im);
- }
- }
- if (!data.CmSp_ITK.empty()) {
- newOneSeedWithCurvaturesComparison(data, data.SP_ITK[b], (*r0), Zob);
+ /// add SP to the list
+ data.SP_ITK[Nb] = (*iter_otherSP);
+ data.R[Nb] = dr; ///< inverse distance to central SP
+ data.U[Nb] = u; ///< transformed U coordinate
+ data.V[Nb] = v; ///< transformed V coordinate
+ data.Er[Nb] = ((covz0 + (*iter_otherSP)->covz()) + (tz * tz) * (covr0 + (*iter_otherSP)->covr())) * r2; ///<squared Error on 1/tan theta coming from the space-point position errors
+ data.SP_ITK[Nb]->setDR(std::sqrt(dxy + dz * dz));
+ data.SP_ITK[Nb]->setDZDR(dZdR);
+ data.Tn[Nb].Fl = tz;
+ data.Tn[Nb].In = Nb;
+
+ /// if we are exceeding the SP capacity of our data object,
+ /// make it resize its vectors. Will add 50 slots by default,
+ /// so rarely should happen more than once per event.
+ if (++Nb == SPcapacity)
+ {
+ data.resizeSPCont();
+ SPcapacity = data.SP_ITK.size();
+ }
+
+ } ///< end of loop over SP in bottom candidate cell
+ } ///< end of loop over bottom candidate cells
+
+ /// if we found no bottom candidates (remember, Nb starts counting at Nt), abort
+ if (!(Nb - Nt))
+ continue;
+
+ sort(data.Tn,0,Nt);
+ sort(data.Tn,Nt,Nb-Nt);
+
+ data.nOneSeeds = 0;
+ data.mapOneSeeds_ITK.clear();
+
+ /// Three space points comparison
+ /// first, loop over the bottom point candidates
+ size_t it0 = 0;
+ for (size_t ib = Nt; ib < Nb; ++ib)
+ {
+
+ if (it0 == Nt)
+ break;
+
+ /// retrieve the geometrical paranmeters w.r.t the central SP for this candidate
+ float Tzb = data.Tn[ib].Fl; ///< 1/tanTheta estimate from central+bottom SP
+ int b = data.Tn[ib].In; /// bottom seed index after sorting
+
+ float Rb2r = data.R[b] * covr0;
+ float Rb2z = data.R[b] * covz0;
+ float Erb = data.Er[b]; ///< this is the uncertainty in 1/tanTheta on the bottom segment resulting from the position errors in the 2 SP
+ float Vb = data.V[b]; ///< v-coordinate of bottom SP
+ float Ub = data.U[b]; ///< u-coordinate of bottom SP
+ float Tzb2 = (1. + Tzb * Tzb); ///< 1+1/tan²theta - converts transverse to total squared pt
+ float sTzb2 = std::sqrt(Tzb2); ///< sqrt (1+1/tan²theta) - used to convert pt to |p|
+
+ float sigmaSquaredScatteringPtDependent = Tzb2 * COFK; ///< this, when divided by the 2R², yields an approximated multiple scattering term assuming the measured pt.
+ float sigmaSquaredScatteringMinPt = Tzb2 * ipt2C; ///< this is an approximate worst case multiple scattering term assuming the lowest
+ /// pt we allow and the estimated theta angle
+ /// max IP
+ float d0max = maxd0cut;
+
+ size_t Nc = 1;
+ if (data.SP_ITK[b]->radius() > m_rmaxPPP)
+ Nc = 0;
+ if (data.nOneSeeds)
+ ++Nc;
+
+ /// inner loop over the top point candidates
+ for (size_t it = it0; it < Nt; ++it)
+ {
+
+ int t = data.Tn[it].In; // index of top seed after sorting
+ float Tzt = data.Tn[it].Fl;
+
+ /// Apply a cut on the compatibility between the r-z slope of the two seed segments.
+ /// This is done by comparing the squared difference between slopes, and comparing
+ /// to the squared uncertainty in this difference - we keep a seed if the difference
+ /// is compatible within the assumed uncertainties.
+
+ /// average value of 1/tan(theta), approximate the slope at the location of the central space point
+ float meanOneOverTanThetaSquare = Tzb * Tzt; // SSS uses arithmetic average, PPP geometric average
+
+ /// squared error on the difference in tan(theta) due to space point position errors.
+ float sigmaSquaredSpacePointErrors = Erb + data.Er[t] /// pre-computed individual squared errors on 1/tan(theta) for the two segments
+ + 2 * Rb2z * data.R[t] /// mixed term with z-uncertainty on central SP
+ + 2 * Rb2r * data.R[t] * meanOneOverTanThetaSquare; // mixed term with r-uncertainy on central SP
+
+ /// start out by subtracting from the squared difference in 1/tanTheta the space-point-related squared error
+ float remainingSquaredDelta = (Tzb - Tzt) * (Tzb - Tzt) - sigmaSquaredSpacePointErrors;
+
+ /// First, we test using a generous scattering term calculated assuming the minimum pt we expect
+ /// to reconstruct.
+ if (remainingSquaredDelta - sigmaSquaredScatteringMinPt > 0)
+ {
+ if (Tzb - Tzt < 0.)
+ break;
+ it0 = it + 1;
+ continue;
+ }
+
+ /**
+ * The following exploits the transformation u:=x/(x²+y²); v:=y/(x²+y²);
+ * This is applied on the x,y coordinates in the frame described above, where the
+ * origin is put in the central SP and the x axis defined to point directly away from the IP.
+ *
+ * In this transformed u,v frame, what would be our circle in x-y space takes the form of
+ * a linear function V = (-x0/y0) x U + 1/(2y0) =: A x U + B/2.
+ * Here, x0 and y0 describe the center point of the circle in the x-y frame.
+ * As the origin of the x-y frame (the middle space point of our seed) is on the circle,
+ * we have x0²+y0²=R² with circle radius R.
+ *
+ * For our seed, we can experimentally obtain A as the slope of the linear function,
+ * delta V / delta U,
+ * estimated using the delta U and delta V between the top and bottom space point.
+ *
+ * B is then obtained by inserting the obtained A into the
+ * linear equation for the bottom SP, A x U + B/2 = V --> B = 2(V - A x U0
+ *
+ * With x0²+y0²=R², and x0=-A/B and y0=1/B, the radius of the circle is
+ * then obtained as R²=(1+A²)/B².
+ **/
+
+ float dU = data.U[t] - Ub;
+ if (dU == 0.)
+ continue;
+ float A = (data.V[t] - Vb) / dU;
+ float onePlusAsquare = 1. + A * A;
+ float B = Vb - A * Ub;
+ float BSquare = B * B;
+
+ /** With this radius (and pT) estimate, we can apply our pt cut.
+ * Reminder, ipt2K is 1 / (K x 0.9 x pt-cut)², where K translates pt into 2R.
+ * So here we can apply the pt cut directly on the (2R)² estimate without
+ * the extra overhead of conversion / division.
+ * The second check is a refinement of the above Tz compatibility cut,
+ * replacing the sigmaSquaredScatteringMinPt scattering contribution which assumes the lowest pt
+ * by one based on the actual estimated pt.
+ *
+ * The second term in this if-statement applies a second version of the
+ * 1/tan(theta) compatibility, this time using a scattering term scaled by the actual measured
+ * pt. This refines the cut applied above, following the same logic ("delta² - sigma² ?<=0")
+ **/
+ if (BSquare > ipt2K * onePlusAsquare)
+ continue;
+ if (remainingSquaredDelta * onePlusAsquare > BSquare * sigmaSquaredScatteringPtDependent)
+ {
+ if (Tzb - Tzt < 0.)
+ break;
+ it0 = it + 1;
+ continue;
+ }
+
+ /** This is an estimate of the transverse impact parameter.
+ * The reasoning is that, in the x-y frame with the central SP as origin and
+ * the x axis pointing away from the IP, we have for the distance between
+ * the IP and the middle of the circle:
+ * (x0 - r_central)²+y0² = (R + d0)²,
+ * with R being the circle radius and r_central
+ * the radial location of the central SP, placing the IP at IP at (-r_central, 0).
+ *
+ * First simplify using R² =x0²+y0², then apply the approximation d0²/R² ~ 0.
+ *
+ * Finally, consider that locally close to the central SP, the circle is parallel to the x axis,
+ * so A = 0 --> expand (2R)²=(1+A²)/B² around this point to obtain
+ * d0 = r_central x (r_central x B - A).
+ * Note that below, the variable R is the radial coordinate fo the central SP,
+ * corresponding to r_central in the notation above.
+ **/
+ float d0 = std::abs((A - B * R) * R);
+
+ /// apply d0 cut to seed
+ if (d0 <= d0max)
+ {
+ /// evaluate distance the two closest-by SP in this seed candidate
+ float dr = data.R[b];
+ if (data.R[t] < data.R[b])
+ dr = data.R[t];
+ /// obtain a quality score - start from the d0 estimate, and add
+ /// a penalty term corresponding to how far the seed segments
+ /// deviate from a straight line in r-z
+ data.SP_ITK[t]->setScorePenalty(std::abs((Tzb - Tzt) / (dr * sTzb2)));
+ data.SP_ITK[t]->setParam(d0);
+
+ float meanOneOverTanTheta = std::sqrt(meanOneOverTanThetaSquare);
+ if (meanOneOverTanTheta < 1e-8)
+ meanOneOverTanTheta = 1e-8;
+ if (BSquare < 1e-8)
+ BSquare = 1e-8;
+ float theta = std::atan(1. / std::sqrt(meanOneOverTanThetaSquare));
+ data.SP_ITK[t]->setEta(-std::log(std::tan(0.5 * theta)));
+ data.SP_ITK[t]->setPt(std::sqrt(onePlusAsquare / BSquare) / (1000 * data.K));
+
+ /// record one possible seed candidate, sort by the curvature
+ data.CmSp_ITK.emplace_back(std::make_pair(B / std::sqrt(onePlusAsquare), data.SP_ITK[t]));
+ /// store the transverse IP, will later be used as a quality estimator
+ if (data.CmSp_ITK.size() == 500)
+ break;
+ }
+
+ } ///< end loop over top space point candidates
+ /// now apply further cleaning on the seed candidates for this central+bottom pair.
+ if (data.CmSp_ITK.size() > Nc)
+ {
+ newOneSeedWithCurvaturesComparisonPPP(data, data.SP_ITK[b], (*iter_centralSP), Z - R * Tzb);
}
- }
+ data.CmSp_ITK.clear(); /// cleared in newOneSeedWithCurvaturesComparisonPPP but need to also be cleared in case previous conditional statement isn't fulfilled
+ } ///< end loop over bottom space points
+ ///record seeds found in this run
fillSeeds(data);
nseed += data.fillOneSeeds;
- if (nseed>=m_maxsize) {
- data.endlist=false;
- ++r0;
- data.rMin_ITK = r0;
- return;
- }
- }
+
+ } ///< end loop over central SP
}
///////////////////////////////////////////////////////////////////
-// Production 3 SCT space points seeds for full scan
+/// Production 3 space points seeds for full scan
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::production3SpSSS
-(EventData& data,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rb ,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rbe,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rt ,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rte,
- int NB, int NT, int& nseed) const
+void InDet::SiSpacePointsSeedMaker_ITK::production3SpSSS(EventData &data,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_bottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_endBottomCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_topCands,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &iter_endTopCands,
+ const int numberBottomCells, const int numberTopCells, int &nseed) const
{
- std::list<InDet::SiSpacePointForSeedITK*>::iterator r0=rb[0], r;
- if (!data.endlist) {
- r0 = data.rMin_ITK;
- data.endlist = true;
- }
- float ipt2K = data.ipt2K;
- float ipt2C = data.ipt2C;
- float COFK = data.COFK;
- float imaxs = m_divermax;
+ /**
+ * This method implements the seed search for a single phi-Z region of the detector.
+ * The central SP is taken from the region, while the top and bottom SP are allowed
+ * to come from either the same or a range of neighbouring cells.
+ **/
+
+ /// iterator across the candidates for the central space point.
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator iter_centralSP = iter_bottomCands[0];
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator iter_otherSP; ///< will be used for iterating over top/bottom SP
+
+ /**
+ * Next, we work out where we are within the ATLAS geometry.
+ * This will help us identify which space-points we need to
+ * consider as potential central points of a seed.
+ **/
+
+ /// find the first central SP candidate above the minimum radius.
+ for (; iter_centralSP != iter_endBottomCands[0]; ++iter_centralSP)
+ {
+ if((*iter_centralSP)->radius() > data.RTmin) break;
+ }
+ /// for the top candidates in the central phi-Z bin, we do not need to start at a smaller
+ /// radius than the lowest-r valid central SP candidate
+ iter_topCands[0] = iter_centralSP;
+ ++iter_topCands[0];
+
+ /// prepare cut values
+ const float &ipt2K = data.ipt2K;
+ const float &ipt2C = data.ipt2C;
+ const float &COFK = data.COFK;
+ const float &maxd0cut = m_maxdImpactSSS;
+ const float &zmax = data.zmaxU;
data.CmSp_ITK.clear();
- // Loop through all trigger space points
- //
- for (; r0!=rbe[0]; ++r0) {
- data.nOneSeeds = 0;
- data.mapOneSeeds_ITK.clear();
+ /// keep track of the SP storace capacity.
+ /// Extend it needed (should rarely be the case)
+ size_t SPcapacity = data.SP_ITK.size();
- float R = (*r0)->radius();
+ /// Loop through all central space point candidates
+ for (; iter_centralSP != iter_endBottomCands[0]; ++iter_centralSP)
+ {
+
+ const float &R = (*iter_centralSP)->radius();
+
+ if(R > data.RTmax) break; ///< stop if we have moved outside our radial region of interest.
+
+ /// global coordinates of the central SP
+ const float &X = (*iter_centralSP)->x();
+ const float &Y = (*iter_centralSP)->y();
+ const float &Z = (*iter_centralSP)->z();
+
+ /// for the central SP, we veto locations on the last disk -
+ /// there would be no "outer" hits to complete a seed.
+ double absZ = std::abs(Z);
+ /// veto the last strip disk
+ if (absZ > m_zmaxSSS)
+ continue;
+
+ /// initialise a counter for found bottom links
+ /// This also serves as an index in the data.SP vector
+ size_t Nt = 0;
+
+ /// Top links production
+ /// Loop over all the cells where we expect to find such SP
+ for (int cell = 0; cell < numberTopCells; ++cell)
+ {
+
+ if (iter_otherSP == iter_endTopCands[cell])
+ continue;
+
+ for (iter_otherSP = iter_topCands[cell]; iter_otherSP != iter_endTopCands[cell]; ++iter_otherSP)
+ {
+ /// evaluate the radial distance,
+ float Rt = (*iter_otherSP)->radius();
+ float dR = Rt - R;
+ if (dR >= m_drminSSS)
+ break;
+ }
+ iter_topCands[cell] = iter_otherSP;
- const Trk::Surface* sur0 = (*r0)->sur();
- const Trk::Surface* surn = (*r0)->sun();
- float X = (*r0)->x();
- float Y = (*r0)->y();
- float Z = (*r0)->z();
- int Nb = 0;
+ /// loop over each SP in each cell
+ for (iter_otherSP = iter_topCands[cell]; iter_otherSP != iter_endTopCands[cell]; ++iter_otherSP)
+ {
- // Bottom links production
- //
- for (int i=0; i<NB; ++i) {
- for (r=rb[i]; r!=rbe[i]; ++r) {
- float Rb =(*r)->radius();
- float dR = R-Rb;
- if (dR > m_drmax) {
- rb[i]=r;
+ /// evaluate the radial distance,
+ float Rt = (*iter_otherSP)->radius();
+ float dR = Rt - R;
+ /// if we are to far, the next ones will be even farther, so abort
+ if (dR > m_drmaxSSS)
+ break;
+
+ const float dz = (*iter_otherSP)->z() - Z;
+ const float dZdR = dz / dR;
+
+ /// Comparison with vertices Z coordinates
+ /// straight line extrapolation to r=0
+ const float z0 = Z - R * dZdR;
+ if (std::abs(dz) > m_dzmaxSSS || std::abs(z0) > zmax)
continue;
+
+ /// add SP to the list
+ data.SP_ITK[Nt] = (*iter_otherSP);
+ data.SP_ITK[Nt]->setDZDR(dZdR);
+ /// if we are exceeding the SP capacity of our data object,
+ /// make it resize its vectors. Will add 50 slots by default,
+ /// so rarely should happen more than once per event.
+ if (++Nt == SPcapacity)
+ {
+ data.resizeSPCont();
+ SPcapacity = data.SP_ITK.size();
}
- if (dR < m_drmin) break;
- if ((*r)->sur()==sur0 || (surn && surn==(*r)->sun())) continue;
- float Tz = (Z-(*r)->z())/dR;
- float aTz =std::abs(Tz);
- if (aTz < data.dzdrmin || aTz > data.dzdrmax) continue;
-
- // Comparison with vertices Z coordinates
- //
- float Zo = Z-R*Tz;
- if (!isZCompatible(data, Zo, Rb, Tz)) continue;
- data.SP_ITK[Nb] = (*r);
- if (++Nb==m_maxsizeSP) goto breakb;
- }
- }
- breakb:
- if (!Nb || Nb==m_maxsizeSP) continue;
- int Nt = Nb;
-
- // Top links production
- //
- for (int i=0; i<NT; ++i) {
- for (r=rt[i]; r!=rte[i]; ++r) {
- float Rt =(*r)->radius();
- float dR = Rt-R;
-
- if (dR<m_drmin) {
- rt[i]=r;
+ } ///< end of loop over SP within top candidate cell
+ } ///< end of loop over top candidate cells
+
+ /// if we did not find ANY top SP, or if we exceed the storage capacity, we abort this seed candidate.
+ if (!Nt)
+ continue;
+
+ /// now continue with the bottom SP search.
+ /// Make the counter start at Nt, as this serves as a running
+ /// index in the SP list for this seed.
+ size_t Nb = Nt;
+
+ /// Bottom links production
+ /// Loop over all the cells where we expect to find such SP
+ for (int cell = 0; cell < numberBottomCells; ++cell)
+ {
+ /// in each cell, loop over the space points
+ for (iter_otherSP = iter_bottomCands[cell]; iter_otherSP != iter_endBottomCands[cell]; ++iter_otherSP)
+ {
+
+ /// evaluate the radial distance between the central and bottom SP
+ const float &Rb = (*iter_otherSP)->radius();
+ float dR = R - Rb;
+
+ /// if the bottom SP is too far, remember this for future iterations and
+ /// don't bother starting from the beginning again
+ if (dR > m_drmaxSSS)
+ {
+ iter_bottomCands[cell] = iter_otherSP;
continue;
}
- if (dR>m_drmax) break;
-
- if ((*r)->sur()==sur0 || (surn && surn==(*r)->sun())) continue;
- float Tz = ((*r)->z()-Z)/dR;
- float aTz = std::abs(Tz);
- if (aTz < data.dzdrmin || aTz > data.dzdrmax) continue;
-
- // Comparison with vertices Z coordinates
- //
- float Zo = Z-R*Tz;
- if (!isZCompatible(data, Zo, R, Tz)) continue;
- data.SP_ITK[Nt] = (*r);
- if (++Nt==m_maxsizeSP) goto breakt;
- }
- }
-
- breakt:
- if (!(Nt-Nb)) continue;
- float covr0 = (*r0)->covr();
- float covz0 = (*r0)->covz();
- float ax = X/R;
- float ay = Y/R;
-
- for (int i=0; i<Nt; ++i) {
- InDet::SiSpacePointForSeedITK* sp = data.SP_ITK[i];
-
- float dx = sp->x()-X;
- float dy = sp->y()-Y;
- float dz = sp->z()-Z;
- float x = dx*ax+dy*ay;
- float y = dy*ax-dx*ay;
- float r2 = 1.f/(x*x+y*y);
- float dr = std::sqrt(r2);
- float tz = dz*dr;
- if (i < Nb) tz = -tz;
-
- data.X [i] = x;
- data.Y [i] = y;
- data.Tz[i] = tz;
- data.Zo[i] = Z-R*tz;
- data.R [i] = dr;
- data.U [i] = x*r2;
- data.V [i] = y*r2;
- data.Er[i] = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;
+ /// if the points are too close in r, abort (future ones will be even closer).
+ if (dR < m_drminSSS)
+ break;
+
+ const float dz = Z - (*iter_otherSP)->z();
+ const float dZdR = dz / dR;
+
+ /// Comparison with vertices Z coordinates
+ /// straight line extrapolation to r=0
+ const float z0 = Z - R * dZdR;
+ if (std::abs(dz) > m_dzmaxSSS || std::abs(z0) > zmax)
+ continue;
+ /// found a bottom SP candidate, write it into the data object
+ data.SP_ITK[Nb] = (*iter_otherSP);
+ data.SP_ITK[Nb]->setDZDR(dZdR);
+ /// if we are exceeding the SP capacity of our data object,
+ /// make it resize its vectors. Will add 50 slots by default,
+ /// so rarely should happen more than once per event.
+ if (++Nb == SPcapacity)
+ {
+ data.resizeSPCont();
+ SPcapacity = data.SP_ITK.size();
+ }
+ } ///< end of loop over SP in bottom candidate cell
+ } ///< end of loop over bottom candidate cells
+
+ /// if we found no bottom candidates (remember, Nb starts counting at Nt), abort
+ if (!(Nb - Nt))
+ continue;
+
+ /// get covariance on r and z for the central SP
+ float covr0 = (*iter_centralSP)->covr();
+ float covz0 = (*iter_centralSP)->covz();
+
+ /// build a unit direction vector pointing from the IP to the central SP
+ float ax = X / R;
+ float ay = Y / R;
+
+ /// check all SP candidates we found during our loop and
+ /// compute geometrical variables w.r.t the central point.
+ for (size_t i = 0; i < Nb; ++i)
+ {
+
+ InDet::SiSpacePointForSeedITK *sp = data.SP_ITK[i];
+
+ /// transform the space point coordinates into a frame centered around the middle SP,
+ /// where the x axis points away from the detector frame origin
+ float dx = sp->x() - X;
+ float dy = sp->y() - Y;
+ float dz = sp->z() - Z;
+ float x = dx * ax + dy * ay;
+ float y = dy * ax - dx * ay;
+
+ /// inverse square distance of the candidate space point to the central point
+ float r2 = 1. / (x * x + y * y);
+ /// inverse distance of the candidate space point to the central point
+ float dr = std::sqrt(r2);
+ /// estimate slope in z - distance traveled in transverse plane vs z direction.
+ /// rough estimate of 1/tan theta from 2 points
+ float tz = dz * dr;
+
+ /// if we are looking at a bottom SP candidate, flip the sign to account for
+ /// different direction of flight (from bottom to central)
+ if (i >= Nt)
+ tz = -tz;
+
+ /// save this into our data object
+ data.X[i] = x;
+ data.Y[i] = y;
+ data.Tz[i] = tz; ///< 1/ tan theta
+ data.Zo[i] = Z - R * tz; ///< z0 estimate.
+ data.R[i] = dr; ///< inverse distance to central SP
+ data.U[i] = x * r2; ///< transformed U coordinate
+ data.V[i] = y * r2; ///< transformed V coordinate
+ data.Er[i] = ((covz0 + sp->covz()) + (tz * tz) * (covr0 + sp->covr())) * r2; ///<squared Error on 1/tan theta coming from the space-point position errors
}
- covr0 *= .5;
- covz0 *= 2.;
-
- // Three space points comparison
- //
- for (int b=0; b<Nb; ++b) {
- float Zob = data.Zo[b];
- float Tzb = data.Tz[b];
- float Rb2r = data.R [b]*covr0;
- float Rb2z = data.R [b]*covz0;
- float Erb = data.Er[b];
- float Vb = data.V [b];
- float Ub = data.U [b];
- float Tzb2 = (1.f+Tzb*Tzb);
- float sTzb2 = std::sqrt(Tzb2);
- float CSA = Tzb2*COFK;
- float ICSA = Tzb2*ipt2C;
- float imax = imaxs;
-
- float Se = 1.f/sqrt(1.f+Tzb*Tzb);
- float Ce = Se*Tzb;
- float Sx = Se*ax;
- float Sy = Se*ay;
-
- for (int t=Nb; t<Nt; ++t) {
- // Trigger point
- //
- float dU0 = data.U[t]-Ub;
- if (dU0 == 0.) continue;
- float A0 = (data.V[t]-Vb)/dU0;
- float C0 = 1.f/std::sqrt(1.f+A0*A0);
- float S0 = A0*C0;
- float d0[3] = {Sx*C0-Sy*S0, Sx*S0+Sy*C0, Ce};
- float rn[3];
- if (!(*r0)->coordinates(d0,rn)) continue;
-
- // Bottom point
- //
- float B0 = 2.*(Vb-A0*Ub);
- float Cb = (1.-B0*data.Y[b])*C0;
- float Sb = (A0+B0*data.X[b])*C0;
- float db[3] = {Sx*Cb-Sy*Sb,Sx*Sb+Sy*Cb,Ce};
- float rbDup[3]; //a new and different rb
- if (!data.SP_ITK[b]->coordinates(db,rbDup)) continue;
-
- // Top point
- //
- float Ct = (1.f-B0*data.Y[t])*C0;
- float St = (A0+B0*data.X[t])*C0;
- float dt[3] = {Sx*Ct-Sy*St,Sx*St+Sy*Ct,Ce};
- float rtDup[3]; //doesnt hide previous declaration of rt
- if (!data.SP_ITK[t]->coordinates(dt,rtDup)) continue;
-
- float xb = rbDup[0]-rn[0];
- float yb = rbDup[1]-rn[1];
- float xt = rtDup[0]-rn[0];
- float yt = rtDup[1]-rn[1];
-
- float rb2 = 1.f/(xb*xb+yb*yb);
- float rt2 = 1.f/(xt*xt+yt*yt);
-
- float tb = (rn[2]-rbDup[2])*std::sqrt(rb2);
- float tz = (rtDup[2]-rn[2])*std::sqrt(rt2);
-
- float dT = ((tb-tz)*(tb-tz)-data.R[t]*Rb2z-(Erb+data.Er[t]))-(data.R[t]*Rb2r)*((tb+tz)*(tb+tz));
- if ( dT > ICSA) continue;
-
- float Rn = std::sqrt(rn[0]*rn[0]+rn[1]*rn[1]);
- float Ax = rn[0]/Rn;
- float Ay = rn[1]/Rn;
-
- float ub = (xb*Ax+yb*Ay)*rb2;
- float vb = (yb*Ax-xb*Ay)*rb2;
- float ut = (xt*Ax+yt*Ay)*rt2;
- float vt = (yt*Ax-xt*Ay)*rt2;
-
- float dU = ut-ub;
- if (dU == 0.) continue;
- float A = (vt-vb)/dU;
- float S2 = 1.f+A*A;
- float B = vb-A*ub;
- float B2 = B*B;
- if (B2 > ipt2K*S2 || dT*S2 > B2*CSA) continue;
-
- float Im = std::abs((A-B*Rn)*Rn);
-
- if (Im <= imax) {
- float dr;
- data.R[t] < data.R[b] ? dr = data.R[t] : dr = data.R[b];
- Im+=std::abs((Tzb-data.Tz[t])/(dr*sTzb2));
- data.CmSp_ITK.emplace_back(std::make_pair(B/std::sqrt(S2), data.SP_ITK[t]));
- data.SP_ITK[t]->setParam(Im);
- }
-
- }
- if (!data.CmSp_ITK.empty()) {
- newOneSeedWithCurvaturesComparison(data, data.SP_ITK[b], (*r0), Zob);
+
+ data.nOneSeeds = 0;
+ data.mapOneSeeds_ITK.clear();
+
+ /// Three space points comparison
+ /// first, loop over the bottom point candidates
+ for (size_t b = Nt; b < Nb; ++b)
+ {
+
+ /// retrieve the geometrical paranmeters w.r.t the central SP for this candidate
+ float Zob = data.Zo[b]; ///< z0 estimate from central+bottom SP
+ float Tzb = data.Tz[b]; ///< 1/tanTheta estimate from central+bottom SP
+ float Rb2r = data.R[b] * covr0;
+ float Rb2z = data.R[b] * covz0;
+ float Erb = data.Er[b]; ///< this is the uncertainty in 1/tanTheta on the bottom segment resulting from the position errors in the 2 SP
+ float Vb = data.V[b]; ///< v-coordinate of bottom SP
+ float Ub = data.U[b]; ///< u-coordinate of bottom SP
+ float Tzb2 = (1. + Tzb * Tzb); ///< 1+1/tan²theta - converts transverse to total squared pt
+ float sTzb2 = std::sqrt(Tzb2); ///< sqrt (1+1/tan²theta) - used to convert pt to |p|
+ float Se = 1. / std::sqrt(Tzb2);
+ float Ce = Se * Tzb;
+ float Sx = Se * ax;
+ float Sy = Se * ay;
+
+ float sigmaSquaredScatteringPtDependent = Tzb2 * COFK; ///< this, when divided by the 2R², yields an approximated multiple scattering term assuming the measured pt.
+ float sigmaSquaredScatteringMinPt = Tzb2 * ipt2C; ///< this is an approximate worst case multiple scattering term assuming the lowest
+ /// pt we allow and the estimated theta angle
+ /// max IP
+ float d0max = maxd0cut;
+
+ /// inner loop over the top point candidates
+ for (size_t t = 0; t < Nt; ++t)
+ {
+
+ /**
+ * The following exploits the transformation u:=x/(x²+y²); v:=y/(x²+y²);
+ * This is applied on the x,y coordinates in the frame described above, where the
+ * origin is put in the central SP and the x axis defined to point directly away from the IP.
+ *
+ * In this transformed u,v frame, what would be our circle in x-y space takes the form of
+ * a linear function V = (-x0/y0) x U + 1/(2y0) =: A x U + B/2.
+ * Here, x0 and y0 describe the center point of the circle in the x-y frame.
+ * As the origin of the x-y frame (the middle space point of our seed) is on the circle,
+ * we have x0²+y0²=R² with circle radius R.
+ *
+ * For our seed, we can experimentally obtain A as the slope of the linear function,
+ * delta V / delta U,
+ * estimated using the delta U and delta V between the top and bottom space point.
+ *
+ * B is then obtained by inserting the obtained A into the
+ * linear equation for the bottom SP, A x U + B/2 = V --> B = 2(V - A x U0
+ *
+ * With x0²+y0²=R², and x0=-A/B and y0=1/B, the radius of the circle is
+ * then obtained as R²=(1+A²)/B².
+ **/
+
+ // Trigger point
+ //
+ float dU0 = data.U[t] - Ub;
+ if (dU0 == 0.)
+ continue;
+ float A0 = (data.V[t] - Vb) / dU0;
+
+ float Cn = Ce * std::sqrt(1. + A0 * A0);
+
+ float dn[3] = {Sx - Sy * A0, Sx * A0 + Sy, Cn};
+ float rn[3];
+ if (!(*iter_centralSP)->coordinates(dn, rn))
+ continue;
+
+ // Bottom point
+ //
+ float B0 = 2. * (Vb - A0 * Ub);
+ float Cb = 1. - B0 * data.Y[b];
+ float Sb = A0 + B0 * data.X[b];
+ float db[3] = {Sx * Cb - Sy * Sb, Sx * Sb + Sy * Cb, Cn};
+ float rb[3];
+ if (!data.SP_ITK[b]->coordinates(db, rb))
+ continue;
+
+ // Top point
+ //
+ float Ct = 1. - B0 * data.Y[t];
+ float St = A0 + B0 * data.X[t];
+ float dt[3] = {Sx * Ct - Sy * St, Sx * St + Sy * Ct, Cn};
+ float rt[3];
+ if (!data.SP_ITK[t]->coordinates(dt, rt))
+ continue;
+
+ float xb = rb[0] - rn[0];
+ float yb = rb[1] - rn[1];
+ float zb = rb[2] - rn[2];
+ float xt = rt[0] - rn[0];
+ float yt = rt[1] - rn[1];
+ float zt = rt[2] - rn[2];
+
+ float rb2 = 1. / (xb * xb + yb * yb);
+ float rt2 = 1. / (xt * xt + yt * yt);
+
+ float tb = -zb * std::sqrt(rb2);
+ float tz = zt * std::sqrt(rt2);
+
+ /// Apply a cut on the compatibility between the r-z slope of the two seed segments.
+ /// This is done by comparing the squared difference between slopes, and comparing
+ /// to the squared uncertainty in this difference - we keep a seed if the difference
+ /// is compatible within the assumed uncertainties.
+
+ /// average value of 1/tan(theta), approximate the slope at the location of the central space point
+ float meanOneOverTanTheta = (tb + tz) / 2.;
+
+ /// squared error on the difference in tan(theta) due to space point position errors.
+ float sigmaSquaredSpacePointErrors = Erb + data.Er[t] /// pre-computed individual squared errors on 1/tan(theta) for the two segments
+ + 2 * Rb2z * data.R[t] /// mixed term with z-uncertainty on central SP
+ + 2 * Rb2r * data.R[t] * meanOneOverTanTheta * meanOneOverTanTheta; // mixed term with r-uncertainy on central SP
+
+ /// start out by subtracting from the squared difference in 1/tanTheta the space-point-related squared error
+ float remainingSquaredDelta = (tb - tz) * (tb - tz) - sigmaSquaredSpacePointErrors;
+
+ /// First, we test using a generous scattering term calculated assuming the minimum pt we expect
+ /// to reconstruct.
+ if (remainingSquaredDelta - sigmaSquaredScatteringMinPt > 0)
+ continue;
+
+ float Rn = std::sqrt(rn[0] * rn[0] + rn[1] * rn[1]);
+ float Ax = rn[0] / Rn;
+ float Ay = rn[1] / Rn;
+
+ float ub = (xb * Ax + yb * Ay) * rb2;
+ float vb = (yb * Ax - xb * Ay) * rb2;
+ float ut = (xt * Ax + yt * Ay) * rt2;
+ float vt = (yt * Ax - xt * Ay) * rt2;
+
+ float dU = ut - ub;
+ if (dU == 0.)
+ continue;
+ float A = (vt - vb) / dU;
+ float onePlusAsquare = 1. + A * A;
+ float B = vb - A * ub;
+ float BSquare = B * B;
+
+ /** With this radius (and pT) estimate, we can apply our pt cut.
+ * Reminder, ipt2K is 1 / (K x 0.9 x pt-cut)², where K translates pt into 2R.
+ * So here we can apply the pt cut directly on the (2R)² estimate without
+ * the extra overhead of conversion / division.
+ * The second check is a refinement of the above Tz compatibility cut,
+ * replacing the sigmaSquaredScatteringMinPt scattering contribution which assumes the lowest pt
+ * by one based on the actual estimated pt.
+ *
+ * The second term in this if-statement applies a second version of the
+ * 1/tan(theta) compatibility, this time using a scattering term scaled by the actual measured
+ * pt. This refines the cut applied above, following the same logic ("delta² - sigma² ?<=0")
+ **/
+ if (BSquare > ipt2K * onePlusAsquare || remainingSquaredDelta * onePlusAsquare > BSquare * sigmaSquaredScatteringPtDependent)
+ continue;
+
+ /** This is an estimate of the transverse impact parameter.
+ * The reasoning is that, in the x-y frame with the central SP as origin and
+ * the x axis pointing away from the IP, we have for the distance between
+ * the IP and the middle of the circle:
+ * (x0 - r_central)²+y0² = (R + d0)²,
+ * with R being the circle radius and r_central
+ * the radial location of the central SP, placing the IP at IP at (-r_central, 0).
+ *
+ * First simplify using R² =x0²+y0², then apply the approximation d0²/R² ~ 0.
+ *
+ * Finally, consider that locally close to the central SP, the circle is parallel to the x axis,
+ * so A = 0 --> expand (2R)²=(1+A²)/B² around this point to obtain
+ * d0 = r_central x (r_central x B - A).
+ * Note that below, the variable R is the radial coordinate fo the central SP,
+ * corresponding to r_central in the notation above.
+ **/
+ float d0 = std::abs((A - B * Rn) * Rn);
+
+ /// apply d0 cut to seed
+ if (d0 <= d0max)
+ {
+ /// evaluate distance the two closest-by SP in this seed candidate
+ float dr = std::sqrt(1 / rb2);
+ if (data.R[t] < data.R[b])
+ dr = std::sqrt(1 / rt2);
+ /// obtain a quality score - start from the d0 estimate, and add
+ /// a penalty term corresponding to how far the seed segments
+ /// deviate from a straight line in r-z
+ data.SP_ITK[t]->setScorePenalty(std::abs((tb - tz) / (dr * sTzb2)));
+ data.SP_ITK[t]->setParam(d0);
+ float DR = xt * xt + yt * yt + zt * zt; // distance between top and central SP
+ data.SP_ITK[t]->setDR(DR);
+
+ if (meanOneOverTanTheta < 1e-8)
+ meanOneOverTanTheta = 1e-8;
+ if (BSquare < 1e-8)
+ BSquare = 1e-8;
+ float theta = std::atan(1. / meanOneOverTanTheta);
+ data.SP_ITK[t]->setEta(-std::log(std::tan(0.5 * theta)));
+ data.SP_ITK[t]->setPt(std::sqrt(onePlusAsquare / BSquare) / (1000 * data.K));
+
+ /// record one possible seed candidate, sort by the curvature
+ data.CmSp_ITK.emplace_back(std::make_pair(B / std::sqrt(onePlusAsquare), data.SP_ITK[t]));
+ /// store the transverse IP, will later be used as a quality estimator
+ if (data.CmSp_ITK.size() == 500)
+ break;
+ }
+
+ } ///< end loop over top space point candidates
+ /// now apply further cleaning on the seed candidates for this central+bottom pair.
+ if (!data.CmSp_ITK.empty())
+ {
+ newOneSeedWithCurvaturesComparisonSSS(data, data.SP_ITK[b], (*iter_centralSP), Zob);
}
- }
+ } ///< end loop over bottom space points
+ ///record seeds found in this run
fillSeeds(data);
nseed += data.fillOneSeeds;
- if (nseed>=m_maxsize) {
- data.endlist=false;
- ++r0;
- data.rMin_ITK = r0;
- return;
- }
- }
+
+ } ///< end loop over central SP
}
///////////////////////////////////////////////////////////////////
// Production 3 space points seeds in ROI
///////////////////////////////////////////////////////////////////
-
-void InDet::SiSpacePointsSeedMaker_ITK::production3SpTrigger
-(EventData& data,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rb ,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rbe,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rt ,
- std::list<InDet::SiSpacePointForSeedITK*>::iterator* rte,
- int NB, int NT, int& nseed) const
+void InDet::SiSpacePointsSeedMaker_ITK::production3SpTrigger(EventData &data,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &rb,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &rbe,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &rt,
+ std::array<std::list<InDet::SiSpacePointForSeedITK *>::iterator, arraySizeNeighbourBins> &rte,
+ const int numberBottomCells, const int numberTopCells, int &nseed) const
{
- std::list<InDet::SiSpacePointForSeedITK*>::iterator r0=rb[0], r;
- if (!data.endlist) {
+ std::list<InDet::SiSpacePointForSeedITK *>::iterator r0 = rb[0], r;
+ if (!data.endlist)
+ {
r0 = data.rMin_ITK;
data.endlist = true;
}
- constexpr float pi2 = 2.*M_PI;
+ constexpr float pi2 = 2. * M_PI;
float ipt2K = data.ipt2K;
float ipt2C = data.ipt2C;
- float COFK = data.COFK;
- float imaxp = m_diver;
- float imaxs = m_diversss;
+ float COFK = data.COFK;
+ float imaxp = m_maxdImpact;
+ float imaxs = m_maxdImpactSSS;
data.CmSp_ITK.clear();
// Loop through all trigger space points
//
- for (; r0!=rbe[0]; ++r0) {
+ for (; r0 != rbe[0]; ++r0)
+ {
data.nOneSeeds = 0;
data.mapOneSeeds_ITK.clear();
-
+
float R = (*r0)->radius();
- const Trk::Surface* sur0 = (*r0)->sur();
- float X = (*r0)->x();
- float Y = (*r0)->y();
- float Z = (*r0)->z();
- int Nb = 0;
+ const Trk::Surface *sur0 = (*r0)->sur();
+ float X = (*r0)->x();
+ float Y = (*r0)->y();
+ float Z = (*r0)->z();
+ int Nb = 0;
// Bottom links production
//
- for (int i=0; i<NB; ++i) {
- for (r=rb[i]; r!=rbe[i]; ++r) {
- float Rb =(*r)->radius();
-
- float dR = R-Rb;
- if (dR < m_drmin || (data.iteration && (*r)->spacepoint->clusterList().second)) break;
- if (dR > m_drmax || (*r)->sur()==sur0) continue;
-
- // Comparison with bottom and top Z
- //
- float Tz = (Z-(*r)->z())/dR;
- float Zo = Z-R*Tz;
- if (Zo < data.zminB || Zo > data.zmaxB) continue;
- float Zu = Z+(550.f-R)*Tz;
- if (Zu < data.zminU || Zu > data.zmaxU) continue;
- data.SP_ITK[Nb] = (*r);
- if (++Nb==m_maxsizeSP) goto breakb;
+ for (int i = 0; i < numberBottomCells; ++i)
+ {
+ for (r = rb[i]; r != rbe[i]; ++r)
+ {
+ float Rb = (*r)->radius();
+
+ float dR = R - Rb;
+ if (dR < m_drmin || (data.iteration && (*r)->spacepoint->clusterList().second))
+ break;
+ if (dR > m_drmax || (*r)->sur() == sur0)
+ continue;
+
+ // Comparison with bottom and top Z
+ //
+ float Tz = (Z - (*r)->z()) / dR;
+ float Zo = Z - R * Tz;
+ if (Zo < data.zminB || Zo > data.zmaxB)
+ continue;
+ float Zu = Z + (550.f - R) * Tz;
+ if (Zu < data.zminU || Zu > data.zmaxU)
+ continue;
+ data.SP_ITK[Nb] = (*r);
+ if (++Nb == m_maxsizeSP)
+ goto breakb;
}
}
breakb:
- if (!Nb || Nb==m_maxsizeSP) continue;
+ if (!Nb || Nb == m_maxsizeSP)
+ continue;
int Nt = Nb;
-
+
// Top links production
//
- for (int i=0; i<NT; ++i) {
- for (r=rt[i]; r!=rte[i]; ++r) {
- float Rt =(*r)->radius();
- float dR = Rt-R;
-
- if (dR<m_drmin) {
- rt[i]=r;
+ for (int i = 0; i < numberTopCells; ++i)
+ {
+ for (r = rt[i]; r != rte[i]; ++r)
+ {
+ float Rt = (*r)->radius();
+ float dR = Rt - R;
+
+ if (dR < m_drmin)
+ {
+ rt[i] = r;
continue;
}
- if (dR>m_drmax) break;
-
- if ((*r)->sur()==sur0) continue;
-
- // Comparison with bottom and top Z
- //
- float Tz = ((*r)->z()-Z)/dR;
- float Zo = Z-R*Tz;
- if (Zo < data.zminB || Zo > data.zmaxB) continue;
- float Zu = Z+(550.f-R)*Tz;
- if (Zu < data.zminU || Zu > data.zmaxU) continue;
- data.SP_ITK[Nt] = (*r);
- if (++Nt==m_maxsizeSP) goto breakt;
+ if (dR > m_drmax)
+ break;
+
+ if ((*r)->sur() == sur0)
+ continue;
+
+ // Comparison with bottom and top Z
+ //
+ float Tz = ((*r)->z() - Z) / dR;
+ float Zo = Z - R * Tz;
+ if (Zo < data.zminB || Zo > data.zmaxB)
+ continue;
+ float Zu = Z + (550.f - R) * Tz;
+ if (Zu < data.zminU || Zu > data.zmaxU)
+ continue;
+ data.SP_ITK[Nt] = (*r);
+ if (++Nt == m_maxsizeSP)
+ goto breakt;
}
}
-
+
breakt:
- if (!(Nt-Nb)) continue;
- float covr0 = (*r0)->covr ();
- float covz0 = (*r0)->covz ();
+ if (!(Nt - Nb))
+ continue;
+ float covr0 = (*r0)->covr();
+ float covz0 = (*r0)->covz();
- float ax = X/R;
- float ay = Y/R;
-
- for (int i=0; i<Nt; ++i) {
- InDet::SiSpacePointForSeedITK* sp = data.SP_ITK[i];
-
- float dx = sp->x()-X;
- float dy = sp->y()-Y;
- float dz = sp->z()-Z;
- float x = dx*ax+dy*ay;
- float y = dy*ax-dx*ay;
- float r2 = 1.f/(x*x+y*y);
- float dr = std::sqrt(r2);
- float tz = dz*dr;
- if (i < Nb) tz = -tz;
-
- data.X [i] = x;
- data.Y [i] = y;
- data.Tz[i] = tz;
- data.Zo[i] = Z-R*tz;
- data.R [i] = dr;
- data.U [i] = x*r2;
- data.V [i] = y*r2;
- data.Er[i] = ((covz0+sp->covz())+(tz*tz)*(covr0+sp->covr()))*r2;
+ float ax = X / R;
+ float ay = Y / R;
+
+ for (int i = 0; i < Nt; ++i)
+ {
+ InDet::SiSpacePointForSeedITK *sp = data.SP_ITK[i];
+
+ float dx = sp->x() - X;
+ float dy = sp->y() - Y;
+ float dz = sp->z() - Z;
+ float x = dx * ax + dy * ay;
+ float y = dy * ax - dx * ay;
+ float r2 = 1.f / (x * x + y * y);
+ float dr = std::sqrt(r2);
+ float tz = dz * dr;
+ if (i < Nb)
+ tz = -tz;
+
+ data.X[i] = x;
+ data.Y[i] = y;
+ data.Tz[i] = tz;
+ data.Zo[i] = Z - R * tz;
+ data.R[i] = dr;
+ data.U[i] = x * r2;
+ data.V[i] = y * r2;
+ data.Er[i] = ((covz0 + sp->covz()) + (tz * tz) * (covr0 + sp->covr())) * r2;
}
covr0 *= .5;
covz0 *= 2.;
-
+
// Three space points comparison
//
- for (int b=0; b<Nb; ++b) {
- float Zob = data.Zo[b];
- float Tzb = data.Tz[b];
- float Rb2r = data.R [b]*covr0;
- float Rb2z = data.R [b]*covz0;
- float Erb = data.Er[b];
- float Vb = data.V [b];
- float Ub = data.U [b];
- float Tzb2 = (1.f+Tzb*Tzb);
- float CSA = Tzb2*COFK;
- float ICSA = Tzb2*ipt2C;
- float imax = imaxp;
- if (data.SP_ITK[b]->spacepoint->clusterList().second) imax = imaxs;
-
- for (int t=Nb; t!=Nt; ++t) {
- float dT = ((Tzb-data.Tz[t])*(Tzb-data.Tz[t])-data.R[t]*Rb2z-(Erb+data.Er[t]))-(data.R[t]*Rb2r)*((Tzb+data.Tz[t])*(Tzb+data.Tz[t]));
- if ( dT > ICSA) continue;
- float dU = data.U[t]-Ub;
- if (dU == 0.) continue;
- float A = (data.V[t]-Vb)/dU;
- float S2 = 1.f+A*A;
- float B = Vb-A*Ub;
- float B2 = B*B;
- if (B2 > ipt2K*S2 || dT*S2 > B2*CSA) continue;
-
- float Im = std::abs((A-B*R)*R);
- if (Im > imax) continue;
-
- // Azimuthal angle test
- //
- float y = 1.;
- float x = 2.f*B*R-A;
- float df = std::abs(std::atan2(ay*y-ax*x,ax*y+ay*x)-data.ftrig);
- if (df > M_PI) df=pi2-df;
- if (df > data.ftrigW) continue;
- data.CmSp_ITK.emplace_back(std::make_pair(B/std::sqrt(S2), data.SP_ITK[t]));
+ for (int b = 0; b < Nb; ++b)
+ {
+ float Zob = data.Zo[b];
+ float Tzb = data.Tz[b];
+ float Rb2r = data.R[b] * covr0;
+ float Rb2z = data.R[b] * covz0;
+ float Erb = data.Er[b];
+ float Vb = data.V[b];
+ float Ub = data.U[b];
+ float Tzb2 = (1.f + Tzb * Tzb);
+ float CSA = Tzb2 * COFK;
+ float ICSA = Tzb2 * ipt2C;
+ float imax = imaxp;
+ if (data.SP_ITK[b]->spacepoint->clusterList().second)
+ imax = imaxs;
+
+ for (int t = Nb; t != Nt; ++t)
+ {
+ float dT = ((Tzb - data.Tz[t]) * (Tzb - data.Tz[t]) - data.R[t] * Rb2z - (Erb + data.Er[t])) - (data.R[t] * Rb2r) * ((Tzb + data.Tz[t]) * (Tzb + data.Tz[t]));
+ if (dT > ICSA)
+ continue;
+ float dU = data.U[t] - Ub;
+ if (dU == 0.)
+ continue;
+ float A = (data.V[t] - Vb) / dU;
+ float S2 = 1.f + A * A;
+ float B = Vb - A * Ub;
+ float B2 = B * B;
+ if (B2 > ipt2K * S2 || dT * S2 > B2 * CSA)
+ continue;
+
+ float Im = std::abs((A - B * R) * R);
+ if (Im > imax)
+ continue;
+
+ // Azimuthal angle test
+ //
+ float y = 1.;
+ float x = 2.f * B * R - A;
+ float df = std::abs(std::atan2(ay * y - ax * x, ax * y + ay * x) - data.ftrig);
+ if (df > M_PI)
+ df = pi2 - df;
+ if (df > data.ftrigW)
+ continue;
+ data.CmSp_ITK.emplace_back(std::make_pair(B / std::sqrt(S2), data.SP_ITK[t]));
data.SP_ITK[t]->setParam(Im);
}
- if (!data.CmSp_ITK.empty()) {
+ if (!data.CmSp_ITK.empty())
+ {
newOneSeedWithCurvaturesComparison(data, data.SP_ITK[b], (*r0), Zob);
}
}
fillSeeds(data);
nseed += data.fillOneSeeds;
- if (nseed>=m_maxsize) {
- data.endlist=false;
+ if (nseed >= m_maxsize)
+ {
+ data.endlist = false;
++r0;
data.rMin_ITK = r0;
return;
- }
+ }
}
}
///////////////////////////////////////////////////////////////////
-// New 3 space points pro seeds
+// New 3 space points pro seeds
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::newOneSeed
-(EventData& data,
- InDet::SiSpacePointForSeedITK*& p1, InDet::SiSpacePointForSeedITK*& p2,
- InDet::SiSpacePointForSeedITK*& p3, float z, float q) const
+void InDet::SiSpacePointsSeedMaker_ITK::newOneSeed(EventData &data,
+ InDet::SiSpacePointForSeedITK *&p1, InDet::SiSpacePointForSeedITK *&p2,
+ InDet::SiSpacePointForSeedITK *&p3, float z, float seedCandidateQuality) const
{
- if (data.nOneSeeds < m_maxOneSize) {
- data.OneSeeds_ITK[data.nOneSeeds].set(p1,p2,p3,z);
- data.mapOneSeeds_ITK.insert(std::make_pair(q, &(data.OneSeeds_ITK[data.nOneSeeds])));
+ /// get the worst seed so far
+ float worstQualityInMap = std::numeric_limits<float>::min();
+ InDet::SiSpacePointsProSeedITK *worstSeedSoFar = nullptr;
+ if (!data.mapOneSeeds_ITK.empty())
+ {
+ std::multimap<float, InDet::SiSpacePointsProSeedITK *>::reverse_iterator l = data.mapOneSeeds_ITK.rbegin();
+ worstQualityInMap = (*l).first;
+ worstSeedSoFar = (*l).second;
+ }
+ /// There are three cases where we simply add our new seed to the list and push it into the map:
+ /// a) we have not yet reached our max number of seeds
+ if (data.nOneSeeds < data.maxSeedsPerSP
+ /// b) we have reached the max number but always want to keep confirmed seeds
+ /// and the new seed is a confirmed one, with worse quality than the worst one so far
+ || (m_useSeedConfirmation && data.keepAllConfirmedSeeds && worstQualityInMap <= seedCandidateQuality && isConfirmedSeed(p1, p3, seedCandidateQuality) && data.nOneSeeds < data.seedPerSpCapacity)
+ /// c) we have reached the max number but always want to keep confirmed seeds
+ ///and the new seed of higher quality than the worst one so far, with the latter however being confirmed
+ || (m_useSeedConfirmation && data.keepAllConfirmedSeeds && worstQualityInMap > seedCandidateQuality && isConfirmedSeed(worstSeedSoFar->spacepoint0(), worstSeedSoFar->spacepoint2(), worstQualityInMap) && data.nOneSeeds < data.seedPerSpCapacity))
+ {
+ data.OneSeeds_ITK[data.nOneSeeds].set(p1, p2, p3, z);
+ data.mapOneSeeds_ITK.insert(std::make_pair(seedCandidateQuality, &data.OneSeeds_ITK[data.nOneSeeds]));
++data.nOneSeeds;
- } else {
- std::multimap<float,InDet::SiSpacePointsProSeedITK*>::reverse_iterator
- l = data.mapOneSeeds_ITK.rbegin();
+ }
- if ((*l).first <= q) return;
-
- InDet::SiSpacePointsProSeedITK* s = (*l).second;
- s->set(p1,p2,p3,z);
-
- std::multimap<float,InDet::SiSpacePointsProSeedITK*>::iterator
- i = data.mapOneSeeds_ITK.insert(std::make_pair(q,s));
-
- for (++i; i!=data.mapOneSeeds_ITK.end(); ++i) {
- if ((*i).second==s) {
+ /// otherwise, we check if there is a poorer-quality seed that we can kick out
+ else if (worstQualityInMap > seedCandidateQuality)
+ {
+ /// Overwrite the parameters of the worst seed with the new one
+ worstSeedSoFar->set(p1, p2, p3, z);
+ /// re-insert it with its proper quality to make sure it ends up in the right place
+ std::multimap<float, InDet::SiSpacePointsProSeedITK *>::iterator
+ i = data.mapOneSeeds_ITK.insert(std::make_pair(seedCandidateQuality, worstSeedSoFar));
+ /// and remove the entry with the old quality to avoid duplicates
+ for (++i; i != data.mapOneSeeds_ITK.end(); ++i)
+ {
+ if ((*i).second == worstSeedSoFar)
+ {
data.mapOneSeeds_ITK.erase(i);
return;
}
@@ -1880,57 +3194,78 @@ void InDet::SiSpacePointsSeedMaker_ITK::newOneSeed
// New 3 space points pro seeds production
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::newOneSeedWithCurvaturesComparison
-(EventData& data, SiSpacePointForSeedITK*& SPb, SiSpacePointForSeedITK*& SP0, float Zob) const
+void InDet::SiSpacePointsSeedMaker_ITK::newOneSeedWithCurvaturesComparison(EventData &data, SiSpacePointForSeedITK *&SPb, SiSpacePointForSeedITK *&SP0, float Zob) const
{
const float dC = .00003;
- bool pixb = !SPb->spacepoint->clusterList().second;
+ bool pixb = !SPb->spacepoint->clusterList().second;
std::sort(data.CmSp_ITK.begin(), data.CmSp_ITK.end(), comCurvatureITK());
- std::vector<std::pair<float,InDet::SiSpacePointForSeedITK*>>::iterator j,jn,i = data.CmSp_ITK.begin(),ie = data.CmSp_ITK.end();
- jn=i;
-
- for (; i!=ie; ++i) {
- float u = (*i).second->param();
- bool pixt = !(*i).second->spacepoint->clusterList().second;
- if (pixt && std::abs(SPb->z() -(*i).second->z()) > m_dzmaxPPP) continue;
-
- const Trk::Surface* Sui = (*i).second->sur ();
- float Ri = (*i).second->radius();
- float Ci1 =(*i).first-dC;
- float Ci2 =(*i).first+dC;
- float Rmi = 0.;
- float Rma = 0.;
- bool in = false;
-
- if (!pixb) u-=400.;
- else if ( pixt) u-=200.;
-
- for (j=jn; j!=ie; ++j) {
- if ( j == i ) continue;
- if ( (*j).first < Ci1 ) {jn=j; ++jn; continue;}
- if ( (*j).first > Ci2 ) break;
- if ( (*j).second->sur()==Sui) continue;
-
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator j, jn, i = data.CmSp_ITK.begin(), ie = data.CmSp_ITK.end();
+ jn = i;
+
+ for (; i != ie; ++i)
+ {
+ float u = (*i).second->param();
+ bool pixt = !(*i).second->spacepoint->clusterList().second;
+ if (pixt && std::abs(SPb->z() - (*i).second->z()) > m_dzmaxPPP)
+ continue;
+
+ const Trk::Surface *Sui = (*i).second->sur();
+ float Ri = (*i).second->radius();
+ float Ci1 = (*i).first - dC;
+ float Ci2 = (*i).first + dC;
+ float Rmi = 0.;
+ float Rma = 0.;
+ bool in = false;
+
+ if (!pixb)
+ u -= 400.;
+ else if (pixt)
+ u -= 200.;
+
+ for (j = jn; j != ie; ++j)
+ {
+ if (j == i)
+ continue;
+ if ((*j).first < Ci1)
+ {
+ jn = j;
+ ++jn;
+ continue;
+ }
+ if ((*j).first > Ci2)
+ break;
+ if ((*j).second->sur() == Sui)
+ continue;
+
float Rj = (*j).second->radius();
- if (std::abs(Rj-Ri) < m_drmin) continue;
-
- if (in) {
- if (Rj > Rma) Rma = Rj;
- else if (Rj < Rmi) Rmi = Rj;
- else continue;
- if ( (Rma-Rmi) > 20.) {
- u-=200.;
+ if (std::abs(Rj - Ri) < m_drmin)
+ continue;
+
+ if (in)
+ {
+ if (Rj > Rma)
+ Rma = Rj;
+ else if (Rj < Rmi)
+ Rmi = Rj;
+ else
+ continue;
+ if ((Rma - Rmi) > 20.)
+ {
+ u -= 200.;
break;
}
- } else {
- in=true;
- Rma=Rmi=Rj;
- u-=200.;
+ }
+ else
+ {
+ in = true;
+ Rma = Rmi = Rj;
+ u -= 200.;
}
}
- if (u > m_umax) continue;
+ if (u > m_umax)
+ continue;
newOneSeed(data, SPb, SP0, (*i).second, Zob, u);
}
@@ -1941,148 +3276,615 @@ void InDet::SiSpacePointsSeedMaker_ITK::newOneSeedWithCurvaturesComparison
// Fill seeds
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::fillSeeds(EventData& data) const
+void InDet::SiSpacePointsSeedMaker_ITK::fillSeeds(EventData &data) const
{
data.fillOneSeeds = 0;
- std::multimap<float,InDet::SiSpacePointsProSeedITK*>::iterator
- lf = data.mapOneSeeds_ITK.begin(),
- l = data.mapOneSeeds_ITK.begin(),
- le = data.mapOneSeeds_ITK.end ();
-
- if (l==le) return;
+ std::multimap<float, InDet::SiSpacePointsProSeedITK *>::iterator it_seedCandidate = data.mapOneSeeds_ITK.begin();
+ std::multimap<float, InDet::SiSpacePointsProSeedITK *>::iterator it_endSeedCandidates = data.mapOneSeeds_ITK.end();
- SiSpacePointsProSeedITK* s = nullptr;
+ /// no seeds - nothing to do.
+ if (it_seedCandidate == it_endSeedCandidates)
+ return;
- for (; l!=le; ++l) {
- float w = (*l).first;
- s = (*l).second;
- if (l!=lf && s->spacepoint0()->radius() < 43. && w > -200.) continue;
- if (!s->setQuality(w)) continue;
-
- if (data.i_seede_ITK!=data.l_seeds_ITK.end()) {
- s = &(*data.i_seede_ITK++);
- *s = *(*l).second;
- } else {
- data.l_seeds_ITK.emplace_back(SiSpacePointsProSeedITK(*(*l).second));
- s = &(data.l_seeds_ITK.back());
+ SiSpacePointsProSeedITK *theSeed{nullptr};
+
+ /// loop over the seed candidates we have stored in the event data
+ for (; it_seedCandidate != it_endSeedCandidates; ++it_seedCandidate)
+ {
+
+ /// quality score of the seed, lower = better, list is sorted by quality
+ float quality = (*it_seedCandidate).first;
+ theSeed = (*it_seedCandidate).second;
+
+ /// this will set the quality member of all points on the seed to the quality score of this candidate
+ if (!theSeed->setQuality(quality))
+ continue;
+
+ /// if we have space, write the seed directly into an existing slot
+ if (data.i_seede_ITK != data.l_seeds_ITK.end())
+ {
+ theSeed = &(*data.i_seede_ITK++);
+ *theSeed = *(*it_seedCandidate).second;
+ }
+ else
+ {
+ /// otherwise, extend the seed list and update the iterators
+ data.l_seeds_ITK.emplace_back(SiSpacePointsProSeedITK(*(*it_seedCandidate).second));
+ theSeed = &(data.l_seeds_ITK.back());
data.i_seede_ITK = data.l_seeds_ITK.end();
}
-
- if (s->spacepoint0()->spacepoint->clusterList().second) w-=3000.;
- else if (s->spacepoint1()->spacepoint->clusterList().second) w-=2000.;
- else if (s->spacepoint2()->spacepoint->clusterList().second) w-=1000.;
- data.seeds_ITK.insert(std::make_pair(w,s));
++data.fillOneSeeds;
- }
+ } ///< end loop over seed candidates
}
-const InDet::SiSpacePointsSeed* InDet::SiSpacePointsSeedMaker_ITK::next(const EventContext&, EventData& data) const
+const InDet::SiSpacePointsSeed *InDet::SiSpacePointsSeedMaker_ITK::next(const EventContext &, EventData &data) const
{
- if (not data.initialized) initializeEventData(data);
+ /// This only holds if we call next() without manually calling newEvent/find3Sp
+ if (not data.initialized)
+ initializeEventData(data);
- if (data.nspoint==3) {
- do {
- if (data.i_seed_ITK==data.i_seede_ITK) {
+ if (data.nspoint == 3)
+ {
+ do
+ {
+ /// If we are out of seeds, call findNext to see if we can find more.
+ if (data.i_seed_ITK == data.i_seede_ITK)
+ {
+ /// findNext will call production3Sp again IF data.endlist is false,
+ /// which is only the case if the last run of production3Sp did not run to the end
+ /// or if we did not run seed finding before
+ /// For run-3 offline, this will not do anything.
findNext(data);
- if (data.i_seed_ITK==data.i_seede_ITK) return nullptr;
+ /// if no new seeds were found, exit
+ if (data.i_seed_ITK == data.i_seede_ITK)
+ return nullptr;
}
- ++data.i_seed_ITK;
- } while (!(*data.seed_ITK++).second->set3(data.seedOutput));
+ /// iterate until we find a valid seed satisfying certain quality cuts in set3
+ } while (!(*data.i_seed_ITK++).set3(data.seedOutput));
+ /// then return this next seed candidate
return &data.seedOutput;
- } else {
- if (data.i_seed_ITK==data.i_seede_ITK) {
+ }
+ else
+ {
+ /// same as above for 2SP
+ if (data.i_seed_ITK == data.i_seede_ITK)
+ {
findNext(data);
- if (data.i_seed_ITK==data.i_seede_ITK) return nullptr;
- }
+ if (data.i_seed_ITK == data.i_seede_ITK)
+ return nullptr;
+ }
(*data.i_seed_ITK++).set2(data.seedOutput);
return &data.seedOutput;
}
return nullptr;
}
-
-bool InDet::SiSpacePointsSeedMaker_ITK::isZCompatible
-(EventData& data, float& Zv, float& R, float& T) const
+bool InDet::SiSpacePointsSeedMaker_ITK::isZCompatible(EventData &data, float &Zv, float &R, float &T) const
{
- if (Zv < data.zminU || Zv > data.zmaxU) return false;
- if (!data.isvertex) return true;
+ if (Zv < data.zminU || Zv > data.zmaxU)
+ return false;
+ if (!data.isvertex)
+ return true;
float dZmin = std::numeric_limits<float>::max();
- for (const float& v: data.l_vertex) {
- float dZ = std::abs(v-Zv);
- if (dZ >= dZmin) break;
- dZmin=dZ;
+ for (const float &v : data.l_vertex)
+ {
+ float dZ = std::abs(v - Zv);
+ if (dZ >= dZmin)
+ break;
+ dZmin = dZ;
}
- return dZmin < (m_dzver+m_dzdrver*R)*sqrt(1.+T*T);
+ return dZmin < (m_dzver + m_dzdrver * R) * sqrt(1. + T * T);
}
///////////////////////////////////////////////////////////////////
-// New space point for seeds
+// New space point for seeds
///////////////////////////////////////////////////////////////////
-InDet::SiSpacePointForSeedITK* InDet::SiSpacePointsSeedMaker_ITK::newSpacePoint
-(EventData& data, const Trk::SpacePoint*const& sp) const
+InDet::SiSpacePointForSeedITK *InDet::SiSpacePointsSeedMaker_ITK::newSpacePoint(EventData &data, const Trk::SpacePoint *const &sp) const
{
- InDet::SiSpacePointForSeedITK* sps = nullptr;
-
float r[15];
+ return newSpacePoint(data, sp, r, true);
+}
+
+InDet::SiSpacePointForSeedITK *InDet::SiSpacePointsSeedMaker_ITK::newSpacePoint(EventData &data, const Trk::SpacePoint *const &sp, float *r, bool usePixSctInform) const
+{
+
+ InDet::SiSpacePointForSeedITK *sps = nullptr;
+
+ /// r will store the coordinates of the space point relative
+ /// to the beam spot
convertToBeamFrameWork(data, sp, r);
- if (data.checketa) {
- float z = (std::abs(r[2])+m_zmax);
- float x = r[0]*data.dzdrmin;
- float y = r[1]*data.dzdrmin;
- if ((z*z )<(x*x+y*y)) return sps;
+ /// if needed, apply eta criterion
+ if (data.checketa)
+ {
+ float z = (std::abs(r[2]) + m_zmax);
+ float x = r[0] * data.dzdrmin;
+ float y = r[1] * data.dzdrmin;
+ if ((z * z) < (x * x + y * y))
+ return sps;
}
- if (data.i_spforseed_ITK!=data.l_spforseed_ITK.end()) {
+ if (m_fastTracking)
+ {
+ float R2 = r[0] * r[0] + r[1] * r[1];
+ if (std::abs(r[2]) > m_dzMaxFast && R2 < m_R2MaxFast)
+ return 0;
+ if (std::abs(r[2]) - m_zmax > data.dzdrmax * std::sqrt(R2))
+ return 0;
+ }
+
+ if (usePixSctInform)
+ {
+ if (!sp->clusterList().second)
+ pixInform(sp, r);
+ else
+ sctInform(data, sp, r);
+ }
+
+ /// If we have previously populated the list and just reset
+ /// the iterator when re-initialising the data object,
+ /// then we re-use existing entries
+ if (data.i_spforseed_ITK != data.l_spforseed_ITK.end())
+ {
+ /// re-use existing entry at the current location
sps = &(*data.i_spforseed_ITK++);
+ /// and then update the existing entry with the new SP and location.
+ /// Unfortunately, set still relies on C-arrays...
sps->set(sp, r);
- } else {
- data.l_spforseed_ITK.emplace_back(InDet::SiSpacePointForSeedITK(sp, r));
+ }
+ else
+ {
+ /// otherwise, the list needs to grow
+ data.l_spforseed_ITK.emplace_back(InDet::SiSpacePointForSeedITK(sp, &(r[0])));
+ /// set our return pointer
sps = &(data.l_spforseed_ITK.back());
+ /// and make sure to update the iterator
data.i_spforseed_ITK = data.l_spforseed_ITK.end();
}
-
+
return sps;
}
///////////////////////////////////////////////////////////////////
-// New 2 space points seeds
+// New 2 space points seeds
///////////////////////////////////////////////////////////////////
-void InDet::SiSpacePointsSeedMaker_ITK::newSeed
-(EventData& data,
- InDet::SiSpacePointForSeedITK*& p1, InDet::SiSpacePointForSeedITK*& p2, float z) const
+void InDet::SiSpacePointsSeedMaker_ITK::newSeed(EventData &data,
+ InDet::SiSpacePointForSeedITK *&p1, InDet::SiSpacePointForSeedITK *&p2, float z) const
{
- InDet::SiSpacePointForSeedITK* p3 = nullptr;
+ InDet::SiSpacePointForSeedITK *p3 = nullptr;
- if (data.i_seede_ITK!=data.l_seeds_ITK.end()) {
- InDet::SiSpacePointsProSeedITK* s = &(*data.i_seede_ITK++);
+ if (data.i_seede_ITK != data.l_seeds_ITK.end())
+ {
+ InDet::SiSpacePointsProSeedITK *s = &(*data.i_seede_ITK++);
s->set(p1, p2, p3, z);
- } else {
+ }
+ else
+ {
data.l_seeds_ITK.emplace_back(InDet::SiSpacePointsProSeedITK(p1, p2, p3, z));
data.i_seede_ITK = data.l_seeds_ITK.end();
}
}
-
-void InDet::SiSpacePointsSeedMaker_ITK::initializeEventData(EventData& data) const {
+
+void InDet::SiSpacePointsSeedMaker_ITK::initializeEventData(EventData &data) const
+{
+ int seedArrayPerSPSize = (m_maxOneSizePPP > m_maxOneSizeSSS ? m_maxOneSizePPP : m_maxOneSizeSSS);
+ if (m_alwaysKeepConfirmedStripSeeds || m_alwaysKeepConfirmedPixelSeeds)
+ seedArrayPerSPSize = 50;
data.initialize(EventData::ITK,
m_maxsizeSP,
- m_maxOneSize,
- 0, // maxsize not used
- m_r_size,
- 0, // sizeRF not used
- SizeRFZ,
- SizeRFZV,
+ seedArrayPerSPSize,
+ 0, /// maxsize not used
+ m_nBinsR,
+ 0, /// sizeRF not used
+ arraySizePhiZ,
+ arraySizePhiZV,
m_checketa);
}
-void InDet::SiSpacePointsSeedMaker_ITK::writeNtuple(const SiSpacePointsSeed*, const Trk::Track*, int, long) const{
+///////////////////////////////////////////////////////////////////
+// New 3 space points pro seeds production
+///////////////////////////////////////////////////////////////////
+
+void InDet::SiSpacePointsSeedMaker_ITK::newOneSeedWithCurvaturesComparisonSSS(EventData &data,
+ SiSpacePointForSeedITK *&SPb, SiSpacePointForSeedITK *&SP0, float Zob) const
+{
+
+ if (m_useSeedConfirmation)
+ {
+ newOneSeedWithCurvaturesComparisonSeedConfirmation(data, SPb, SP0, Zob);
+ }
+
+ else
+ {
+
+ static const float curvatureInterval = .00003;
+
+ /// sort common SP by curvature
+ if (data.CmSp_ITK.size() > 2)
+ std::sort(data.CmSp_ITK.begin(), data.CmSp_ITK.end(), comCurvatureITK());
+
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_otherSP;
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_commonTopSP = data.CmSp_ITK.begin(), ie = data.CmSp_ITK.end();
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_startInnerLoop = it_commonTopSP;
+
+ float Lt[4];
+
+ /// check all possible common top SP
+ for (; it_commonTopSP != ie; ++it_commonTopSP)
+ {
+
+ SiSpacePointForSeedITK *SPt = (*it_commonTopSP).second;
+ int NT = 1;
+ Lt[0] = SPt->dR();
+ float seedIP = SPt->param();
+
+ /// form a curvature interval cut
+ float minCurvature = (*it_commonTopSP).first - curvatureInterval;
+ float maxCurvature = (*it_commonTopSP).first + curvatureInterval;
+
+ /**
+ * Now we look at the other SP candidates and try to find a confirmation seed,
+ * including the same centre/lower SP and giving a compatible curvature,
+ * but with the top SP in a different layer
+ **/
+
+ for (it_otherSP = it_startInnerLoop; it_otherSP != ie; ++it_otherSP)
+ {
+ /// if we are looking at the same SP, skip it
+ if (it_otherSP == it_commonTopSP)
+ continue;
+ /// if we have a lower curvature than the minimum, skip - and remember to
+ /// not bother with this candidate again later, as the vectors are curvature-sorted
+ if ((*it_otherSP).first < minCurvature)
+ {
+ it_startInnerLoop = it_otherSP;
+ ++it_startInnerLoop;
+ continue;
+ }
+ /// abort once the the curvature gets too large
+ if ((*it_otherSP).first > maxCurvature)
+ break;
+
+ float L = (*it_otherSP).second->dR();
+
+ int k = 0;
+ for (; k != NT; ++k)
+ {
+ if (std::abs(L - Lt[k]) < 20.)
+ break;
+ }
+ if (k == NT)
+ {
+ Lt[NT] = L;
+ if (++NT == 4)
+ break;
+ }
+ }
+
+ // ITK seed quality used so far
+ float Q = seedIP - float(NT) * 100.;
+ if (NT > 2)
+ Q -= 100000.;
+ /// this is a good seed, save it (unless we have too many seeds per SP)
+ newOneSeed(data, SPb, SP0, SPt, Zob, Q);
+ } ///< end of loop over top SP candidates
+ data.CmSp_ITK.clear();
+ }
+}
+
+void InDet::SiSpacePointsSeedMaker_ITK::newOneSeedWithCurvaturesComparisonPPP(EventData &data,
+ SiSpacePointForSeedITK *&SPb, SiSpacePointForSeedITK *&SP0, float Zob) const
+{
+
+ if (m_useSeedConfirmation)
+ {
+ newOneSeedWithCurvaturesComparisonSeedConfirmation(data, SPb, SP0, Zob);
+ }
+
+ else
+ {
+
+ static const float curvatureInterval = .00003;
+
+ /// sort common SP by curvature
+ if (data.CmSp_ITK.size() > 2)
+ std::sort(data.CmSp_ITK.begin(), data.CmSp_ITK.end(), comCurvatureITK());
+
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_otherSP;
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_commonTopSP = data.CmSp_ITK.begin(), ie = data.CmSp_ITK.end();
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_startInnerLoop = it_commonTopSP;
+
+ float Lt[4];
+
+ float Qmin = 1.e20;
+ float Rb = 2. * SPb->radius();
+ int NTc(2);
+ if (Rb > 280.)
+ NTc = 1;
+ SiSpacePointForSeedITK *SPmin = nullptr;
+ bool Qm = Rb < 120. || std::abs(Zob) > 150.;
+
+ /// check all possible common top SP
+ for (; it_commonTopSP != ie; ++it_commonTopSP)
+ {
+
+ SiSpacePointForSeedITK *SPt = (*it_commonTopSP).second;
+ int NT = 1;
+ Lt[0] = SPt->dR();
+ float seedIP = SPt->param();
+
+ /// form a curvature interval cut
+ float minCurvature = (*it_commonTopSP).first - curvatureInterval;
+ float maxCurvature = (*it_commonTopSP).first + curvatureInterval;
+
+ /**
+ * Now we look at the other SP candidates and try to find a confirmation seed,
+ * including the same centre/lower SP and giving a compatible curvature,
+ * but with the top SP in a different layer
+ **/
+
+ for (it_otherSP = it_startInnerLoop; it_otherSP != ie; ++it_otherSP)
+ {
+ /// if we are looking at the same SP, skip it
+ if (it_otherSP == it_commonTopSP)
+ continue;
+ /// if we have a lower curvature than the minimum, skip - and remember to
+ /// not bother with this candidate again later, as the vectors are curvature-sorted
+ if ((*it_otherSP).first < minCurvature)
+ {
+ it_startInnerLoop = it_otherSP;
+ ++it_startInnerLoop;
+ continue;
+ }
+ /// abort once the the curvature gets too large
+ if ((*it_otherSP).first > maxCurvature)
+ break;
+
+ float L = (*it_otherSP).second->dR();
+
+ int k = 0;
+ for (; k != NT; ++k)
+ {
+ if (std::abs(L - Lt[k]) < 20.)
+ break;
+ }
+ if (k == NT)
+ {
+ Lt[NT] = L;
+ if (++NT == 4)
+ break;
+ }
+ }
+
+ int dN = NT - NTc;
+ if (dN < 0 || (data.nOneSeeds && !dN))
+ continue;
+ if (Qm && !dN && seedIP > 1.)
+ continue;
+
+ // ITK seed quality used so far
+ float Q = 100. * seedIP + (std::abs(Zob) - float(NT) * 100.);
+ if (Q > SPb->quality() && Q > SP0->quality() && Q > SPt->quality())
+ continue;
+
+ if (dN)
+ newOneSeed(data, SPb, SP0, SPt, Zob, Q);
+ else if (Q < Qmin)
+ {
+ Qmin = Q;
+ SPmin = SPt;
+ }
+ } ///< end of loop over top SP candidates
+ if (SPmin && !data.nOneSeeds)
+ newOneSeed(data, SPb, SP0, SPmin, Zob, Qmin);
+ data.CmSp_ITK.clear();
+ }
}
-bool InDet::SiSpacePointsSeedMaker_ITK::getWriteNtupleBoolProperty() const{
- return false;
+void InDet::SiSpacePointsSeedMaker_ITK::newOneSeedWithCurvaturesComparisonSeedConfirmation(EventData &data,
+ SiSpacePointForSeedITK *&SPb, SiSpacePointForSeedITK *&SP0, float Zob) const
+{
+ static const float curvatureInterval = .00003;
+ bool bottomSPisPixel = !SPb->spacepoint->clusterList().second;
+ float bottomSPQuality = SPb->quality();
+ float centralSPQuality = SP0->quality();
+
+ /// sort common SP by curvature
+ if (data.CmSp_ITK.size() > 2)
+ std::sort(data.CmSp_ITK.begin(), data.CmSp_ITK.end(), comCurvatureITK());
+
+ float bottomR = SPb->radius();
+ float bottomZ = SPb->z();
+
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_otherSP;
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_commonTopSP = data.CmSp_ITK.begin(), ie = data.CmSp_ITK.end();
+ std::vector<std::pair<float, InDet::SiSpacePointForSeedITK *>>::iterator it_startInnerLoop = it_commonTopSP;
+
+ /// check all possible common top SP
+ for (; it_commonTopSP != ie; ++it_commonTopSP)
+ {
+
+ SiSpacePointForSeedITK *SPt = (*it_commonTopSP).second;
+ /// the seed quality is set to d0 initially
+ float seedIP = SPt->param();
+ float seedQuality = seedIP + SPt->scorePenalty();
+ float originalSeedQuality = seedQuality;
+
+ if (m_maxdImpact > 50)
+ { //This only applies to LRT
+
+ float topR = SPt->radius();
+ float topZ = SPt->z();
+
+ float Zot = std::abs(topR - bottomR) > 10e-9 ? bottomZ - (bottomR - originalSeedQuality) * ((topZ - bottomZ) / (topR - bottomR)) : bottomZ;
+
+ float theta1 = std::abs(topR - bottomR) > 10e-9 ? std::atan2(topR - bottomR, topZ - bottomZ) : 0.;
+ float eta1 = theta1 > 0 ? -std::log(std::tan(.5 * theta1)) : 0.;
+
+ float theta0 = seedIP > 0 ? std::atan2(seedIP, Zot) : 0;
+ float eta0 = theta0 > 0 ? -std::log(std::tan(.5 * theta0)) : 0.;
+
+ float deltaEta = std::abs(eta1 - eta0); //For LLP daughters, the direction of the track is correlated with the direction of the LLP (which is correlated with the direction of the point of closest approach
+ //calculate weighted average of d0 and deltaEta, normalized by their maximum values
+ float f = std::min(0.5, originalSeedQuality / 200.); //0.5 and 200 are parameters chosen from a grid scan to optimize efficiency
+ seedQuality *= (1 - f) / 300.;
+ seedQuality += f * deltaEta / 2.5;
+ }
+
+ bool topSPisPixel = !SPt->spacepoint->clusterList().second;
+
+ /// check the surface the hit is on
+ const Trk::Surface *surfaceTopSP = SPt->sur();
+ float radiusTopSP = SPt->radius();
+ /// form a curvature interval cut
+ float minCurvature = (*it_commonTopSP).first - curvatureInterval;
+ float maxCurvature = (*it_commonTopSP).first + curvatureInterval;
+
+ /** Note: The score modifiers used here have the purpose of separating the candidates into
+ * classes / groups disjoint from each other.
+ * So the score increment (200 by default) should exceed the maximum
+ * |d0| (base score) we expect to encounter to avoid overlap.
+ * For LRT, we may want to tune this!
+ **/
+
+ /// if we have a SSS seed, boost the quality score by 400
+ if (!bottomSPisPixel)
+ seedQuality += m_seedScoreBonusSSS;
+ /// if we have a PPP, boost the quality by 200
+ else if (topSPisPixel)
+ seedQuality += m_seedScoreBonusPPP;
+
+ /**
+ * Now we look at the other SP candidates and try to find a confirmation seed,
+ * including the same centre/lower SP and giving a compatible curvature,
+ * but with the top SP in a different layer
+ **/
+
+ for (it_otherSP = it_startInnerLoop; it_otherSP != ie; ++it_otherSP)
+ {
+ /// if we are looking at the same SP, skip it
+ if (it_otherSP == it_commonTopSP)
+ continue;
+ /// if we have a lower curvature than the minimum, skip - and remember to
+ /// not bother with this candidate again later, as the vectors are curvature-sorted
+ if ((*it_otherSP).first < minCurvature)
+ {
+ it_startInnerLoop = it_otherSP;
+ ++it_startInnerLoop;
+ continue;
+ }
+ /// abort once the the curvature gets too large
+ if ((*it_otherSP).first > maxCurvature)
+ break;
+ /// if both SP are on the surface, skip it
+ if ((*it_otherSP).second->sur() == surfaceTopSP)
+ continue;
+ /// if the other SP is too close to the current top one, skip
+ float radiusOtherSP = (*it_otherSP).second->radius();
+ if (std::abs(radiusOtherSP - radiusTopSP) < m_drminSeedConf)
+ continue;
+ // if we have a confirmation seed, we improve the score of the seed.
+ seedQuality += m_seedScoreBonusConfirmationSeed;
+ // apply confirmation bonus only once
+ break;
+ }
+
+ /// kick this seed candidate if the score is too high (lower values = better)
+ if (seedQuality > data.maxScore)
+ continue;
+
+ /// if we have PPS seeds and no confirmation SP exists (which would give the -200 bonus)
+ /// or the hits on this seed were already used on a higher quality PPP/SSS seed, kick this one
+ if (bottomSPisPixel != topSPisPixel)
+ {
+ if (seedQuality > 0. ||
+ (seedQuality > bottomSPQuality && seedQuality > centralSPQuality && seedQuality > SPt->quality()))
+ continue;
+ }
+ /// If we have a non-confirmed seed, apply a stricter d0 cut.
+ /// This, is determined using the original cut and the score penalty modifier.
+ if (!isConfirmedSeed(SPb, SPt, seedQuality))
+ {
+ /// PPP seeds
+ double maxdImpact = m_maxdImpact - (m_dImpactCutSlopeUnconfirmedPPP * SPt->scorePenalty());
+ /// SSS seeds
+ if (!bottomSPisPixel)
+ maxdImpact = m_maxdImpactSSS - (m_dImpactCutSlopeUnconfirmedSSS * SPt->scorePenalty());
+ if (seedIP > maxdImpact)
+ continue;
+ }
+ /// this is a good seed, save it (unless we have too many seeds per SP)
+ newOneSeed(data, SPb, SP0, SPt, Zob, seedQuality);
+ } ///< end of loop over top SP candidates
+ data.CmSp_ITK.clear();
+}
+
+bool InDet::SiSpacePointsSeedMaker_ITK::isConfirmedSeed(const InDet::SiSpacePointForSeedITK *bottomSP,
+ const InDet::SiSpacePointForSeedITK *topSP, float quality) const
+{
+
+ /// SSS seeds
+ if (bottomSP->spacepoint->clusterList().second)
+ {
+ return (quality < m_seedScoreThresholdSSSConfirmationSeed);
+ }
+ /// PPP seeds
+ else if (!topSP->spacepoint->clusterList().second)
+ {
+ return (quality < m_seedScoreThresholdPPPConfirmationSeed);
+ }
+ /// PPS: the confirmation is the only quality modifier applied
+ else
+ return (quality < 0.);
+}
+
+void InDet::SiSpacePointsSeedMaker_ITK::writeNtuple(const SiSpacePointsSeed* seed, const Trk::Track* track, int seedType, long eventNumber) const
+{
+ if(m_writeNtuple) {
+ std::lock_guard<std::mutex> lock(m_mutex);
+
+ if(track != nullptr) {
+ m_trackPt = (track->trackParameters()->front()->pT())/1000.f;
+ m_trackEta = std::abs(track->trackParameters()->front()->eta());
+ }
+ else {
+ m_trackPt = -1.;
+ m_trackEta = -1.;
+ }
+ m_d0 = seed->d0();
+ m_z0 = seed->zVertex();
+ m_eta = seed->eta();
+ m_x1 = seed->x1();
+ m_x2 = seed->x2();
+ m_x3 = seed->x3();
+ m_y1 = seed->y1();
+ m_y2 = seed->y2();
+ m_y3 = seed->y3();
+ m_z1 = seed->z1();
+ m_z2 = seed->z2();
+ m_z3 = seed->z3();
+ m_r1 = seed->r1();
+ m_r2 = seed->r2();
+ m_r3 = seed->r3();
+ m_type = seedType;
+ m_dzdr_b = seed->dzdr_b();
+ m_dzdr_t = seed->dzdr_t();
+ m_pt = seed->pt();
+ m_givesTrack = !(track == nullptr);
+ m_eventNumber = eventNumber;
+
+ m_outputTree->Fill();
+
+ }
+
+}
+
+bool InDet::SiSpacePointsSeedMaker_ITK::getWriteNtupleBoolProperty() const
+{
+ return m_writeNtuple;
}