ATLAS Offline Software
ConversionPostSelector.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 /***************************************************************************
6  ConversionPostSelector.cxx - Description
7  -------------------
8  begin : 01-01-2008
9  authors : Tatjana Lenz, Thomas Koffas
10  email : tatjana.lenz@cern.ch, Thomas.Koffas@cern.ch
11  changes :
12 ***************************************************************************/
13 #include "InDetConversionFinderTools/ConversionPostSelector.h"
14 #include "xAODTracking/TrackParticle.h"
15 #include "xAODTracking/Vertex.h"
16 #include "CLHEP/Vector/LorentzVector.h"
17 #include "TrkEventPrimitives/ParticleHypothesis.h" //ParticleMasses struct
18 
19 #include <cmath>
20 
21 namespace{
22  constexpr double twopi{2.*M_PI};
23 }
24 
25 namespace InDet {
26 
27 
28  static const InterfaceID IID_IConversionPostSelector("InDet::ConversionPostSelector", 1, 0);
29 
30  ConversionPostSelector::ConversionPostSelector(const std::string& type,
31  const std::string& name,
32  const IInterface* parent)
33  : AthAlgTool(type, name, parent),
34  m_minPt{0},
35  m_maxdR{-10000.},
36  m_maxPhiVtxTrk{0.2},
37  m_decorateVertices{true},
38  m_massK0{497.672},
39  m_sigmaK0{8.5},
40  m_massLambda{1115.683},
41  m_sigmaLambda{3.5},
42  m_nsig{5} {
43 
44  m_maxChi2.push_back(35.);
45  m_maxChi2.push_back(25.);
46  m_maxChi2.push_back(20.);
47 
48  m_invMassCut.push_back(10000.);
49  m_invMassCut.push_back(10000.);
50  m_invMassCut.push_back(10000.);
51 
52  m_fitMomentum.push_back(0.);
53  m_fitMomentum.push_back(0.);
54  m_fitMomentum.push_back(0.);
55 
56  m_minRadius.push_back(-10000.);
57  m_minRadius.push_back(-10000.);
58  m_minRadius.push_back(-10000.);
59 
60  declareInterface<ConversionPostSelector>(this);
61  declareProperty("MaxChi2Vtx", m_maxChi2);
62  declareProperty("MaxInvariantMass", m_invMassCut);
63  declareProperty("MinFitMomentum", m_fitMomentum);
64  declareProperty("MinRadius", m_minRadius);
65  declareProperty("MinPt", m_minPt);
66  declareProperty("MaxdR", m_maxdR);
67  declareProperty("MaxPhiVtxTrk", m_maxPhiVtxTrk);
68  declareProperty("NSigma", m_nsig);
69  declareProperty("DecorateVertices", m_decorateVertices);
70  }
71 
72  ConversionPostSelector::~ConversionPostSelector() = default;
73 
74  const InterfaceID& ConversionPostSelector::interfaceID() {
75  return IID_IConversionPostSelector;
76  }
77 
78  StatusCode ConversionPostSelector::initialize() {
79  return StatusCode::SUCCESS;
80  }
81 
82  StatusCode ConversionPostSelector::finalize() {
83  return StatusCode::SUCCESS;
84  }
85 
86  bool ConversionPostSelector::selectConversionCandidate(xAOD::Vertex * vertex, int flag,
87  std::vector<Amg::Vector3D>& trkL) const{
88  bool pass = true;
89 
90  //Determine the cuts
91  double maxChi2 = 1000.;
92  double invMassCut = 1000.;
93  double fitMomentum = 0.;
94  double radius = 1000.;
95 
96  if(flag==0) { // Si+Si
97  maxChi2 = m_maxChi2[0] ; //Vertex fit chi2 cut
98  invMassCut = m_invMassCut[0] ; //Fit invariant mass cut
99  fitMomentum = m_fitMomentum[0]; //Photon fitted momentum
100  radius = m_minRadius[0] ; //Minimum acceptable radius of conversion vertex
101  }
102  if(flag==1) { // Si+TRT
103  maxChi2 = m_maxChi2[1] ; //Vertex fit chi2 cut
104  invMassCut = m_invMassCut[1] ; //Fit invariant mass cut
105  fitMomentum = m_fitMomentum[1]; //Photon fitted momentum
106  radius = m_minRadius[1] ; //Minimum acceptable radius of conversion vertex
107  }
108  if(flag==2) { // TRT+TRT
109  maxChi2 = m_maxChi2[2] ; //Vertex fit chi2 cut
110  invMassCut = m_invMassCut[2] ; //Fit invariant mass cut
111  fitMomentum = m_fitMomentum[2]; //Photon fitted momentum
112  radius = m_minRadius[2] ; //Minimum acceptable radius of conversion vertex
113  }
114 
115  //chi2 cut
116  if (vertex->nTrackParticles() != 2) {
117  ATH_MSG_DEBUG("Incorrect number of tracks used in conversion fit.");
118  pass = false;
119  } else {
120 
121  float reducedChi2 = vertex->chiSquared()/vertex->numberDoF();
122  if (reducedChi2 > maxChi2) pass = false;
123  if (reducedChi2 > maxChi2) pass = false;
124 
125  //Minimum radius cut
126  double vtxR = vertex->position().perp();
127  if(vtxR < radius) pass = false;
128 
129  // Parameters at vertex
130  std::vector< Trk::VxTrackAtVertex >& trkAtVx = vertex->vxTrackAtVertex();
131  if (trkAtVx.size() != 2 ||
132  !trkAtVx[0].perigeeAtVertex() ||
133  !trkAtVx[1].perigeeAtVertex())
134  {
135  ATH_MSG_DEBUG("VxTrackAtVertex or perigeeAtVertex not available");
136  return false;
137  }
138  const Trk::TrackParameters& perigee1 = *(trkAtVx[0].perigeeAtVertex());
139  const Trk::TrackParameters& perigee2 = *(trkAtVx[1].perigeeAtVertex());
140 
141  //invariant mass
142  CLHEP::HepLorentzVector momentum;
143  Amg::Vector3D sum_mom = perigee1.momentum() + perigee2.momentum();
144  double m2 = std::pow(Trk::ParticleMasses::mass[Trk::electron],2);
145  double ee = std::sqrt(m2 + perigee1.momentum().mag2()) + std::sqrt(m2 + perigee2.momentum().mag2());
146  momentum.setPx(sum_mom.x()); momentum.setPy(sum_mom.y()); momentum.setPz(sum_mom.z()); momentum.setE(ee);
147  double inv_mass = momentum.m();
148  double photonP = std::sqrt(momentum.x()*momentum.x() + momentum.y()*momentum.y());
149  double pt1 = perigee1.pT(); double pt2 = perigee2.pT();
150 
151  if (pt1<m_minPt || pt2<m_minPt) pass = false;
152  if (std::fabs(inv_mass) > invMassCut) pass = false;
153  if (photonP < fitMomentum) pass = false;
154 
155  double fR = 1000.;
156  std::vector<Amg::Vector3D>::const_iterator ipb=trkL.begin();
157  std::vector<Amg::Vector3D>::const_iterator ipbe=trkL.end();
158  for(; ipb!=ipbe;++ipb){
159  double tmpfR = (*ipb).perp();
160  if(tmpfR<fR) fR = tmpfR;
161  }
162  if(flag==1 && fR-vtxR<m_maxdR) pass = false;
163 
164  double PhiVtxTrk = vertex->position().phi() - perigee1.parameters()[Trk::phi0];
165  if (PhiVtxTrk < -M_PI) PhiVtxTrk += twopi;
166  if (PhiVtxTrk > M_PI) PhiVtxTrk -= twopi;
167  if (std::fabs(PhiVtxTrk)>m_maxPhiVtxTrk) pass = false;
168 
169  if (pass && m_decorateVertices)
170  {
171  ATH_MSG_DEBUG("Decorating vertex with values used in post selector");
172  decorateVertex(*vertex, inv_mass, pt1, pt2, fR, std::fabs(PhiVtxTrk) );
173  }
174 
175  }
176  return pass;
177  }
178 
179  bool
180  ConversionPostSelector::selectSecVtxCandidate(
181  xAOD::Vertex* vertex,
182  int flag,
183  std::vector<Amg::Vector3D>& trkL,
184  int& type) const
185  {
186  bool pass = true;
187  bool isK0 = false;
188  bool isLambda = false;
189  bool isLambdaBar = false;
190  int kind = -1;
191 
192  //Determine the cuts
193  double maxChi2 = 1000.;
194  double fitMomentum = 0.;
195  double radius = 1000.;
196 
197  if(flag==0) {
198  maxChi2 = m_maxChi2[0] ; //Vertex fit chi2 cut
199  fitMomentum = m_fitMomentum[0]; //Photon fitted momentum
200  radius = m_minRadius[0] ; //Minimum acceptable radius of conversion vertex
201  }
202  if(flag==1) {
203  maxChi2 = m_maxChi2[1] ; //Vertex fit chi2 cut
204  fitMomentum = m_fitMomentum[1]; //Photon fitted momentum
205  radius = m_minRadius[1] ; //Minimum acceptable radius of conversion vertex
206  }
207  if(flag==2) {
208  maxChi2 = m_maxChi2[2] ; //Vertex fit chi2 cut
209  fitMomentum = m_fitMomentum[2]; //Photon fitted momentum
210  radius = m_minRadius[2] ; //Minimum acceptable radius of conversion vertex
211  }
212 
213  //chi2 cut
214 
215  if (vertex->nTrackParticles() != 2) {
216  ATH_MSG_DEBUG("Incorrect number of tracks used in conversion fit.");
217  pass = false;
218  } else {
219  float reducedChi2 = vertex->chiSquared()/vertex->numberDoF();
220  if (reducedChi2 > maxChi2) pass = false;
221 
222  //Minimum radius cut
223  double vtxR = vertex->position().perp();
224  if(vtxR < radius) pass = false;
225 
226  // Parameters at vertex
227  std::vector< Trk::VxTrackAtVertex >& trkAtVx = vertex->vxTrackAtVertex();
228  if (trkAtVx.size() != 2 ||
229  !trkAtVx[0].perigeeAtVertex() ||
230  !trkAtVx[1].perigeeAtVertex())
231  {
232  ATH_MSG_DEBUG("VxTrackAtVertex or perigeeAtVertex not available");
233  return false;
234  }
235  const Trk::TrackParameters& perigee1 = *(trkAtVx[0].perigeeAtVertex());
236  const Trk::TrackParameters& perigee2 = *(trkAtVx[1].perigeeAtVertex());
237 
238  double pt1 = perigee1.pT(); double pt2 = perigee2.pT();
239  if (pt1<m_minPt || pt2<m_minPt) pass = false;
240 
241  double fR = 1000.;
242  std::vector<Amg::Vector3D>::const_iterator ipb=trkL.begin();
243  std::vector<Amg::Vector3D>::const_iterator ipbe=trkL.end();
244  for(; ipb!=ipbe;++ipb){
245  double tmpfR = (*ipb).perp();
246  if(tmpfR<fR) fR = tmpfR;
247  }
248  if(flag==1 && fR-vtxR<m_maxdR) pass = false;
249 
250  //invariant mass. First assume K0, if failed assume Lambda
251  CLHEP::HepLorentzVector momentumK0 = fourP(perigee1,perigee2,m_massK0,false);
252  double inv_massK0 = momentumK0.m();
253  if (std::fabs(inv_massK0-m_massK0) <= m_nsig*m_sigmaK0) isK0 = true;
254  CLHEP::HepLorentzVector momentumL = fourP(perigee1,perigee2,m_massLambda,false);
255  double inv_massL = momentumL.m();
256  if (std::fabs(inv_massL-m_massLambda) <= m_nsig*m_sigmaLambda) isLambda = true;
257  CLHEP::HepLorentzVector momentumLb = fourP(perigee1,perigee2,m_massLambda,true);
258  double inv_massLb = momentumLb.m();
259  if (std::fabs(inv_massLb-m_massLambda) <= m_nsig*m_sigmaLambda) isLambdaBar = true;
260  if (!isLambdaBar && !isLambda && !isK0) pass = false;
261  CLHEP::HepLorentzVector momentum;
262  if(isK0 && isLambda && !isLambdaBar) {momentum = momentumK0; kind = 110;}
263  if(isK0 && isLambdaBar && !isLambda) {momentum = momentumK0; kind = 101;}
264  if(isK0 && !isLambda && !isLambdaBar) {momentum = momentumK0; kind = 100;}
265  if(!isK0 && isLambda && !isLambdaBar) {momentum = momentumL; kind = 10;}
266  if(!isK0 && isLambdaBar && !isLambda) {momentum = momentumLb; kind = 1;}
267  if(!isK0 && isLambda && isLambdaBar) {momentum = momentumL; kind = 11;}
268  double particleP = std::sqrt(momentum.x()*momentum.x() + momentum.y()*momentum.y());
269  if (particleP < fitMomentum) pass = false;
270  }
271  type = kind;
272  return pass;
273  }
274 
275  CLHEP::HepLorentzVector
276  ConversionPostSelector::fourP(const Trk::TrackParameters& per1,
277  const Trk::TrackParameters& per2,
278  double mass,
279  bool isBar) const
280  {
281  CLHEP::HepLorentzVector momentum;
282  Amg::Vector3D sum_mom = per1.momentum() + per2.momentum();
283  double mp1 = 0.; double mp2 = 0.;
284  if(mass==m_massK0) {
285  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
286  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
287  }else{
288  if(!isBar){
289  if(per1.charge()>0) {
290  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
291  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
292  } else {
293  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
294  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
295  }
296  }else{
297  if(per1.charge()>0) {
298  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
299  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
300  } else {
301  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
302  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
303  }
304  }
305  }
306  double ee = std::sqrt(mp1 + per1.momentum().mag2()) + std::sqrt(mp2 + per2.momentum().mag2());
307  momentum.setPx(sum_mom.x()); momentum.setPy(sum_mom.y()); momentum.setPz(sum_mom.z()); momentum.setE(ee);
308  return momentum;
309  }
310 
311  void
312  ConversionPostSelector::decorateVertex(xAOD::Vertex& vertex,
313  float inv_mass,
314  float pt1,
315  float pt2,
316  float fR,
317  float deltaPhiVtxTrk)
318  {
319  static const SG::AuxElement::Accessor<float> accMass("mass");
320  accMass(vertex) = inv_mass;
321  static const SG::AuxElement::Accessor<float> accPt1("pt1");
322  accPt1(vertex) = pt1;
323  static const SG::AuxElement::Accessor<float> accPt2("pt2");
324  accPt2(vertex) = pt2;
325  static const SG::AuxElement::Accessor<float> accMinRfirstHit("minRfirstHit");
326  accMinRfirstHit(vertex) = fR;
327  static const SG::AuxElement::Accessor<float> accDeltaPhiVtxTrk("deltaPhiVtxTrk");
328  accDeltaPhiVtxTrk(vertex) = deltaPhiVtxTrk;
329  }
330 
331 } // namespace InDet
332 
Trk::proton
@ proton
Definition: ParticleHypothesis.h:31
python.SystemOfUnits.m2
int m2
Definition: SystemOfUnits.py:92
InDet::ConversionPostSelector::fourP
CLHEP::HepLorentzVector fourP(const Trk::TrackParameters &, const Trk::TrackParameters &, double, bool) const
Compute the four-momentum of a particle according to a mass hypothesis.
Definition: ConversionPostSelector.cxx:276
InDet::ConversionPostSelector::m_invMassCut
std::vector< double > m_invMassCut
Chi2 cut.
Definition: ConversionPostSelector.h:60
InDet::ConversionPostSelector::selectSecVtxCandidate
bool selectSecVtxCandidate(xAOD::Vertex *myCandidate, int flag, std::vector< Amg::Vector3D > &trkL, int &) const
Definition: ConversionPostSelector.cxx:180
Trk::ParametersBase::charge
double charge() const
Returns the charge.
InDet::ConversionPostSelector::m_minRadius
std::vector< double > m_minRadius
Converted photon reconstructed momentum at vertex cut.
Definition: ConversionPostSelector.h:62
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:66
InDet
DUMMY Primary Vertex Finder.
Definition: VP1ErrorUtils.h:36
conifer::pow
constexpr int pow(int x)
Definition: conifer.h:20
M_PI
#define M_PI
Definition: ActiveFraction.h:11
InDet::ConversionPostSelector::m_minPt
double m_minPt
Converted photon reconstructed vertex radial position cut.
Definition: ConversionPostSelector.h:63
dqt_zlumi_pandas.mass
mass
Definition: dqt_zlumi_pandas.py:170
InDet::ConversionPostSelector::m_maxChi2
std::vector< double > m_maxChi2
Properties for track selection: all cuts are ANDed.
Definition: ConversionPostSelector.h:59
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
InDet::ConversionPostSelector::m_sigmaLambda
double m_sigmaLambda
Definition: ConversionPostSelector.h:73
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Trk::electron
@ electron
Definition: ParticleHypothesis.h:27
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::pion
@ pion
Definition: ParticleHypothesis.h:29
master.flag
bool flag
Definition: master.py:29
InDet::ConversionPostSelector::decorateVertex
static void decorateVertex(xAOD::Vertex &vertex, float inv_mass, float pt1, float pt2, float fR, float deltaPhiVtxTrk)
Decorate vertices with values used in post selector.
Definition: ConversionPostSelector.cxx:312
test_pyathena.parent
parent
Definition: test_pyathena.py:15
InDet::ConversionPostSelector::m_decorateVertices
bool m_decorateVertices
Maximum difference in phi between reconstructed vertex and track at vertex.
Definition: ConversionPostSelector.h:66
Trk::ParametersBase
Definition: ParametersBase.h:55
InDet::ConversionPostSelector::ConversionPostSelector
ConversionPostSelector(const std::string &type, const std::string &name, const IInterface *parent)
Definition: ConversionPostSelector.cxx:30
InDet::ConversionPostSelector::initialize
virtual StatusCode initialize() override
Definition: ConversionPostSelector.cxx:78
Vertex.h
ParticleHypothesis.h
Trk::ParticleMasses::mass
constexpr double mass[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:53
twopi
constexpr double twopi
Definition: VertexPointEstimator.cxx:16
InDet::ConversionPostSelector::interfaceID
static const InterfaceID & interfaceID()
Definition: ConversionPostSelector.cxx:74
Trk::ParametersBase::pT
double pT() const
Access method for transverse momentum.
InDet::ConversionPostSelector::m_fitMomentum
std::vector< double > m_fitMomentum
Invariant mass cut.
Definition: ConversionPostSelector.h:61
InDet::ConversionPostSelector::selectConversionCandidate
bool selectConversionCandidate(xAOD::Vertex *myCandidate, int flag, std::vector< Amg::Vector3D > &trkL) const
Conversion candidate post-fit selectors.
Definition: ConversionPostSelector.cxx:86
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
InDet::ConversionPostSelector::m_maxdR
double m_maxdR
Pt of the two participating tracks at the vertex.
Definition: ConversionPostSelector.h:64
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
InDet::ConversionPostSelector::m_massLambda
double m_massLambda
Definition: ConversionPostSelector.h:72
InDet::ConversionPostSelector::finalize
virtual StatusCode finalize() override
Definition: ConversionPostSelector.cxx:82
TrackParticle.h
Trk::ParametersBase::momentum
const Amg::Vector3D & momentum() const
Access method for the momentum.
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
InDet::ConversionPostSelector::m_massK0
double m_massK0
Decorate vertices with values used for vertex selection.
Definition: ConversionPostSelector.h:70
declareProperty
#define declareProperty(n, p, h)
Definition: BaseFakeBkgTool.cxx:15
InDet::ConversionPostSelector::~ConversionPostSelector
virtual ~ConversionPostSelector()
AthAlgTool
Definition: AthAlgTool.h:26
ConversionPostSelector.h
InDet::ConversionPostSelector::m_sigmaK0
double m_sigmaK0
Definition: ConversionPostSelector.h:71
InDet::ConversionPostSelector::m_maxPhiVtxTrk
double m_maxPhiVtxTrk
Distance of first track hit- reconstructed vertex radial position.
Definition: ConversionPostSelector.h:65
InDet::ConversionPostSelector::m_nsig
int m_nsig
Definition: ConversionPostSelector.h:74
Trk::phi0
@ phi0
Definition: ParamDefs.h:71
Generated on Sun Jun 30 2024 21:12:57 for ATLAS Offline Software by   doxygen 1.8.18