ATLAS Offline Software
ConversionPostSelector.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 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  declareInterface<ConversionPostSelector>(this);
35  }
36 
37  const InterfaceID& ConversionPostSelector::interfaceID() {
38  return IID_IConversionPostSelector;
39  }
40 
41  StatusCode ConversionPostSelector::initialize() {
42  return StatusCode::SUCCESS;
43  }
44 
45  StatusCode ConversionPostSelector::finalize() {
46  return StatusCode::SUCCESS;
47  }
48 
49  bool ConversionPostSelector::selectConversionCandidate(xAOD::Vertex * vertex, int flag,
50  std::vector<Amg::Vector3D>& trkL) const{
51  bool pass = true;
52 
53  //Determine the cuts
54  double maxChi2 = 1000.;
55  double invMassCut = 1000.;
56  double fitMomentum = 0.;
57  double radius = 1000.;
58 
59  if(flag==0) { // Si+Si
60  maxChi2 = m_maxChi2[0] ; //Vertex fit chi2 cut
61  invMassCut = m_invMassCut[0] ; //Fit invariant mass cut
62  fitMomentum = m_fitMomentum[0]; //Photon fitted momentum
63  radius = m_minRadius[0] ; //Minimum acceptable radius of conversion vertex
64  }
65  if(flag==1) { // Si+TRT
66  maxChi2 = m_maxChi2[1] ; //Vertex fit chi2 cut
67  invMassCut = m_invMassCut[1] ; //Fit invariant mass cut
68  fitMomentum = m_fitMomentum[1]; //Photon fitted momentum
69  radius = m_minRadius[1] ; //Minimum acceptable radius of conversion vertex
70  }
71  if(flag==2) { // TRT+TRT
72  maxChi2 = m_maxChi2[2] ; //Vertex fit chi2 cut
73  invMassCut = m_invMassCut[2] ; //Fit invariant mass cut
74  fitMomentum = m_fitMomentum[2]; //Photon fitted momentum
75  radius = m_minRadius[2] ; //Minimum acceptable radius of conversion vertex
76  }
77 
78  //chi2 cut
79  if (vertex->nTrackParticles() != 2) {
80  ATH_MSG_DEBUG("Incorrect number of tracks used in conversion fit.");
81  pass = false;
82  } else {
83 
84  float reducedChi2 = vertex->chiSquared()/vertex->numberDoF();
85  if (reducedChi2 > maxChi2) pass = false;
86  if (reducedChi2 > maxChi2) pass = false;
87 
88  //Minimum radius cut
89  double vtxR = vertex->position().perp();
90  if(vtxR < radius) pass = false;
91 
92  // Parameters at vertex
93  std::vector< Trk::VxTrackAtVertex >& trkAtVx = vertex->vxTrackAtVertex();
94  if (trkAtVx.size() != 2 ||
95  !trkAtVx[0].perigeeAtVertex() ||
96  !trkAtVx[1].perigeeAtVertex())
97  {
98  ATH_MSG_DEBUG("VxTrackAtVertex or perigeeAtVertex not available");
99  return false;
100  }
101  const Trk::TrackParameters& perigee1 = *(trkAtVx[0].perigeeAtVertex());
102  const Trk::TrackParameters& perigee2 = *(trkAtVx[1].perigeeAtVertex());
103 
104  //invariant mass
105  CLHEP::HepLorentzVector momentum;
106  Amg::Vector3D sum_mom = perigee1.momentum() + perigee2.momentum();
107  double m2 = std::pow(Trk::ParticleMasses::mass[Trk::electron],2);
108  double ee = std::sqrt(m2 + perigee1.momentum().mag2()) + std::sqrt(m2 + perigee2.momentum().mag2());
109  momentum.setPx(sum_mom.x()); momentum.setPy(sum_mom.y()); momentum.setPz(sum_mom.z()); momentum.setE(ee);
110  double inv_mass = momentum.m();
111  double photonP = std::sqrt(momentum.x()*momentum.x() + momentum.y()*momentum.y());
112  double pt1 = perigee1.pT(); double pt2 = perigee2.pT();
113 
114  if (pt1<m_minPt || pt2<m_minPt) pass = false;
115  if (std::fabs(inv_mass) > invMassCut) pass = false;
116  if (photonP < fitMomentum) pass = false;
117 
118  double fR = 1000.;
119  std::vector<Amg::Vector3D>::const_iterator ipb=trkL.begin();
120  std::vector<Amg::Vector3D>::const_iterator ipbe=trkL.end();
121  for(; ipb!=ipbe;++ipb){
122  double tmpfR = (*ipb).perp();
123  if(tmpfR<fR) fR = tmpfR;
124  }
125  if(flag==1 && fR-vtxR<m_maxdR) pass = false;
126 
127  double PhiVtxTrk = vertex->position().phi() - perigee1.parameters()[Trk::phi0];
128  if (PhiVtxTrk < -M_PI) PhiVtxTrk += twopi;
129  if (PhiVtxTrk > M_PI) PhiVtxTrk -= twopi;
130  if (std::fabs(PhiVtxTrk)>m_maxPhiVtxTrk) pass = false;
131 
132  if (pass && m_decorateVertices)
133  {
134  ATH_MSG_DEBUG("Decorating vertex with values used in post selector");
135  decorateVertex(*vertex, inv_mass, pt1, pt2, fR, std::fabs(PhiVtxTrk) );
136  }
137 
138  }
139  return pass;
140  }
141 
142  bool
143  ConversionPostSelector::selectSecVtxCandidate(
144  xAOD::Vertex* vertex,
145  int flag,
146  std::vector<Amg::Vector3D>& trkL,
147  int& type) const
148  {
149  bool pass = true;
150  bool isK0 = false;
151  bool isLambda = false;
152  bool isLambdaBar = false;
153  int kind = -1;
154 
155  //Determine the cuts
156  double maxChi2 = 1000.;
157  double fitMomentum = 0.;
158  double radius = 1000.;
159 
160  if(flag==0) {
161  maxChi2 = m_maxChi2[0] ; //Vertex fit chi2 cut
162  fitMomentum = m_fitMomentum[0]; //Photon fitted momentum
163  radius = m_minRadius[0] ; //Minimum acceptable radius of conversion vertex
164  }
165  if(flag==1) {
166  maxChi2 = m_maxChi2[1] ; //Vertex fit chi2 cut
167  fitMomentum = m_fitMomentum[1]; //Photon fitted momentum
168  radius = m_minRadius[1] ; //Minimum acceptable radius of conversion vertex
169  }
170  if(flag==2) {
171  maxChi2 = m_maxChi2[2] ; //Vertex fit chi2 cut
172  fitMomentum = m_fitMomentum[2]; //Photon fitted momentum
173  radius = m_minRadius[2] ; //Minimum acceptable radius of conversion vertex
174  }
175 
176  //chi2 cut
177 
178  if (vertex->nTrackParticles() != 2) {
179  ATH_MSG_DEBUG("Incorrect number of tracks used in conversion fit.");
180  pass = false;
181  } else {
182  float reducedChi2 = vertex->chiSquared()/vertex->numberDoF();
183  if (reducedChi2 > maxChi2) pass = false;
184 
185  //Minimum radius cut
186  double vtxR = vertex->position().perp();
187  if(vtxR < radius) pass = false;
188 
189  // Parameters at vertex
190  std::vector< Trk::VxTrackAtVertex >& trkAtVx = vertex->vxTrackAtVertex();
191  if (trkAtVx.size() != 2 ||
192  !trkAtVx[0].perigeeAtVertex() ||
193  !trkAtVx[1].perigeeAtVertex())
194  {
195  ATH_MSG_DEBUG("VxTrackAtVertex or perigeeAtVertex not available");
196  return false;
197  }
198  const Trk::TrackParameters& perigee1 = *(trkAtVx[0].perigeeAtVertex());
199  const Trk::TrackParameters& perigee2 = *(trkAtVx[1].perigeeAtVertex());
200 
201  double pt1 = perigee1.pT(); double pt2 = perigee2.pT();
202  if (pt1<m_minPt || pt2<m_minPt) pass = false;
203 
204  double fR = 1000.;
205  std::vector<Amg::Vector3D>::const_iterator ipb=trkL.begin();
206  std::vector<Amg::Vector3D>::const_iterator ipbe=trkL.end();
207  for(; ipb!=ipbe;++ipb){
208  double tmpfR = (*ipb).perp();
209  if(tmpfR<fR) fR = tmpfR;
210  }
211  if(flag==1 && fR-vtxR<m_maxdR) pass = false;
212 
213  //invariant mass. First assume K0, if failed assume Lambda
214  CLHEP::HepLorentzVector momentumK0 = fourP(perigee1,perigee2,m_massK0,false);
215  double inv_massK0 = momentumK0.m();
216  if (std::fabs(inv_massK0-m_massK0) <= m_nsig*m_sigmaK0) isK0 = true;
217  CLHEP::HepLorentzVector momentumL = fourP(perigee1,perigee2,m_massLambda,false);
218  double inv_massL = momentumL.m();
219  if (std::fabs(inv_massL-m_massLambda) <= m_nsig*m_sigmaLambda) isLambda = true;
220  CLHEP::HepLorentzVector momentumLb = fourP(perigee1,perigee2,m_massLambda,true);
221  double inv_massLb = momentumLb.m();
222  if (std::fabs(inv_massLb-m_massLambda) <= m_nsig*m_sigmaLambda) isLambdaBar = true;
223  if (!isLambdaBar && !isLambda && !isK0) pass = false;
224  CLHEP::HepLorentzVector momentum;
225  if(isK0 && isLambda && !isLambdaBar) {momentum = momentumK0; kind = 110;}
226  if(isK0 && isLambdaBar && !isLambda) {momentum = momentumK0; kind = 101;}
227  if(isK0 && !isLambda && !isLambdaBar) {momentum = momentumK0; kind = 100;}
228  if(!isK0 && isLambda && !isLambdaBar) {momentum = momentumL; kind = 10;}
229  if(!isK0 && isLambdaBar && !isLambda) {momentum = momentumLb; kind = 1;}
230  if(!isK0 && isLambda && isLambdaBar) {momentum = momentumL; kind = 11;}
231  double particleP = std::sqrt(momentum.x()*momentum.x() + momentum.y()*momentum.y());
232  if (particleP < fitMomentum) pass = false;
233  }
234  type = kind;
235  return pass;
236  }
237 
238  CLHEP::HepLorentzVector
239  ConversionPostSelector::fourP(const Trk::TrackParameters& per1,
240  const Trk::TrackParameters& per2,
241  double mass,
242  bool isBar) const
243  {
244  CLHEP::HepLorentzVector momentum;
245  Amg::Vector3D sum_mom = per1.momentum() + per2.momentum();
246  double mp1 = 0.; double mp2 = 0.;
247  if(mass==m_massK0) {
248  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
249  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
250  }else{
251  if(!isBar){
252  if(per1.charge()>0) {
253  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
254  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
255  } else {
256  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
257  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
258  }
259  }else{
260  if(per1.charge()>0) {
261  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
262  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
263  } else {
264  mp2 = std::pow(Trk::ParticleMasses::mass[Trk::pion],2);
265  mp1 = std::pow(Trk::ParticleMasses::mass[Trk::proton],2);
266  }
267  }
268  }
269  double ee = std::sqrt(mp1 + per1.momentum().mag2()) + std::sqrt(mp2 + per2.momentum().mag2());
270  momentum.setPx(sum_mom.x()); momentum.setPy(sum_mom.y()); momentum.setPz(sum_mom.z()); momentum.setE(ee);
271  return momentum;
272  }
273 
274  void
275  ConversionPostSelector::decorateVertex(xAOD::Vertex& vertex,
276  float inv_mass,
277  float pt1,
278  float pt2,
279  float fR,
280  float deltaPhiVtxTrk)
281  {
282  static const SG::AuxElement::Accessor<float> accMass("mass");
283  accMass(vertex) = inv_mass;
284  static const SG::AuxElement::Accessor<float> accPt1("pt1");
285  accPt1(vertex) = pt1;
286  static const SG::AuxElement::Accessor<float> accPt2("pt2");
287  accPt2(vertex) = pt2;
288  static const SG::AuxElement::Accessor<float> accMinRfirstHit("minRfirstHit");
289  accMinRfirstHit(vertex) = fR;
290  static const SG::AuxElement::Accessor<float> accDeltaPhiVtxTrk("deltaPhiVtxTrk");
291  accDeltaPhiVtxTrk(vertex) = deltaPhiVtxTrk;
292  }
293 
294 } // namespace InDet
295 
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:239
InDet::ConversionPostSelector::selectSecVtxCandidate
bool selectSecVtxCandidate(xAOD::Vertex *myCandidate, int flag, std::vector< Amg::Vector3D > &trkL, int &) const
Definition: ConversionPostSelector.cxx:143
InDet::ConversionPostSelector::m_maxPhiVtxTrk
DoubleProperty m_maxPhiVtxTrk
Definition: ConversionPostSelector.h:73
Trk::ParametersBase::charge
double charge() const
Returns the charge.
InDet::ConversionPostSelector::m_massK0
static constexpr double m_massK0
Masses and mass ranges for different V0 hypotheses.
Definition: ConversionPostSelector.h:80
Base_Fragment.mass
mass
Definition: Sherpa_i/share/common/Base_Fragment.py:59
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition: Control/AthContainers/AthContainers/Accessor.h:68
InDet::ConversionPostSelector::m_decorateVertices
BooleanProperty m_decorateVertices
Definition: ConversionPostSelector.h:75
InDet
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
DoubleProperty m_minPt
Definition: ConversionPostSelector.h:69
InDet::ConversionPostSelector::m_sigmaLambda
static constexpr double m_sigmaLambda
Definition: ConversionPostSelector.h:83
InDet::ConversionPostSelector::m_sigmaK0
static constexpr double m_sigmaK0
Definition: ConversionPostSelector.h:81
ParticleGun_EoverP_Config.momentum
momentum
Definition: ParticleGun_EoverP_Config.py:63
InDet::ConversionPostSelector::m_maxdR
DoubleProperty m_maxdR
Definition: ConversionPostSelector.h:71
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:275
test_pyathena.parent
parent
Definition: test_pyathena.py:15
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:41
InDet::ConversionPostSelector::m_maxChi2
DoubleArrayProperty m_maxChi2
Properties for track selection: all cuts are ANDed.
Definition: ConversionPostSelector.h:60
Vertex.h
ParticleHypothesis.h
InDet::ConversionPostSelector::m_minRadius
DoubleArrayProperty m_minRadius
Definition: ConversionPostSelector.h:67
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:37
Trk::ParametersBase::pT
double pT() const
Access method for transverse momentum.
InDet::ConversionPostSelector::selectConversionCandidate
bool selectConversionCandidate(xAOD::Vertex *myCandidate, int flag, std::vector< Amg::Vector3D > &trkL) const
Conversion candidate post-fit selectors.
Definition: ConversionPostSelector.cxx:49
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
InDet::ConversionPostSelector::m_nsig
IntegerProperty m_nsig
Definition: ConversionPostSelector.h:84
InDet::ConversionPostSelector::m_massLambda
static constexpr double m_massLambda
Definition: ConversionPostSelector.h:82
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
ParticleGun_SamplingFraction.radius
radius
Definition: ParticleGun_SamplingFraction.py:96
InDet::ConversionPostSelector::m_fitMomentum
DoubleArrayProperty m_fitMomentum
Definition: ConversionPostSelector.h:64
InDet::ConversionPostSelector::m_invMassCut
DoubleArrayProperty m_invMassCut
Definition: ConversionPostSelector.h:62
InDet::ConversionPostSelector::finalize
virtual StatusCode finalize() override
Definition: ConversionPostSelector.cxx:45
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
AthAlgTool
Definition: AthAlgTool.h:26
ConversionPostSelector.h
Trk::phi0
@ phi0
Definition: ParamDefs.h:65
Generated on Thu Nov 7 2024 21:12:29 for ATLAS Offline Software by   doxygen 1.8.18