ATLAS Offline Software
Loading...
Searching...
No Matches
MuonPhiHitSelector Class Reference

#include <MuonPhiHitSelector.h>

Inheritance diagram for MuonPhiHitSelector:
Collaboration diagram for MuonPhiHitSelector:

Public Member Functions

 MuonPhiHitSelector (const std::string &, const std::string &, const IInterface *)
virtual ~MuonPhiHitSelector ()=default
StatusCode initialize () override
std::vector< std::unique_ptr< const Trk::MeasurementBase > > select_rio (const double pmom, const std::vector< const Trk::RIO_OnTrack * > &associatedHits, const std::vector< const Trk::PrepRawData * > &unassociatedHits) const override
 Selects and builds a cleaned vector of RIO fits the associatedHits and build new RIOs, if m_competingRios true then for ambiguous hits competing rios are built.
ServiceHandle< StoreGateSvc > & evtStore ()
 The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.
const ServiceHandle< StoreGateSvc > & detStore () const
 The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
virtual StatusCode sysInitialize () override
 Perform system initialization for an algorithm.
virtual StatusCode sysStart () override
 Handle START transition.
virtual std::vector< Gaudi::DataHandle * > inputHandles () const override
 Return this algorithm's input handles.
virtual std::vector< Gaudi::DataHandle * > outputHandles () const override
 Return this algorithm's output handles.
Gaudi::Details::PropertyBase & declareProperty (Gaudi::Property< T, V, H > &t)
void updateVHKA (Gaudi::Details::PropertyBase &)
MsgStream & msg () const
bool msgLvl (const MSG::Level lvl) const

Static Public Member Functions

static const InterfaceID & interfaceID ()

Protected Member Functions

void renounceArray (SG::VarHandleKeyArray &handlesArray)
 remove all handles from I/O resolution
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce (T &h)
void extraDeps_update_handler (Gaudi::Details::PropertyBase &ExtraDeps)
 Add StoreName to extra input/output deps as needed.

Private Types

typedef ServiceHandle< StoreGateSvcStoreGateSvc_t

Private Member Functions

void fitRecPhi (const double pmom, const std::vector< Identifier > &phiId, const std::vector< double > &phiHitx, const std::vector< double > &phiHity, const std::vector< double > &phiHitz, const std::vector< double > &phiError, std::vector< int > &quality, const int nphi, std::vector< double > &phiPull, std::vector< int > &phiMult, std::vector< int > &phiSelect, double &chi2, double &r0, double &phi, std::vector< double > &errorM, int &nfit) const
 fit method curved track model
void fitPhiSL (const double pmom, const std::vector< Identifier > &id, const std::vector< double > &hitx, const std::vector< double > &hity, const std::vector< double > &hitz, const std::vector< double > &error, std::vector< int > &select, const int n, std::vector< double > &pull, int &imax, double &chi2, double &r0, double &phi, std::vector< double > &errorM, bool fast) const
 fit method straight line model
void clusterPhi (const std::vector< Identifier > &id, const std::vector< double > &hitx, const std::vector< double > &hity, const std::vector< double > &hitz, const std::vector< double > &error, const std::vector< double > &pull, std::vector< int > &select, const int n, std::vector< double > &clusterX, std::vector< double > &clusterY, std::vector< double > &clusterZ, std::vector< double > &clusterError, std::vector< Identifier > &clusterId, std::vector< int > &clusterHits, std::vector< int > &clusterSelect, std::vector< int > &clusterInt, int &ncl) const
 clusterization method
Gaudi::Details::PropertyBase & declareGaudiProperty (Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
 specialization for handling Gaudi::Property<SG::VarHandleKey>

Private Attributes

ServiceHandle< Muon::IMuonIdHelperSvcm_idHelperSvc
ToolHandle< Muon::IMuonCompetingClustersOnTrackCreatorm_competingRIOsOnTrackTool
 Toolhandle to CompetingRIOsOnTrackTool creator.
ToolHandle< Muon::IMuonClusterOnTrackCreatorm_clusterCreator {this, "MuonClusterOnTrackCreator",""}
 Toolhandle to ClusterOnTrackTool creator.
ToolHandle< Muon::IMuonClusterOnTrackCreatorm_cscRotCreator {this, "CscRotCreator",""}
Gaudi::Property< bool > m_summary {this, "DoSummary", false}
 flag to print out a summary of what comes in and what comes out
Gaudi::Property< bool > m_cosmics {this, "DoCosmics", false}
 flag for use of cosmics, straight line model will be used, no interaction point constraint
Gaudi::Property< bool > m_makeClusters {this, "MakeClusters", false}
 flag that performs a clusterization and return clusters (default: false)
Gaudi::Property< bool > m_competingRios {this, "CompetingRios", false}
 flag that build competing rios on track for amibguous trigger hits (default: false)
StoreGateSvc_t m_evtStore
 Pointer to StoreGate (event store by default)
StoreGateSvc_t m_detStore
 Pointer to StoreGate (detector store by default)
std::vector< SG::VarHandleKeyArray * > m_vhka
bool m_varHandleArraysDeclared

Detailed Description

Definition at line 21 of file MuonPhiHitSelector.h.

Member Typedef Documentation

◆ StoreGateSvc_t

typedef ServiceHandle<StoreGateSvc> AthCommonDataStore< AthCommonMsg< AlgTool > >::StoreGateSvc_t
privateinherited

Definition at line 388 of file AthCommonDataStore.h.

Constructor & Destructor Documentation

◆ MuonPhiHitSelector()

MuonPhiHitSelector::MuonPhiHitSelector ( const std::string & type,
const std::string & name,
const IInterface * parent )

Definition at line 29 of file MuonPhiHitSelector.cxx.

30 : AthAlgTool(type, name, parent)
31{
32 declareInterface<IMuonHitSelector>(this);
33
34
35}
AthAlgTool()
Default constructor:

◆ ~MuonPhiHitSelector()

virtual MuonPhiHitSelector::~MuonPhiHitSelector ( )
virtualdefault

Member Function Documentation

◆ clusterPhi()

void MuonPhiHitSelector::clusterPhi ( const std::vector< Identifier > & id,
const std::vector< double > & hitx,
const std::vector< double > & hity,
const std::vector< double > & hitz,
const std::vector< double > & error,
const std::vector< double > & pull,
std::vector< int > & select,
const int n,
std::vector< double > & clusterX,
std::vector< double > & clusterY,
std::vector< double > & clusterZ,
std::vector< double > & clusterError,
std::vector< Identifier > & clusterId,
std::vector< int > & clusterHits,
std::vector< int > & clusterSelect,
std::vector< int > & clusterInt,
int & ncl ) const
private

clusterization method

Use hits (select > 0) and pulls to make clusters

Inputs id = identifiers hits hitx hity hitz = position in space error = associated error (in x-y plane) pull (from fit)= residual (hit -fit) /error select = > 0 for selected hits n = total number of hits fast = true = fast fit without scattering centres and no error propagation false = fit with scattering centres and error propagation

Outputs clusterX Y Z = cluster positions clusterError = errors clusterId = cluster identifier (smallest pull) clusterHits = number of hits per cluster ncl = number of clusters chi2 = total chi2 r0 = perigee parameter of fit (0,0) phi = azimuthal angle of fit at perigee

Definition at line 302 of file MuonPhiHitSelector.cxx.

310{
311
312 //
313 // Use hits (select > 0) and pulls to make clusters
314 //
315 //
316 // Inputs
317 // id = identifiers hits
318 // hitx hity hitz = position in space
319 // error = associated error (in x-y plane)
320 // pull (from fit)= residual (hit -fit) /error
321 // select = > 0 for selected hits
322 // n = total number of hits
323
324 // Outputs
325 // clusterX Y Z = cluster positions
326 // clusterError = errors
327 // clusterId = cluster identifier (smallest pull)
328 // clusterHits = number of hits per cluster
329 // ncl = number of clusters
330 // chi2 = total chi2
331 // r0 = perigee parameter of fit (0,0)
332 // phi = azimuthal angle of fit at perigee
333
334
335 ATH_MSG_DEBUG("Start phi clustering");
336
337 ncl = 0;
338 if (n == 0) return;
339
340 std::vector<int> scode(n);
341
342 for (int i = 0; i < n; ++i) {
343 Identifier idi = id[i];
344 int code = 0;
345 if (m_idHelperSvc->isRpc(idi)) {
346 int doubZ = m_idHelperSvc->rpcIdHelper().doubletZ(idi);
347 int doubPhi = m_idHelperSvc->rpcIdHelper().doubletPhi(idi);
348 code = 100000000 * (m_idHelperSvc->rpcIdHelper().stationName(idi))
349 + 1000000 * (m_idHelperSvc->rpcIdHelper().stationPhi(idi))
350 + 10000 * ((m_idHelperSvc->rpcIdHelper().stationEta(idi)) + 1000);
351 code += 1000 * (doubZ - 1) + 100 * (doubPhi - 1);
352 code += 2 * ((m_idHelperSvc->rpcIdHelper().doubletR(idi)) - 1)
353 + 16 * ((m_idHelperSvc->rpcIdHelper().gasGap(idi)) - 1);
354 } else if (m_idHelperSvc->isTgc(idi)) {
355 code = 1000000 * (m_idHelperSvc->tgcIdHelper().stationName(idi))
356 + 10000 * (m_idHelperSvc->tgcIdHelper().stationPhi(idi))
357 + 100 * ((m_idHelperSvc->tgcIdHelper().stationEta(idi)) + 10);
358 code = code + m_idHelperSvc->tgcIdHelper().gasGap(idi);
359 } else if (m_idHelperSvc->isCsc(idi)) {
360 code = 1000000 * (m_idHelperSvc->cscIdHelper().stationName(idi))
361 + 10000 * (m_idHelperSvc->cscIdHelper().stationPhi(idi))
362 + 100 * ((m_idHelperSvc->cscIdHelper().stationEta(idi)) + 10);
363 code = code + m_idHelperSvc->cscIdHelper().wireLayer(idi);
364 }
365 scode[i] = code;
366 }
367
368 // std::vector<int> clusterInt(n);
369
370 for (int i = 0; i < n; ++i) {
371 clusterInt[i] = -1;
372 }
373
374 int icl = -1;
375 for (int i = 0; i < n; ++i) {
376 if (error[i] != 0 && select[i] > 0 && clusterInt[i] == -1) {
377 icl++;
378 clusterInt[i] = icl;
379 for (int j = i + 1; j < n; ++j) {
380 if (clusterInt[j] == -1) {
381 if (error[j] != 0 && select[j] > 0) {
382 if (scode[i] == scode[j]) clusterInt[j] = icl;
383 }
384 }
385 }
386 }
387 }
388
389 std::vector<double> clusterCommon2Error(icl + 1);
390 std::vector<int> clusterCode(icl + 1);
391 ncl = icl + 1;
392 for (int ic = 0; ic < icl + 1; ++ic) {
393 clusterX[ic] = 0.;
394 clusterY[ic] = 0.;
395 clusterZ[ic] = 0.;
396 clusterError[ic] = 0.;
397 clusterCommon2Error[ic] = 0.;
398 clusterHits[ic] = 0;
399 clusterCode[ic] = 0;
400 clusterSelect[ic] = 0;
401 double pullMax = 10.;
402 for (int i = 0; i < n; ++i) {
403 if (select[i] > 0) {
404 if (ic == clusterInt[i]) {
405 clusterSelect[ic] = select[i];
406 double w = 1. / (error[i] * error[i]);
407 clusterX[ic] += hitx[i] * w;
408 clusterY[ic] += hity[i] * w;
409 clusterZ[ic] += hitz[i] * w;
410 clusterError[ic] += w;
411 if (std::abs(pull[i]) < std::abs(pullMax)) {
412 pullMax = pull[i];
413 clusterId[ic] = id[i];
414 clusterCode[ic] = scode[i];
415 clusterSelect[ic] = select[i];
416 }
417 clusterHits[ic]++;
418 if (clusterHits[ic] == 1) clusterCommon2Error[ic] = 0.;
419 }
420 }
421 }
422 clusterX[ic] = clusterX[ic] / clusterError[ic];
423 clusterY[ic] = clusterY[ic] / clusterError[ic];
424 clusterZ[ic] = clusterZ[ic] / clusterError[ic];
425 // Don't assume improvement on errors due to clustering
426 clusterError[ic] = std::sqrt(clusterHits[ic]) / std::sqrt(clusterError[ic]);
427 {
428 ATH_MSG_DEBUG("cluster phi " << ic << " x " << clusterX[ic] << " y " << clusterY[ic] << " z "
429 << clusterZ[ic] << " error " << clusterError[ic] << " hits " << clusterHits[ic]
430 << " select " << clusterSelect[ic] << " Code " << clusterCode[ic]);
431 }
432 }
433}
#define ATH_MSG_DEBUG(x)
ServiceHandle< Muon::IMuonIdHelperSvc > m_idHelperSvc
void select(const xAOD::IParticle *particle, const float coneSize, const xAOD::CaloClusterContainer *clusters, std::vector< bool > &mask)
constexpr uint8_t stationPhi
station Phi 1 to 8
int ic
Definition grepfile.py:33

◆ declareGaudiProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareGaudiProperty ( Gaudi::Property< T, V, H > & hndl,
const SG::VarHandleKeyType &  )
inlineprivateinherited

specialization for handling Gaudi::Property<SG::VarHandleKey>

Definition at line 156 of file AthCommonDataStore.h.

158 {
160 hndl.value(),
161 hndl.documentation());
162
163 }
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)

◆ declareProperty()

Gaudi::Details::PropertyBase & AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty ( Gaudi::Property< T, V, H > & t)
inlineinherited

Definition at line 145 of file AthCommonDataStore.h.

145 {
146 typedef typename SG::HandleClassifier<T>::type htype;
148 }
Gaudi::Details::PropertyBase & declareGaudiProperty(Gaudi::Property< T, V, H > &hndl, const SG::VarHandleKeyType &)
specialization for handling Gaudi::Property<SG::VarHandleKey>

◆ detStore()

const ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::detStore ( ) const
inlineinherited

The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 95 of file AthCommonDataStore.h.

◆ evtStore()

ServiceHandle< StoreGateSvc > & AthCommonDataStore< AthCommonMsg< AlgTool > >::evtStore ( )
inlineinherited

The standard StoreGateSvc (event store) Returns (kind of) a pointer to the StoreGateSvc.

Definition at line 85 of file AthCommonDataStore.h.

◆ extraDeps_update_handler()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::extraDeps_update_handler ( Gaudi::Details::PropertyBase & ExtraDeps)
protectedinherited

Add StoreName to extra input/output deps as needed.

use the logic of the VarHandleKey to parse the DataObjID keys supplied via the ExtraInputs and ExtraOuputs Properties to add the StoreName if it's not explicitly given

◆ fitPhiSL()

void MuonPhiHitSelector::fitPhiSL ( const double pmom,
const std::vector< Identifier > & id,
const std::vector< double > & hitx,
const std::vector< double > & hity,
const std::vector< double > & hitz,
const std::vector< double > & error,
std::vector< int > & select,
const int n,
std::vector< double > & pull,
int & imax,
double & chi2,
double & r0,
double & phi,
std::vector< double > & errorM,
bool fast ) const
private

fit method straight line model

Definition at line 812 of file MuonPhiHitSelector.cxx.

817{
818
819 // Perform straight line fit to hits: good hits have select > 0
820 // in the fit scattering centres are added for nfit-1 angles
821 // WITH beamspot constraint
822 // degrees of freedom = 2*nfit
823
824 // Fit is based on matrix inversions formulae
825
826 // Inputs pmom = estimate of momentum
827 // id = identifiers hits
828 // hitx hity hitz = position in space
829 // error = associated error (in x-y plane)
830 // select = > 0 for selected hits
831 // n = total number of hits
832
833
834 // Outputs pull = residual (hit -fit) /error
835 // imax = index for hit with maximum pull
836 // chi2 = total chi2
837 // r0 = perigee parameter of fit (0,0)
838 // phi = azimuthal angle of fit at perigee
839
840 double pest = pmom;
841 if (pest > 20000.) pest = 20000.;
842
843 r0 = 0.;
844 phi = 0.;
845 chi2 = 0.;
846 imax = 0;
847
848 // Calculate mean position
849 double xm = 0.;
850 double ym = 0.;
851 double dtot = 0.;
852 double em = 0.;
853 for (int i = 0; i < n; ++i) {
854 if (error[i] != 0 && select[i] > 0) {
855 double inver2 = 1. / (error[i] * error[i]);
856 xm += hitx[i] * inver2;
857 ym += hity[i] * inver2;
858 dtot += std::sqrt(hitx[i] * hitx[i] + hity[i] * hity[i] + hitz[i] * hitz[i]) * inver2;
859 em += inver2;
860 }
861 }
862
863 if (em == 0) return;
864
865 dtot = dtot / em;
866
867 // Beamspot error 10 mm for cosmics 10000
868
869 double ebs = 0.1;
870 if (m_cosmics) ebs = 10000.;
871
872 ATH_MSG_DEBUG("pmom " << pmom << " error beam " << ebs);
873 double ebs2 = ebs * ebs;
874 double invebs2 = 1. / ebs2;
875 double xmc = xm / (em + invebs2);
876 double ymc = ym / (em + invebs2);
877 xm = xm / em;
878 ym = ym / em;
879
880 // Constraint on beam spot
881
882 double len2 = xmc * xmc + ymc * ymc;
883 double xcc = len2 * xmc * invebs2;
884 double ycc = len2 * ymc * invebs2;
885
886 for (int i = 0; i < n; ++i) {
887 if (error[i] != 0 && select[i] > 0) {
888 double inver2 = 1. / (error[i] * error[i]);
889 double xdiff = hitx[i] - xmc;
890 double ydiff = hity[i] - ymc;
891 double xdiff2 = xdiff * xdiff;
892 double ydiff2 = ydiff * ydiff;
893 len2 = xdiff2 + ydiff2;
894 double sign = 1.;
895 // Non Cosmics assume IP at 0 0
896 if (xdiff * hitx[i] + ydiff * hity[i] < 0 && !m_cosmics) sign = -1;
897 // Cosmics assume down going
898 if (ydiff < 0 && m_cosmics) sign = -1;
899 xcc += len2 * sign * xdiff * inver2;
900 ycc += len2 * sign * ydiff * inver2;
901 }
902 }
903
904 if (em > 0) phi = std::atan2(ycc, xcc);
905 CxxUtils::sincos scphi(phi);
906
907 r0 = xmc * scphi.sn - ymc * scphi.cs;
908 double x0 = r0 * scphi.sn;
909 double y0 = -r0 * scphi.cs;
910
911 if (msgLvl(MSG::DEBUG))
912 ATH_MSG_DEBUG("Constraint r0 " << r0 << " xpos " << xmc << " ypos " << ymc << " phi " << phi);
913 // assume 0,0
914 std::vector<double> d(n);
915 std::vector<double> dist(n);
916 std::map<double, int> distanceSort;
917 double pullmax = 0.;
918 for (int i = 0; i < n; ++i) {
919 if (error[i] != 0 && select[i] > 0) {
920 double xdiff = hitx[i] - x0;
921 double ydiff = hity[i] - y0;
922 double xdiff2 = xdiff * xdiff;
923 double ydiff2 = ydiff * ydiff;
924 d[i] = std::sqrt(xdiff2 + ydiff2);
925 dist[i] = std::sqrt(xdiff2 + ydiff2 + hitz[i] * hitz[i]);
926 distanceSort[dist[i]] = i;
927 pull[i] = hitx[i] * scphi.sn - hity[i] * scphi.cs - r0;
928 if (std::abs(pull[i]) > std::abs(pullmax)) {
929 pullmax = pull[i];
930 imax = i;
931 }
932 }
933 }
934
935 if (fast) return;
936
937 std::map<double, int>::iterator it = distanceSort.begin();
938 std::map<double, int>::iterator it_end = distanceSort.end();
939
940 int nfit = 0;
941 std::vector<double> xf(2 * n);
942 std::vector<double> lf(2 * n);
943 std::vector<double> yf(2 * n);
944 std::vector<double> ef(2 * n);
945 std::vector<int> indexf(n);
946 //
947 // measurements yf error ef at distance xf (0:nfit)
948 // scattering centra angle zero yf error ef at distance xf(nfit+1..2nfit-1)
949 // beamspot at yf(2 nfit) = 0 error ebs2 at distance xf(2 nfit)
950
951 for (; it != it_end; ++it) {
952 int index = it->second;
953 xf[nfit] = d[index];
954 lf[nfit] = dist[index];
955 yf[nfit] = (hitx[index] - xmc) * scphi.sn - (hity[index] - ymc) * scphi.cs;
956 ef[nfit] = error[index];
957 indexf[nfit] = index;
958 nfit++;
959 }
960
961 // NB start at 1 to add scattering centra
962
963 double erang = 0.030 * 5000. / (pest + 1000.);
964 for (int i = 1; i < nfit; ++i) {
965 xf[nfit + i - 1] = xf[i - 1];
966 yf[nfit + i - 1] = 0.;
967 double scale = 1.;
968 if (select[i] == 1)
969 scale = 1.;
970 else if (select[i] == 3)
971 scale = 0.5;
972 else if (select[i] == 2)
973 scale = 2.5;
974 ef[nfit + i - 1] = scale * erang;
975 }
976 // Beamspot
977 yf[2 * nfit - 1] = 0.;
978 xf[2 * nfit - 1] = 0.;
979 ef[2 * nfit - 1] = ebs;
980
981 Amg::MatrixX v(nfit + 1, 1);
982 v.setIdentity();
983
984 if (msgLvl(MSG::DEBUG))
985 ATH_MSG_DEBUG("fitPhiSL "
986 << " nfit " << nfit);
987
988 for (int i = 0; i < nfit + 1; ++i) {
989 v(i, 0) = 0.;
990 for (int j = 0; j < nfit; ++j) {
991 double inver2 = 1. / (ef[j] * ef[j]);
992 if (i == 0)
993 v(i, 0) += yf[j] * inver2;
994 else if (i == 1)
995 v(i, 0) += yf[j] * xf[j] * inver2;
996 else if (i > 1 && j > i - 2) {
997 v(i, 0) += yf[j] * (lf[j] - lf[i - 2]) * inver2;
998 }
999 }
1000 }
1001
1002 // Track Model Matrix
1003
1004 Amg::MatrixX model(nfit + 1, 2 * nfit);
1005 model.setIdentity();
1006 // Measurements related to position and slope
1007
1008 for (int i = 0; i < nfit + 1; ++i) {
1009 for (int j = 0; j < nfit; ++j) {
1010 model(i, j) = 0.;
1011 if (i == 0)
1012 model(i, j) = 1.;
1013 else if (i == 1)
1014 model(i, j) = xf[j];
1015 // scattering angle
1016 else if (i > 1 && j > i - 2)
1017 model(i, j) = lf[j] - lf[i - 2];
1018 }
1019 }
1020
1021 // Constraints on scattering angles and beamspot
1022
1023 for (int i = 0; i < nfit + 1; ++i) {
1024 for (int j = nfit; j < 2 * nfit; ++j) {
1025 model(i, j) = 0.;
1026 // scattering angle
1027 if (i == j - nfit + 2) model(i, j) = 1.;
1028 // Beam spot
1029 if (i == 0 && j == 2 * nfit - 1) model(i, j) = 1.;
1030 }
1031 }
1032
1033 // Covariance Inverse of Track parameters
1034
1035 Amg::MatrixX covT(nfit + 1, nfit + 1);
1036 for (int i = 0; i < nfit + 1; ++i) {
1037 for (int j = 0; j < nfit + 1; ++j) {
1038 covT(i, j) = 0.;
1039 for (int k = 0; k < 2 * nfit; ++k) {
1040 double er2 = ef[k] * ef[k];
1041 covT(i, j) += model(i, k) * model(j, k) / er2;
1042 }
1043 }
1044 }
1045
1046 // Invert covariance matrix and replace it (should be CovT)
1047 Amg::MatrixX covTI = covT.inverse();
1048
1049 Amg::MatrixX t(nfit + 1, 1);
1050 // Solution for Track parameters
1051 t = covTI * v;
1052
1053 if (msgLvl(MSG::DEBUG) && std::abs(t(1, 0)) > 0.2) {
1054 ATH_MSG_DEBUG("Don't trust fit result " << t(1, 0) << " Keep Old result");
1055 }
1056 if (std::abs(t(1, 0)) > 0.2) return;
1057
1058 // calculate residuals and chi2
1059 std::vector<double> resi(2 * nfit);
1060 std::vector<double> pulli(2 * nfit);
1061 std::vector<double> errf(2 * nfit);
1062 std::vector<double> pullf(2 * nfit);
1063 std::vector<double> resiOut(2 * nfit);
1064 std::vector<double> pullOut(2 * nfit);
1065 pullmax = 0.;
1066 int jmax = 0;
1067
1068 errorM[0] = covTI(0, 0);
1069 errorM[1] = covTI(0, 1);
1070 errorM[2] = covTI(1, 1);
1071 errorM[3] = 0.;
1072 if (nfit > 2) {
1073 double invlt = 1. / (lf[nfit - 1] - lf[1]);
1074 for (int i = 1; i < nfit - 1; ++i) {
1075 double w = (lf[nfit - 1] - lf[i]) * invlt;
1076 errorM[3] += covTI(i + 1, i + 1) * w * w;
1077 }
1078 }
1079
1080 {
1081 if (nfit >= 3) {
1082 ATH_MSG_DEBUG("Error angle " << covTI(3, 3));
1083 } // covTI has dim nfit+1
1084 ATH_MSG_DEBUG("errorM[3] " << errorM[3]);
1085 }
1086
1087 for (int i = 0; i < 2 * nfit; ++i) {
1088
1089 // Calculate prediction at each measurement i
1090 // propagate error of track parameters to measurement i
1091 double error2 = 0.;
1092 double ypred = 0.;
1093 for (int j = 0; j < nfit + 1; ++j) {
1094 if (msgLvl(MSG::DEBUG) && i == 0) ATH_MSG_DEBUG("Parameter j " << j << " t(j,0) " << t(j, 0));
1095 if (msgLvl(MSG::DEBUG) && model(j, i) != 0) ATH_MSG_DEBUG("i " << i << " model ij " << model(j, i));
1096 ypred += model(j, i) * t(j, 0);
1097 for (int k = 0; k < nfit + 1; ++k) {
1098 error2 += model(j, i) * covTI(j, k) * model(k, i);
1099 }
1100 }
1101 double ef_i2 = ef[i] * ef[i];
1102 double inv_ef_i2 = 1. / ef_i2;
1103
1104 resi[i] = ypred - yf[i];
1105 pulli[i] = resi[i] / ef[i];
1106
1107 // errf propagated error and pullf
1108 errf[i] = std::sqrt(error2);
1109 pullf[i] = resi[i] / errf[i];
1110
1111 // calculate residual without hit and error without hit
1112 // Think of Kalmanm method to exclude hit and error
1113 double err2invOut = 1. / error2 - inv_ef_i2;
1114 if (err2invOut > 0) {
1115 resiOut[i] = (ypred / error2 - yf[i] * inv_ef_i2) / err2invOut - yf[i];
1116 pullOut[i] = resiOut[i] / std::sqrt(1. / err2invOut + ef_i2);
1117 }
1118
1119 if ((i < nfit) and (std::abs(pullOut[i]) > std::abs(pullmax)) ) {
1120 imax = indexf[i];
1121 jmax = i;
1122 pullmax = pullOut[i];
1123 }
1124 chi2 += resi[i] * resi[i] * inv_ef_i2;
1125
1126 if (i < nfit) {
1127 pull[indexf[i]] = pullOut[i];
1128 }
1129 if (msgLvl(MSG::DEBUG) && i < nfit)
1130 ATH_MSG_DEBUG("i " << i << " index " << indexf[i] << " det " << select[indexf[i]] << " ypred " << ypred
1131 << " mst " << yf[i] << " residual " << resi[i] << " error " << ef[i] << " dist "
1132 << dist[i] << " hitz " << hitz[i] << " Pull " << pulli[i] << " Pullf " << pullf[i]
1133 << " resi out " << resiOut[i] << " pull out " << pullOut[i]);
1134 if (msgLvl(MSG::DEBUG) && i > nfit)
1135 ATH_MSG_DEBUG("i " << i << " ypred " << ypred << " mst " << yf[i] << " residual " << resi[i] << " error "
1136 << ef[i]);
1137 }
1138 r0 = r0 + t(0, 0);
1139 phi = phi + t(1, 0);
1140
1141 ATH_MSG_DEBUG("delta phi " << t(1, 0));
1142 if (msgLvl(MSG::DEBUG) && std::abs(t(1, 0)) > 0.1) ATH_MSG_DEBUG("ALARM delta phi " << t(1, 0));
1143
1144 if (msgLvl(MSG::DEBUG))
1145 ATH_MSG_DEBUG("Track parameters r0 " << r0 << " phi " << phi << " chi2 " << chi2 << " jmax " << jmax << " imax "
1146 << imax << " pullmax " << pullmax);
1147}
Scalar phi() const
phi method
int sign(int a)
int imax(int i, int j)
bool msgLvl(const MSG::Level lvl) const
Gaudi::Property< bool > m_cosmics
flag for use of cosmics, straight line model will be used, no interaction point constraint
double chi2(TH1 *h0, TH1 *h1)
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
str index
Definition DeMoScan.py:362
const double r0
electron radius{cm}

◆ fitRecPhi()

void MuonPhiHitSelector::fitRecPhi ( const double pmom,
const std::vector< Identifier > & phiId,
const std::vector< double > & phiHitx,
const std::vector< double > & phiHity,
const std::vector< double > & phiHitz,
const std::vector< double > & phiError,
std::vector< int > & quality,
const int nphi,
std::vector< double > & phiPull,
std::vector< int > & phiMult,
std::vector< int > & phiSelect,
double & chi2,
double & r0,
double & phi,
std::vector< double > & errorM,
int & nfit ) const
private

fit method curved track model

Definition at line 435 of file MuonPhiHitSelector.cxx.

441{
442
443 //
444 // Use reconstructed hits to perform fit for phi
445 //
446
447 ATH_MSG_DEBUG("Start phi fit reconstruction");
448
449 chi2 = 0.;
450 r0 = 0.;
451 nfit = 0;
452 phi = 0.;
453
454 int ncsc = 0;
455 int ntgc = 0;
456 int nrpc = 0;
457
458
459 if (nphi == 0) return;
460
461 std::vector<double> error0(nphi);
462 std::vector<double> error(nphi);
463 std::vector<double> errorf(nphi);
464 std::vector<int> scode(nphi);
465 std::vector<int> srcode(nphi);
466 std::vector<int> phiSelectKeep(nphi);
467 std::map<int, int> clusters;
468 std::map<int, int> clustersr;
469 std::map<int, int> clusterspat;
470
471 for (int i = 0; i < nphi; ++i) {
472
473 Identifier idi = phiId[i];
474 int code = 0;
475 int rcode = 0;
476 if (m_idHelperSvc->isRpc(idi)) {
477 code = 1000000 * (m_idHelperSvc->rpcIdHelper().stationName(idi))
478 + 10000 * (m_idHelperSvc->rpcIdHelper().stationPhi(idi))
479 + 100 * ((m_idHelperSvc->rpcIdHelper().stationEta(idi)) + 10);
480 code = code + 2 * ((m_idHelperSvc->rpcIdHelper().doubletR(idi)) - 1)
481 + 16 * ((m_idHelperSvc->rpcIdHelper().gasGap(idi)) - 1);
482 rcode = 1000000 * (m_idHelperSvc->rpcIdHelper().stationName(idi))
483 + 10000 * (m_idHelperSvc->rpcIdHelper().stationPhi(idi))
484 + 0 * ((m_idHelperSvc->rpcIdHelper().stationEta(idi)) + 10);
485 rcode = rcode + 2 * ((m_idHelperSvc->rpcIdHelper().doubletR(idi)) - 1)
486 + 16 * ((m_idHelperSvc->rpcIdHelper().gasGap(idi)) - 1);
487 } else if (m_idHelperSvc->isTgc(idi)) {
488 code = 1000000 * (m_idHelperSvc->tgcIdHelper().stationName(idi))
489 + 10000 * (m_idHelperSvc->tgcIdHelper().stationPhi(idi))
490 + 100 * ((m_idHelperSvc->tgcIdHelper().stationEta(idi)) + 10);
491 code = code + m_idHelperSvc->tgcIdHelper().gasGap(idi);
492 rcode = 1000000 * (m_idHelperSvc->tgcIdHelper().stationName(idi))
493 + 10000 * (m_idHelperSvc->tgcIdHelper().stationPhi(idi))
494 + 0 * ((m_idHelperSvc->tgcIdHelper().stationEta(idi)) + 10);
495 rcode = rcode + m_idHelperSvc->tgcIdHelper().gasGap(idi);
496 } else if (m_idHelperSvc->isCsc(idi)) {
497 code = 1000000 * (m_idHelperSvc->cscIdHelper().stationName(idi))
498 + 10000 * (m_idHelperSvc->cscIdHelper().stationPhi(idi))
499 + 100 * ((m_idHelperSvc->cscIdHelper().stationEta(idi)) + 10);
500 code = code + m_idHelperSvc->cscIdHelper().wireLayer(idi);
501 rcode = 1000000 * (m_idHelperSvc->cscIdHelper().stationName(idi))
502 + 10000 * (m_idHelperSvc->cscIdHelper().stationPhi(idi))
503 + 0 * ((m_idHelperSvc->cscIdHelper().stationEta(idi)) + 10);
504 rcode = rcode + m_idHelperSvc->cscIdHelper().wireLayer(idi);
505 }
506
507 scode[i] = code;
508 srcode[i] = rcode;
509 int idet = 0;
510 if (m_idHelperSvc->isRpc(idi))
511 idet = 1;
512 else if (m_idHelperSvc->isTgc(idi))
513 idet = 2;
514 else if (m_idHelperSvc->isCsc(idi))
515 idet = 3;
516 phiSelect[i] = idet;
517 phiSelectKeep[i] = idet;
518 }
519 // Hits on segments
520 for (int i = 0; i < nphi; ++i) {
521 if (phiError[i] != 0 && quality[i] > 100) {
522 clusters[scode[i]]++;
523 clustersr[srcode[i]]++;
524 }
525 }
526 // Drop hits on patterns that are in same station and layer as segment hit
527 // Avoid adding again (solved) ambiguous hits
528
529 for (int i = 0; i < nphi; ++i) {
530 if (phiError[i] != 0 && quality[i] > 0 && quality[i] < 100) {
531 if (clustersr.count(srcode[i]) > 0) {
532 quality[i] = 0;
533 } else {
534 clusterspat[scode[i]]++;
535 }
536 }
537 }
538
539 // Assign errors according to multiplicities
540 if (msgLvl(MSG::DEBUG))
541 ATH_MSG_DEBUG("phi hits " << nphi << " segment clusters " << clusters.size() << " pattern clusters "
542 << clusterspat.size());
543
544 for (int i = 0; i < nphi; ++i) {
545 error0[i] = 0;
546 Identifier id = phiId[i];
547 phiMult[i] = 0;
548 if (phiError[i] != 0 && quality[i] > 0) {
549 int n = 0;
550 if (quality[i] > 100) {
551 n = clusters[scode[i]];
552 // Treat phi hits from segment with high multiplicity > 10 as lying on patterm
553 if (n > 10) quality[i] = 10;
554 } else if (quality[i] < 100) {
555 n = clusterspat[scode[i]];
556 // Drop phi hits patterns with high multiplicity
557 if (clusters.count(scode[i]) == 1 || n > 10) {
558 n = 0;
559 // drop phi hits from pattern if already segment hits in same layer
560 quality[i] = 0;
561 phiSelect[i] = 0;
562 phiSelectKeep[i] = 0;
563 continue;
564 }
565 }
566 phiMult[i] = n;
567 double fact = 1.;
568 if (m_idHelperSvc->isRpc(id))
569 fact = 1.2;
570 else if (m_idHelperSvc->isTgc(id))
571 n = 1;
572 else if (m_idHelperSvc->isCsc(id))
573 n = 1;
574
575 error0[i] = phiError[i] * std::sqrt(n) * fact;
576 error[i] = phiError[i] * std::sqrt(n) * fact;
577 double phiHit = std::atan2(phiHity[i], phiHitx[i]);
578 {
579 ATH_MSG_DEBUG(i << " Station " << int(scode[i] / 1000000) << " Hit x " << phiHitx[i] << " Hit y "
580 << phiHity[i] << " Hit z " << phiHitz[i] << " error " << phiError[i] << " phi Hit "
581 << phiHit);
582 ATH_MSG_DEBUG("station " << phiSelect[i]);
583 ATH_MSG_DEBUG("code " << scode[i] << " multiplicity " << n << " error " << error0[i] << " quality "
584 << quality[i]);
585 if (error0[i] < 1.) ATH_MSG_DEBUG("TOO small error ");
586 }
587 }
588 }
589
590 // Count layers hit
591
592 std::map<int, int> layersHit;
593 for (int i = 0; i < nphi; ++i) {
594 if (phiError[i] != 0 && quality[i] > 0) {
595 layersHit[srcode[i]]++;
596 }
597 }
598 int allLayerHits = layersHit.size();
599 int allLayerRecoHits = 0;
600 double pfit = 20000.;
601
602 for (int iqua = 0; iqua < 3; ++iqua) {
603
604 double quacut = 10;
605 if (iqua == 1)
606 quacut = 0;
607 else if (iqua == 2) {
608 quacut = 10;
609 pfit = pmom;
610 }
611
612 ATH_MSG_DEBUG("Quality loop " << iqua << " quality cut " << quacut);
613 int nsel = 0;
614 int nselseg = 0;
615 for (int i = 0; i < nphi; ++i) {
616
617 if (iqua == 1) phiSelect[i] = phiSelectKeep[i];
618
619 if (phiError[i] != 0 && quality[i] > quacut) {
620 nsel++;
621 if (quality[i] > 100) nselseg++;
622 if (quality[i] == 10 && iqua == 1) quality[i] = 11;
623 } else {
624 phiSelect[i] = 0;
625 }
626 if (msgLvl(MSG::DEBUG))
627 ATH_MSG_DEBUG("index i " << i << " phiSelect " << phiSelect[i] << " Quality " << quality[i] << " error "
628 << error[i]);
629 }
630
631 int imax = -1;
632 if (iqua == 1 && nselseg > 0) {
633 // Test and drop pattern Hits if far off
634 double errorScaleFactor = 25.;
635 std::vector<int> phiPatSelect(nphi, 0);
636 for (int i = 0; i < nphi; ++i) {
637 phiPatSelect[i] = 0;
638 if (phiSelect[i] > 0 && quality[i] > 0 && quality[i] < 100) {
639 phiPatSelect[i] = 1;
640 error[i] = errorScaleFactor * error[i];
641 }
642 if (msgLvl(MSG::DEBUG))
643 ATH_MSG_DEBUG("select " << phiSelect[i] << " quality " << quality[i] << " error " << error[i]);
644 }
645 ATH_MSG_DEBUG("performing outlier removal for pattern hits ");
646 fitPhiSL(pfit, phiId, phiHitx, phiHity, phiHitz, error, phiSelect, nphi, phiPull, imax, chi2, r0, phi,
647 errorM, false);
648 for (int i = 0; i < nphi; ++i) {
649 if (phiPatSelect[i] == 1) {
650 error[i] = error[i] / errorScaleFactor;
651 double rescaledPull = phiPull[i] * errorScaleFactor;
652 // 3 sigma cut
653 if (std::abs(rescaledPull) < 3.) {
654 phiSelect[i] = phiSelectKeep[i];
655 } else {
656 phiSelect[i] = 0;
657 phiSelectKeep[i] = 0;
658 if (msgLvl(MSG::DEBUG))
659 ATH_MSG_DEBUG("Drop Pattern Hits with Quality == 1 "
660 << i << " quality " << quality[i] << " Pull " << rescaledPull << " phiSelect "
661 << phiSelect[i]);
662 }
663 }
664 }
665 }
666
667 const double pfitc = pfit;
668 imax = -1;
669
670 if (iqua == 2) {
671 // low momentum fit with scaled error (factor 10) for dropped segment hits
672 std::vector<int> phiReSelect(nphi);
673 for (int i = 0; i < nphi; ++i) {
674 ATH_MSG_DEBUG("select " << phiSelect[i] << " quality " << quality[i]);
675 phiReSelect[i] = 0;
676 if (phiSelect[i] == 0 && quality[i] > 99) {
677 phiReSelect[i] = 1;
678 phiSelect[i] = phiSelectKeep[i];
679 error[i] = 10. * error[i];
680 }
681 }
682 fitPhiSL(pfitc, phiId, phiHitx, phiHity, phiHitz, error, phiSelect, nphi, phiPull, imax, chi2, r0, phi,
683 errorM, false);
684 for (int i = 0; i < nphi; ++i) {
685 if (phiReSelect[i] == 1) {
686 error[i] = error[i] / 10.;
687 // 10 sigma cut (error rescale = 10)
688 if (std::abs(phiPull[i]) < 1) {
689 phiSelect[i] = phiSelectKeep[i];
690 } else {
691 phiSelect[i] = 0;
692 }
693 if (msgLvl(MSG::DEBUG))
694 ATH_MSG_DEBUG("Low momentum Quality == 2 add hit nr "
695 << i << " quality " << quality[i] << " Pull " << phiPull[i] << " phiSelect "
696 << phiSelect[i]);
697 }
698 }
699 }
700 if (iqua == 1 && msgLvl(MSG::DEBUG)) ATH_MSG_DEBUG("Quality loop ");
701 nsel = 0;
702 for (int i = 0; i < nphi; ++i) {
703 errorf[i] = error[i];
704 if (iqua == 1) phiSelect[i] = phiSelectKeep[i];
705
706 if (phiError[i] != 0 && quality[i] > quacut) {
707 nsel++;
708 if (quality[i] == 10 && iqua == 1) quality[i] = 11;
709 } else {
710 phiSelect[i] = 0;
711 }
712 }
713
714 ATH_MSG_DEBUG("Selected PHI hits in fit " << nsel << " iqua " << iqua);
715 if (nsel == 0) continue;
716
717 int niter = -1;
718 // do hit dropping in maximal 10 iterations by putting quality to 0
719
720 for (int iter = 0; iter < 100; ++iter) {
721 niter++;
722 double power = (iter - 10) / 20.;
723 if (power < 0.) power = 0.;
724 chi2 = 0.;
725 nfit = 0;
726 if (iter > 10) {
727 // Shower treatment inflate errors with multiplicity
728 for (int i = 0; i < nphi; ++i) {
729 errorf[i] = error[i] * std::pow(phiMult[i], power);
730 }
731 }
732 fitPhiSL(pfitc, phiId, phiHitx, phiHity, phiHitz, errorf, phiSelect, nphi, phiPull, imax, chi2, r0, phi,
733 errorM, false);
734
735 ncsc = 0;
736 ntgc = 0;
737 nrpc = 0;
738
739 // Count layers hit in Reconstruction
740
741 std::map<int, int> layersRecoHit;
742
743 for (int i = 0; i < nphi; ++i) {
744 Identifier id = phiId[i];
745 if (error[i] == 0 || quality[i] < quacut) phiSelect[i] = 0;
746 if (error[i] != 0 && quality[i] > quacut) {
747 layersRecoHit[srcode[i]]++;
748 {
749 if (m_idHelperSvc->isRpc(id))
750 nrpc++;
751 else if (m_idHelperSvc->isTgc(id))
752 ntgc++;
753 else if (m_idHelperSvc->isCsc(id))
754 ncsc++;
755 }
756 nfit++;
757 }
758 }
759 allLayerRecoHits = layersRecoHit.size();
760 double frac = allLayerRecoHits / (allLayerHits + 0.0001);
761
762 if (nfit == 1) break;
763
764 if (imax < 0 || imax > nphi) {
765 ATH_MSG_DEBUG("Fitphi imax " << imax);
766 break;
767 }
768
769 if (chi2 < 5 * (nfit + 1) || std::abs(phiPull[imax]) < 3.0) {
770
771 ATH_MSG_DEBUG("Final phi " << phi << " frac " << frac << " chi2 " << chi2);
772 break;
773 }
774
775 phiSelect[imax] = 0;
776
777 {
778 ATH_MSG_DEBUG("Start hit dropping " << imax << " pullmax " << phiPull[imax] << " phi " << phi
779 << " chi2 " << chi2);
780 }
781 }
782
783 {
784 ATH_MSG_DEBUG("Fit results phi " << phi << " chi2 " << chi2 << " ndof " << nfit);
785 ATH_MSG_DEBUG("Reco RPC " << nrpc << " TGC " << ntgc << " CSC " << ncsc);
786 }
787
788
789 int nacc = 0;
790 int nshowerdrop = 0;
791 for (int i = 0; i < nphi; ++i) {
792 double power = (niter - 10) / 20.;
793 if (power < 0.) power = 0.;
794 double pull = phiPull[i] * std::pow(phiMult[i], power);
795 if (niter > 10 && std::abs(pull) > 3.0 && phiSelect[i] > 0) {
796 phiSelect[i] = 0;
797 quality[i] = 0;
798 nshowerdrop++;
799 if (msgLvl(MSG::DEBUG))
800 ATH_MSG_DEBUG("Drop shower hit i " << i << " with pull " << pull << " iterations " << niter
801 << " power " << power);
802 }
803 if (phiSelect[i] != 0) nacc++;
804 }
805 if (msgLvl(MSG::DEBUG))
806 ATH_MSG_DEBUG("phi hits " << nphi << " selected for fit " << nfit << " iqua " << iqua << " iterations "
807 << niter << " accepted hits " << nacc << " nshower drop " << nshowerdrop);
808 }
809}
static int layersHit(FPGATrackSimRoad &r)
void fitPhiSL(const double pmom, const std::vector< Identifier > &id, const std::vector< double > &hitx, const std::vector< double > &hity, const std::vector< double > &hitz, const std::vector< double > &error, std::vector< int > &select, const int n, std::vector< double > &pull, int &imax, double &chi2, double &r0, double &phi, std::vector< double > &errorM, bool fast) const
fit method straight line model

◆ initialize()

StatusCode MuonPhiHitSelector::initialize ( )
override

Definition at line 38 of file MuonPhiHitSelector.cxx.

39{
40 ATH_MSG_VERBOSE("MuonPhiHitSelector::Initializing");
42 ATH_CHECK(m_clusterCreator.retrieve());
43 ATH_CHECK(m_idHelperSvc.retrieve());
44 ATH_CHECK(m_cscRotCreator.retrieve(EnableTool{!m_cscRotCreator.empty()}));
45 ATH_MSG_VERBOSE("End of Initializing");
46 return StatusCode::SUCCESS;
47}
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_VERBOSE(x)
ToolHandle< Muon::IMuonCompetingClustersOnTrackCreator > m_competingRIOsOnTrackTool
Toolhandle to CompetingRIOsOnTrackTool creator.
ToolHandle< Muon::IMuonClusterOnTrackCreator > m_cscRotCreator
ToolHandle< Muon::IMuonClusterOnTrackCreator > m_clusterCreator
Toolhandle to ClusterOnTrackTool creator.

◆ inputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::inputHandles ( ) const
overridevirtualinherited

Return this algorithm's input handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ interfaceID()

const InterfaceID & Muon::IMuonHitSelector::interfaceID ( )
inlinestaticinherited

Definition at line 25 of file IMuonHitSelector.h.

25 {
26 static const InterfaceID IID_IMuonHitSelector("Muon::IMuonHitSelector", 1, 0);
27 return IID_IMuonHitSelector;
28 }

◆ msg()

MsgStream & AthCommonMsg< AlgTool >::msg ( ) const
inlineinherited

Definition at line 24 of file AthCommonMsg.h.

24 {
25 return this->msgStream();
26 }

◆ msgLvl()

bool AthCommonMsg< AlgTool >::msgLvl ( const MSG::Level lvl) const
inlineinherited

Definition at line 30 of file AthCommonMsg.h.

30 {
31 return this->msgLevel(lvl);
32 }

◆ outputHandles()

virtual std::vector< Gaudi::DataHandle * > AthCommonDataStore< AthCommonMsg< AlgTool > >::outputHandles ( ) const
overridevirtualinherited

Return this algorithm's output handles.

We override this to include handle instances from key arrays if they have not yet been declared. See comments on updateVHKA.

◆ renounce()

std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > AthCommonDataStore< AthCommonMsg< AlgTool > >::renounce ( T & h)
inlineprotectedinherited

Definition at line 380 of file AthCommonDataStore.h.

381 {
382 h.renounce();
384 }
std::enable_if_t< std::is_void_v< std::result_of_t< decltype(&T::renounce)(T)> > &&!std::is_base_of_v< SG::VarHandleKeyArray, T > &&std::is_base_of_v< Gaudi::DataHandle, T >, void > renounce(T &h)

◆ renounceArray()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::renounceArray ( SG::VarHandleKeyArray & handlesArray)
inlineprotectedinherited

remove all handles from I/O resolution

Definition at line 364 of file AthCommonDataStore.h.

364 {
366 }

◆ select_rio()

std::vector< std::unique_ptr< const Trk::MeasurementBase > > MuonPhiHitSelector::select_rio ( const double pmom,
const std::vector< const Trk::RIO_OnTrack * > & associatedHits,
const std::vector< const Trk::PrepRawData * > & unassociatedHits ) const
overridevirtual

Selects and builds a cleaned vector of RIO fits the associatedHits and build new RIOs, if m_competingRios true then for ambiguous hits competing rios are built.

Implements Muon::IMuonHitSelector.

Definition at line 50 of file MuonPhiHitSelector.cxx.

52{
53 // Make Rios on Track
54
55 int time_start = std::clock() / 1000;
56
57 std::vector<std::unique_ptr<const Trk::MeasurementBase>> selectedHits{}, selectedClusters{};
58
59 ATH_MSG_VERBOSE("Executing MuonPhiHitSelectorTool select_rio ");
60
61 int nhits = associatedHits.size() + unassociatedHits.size();
62
63 ATH_MSG_DEBUG("Executing MuonPhiHitSelectorTool nhits select_rio " << nhits);
64
65 std::vector<double> phiHitx(nhits);
66 std::vector<double> phiHity(nhits);
67 std::vector<double> phiHitz(nhits);
68 std::vector<double> phiError(nhits);
69 std::vector<Identifier> phiId(nhits);
70 std::vector<double> phiPull(nhits);
71 std::vector<int> phiSelect(nhits);
72 std::vector<int> phiMult(nhits);
73 std::vector<int> quality(nhits);
74 std::vector<const Trk::PrepRawData*> phiPrep(nhits);
75
76 std::map<Identifier, int> phiMapId;
77 int nphi = 0;
78
79 for (const Trk::RIO_OnTrack* rot : associatedHits) {
80 Identifier id = rot->identify();
81 phiId[nphi] = id;
82 Amg::Vector3D gHitPos = rot->globalPosition();
83 if (m_idHelperSvc->isRpc(id)) {
84 phiSelect[nphi] = 1;
85 } else if (m_idHelperSvc->isTgc(id)) {
86 phiSelect[nphi] = 2;
87 } else if (m_idHelperSvc->isCsc(id)) {
88 phiSelect[nphi] = 3;
89 }
90 phiHitx[nphi] = gHitPos.x();
91 phiHity[nphi] = gHitPos.y();
92 phiHitz[nphi] = gHitPos.z();
93
94 const Amg::MatrixX& cov = rot->localCovariance();
95 double error = cov(0, 0);
96
97 // for the TGCs diagonalize the error matrix and use the smallest component
98 if (cov.cols() > 1) {
99 AmgSymMatrix(2) Er{AmgSymMatrix(2)::Zero()};
100 Er(0, 0) = cov(0, 0);
101 Er(0, 1) = cov(0, 1);
102 Er(1, 1) = cov(1, 1);
103 Er(1, 0) = Er(0, 1);
104
105 double chi = Er(0, 0) != Er(1, 1) ? std::atan(-2 * Er(0, 1) / (Er(0, 0) - Er(1, 1))) / 2. : 0.;
106
107 CxxUtils::sincos scchi(chi);
108
109 AmgSymMatrix(2) Rot{AmgSymMatrix(2)::Zero()};
110 Rot(0, 0) = scchi.cs;
111 Rot(1, 1) = Rot(0, 0);
112 Rot(0, 1) = scchi.sn;
113 Rot(1, 0) = -Rot(0, 1);
114 AmgMatrix(2, 2) D = Rot.transpose() * Er * Rot;
115 ATH_MSG_DEBUG("Diagonalized error matrix " << D);
116 error = std::min(D(0, 0),D(1, 1));
117 }
118 phiError[nphi] = std::sqrt(error);
119 quality[nphi] = 1000;
120 phiMapId[id] = 1;
121 phiPrep[nphi] = rot->prepRawData();
122 double phipos = std::atan2(phiHity[nphi], phiHitx[nphi]);
123 ATH_MSG_DEBUG("phi Segment Hit " << nphi << " det " << phiSelect[nphi] << " phi " << phipos);
124 nphi++;
125 }
126 int nphiseg = nphi;
127
128 for (const Trk::PrepRawData* prd : unassociatedHits) {
129 Identifier id = prd->identify();
130 phiId[nphi] = id;
131 // Skip phi hits already on segment
132 if (phiMapId.count(id)) continue;
133 const Muon::MuonCluster* clus = dynamic_cast<const Muon::MuonCluster*>(prd);
134 if (!clus) continue;
135 if (m_idHelperSvc->isRpc(id))
136 phiSelect[nphi] = 1;
137 else if (m_idHelperSvc->isTgc(id))
138 phiSelect[nphi] = 2;
139 else if (m_idHelperSvc->isCsc(id))
140 phiSelect[nphi] = 3;
141 Amg::Vector3D gHitPos = clus->globalPosition();
142 phiHitx[nphi] = gHitPos.x();
143 phiHity[nphi] = gHitPos.y();
144 phiHitz[nphi] = gHitPos.z();
145 phiError[nphi] = prd->localCovariance()(Trk::locX);
146 quality[nphi] = 10;
147 phiPrep[nphi] = prd;
148 double phipos = std::atan2(phiHity[nphi], phiHitx[nphi]);
149 ATH_MSG_DEBUG("phi Pattern Hit " << nphi << " phi " << phipos);
150 nphi++;
151 }
152
153 double chi2(0);
154 double r0(0);
155 int nfit;
156 std::vector<double> errorM(4);
157 double phi(DBL_MAX);
158 fitRecPhi(pmom, phiId, phiHitx, phiHity, phiHitz, phiError, quality, nphi, phiPull, phiMult, phiSelect, chi2, r0,
159 phi, errorM, nfit);
160
161 // Define global track parameters (not used 27-8 JS)
162
163 for (int i = 0; i < nphi; ++i) {
164 if (phiSelect[i] > 0) {
165 if (phiSelect[i] == 1) {
166 const Muon::RpcPrepData* prd = dynamic_cast<const Muon::RpcPrepData*>(phiPrep[i]);
167 const Amg::Vector3D globalpos(phiHitx[i], phiHity[i], phiHitz[i]);
168 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
169 if (rio) selectedHits.push_back(std::move(rio));
170 } else if (phiSelect[i] == 2) {
171 const Muon::TgcPrepData* prd = dynamic_cast<const Muon::TgcPrepData*>(phiPrep[i]);
172 const Amg::Vector3D globalpos(phiHitx[i], phiHity[i], phiHitz[i]);
173 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
174 if (rio) selectedHits.push_back(std::move(rio));
175 } else if (phiSelect[i] == 3) {
176 const Muon::CscPrepData* prd = dynamic_cast<const Muon::CscPrepData*>(phiPrep[i]);
177 const Amg::Vector3D globalpos(phiHitx[i], phiHity[i], phiHitz[i]);
178 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_cscRotCreator->createRIO_OnTrack(*prd, globalpos)};
179 if (rio) selectedHits.push_back(std::move(rio));
180 }
181
182 ATH_MSG_DEBUG("Make ONE rio per PrepData");
183 }
184 }
185 ATH_MSG_DEBUG("Fit hit results phi " << phi << " chi2 " << chi2 << " segment hits " << nphiseg
186 << " pattern hits " << nphi - nphiseg << " nfit " << nfit << " rio size "
187 << selectedHits.size());
188
189
190 std::vector<double> clusterX(nphi);
191 std::vector<double> clusterY(nphi);
192 std::vector<double> clusterZ(nphi);
193 std::vector<double> clusterError(nphi);
194 std::vector<Identifier> clusterId(nphi);
195 std::vector<int> clusterHits(nphi);
196 std::vector<double> clusterPull(nphi);
197 std::vector<int> clusterSelect(nphi);
198 // Link from hit to cluster
199 std::vector<int> clusterInt(nphi);
200
201 int ncl, imax;
202 double chi2cl, r0cl, phicl;
203 std::vector<double> errorMcl(4);
204 clusterPhi(phiId, phiHitx, phiHity, phiHitz, phiError, phiPull, phiSelect, nphi, clusterX, clusterY, clusterZ,
205 clusterError, clusterId, clusterHits, clusterSelect, clusterInt, ncl);
206
207
208 for (int ic = 0; ic < ncl; ++ic) {
209 std::list<const Trk::PrepRawData*> prdList;
210 int iic = -1;
211 double avError = 0.;
212 int ip = -1;
213 int np = 0;
214 for (int i = 0; i < nphi; ++i) {
215 if (clusterInt[i] == ic) {
216 ip = i;
217 prdList.push_back(phiPrep[i]);
218 avError += 1. / (phiError[i] * phiError[i]);
219 if (clusterId[ic] == phiId[i]) iic = i;
220 np++;
221 }
222 }
223 if (iic > -1) {
224 ATH_MSG_DEBUG("Phi cluster found np " << np << " ip " << ip);
225 if (np == 1) {
226 // Only one PrepData: create RIO on Track
227 const Amg::Vector3D globalpos(clusterX[ic], clusterY[ic], clusterZ[ic]);
228 if (phiSelect[ip] == 1) {
229 ATH_MSG_DEBUG("Phi RPC rio");
230 const Muon::RpcPrepData* prd = dynamic_cast<const Muon::RpcPrepData*>(phiPrep[ip]);
231 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
232 if (rio) selectedClusters.push_back(std::move(rio));
233 } else if (phiSelect[ip] == 2) {
234 ATH_MSG_DEBUG("Phi TGC rio");
235 const Muon::TgcPrepData* prd = dynamic_cast<const Muon::TgcPrepData*>(phiPrep[ip]);
236 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
237 if (rio) selectedClusters.push_back(std::move(rio));
238 } else if (phiSelect[ip] == 3) {
239 ATH_MSG_DEBUG("Phi CSC rio");
240 const Muon::CscPrepData* prd = dynamic_cast<const Muon::CscPrepData*>(phiPrep[ip]);
241 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_cscRotCreator->createRIO_OnTrack(*prd, globalpos)};
242 if (rio) selectedClusters.push_back(std::move(rio));
243 }
244 } else {
245
246 if (m_competingRios) {
247 // More PrepData's: create Competing RIOs on Track
248 avError = std::sqrt(1. / avError);
249 double scaleFactor = clusterError[ic] / avError;
250 std::unique_ptr<const Trk::CompetingRIOsOnTrack> rio =
251 m_competingRIOsOnTrackTool->createBroadCluster(prdList, scaleFactor);
252 selectedClusters.push_back(std::move(rio));
253 if (msgLvl(MSG::DEBUG))
254 ATH_MSG_DEBUG("Make competing rio/cluster "
255 << " scale factor " << scaleFactor << " number of rios " << prdList.size());
256 } else {
257 // Make one Rio for central cluster
258 ip = iic;
259 const Amg::Vector3D globalpos(clusterX[ic], clusterY[ic], clusterZ[ic]);
260 if (phiSelect[ip] == 1) {
261 ATH_MSG_DEBUG("Phi RPC rio central cluster");
262 const Muon::RpcPrepData* prd = dynamic_cast<const Muon::RpcPrepData*>(phiPrep[ip]);
263 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
264 if (rio) selectedClusters.push_back(std::move(rio));
265 } else if (phiSelect[ip] == 2) {
266 ATH_MSG_DEBUG("Phi TGC rio central cluster");
267 const Muon::TgcPrepData* prd = dynamic_cast<const Muon::TgcPrepData*>(phiPrep[ip]);
268 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{ m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
269 if (rio) selectedClusters.push_back(std::move(rio));
270 } else if (phiSelect[ip] == 3) {
271 ATH_MSG_DEBUG("Phi CSC rio central cluster");
272 const Muon::CscPrepData* prd = dynamic_cast<const Muon::CscPrepData*>(phiPrep[ip]);
273 std::unique_ptr<const Muon::MuonClusterOnTrack> rio{m_clusterCreator->createRIO_OnTrack(*prd, globalpos)};
274 if (rio) selectedClusters.push_back(std::move(rio));
275 }
276 }
277 }
278 } else {
279 ATH_MSG_DEBUG("Phi cluster NOT found ");
280 }
281 }
282
283 fitPhiSL(pmom, clusterId, clusterX, clusterY, clusterZ, clusterError, clusterSelect, ncl, clusterPull, imax, chi2cl,
284 r0cl, phicl, errorMcl, false);
285
286
287 if (msgLvl(MSG::DEBUG) || m_summary) {
288 ATH_MSG_DEBUG("PhiHitSelector Time spent " << std::clock() / 1000 - time_start << " nhits " << nhits
289 << " segment hits " << associatedHits.size() << " nfit " << nfit
290 << " nclusters " << ncl);
291 ATH_MSG_DEBUG("Fit cluster results phi " << phicl << " chi2 " << chi2cl << " number of clusters " << ncl
292 << " size cluster Hits " << selectedClusters.size());
293 }
294 if (m_makeClusters) {
295 return selectedClusters;
296 }
297 return selectedHits;
298}
#define AmgSymMatrix(dim)
#define AmgMatrix(rows, cols)
if(febId1==febId2)
#define min(a, b)
Definition cfImp.cxx:40
Gaudi::Property< bool > m_makeClusters
flag that performs a clusterization and return clusters (default: false)
Gaudi::Property< bool > m_summary
flag to print out a summary of what comes in and what comes out
Gaudi::Property< bool > m_competingRios
flag that build competing rios on track for amibguous trigger hits (default: false)
void clusterPhi(const std::vector< Identifier > &id, const std::vector< double > &hitx, const std::vector< double > &hity, const std::vector< double > &hitz, const std::vector< double > &error, const std::vector< double > &pull, std::vector< int > &select, const int n, std::vector< double > &clusterX, std::vector< double > &clusterY, std::vector< double > &clusterZ, std::vector< double > &clusterError, std::vector< Identifier > &clusterId, std::vector< int > &clusterHits, std::vector< int > &clusterSelect, std::vector< int > &clusterInt, int &ncl) const
clusterization method
void fitRecPhi(const double pmom, const std::vector< Identifier > &phiId, const std::vector< double > &phiHitx, const std::vector< double > &phiHity, const std::vector< double > &phiHitz, const std::vector< double > &phiError, std::vector< int > &quality, const int nphi, std::vector< double > &phiPull, std::vector< int > &phiMult, std::vector< int > &phiSelect, double &chi2, double &r0, double &phi, std::vector< double > &errorM, int &nfit) const
fit method curved track model
Eigen::Matrix< double, 3, 1 > Vector3D
@ locX
Definition ParamDefs.h:37
void Zero(TH1D *hin)
Definition generate.cxx:32

◆ sysInitialize()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysInitialize ( )
overridevirtualinherited

Perform system initialization for an algorithm.

We override this to declare all the elements of handle key arrays at the end of initialization. See comments on updateVHKA.

Reimplemented in asg::AsgMetadataTool, AthCheckedComponent< AthAlgTool >, AthCheckedComponent<::AthAlgTool >, and DerivationFramework::CfAthAlgTool.

◆ sysStart()

virtual StatusCode AthCommonDataStore< AthCommonMsg< AlgTool > >::sysStart ( )
overridevirtualinherited

Handle START transition.

We override this in order to make sure that conditions handle keys can cache a pointer to the conditions container.

◆ updateVHKA()

void AthCommonDataStore< AthCommonMsg< AlgTool > >::updateVHKA ( Gaudi::Details::PropertyBase & )
inlineinherited

Definition at line 308 of file AthCommonDataStore.h.

308 {
309 // debug() << "updateVHKA for property " << p.name() << " " << p.toString()
310 // << " size: " << m_vhka.size() << endmsg;
311 for (auto &a : m_vhka) {
313 for (auto k : keys) {
314 k->setOwner(this);
315 }
316 }
317 }
std::vector< SG::VarHandleKeyArray * > m_vhka

Member Data Documentation

◆ m_clusterCreator

ToolHandle<Muon::IMuonClusterOnTrackCreator> MuonPhiHitSelector::m_clusterCreator {this, "MuonClusterOnTrackCreator",""}
private

Toolhandle to ClusterOnTrackTool creator.

Definition at line 50 of file MuonPhiHitSelector.h.

50{this, "MuonClusterOnTrackCreator",""};

◆ m_competingRios

Gaudi::Property<bool> MuonPhiHitSelector::m_competingRios {this, "CompetingRios", false}
private

flag that build competing rios on track for amibguous trigger hits (default: false)

Definition at line 61 of file MuonPhiHitSelector.h.

61{this, "CompetingRios", false};

◆ m_competingRIOsOnTrackTool

ToolHandle<Muon::IMuonCompetingClustersOnTrackCreator> MuonPhiHitSelector::m_competingRIOsOnTrackTool
private
Initial value:
{
this,
"MuonCompetingClustersOnTrackCreator",
"Muon::MuonCompetingClustersOnTrackCreator/MuonCompetingClustersOnTrackCreator",
}

Toolhandle to CompetingRIOsOnTrackTool creator.

Definition at line 44 of file MuonPhiHitSelector.h.

44 {
45 this,
46 "MuonCompetingClustersOnTrackCreator",
47 "Muon::MuonCompetingClustersOnTrackCreator/MuonCompetingClustersOnTrackCreator",
48 };

◆ m_cosmics

Gaudi::Property<bool> MuonPhiHitSelector::m_cosmics {this, "DoCosmics", false}
private

flag for use of cosmics, straight line model will be used, no interaction point constraint

Definition at line 57 of file MuonPhiHitSelector.h.

57{this, "DoCosmics", false};

◆ m_cscRotCreator

ToolHandle<Muon::IMuonClusterOnTrackCreator> MuonPhiHitSelector::m_cscRotCreator {this, "CscRotCreator",""}
private

Definition at line 52 of file MuonPhiHitSelector.h.

52{this, "CscRotCreator",""};

◆ m_detStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_detStore
privateinherited

Pointer to StoreGate (detector store by default)

Definition at line 393 of file AthCommonDataStore.h.

◆ m_evtStore

StoreGateSvc_t AthCommonDataStore< AthCommonMsg< AlgTool > >::m_evtStore
privateinherited

Pointer to StoreGate (event store by default)

Definition at line 390 of file AthCommonDataStore.h.

◆ m_idHelperSvc

ServiceHandle<Muon::IMuonIdHelperSvc> MuonPhiHitSelector::m_idHelperSvc
private
Initial value:
{
this,
"MuonIdHelperSvc",
"Muon::MuonIdHelperSvc/MuonIdHelperSvc",
}

Definition at line 37 of file MuonPhiHitSelector.h.

37 {
38 this,
39 "MuonIdHelperSvc",
40 "Muon::MuonIdHelperSvc/MuonIdHelperSvc",
41 };

◆ m_makeClusters

Gaudi::Property<bool> MuonPhiHitSelector::m_makeClusters {this, "MakeClusters", false}
private

flag that performs a clusterization and return clusters (default: false)

Definition at line 59 of file MuonPhiHitSelector.h.

59{this, "MakeClusters", false};

◆ m_summary

Gaudi::Property<bool> MuonPhiHitSelector::m_summary {this, "DoSummary", false}
private

flag to print out a summary of what comes in and what comes out

Definition at line 55 of file MuonPhiHitSelector.h.

55{this, "DoSummary", false};

◆ m_varHandleArraysDeclared

bool AthCommonDataStore< AthCommonMsg< AlgTool > >::m_varHandleArraysDeclared
privateinherited

Definition at line 399 of file AthCommonDataStore.h.

◆ m_vhka

std::vector<SG::VarHandleKeyArray*> AthCommonDataStore< AthCommonMsg< AlgTool > >::m_vhka
privateinherited

Definition at line 398 of file AthCommonDataStore.h.


The documentation for this class was generated from the following files: