ATLAS Offline Software
ConstraintFit.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 #include "xAODEgamma/PhotonxAODHelpers.h"
6 #include "xAODTracking/TrackParticle.h"
7 #include "ZMassConstraint/ConstraintFit.h"
8 #include "ZMassConstraint/ConstraintFitInput.h"
9 #include "ZMassConstraint/ConstraintFitOutput.h"
10 
11 namespace ZMassConstraint
12 {
13  ConstraintFit::ConstraintFit( const std::string& name ):
14  asg::AsgTool( name ),
15  m_conHasWidth(true),
16  m_conMass(91187.6),
17  m_conWidth(2495.2),
18  m_resolution(0.01),
19  m_ignoreInputChecks(false),
20  m_parameters(3)
21  {
22  declareProperty( "Z_pdg_mass", m_conMass);
23  declareProperty( "Z_pdg_width", m_conWidth);
24  declareProperty( "hasWidth", m_conHasWidth);
25  declareProperty( "ignoreInputChecks", m_ignoreInputChecks);
26  declareProperty( "EgammaCalibAndSmearingTool", m_energyRescaler = ToolHandle<CP::IEgammaCalibrationAndSmearingTool>("CP::EgammaCalibrationAndSmearingTool"));
27  declareProperty( "MuonCalibrationAndSmearingTool", m_mu_resolSFTool = ToolHandle<CP::IMuonCalibrationAndSmearingTool>("CP::MuonCalibrationAndSmearingTool"));
28  }
29 
30  ConstraintFit::~ConstraintFit(void)
31  {}
32 
33 
34  StatusCode ConstraintFit::initialize()
35  {
36  ATH_MSG_INFO ("initialize: retrieve tools" << m_energyRescaler.name() << " and "
37  << m_mu_resolSFTool.name());
38 
39  if (!m_energyRescaler.retrieve().isSuccess()) {
40  ATH_MSG_ERROR ("initialize: unable to retrieve EgammaCalibrationAndSmearingTool");
41  return StatusCode::FAILURE;
42  }
43 
44  if (!m_mu_resolSFTool.retrieve().isSuccess()) {
45  ATH_MSG_ERROR ("initialize: unable to retrieve MuonCalibrationAndSmearingTool");
46  return StatusCode::FAILURE;
47  }
48 
49  // Return gracefully:
50  return StatusCode::SUCCESS;
51  }
52 
53  StatusCode
54  ConstraintFit::doMassFit (const ConstraintFitInput& input, ConstraintFitOutput& output)
55  {
56 
57  // For the moment, call the 'old' methods
58  if (massFitInterface(input)) {
59  ATH_MSG_ERROR( "Could not add in input" );
60  return StatusCode::FAILURE;
61  }
62 
63  massFitRun(output);
64 
65  // Return gracefully:
66  return StatusCode::SUCCESS;
67  }
68 
69  // void
70  // ConstraintFit::addParticle (const xAOD::IParticle & part, ConstraintFitInput& input)
71  // {
72  // if(part.type() == xAOD::Type::Electron) {
73  // const xAOD::Electron* el = dynamic_cast<const xAOD::Electron*> (&part);
74  // return addParticle(*el, input);
75  // }
76  // else if(part.type()==xAOD::Type::Muon) {
77  // const xAOD::Muon* mu = dynamic_cast<const xAOD::Muon*> (&part);
78  // return addParticle(*mu, input);
79  // }
80  // }
81 
82  void
83  ConstraintFit::addFSRParticle(const xAOD::IParticle& part,
84  const TLorentzVector& fsr4Vec,
85  ConstraintFitInput& input)
86  {
87 
88  ATH_MSG_DEBUG ("addFSRParticle: *** 4vecFsr: " << fsr4Vec.Pt() << "/" << fsr4Vec.Eta() << "/" << fsr4Vec.Phi());
89  // Get energy resolution - done differently for photons and electrons
90  float e_res = 0;
91  PATCore::ParticleType::Type phType = PATCore::ParticleType::Electron;
92  float cl_etaCalo = retrieve_eta_calo(part);
93  if (!m_energyRescaler.empty()) {
94  if (xAOD::Type::Photon == part.type()) {
95  const xAOD::Photon* fsrPf = dynamic_cast<const xAOD::Photon*>(&part);
96  phType = (xAOD::EgammaHelpers::isConvertedPhoton(fsrPf)) ?
97  PATCore::ParticleType::ConvertedPhoton :
98  PATCore::ParticleType::UnconvertedPhoton;
99  }
100  e_res = m_energyRescaler->resolution( fsr4Vec.E(), fsr4Vec.Eta(), cl_etaCalo, phType)*fsr4Vec.E();
101  }
102  else {
103  ATH_MSG_ERROR( "addFSRParticle: - cannot find EgammaCalibrationAndSmearingTool. ");
104  ATH_MSG_ERROR( "addFSRParticle: - please set property 'EgammaCalibAndSmearingToolName' with the name of your EgammaCalibrationAndSmearingTool - default is 'EgammaCalibrationAndSmearingTool'. ");
105  }
106 
107  AmgMatrix(5,5) photonCovarianceMatrix;
108  photonCovarianceMatrix.setZero();
109  photonCovarianceMatrix(d0,d0) = 0.000001;
110  photonCovarianceMatrix(z0,z0) = 0.000001;
111  photonCovarianceMatrix(phi0,phi0) = 0.000001;
112  photonCovarianceMatrix(theta,theta) = 0.000001;
113  photonCovarianceMatrix(qOverP,qOverP) = e_res*e_res;
114  photonCovarianceMatrix(d0,z0) = 0;
115  photonCovarianceMatrix(d0,phi0) = 0;
116  photonCovarianceMatrix(d0,theta) = 0;
117  photonCovarianceMatrix(d0,qOverP) = 0;
118  photonCovarianceMatrix(z0,phi0) = 0;
119  photonCovarianceMatrix(z0,theta) = 0;
120  photonCovarianceMatrix(z0,qOverP) = 0;
121  photonCovarianceMatrix(phi0,theta) = 0;
122  photonCovarianceMatrix(phi0,qOverP) = 0;
123  photonCovarianceMatrix(theta,qOverP) = 0;
124 
125  ATH_MSG_DEBUG ("addFSRParticle: FSR calib type,e, e_res " << phType << " " << fsr4Vec.E() << " " << e_res);
126 
127  // symmetrize it
128  for (int ii = 0; ii < 5; ii++)
129  for (int jj = ii + 1; jj < 5; jj++)
130  photonCovarianceMatrix(jj, ii) = photonCovarianceMatrix(ii, jj);
131 
132  bool isOK = doSanityChecksOnCovariance(fsr4Vec, photonCovarianceMatrix);
133  // Get covariance in cartesian CS
134  AmgMatrix(5,5) covarXYZ;
135  covarXYZ.setZero();
136  convertCovd0z0PhiThetaPToXYZ(fsr4Vec, photonCovarianceMatrix, covarXYZ);
137  // Save in input
138  input.addConstituent_FourVector_d0z0PhiThetaP(fsr4Vec, photonCovarianceMatrix, covarXYZ, isOK);
139 
140  ATH_MSG_DEBUG( "addFSRParticle: - added fsr");
141 
142  }
143 
144  void
145  ConstraintFit::addParticle(const xAOD::Muon& mu,
146  ConstraintFitInput& input,
147  MassConstraintMuonType muonType)
148  {
149 
150  // For muons, we allow the mass constraint to be done in one of three ways:
151  // 1) standard way - use muon momentum and primary track for cov matrix, or
152  // 2,3) ID or MS - use the ID or MS tracks for both the momentum and cov matrix
153 
154  TLorentzVector mu4vec; // to be filled according to muon type
155 
156  // setup accessors for muon decorations - set in MuonCalibrationAndSmearingTool
157  static const SG::AuxElement::ConstAccessor<float> muonSpectrometerPt ("MuonSpectrometerPt");
158  static const SG::AuxElement::ConstAccessor<float> innerDetectorPt ("InnerDetectorPt");
159 
160  // For the momentum,
161  // Get the track particle according to the requested muon type for the covariance matrix
162  const xAOD::TrackParticle* track = mu.primaryTrackParticle();
163  bool set4vec = false;
164  if (((isMS_MCMT == muonType) || (isID_MCMT == muonType)) && xAOD::Muon::Combined == mu.muonType()) {
165  if (isMS_MCMT == muonType) {
166  track = mu.trackParticle(xAOD::Muon::MuonSpectrometerTrackParticle);
167  if (!muonSpectrometerPt.isAvailable(mu))
168  ATH_MSG_ERROR( "addParticle: - could not get muonSpectrometerPt from muon. Please applyCorrection with the MuonCalibAndSmearTool");
169  if (!track) {
170  ATH_MSG_ERROR( "addParticle: - Combined muon is missing MS track particle. Using combined track particle");
171  track = mu.primaryTrackParticle();
172  }
173  mu4vec.SetPtEtaPhiM(muonSpectrometerPt(mu), track->eta(), track->phi(), mu.m());
174  set4vec = true;
175  ATH_MSG_DEBUG ("addParticle: set muon track to MS ");
176  }
177  else {
178  track = mu.trackParticle(xAOD::Muon::InnerDetectorTrackParticle);
179  if (!innerDetectorPt.isAvailable(mu))
180  ATH_MSG_ERROR( "addParticle: - could not get innerDetectorPt from muon. Please applyCorrection with the MuonCalibAndSmearTool");
181  mu4vec.SetPtEtaPhiM(innerDetectorPt(mu), track->eta(), track->phi(), mu.m());
182  set4vec = true;
183  ATH_MSG_DEBUG ("addParticle: set muon track to ID ");
184  }
185  }
186 
187  if (!set4vec) mu4vec.SetPtEtaPhiM(mu.pt(), mu.eta(), mu.phi(), mu.m());
188 
189  if (!track) {
190  ATH_MSG_ERROR( "addParticle: - could not get track particle for muon");
191  return;
192  }
193 
194  ATH_MSG_DEBUG ("addParticle: *** 4vecMu: " << mu4vec.Pt() << "/" << mu4vec.Eta() << "/" << mu4vec.Phi());
195 
196 
197  float muSF = 1;
198  int type = -1;
199 
200  // Comment out use of MuonCalibrationAndSmearingTool until the SF is implemented by Giacomo
201 
202  // if (m_mu_resolSFTool) {
203  // // Get scale factor for the muon covariance qoverp resolution.
204  // // For muonType == isCombMCMT, we get a SF differently according to the type of muon: Comb/ST/CT/SA
205  // // muonType == isID_MCMT, we use the type 2 for ID tracks
206  // // muonType == isMS_MCMT, we use the type 3 for MS tracks
207 
208  // if (muonType == isCombMCMT) {
209  // if (xAOD::Muon::Combined == mu.muonType()) type = 1;
210  // else if (xAOD::Muon::SegmentTagged == mu.muonType()) type = 2;
211  // else if (xAOD::Muon::CaloTagged == mu.muonType()) type = 2;
212  // else if (xAOD::Muon::MuonStandAlone == mu.muonType()) type = 3;
213  // else ATH_MSG_ERROR("addParticle: - **** MUON type is not CB/ST/SA/Calo **** ");
214  // }
215  // else if (muonType == isID_MCMT) {
216  // type = 2;
217  // if (xAOD::Muon::MuonStandAlone == mu.muonType()) type = 3;
218  // }
219  // else if (muonType == isMS_MCMT) {
220  // type = 3;
221  // if (xAOD::Muon::SegmentTagged == mu.muonType() ||
222  // xAOD::Muon::CaloTagged == mu.muonType()) type = 2;
223  // }
224 
225  // if (muonType < isCombMCMT) {
226  // ATH_MSG_ERROR( "addParticle: - Invalid muon type requested, please use .");
227  // }
228  // else {
229  // muSF = m_mu_resolSF->getResolutionScaleFactor(mu4vec, type);
230  // }
231 
232  // ATH_MSG_DEBUG ("addParticle: 4v,type,errSF " << mu4vec.Eta() << " " << type << " " << muSF);
233 
234 
235  // }
236  // else {
237  // ATH_MSG_ERROR( "addParticle: - MuonResSFTool is NOT set. Please call setMuonResSFTool with MuonResolutionAndMomentumScaleFactors tool.");
238  // }
239 
240 
241  AmgMatrix(5,5) covmatrix;
242  covmatrix.setZero();
243  covmatrix(d0,d0) = track->definingParametersCovMatrix()(d0,d0);
244  covmatrix(z0,z0) = track->definingParametersCovMatrix()(z0,z0);
245  covmatrix(phi0,phi0) = track->definingParametersCovMatrix()(phi0,phi0);
246  covmatrix(theta,theta) = track->definingParametersCovMatrix()(theta,theta);
247  covmatrix(qOverP,qOverP) = track->definingParametersCovMatrix()(qOverP,qOverP)*muSF*muSF;
248  covmatrix(d0,z0) = track->definingParametersCovMatrix()(d0,z0);
249  covmatrix(d0,phi0) = track->definingParametersCovMatrix()(d0,phi0);
250  covmatrix(d0,theta) = track->definingParametersCovMatrix()(d0,theta);
251  covmatrix(d0,qOverP) = track->definingParametersCovMatrix()(d0,qOverP)*muSF;
252  covmatrix(z0,phi0) = track->definingParametersCovMatrix()(z0,phi0);
253  covmatrix(z0,theta) = track->definingParametersCovMatrix()(z0,theta);
254  covmatrix(z0,qOverP) = track->definingParametersCovMatrix()(z0,qOverP)*muSF;
255  covmatrix(phi0,theta) = track->definingParametersCovMatrix()(phi0,theta);
256  covmatrix(phi0,qOverP) = track->definingParametersCovMatrix()(phi0,qOverP)*muSF;
257  covmatrix(theta,qOverP) = track->definingParametersCovMatrix()(theta,qOverP)*muSF;
258 
259  ATH_MSG_DEBUG ("addParticle: type,errSF,mu p,cov " << type << " " << muSF
260  << " " << mu4vec.E()
261  << " " << track->definingParametersCovMatrix()(d0, d0)
262  << " " << track->definingParametersCovMatrix()(z0, z0)
263  << " " << track->definingParametersCovMatrix()(phi0, phi0)
264  << " " << track->definingParametersCovMatrix()(theta, theta)
265  << " " << track->definingParametersCovMatrix()(qOverP,qOverP)*muSF*muSF
266  << " " << track->definingParametersCovMatrix()(d0, z0)
267  << " " << track->definingParametersCovMatrix()(d0, phi0)
268  << " " << track->definingParametersCovMatrix()(d0, theta)
269  << " " << track->definingParametersCovMatrix()(d0,qOverP)*muSF
270  << " " << track->definingParametersCovMatrix()(z0, phi0)
271  << " " << track->definingParametersCovMatrix()(z0, theta)
272  << " " << track->definingParametersCovMatrix()(z0,qOverP)*muSF
273  << " " << track->definingParametersCovMatrix()(phi0, theta)
274  << " " << track->definingParametersCovMatrix()(phi0,qOverP)*muSF
275  << " " << track->definingParametersCovMatrix()(theta,qOverP)*muSF);
276 
277  for (int ii = 0; ii < 5; ii++)
278  for (int jj = ii+1; jj < 5; jj++)
279  covmatrix(jj,ii) = covmatrix(ii,jj);
280 
281  ATH_MSG_VERBOSE( "addParticle: cov matrix \n" << covmatrix );
282 
283  // Transform cov of q/p to cov p in covariance matrix
284  double P = mu4vec.P();
285  AmgMatrix(5,5) jacobian0;
286  jacobian0.setZero();
287  jacobian0(0,0) = 1.;
288  jacobian0(1,1) = 1.;
289  jacobian0(2,2) = 1.;
290  jacobian0(3,3) = 1.;
291  jacobian0(4,4) = -P*P;
292 
293  AmgMatrix(5,5) newmatrix = jacobian0.transpose() * covmatrix * jacobian0;
294 
295  ATH_MSG_DEBUG( "addParticle: new matrix \n" << newmatrix );
296 
297  bool isOK = doSanityChecksOnCovariance(mu4vec, newmatrix);
298  // Get covariance in cartesian CS
299  AmgMatrix(5,5) covarXYZ;
300  covarXYZ.setZero();
301  convertCovd0z0PhiThetaPToXYZ(mu4vec, newmatrix, covarXYZ);
302  // Save in input
303  input.addConstituent_FourVector_d0z0PhiThetaP(mu4vec, newmatrix, covarXYZ, isOK);
304 
305  ATH_MSG_VERBOSE( "addParticle: - added muon");
306 
307  }
308 
309  void
310  ConstraintFit::addParticle (const xAOD::Electron& el,
311  float elEnergyRes,
312  ConstraintFitInput& input)
313  {
314  float el_E_res = elEnergyRes;
315  // if (!m_energyRescaler.empty()) {
316  // // Use el.e() which is assumed to be equal to the corrected cluster energy
317  // el_E_res = m_energyRescaler->resolution(el.caloCluster()->e(), el.caloCluster()->eta())*
318  // el.caloCluster()->e();
319  // }
320  // else ATH_MSG_ERROR( "addParticle: - energyRescaler is NOT AVAILABLE. ");
321 
322  const xAOD::TrackParticle* track = el.trackParticle();
323 
324  if (!track) {
325  ATH_MSG_ERROR( "addParticle: - ERROR could not get track particle for electron");
326  return;
327  }
328 
329  ATH_MSG_DEBUG ("addParticle: *** 4vecEl: " << el.pt() << "/" << el.eta() << "/" << el.phi());
330 
331  AmgMatrix(5,5) covmatrix;
332  covmatrix.setZero();
333  covmatrix(d0,d0) = track->definingParametersCovMatrix()(d0,d0);
334  covmatrix(z0,z0) = track->definingParametersCovMatrix()(z0,z0);
335  covmatrix(phi0,phi0) = track->definingParametersCovMatrix()(phi0,phi0);
336  covmatrix(theta,theta) = track->definingParametersCovMatrix()(theta,theta);
337  covmatrix(qOverP,qOverP) = el_E_res*el_E_res;
338  covmatrix(d0,z0) = track->definingParametersCovMatrix()(d0,z0);
339  covmatrix(d0,phi0) = track->definingParametersCovMatrix()(d0,phi0);
340  covmatrix(d0,theta) = track->definingParametersCovMatrix()(d0,theta);
341  covmatrix(d0,qOverP) = 0;
342  covmatrix(z0,phi0) = track->definingParametersCovMatrix()(z0,phi0);
343  covmatrix(z0,theta) = track->definingParametersCovMatrix()(z0,theta);
344  covmatrix(z0,qOverP) = 0;
345  covmatrix(phi0,theta) = track->definingParametersCovMatrix()(phi0,theta);
346  covmatrix(phi0,qOverP) = 0;
347  covmatrix(theta,qOverP) = 0;
348 
349  ATH_MSG_DEBUG( "addParticle: el p,cov"
350  << " " << el_E_res << " " << el.e()
351  << " " << track->definingParametersCovMatrix()(d0, d0)
352  << " " << track->definingParametersCovMatrix()(z0, z0)
353  << " " << track->definingParametersCovMatrix()(phi0, phi0)
354  << " " << track->definingParametersCovMatrix()(theta, theta)
355  << " " << track->definingParametersCovMatrix()(d0, z0)
356  << " " << track->definingParametersCovMatrix()(d0, phi0)
357  << " " << track->definingParametersCovMatrix()(d0, theta)
358  << " " << track->definingParametersCovMatrix()(z0, phi0)
359  << " " << track->definingParametersCovMatrix()(z0, theta)
360  << " " << track->definingParametersCovMatrix()(phi0, theta));
361 
362 
363 
364  for(int ii=0;ii<5;ii++)
365  for(int jj=ii+1;jj<5;jj++)
366  covmatrix(jj,ii) = covmatrix(ii,jj);
367 
368  bool isOK = doSanityChecksOnCovariance(el.p4(), covmatrix);
369  // Get covariance in cartesian CS
370  AmgMatrix(5,5) covarXYZ;
371  covarXYZ.setZero();
372  convertCovd0z0PhiThetaPToXYZ(el.p4(), covmatrix, covarXYZ);
373  // Save in input
374  input.addConstituent_FourVector_d0z0PhiThetaP(el.p4(), covmatrix, covarXYZ, isOK);
375 
376  ATH_MSG_DEBUG( "addParticle: - el \n" << covmatrix);
377 
378  ATH_MSG_VERBOSE( "addParticle: - el\n" << input.getCovarianceCartesian(input.getNConstituents()-1));
379 
380  ATH_MSG_VERBOSE( "addParticle: - added electron");
381  }
382 
383  int
384  ConstraintFit::massFitInterface(const ZMassConstraint::ConstraintFitInput& theInput)
385  {
386  m_objmass.clear();
387  m_theInput = theInput;
388 
389  // the m_parameters first parameters of pin are the fitted parameters
390  // the next one is the mass of the particle
391  if(theInput.getNConstituents() <2 ) {
392  ATH_MSG_ERROR( "massFitInterface: Found number of input particles less than 2");
393  return 1;
394  }
395  else {
396  ATH_MSG_DEBUG( "massFitInterface: nparticles " << theInput.getNConstituents());
397  }
398 
399  m_nobj = theInput.getNConstituents();
400  unsigned int dimension = m_parameters * m_nobj;
401  m_parametersInit = Amg::MatrixX(dimension, 1);
402  m_parametersInit.setZero();
403  m_covarianceInit = Amg::MatrixX(dimension, dimension);
404  m_covarianceInit.setZero();
405  m_parametersFinal= Amg::MatrixX(dimension, 1);
406  m_parametersFinal.setZero();
407  m_covarianceFinal= Amg::MatrixX(dimension, dimension);
408  m_covarianceFinal.setZero();
409 
410  for (unsigned int iobj = 0; iobj < m_nobj; iobj++)
411  {
412  m_objmass.push_back(theInput.getConstituentFourVector(iobj).M());
413  m_parametersInit(0 + iobj*3, 0) = theInput.getConstituentFourVector(iobj).Px();
414  m_parametersInit(1 + iobj*3, 0) = theInput.getConstituentFourVector(iobj).Py();
415  m_parametersInit(2 + iobj*3, 0) = theInput.getConstituentFourVector(iobj).Pz();
416 
417  for (unsigned int i = 0; i < m_parameters; i++)
418  for (unsigned int j = 0; j < m_parameters; j++)
419  m_covarianceInit(i + 3*iobj, j + 3*iobj) = (theInput.getCovarianceCartesian(iobj))(i + 2, j + 2); //keep only the lower right 3x3 part
420  }
421 
422  for(unsigned int i=0;i<dimension;i++)
423  {
424  m_parametersFinal(i,0) = m_parametersInit(i,0);
425  for(unsigned int j=0;j<dimension;j++)
426  m_covarianceFinal(i,j) = m_covarianceInit(i,j);
427  }
428  ATH_MSG_VERBOSE( "massFitInterface: parameters \n" << m_parametersInit );
429  ATH_MSG_VERBOSE( "massFitInterface: cov matrix \n" << m_covarianceInit );
430  return 0;
431  }
432 
433  double
434  ConstraintFit::likelihoodMass(double MassResol)
435  {
436  const int rows =4;
437  const int columns = m_nobj;
438  double **p = 0;
439  p = new double *[rows] ;
440  for( int i = 0 ; i < rows ; i++ )
441  p[i] = new double[columns];
442 
443  for(unsigned int iobj=0; iobj<m_nobj; iobj++)
444  {
445  p[3][iobj] = m_objmass[iobj]*m_objmass[iobj];
446  for(int i=0;i<3;i++)
447  {
448  p[i][iobj] = m_parametersInit(i + 3*iobj, 0);
449  p[3][iobj] += p[i][iobj]*p[i][iobj];
450  }
451  p[3][iobj] = sqrt(p[3][iobj]);
452  }
453  double Etot = 0.;
454  double Xtot = 0.;
455  double Ytot = 0.;
456  double Ztot = 0.;
457  for(unsigned int ij=0;ij<m_nobj;ij++)
458  {
459  Etot += p[3][ij];
460  Xtot += p[0][ij];
461  Ytot += p[1][ij];
462  Ztot += p[2][ij];
463  }
464  for( int i = 0 ; i < rows ; i++ )
465  delete [] p[i];
466  delete [] p;
467 
468  double initMass = Etot * Etot - Xtot * Xtot - Ytot * Ytot - Ztot * Ztot;
469  initMass = sqrt(initMass);
470 
471  double sig = MassResol;
472  double xLeft = initMass;
473  double xRight = m_conMass;
474  if (m_conMass < initMass) {
475  xLeft = m_conMass;
476  xRight = initMass;
477  }
478  double dLinitL = (initMass - xLeft) / sig / sig - 2.*(xLeft - m_conMass) / ((xLeft - m_conMass) * (xLeft - m_conMass) + m_conWidth * m_conWidth);
479  double dLinitR = (initMass - xRight) / sig / sig - 2.*(xRight - m_conMass) / ((xRight - m_conMass) * (xRight - m_conMass) + m_conWidth * m_conWidth);
480  if (dLinitL * dLinitR < 0.) {
481  while (xRight - xLeft > 1.) { //1 MeV
482  double xM = (xRight + xLeft) / 2.;
483  double dL = (initMass - xM) / sig / sig - 2.*(xM - m_conMass) / ((xM - m_conMass) * (xM - m_conMass) + m_conWidth * m_conWidth);
484  if (dL * dLinitL < 0.) {
485  xRight = xM;
486  dLinitR = dL;
487  } else {
488  xLeft = xM;
489  dLinitL = dL;
490  }
491  }
492  return (xLeft + xRight) / 2.;
493  } else {
494  if (dLinitL > dLinitR)
495  return xLeft;
496  else
497  return xRight;
498  }
499  }
500 
501  double
502  ConstraintFit::likelihoodMass2(void)
503  {
504  ATH_MSG_VERBOSE( "likelihoodMass2: entering" );
505  const int rows =4;
506  const int columns = m_nobj;
507  double **p = 0;
508  p = new double *[rows] ;
509  for( int i = 0 ; i < rows ; i++ )
510  p[i] = new double[columns];
511 
512  for(unsigned int iobj=0; iobj<m_nobj; iobj++) {
513  p[3][iobj] = m_objmass[iobj]*m_objmass[iobj];
514  for(int i=0;i<3;i++) {
515  p[i][iobj] = m_parametersInit(i + 3*iobj, 0);
516  p[3][iobj] += p[i][iobj]*p[i][iobj];
517  }
518  p[3][iobj] = sqrt(p[3][iobj]);
519  }
520  double etot = 0.;
521  double xtot = 0.;
522  double ytot = 0.;
523  double ztot = 0.;
524  for(unsigned int ij=0;ij<m_nobj;ij++) {
525  etot += p[3][ij];
526  xtot += p[0][ij];
527  ytot += p[1][ij];
528  ztot += p[2][ij];
529  }
530 
531  double initMass = etot * etot - xtot * xtot - ytot * ytot - ztot * ztot;
532  initMass = sqrt(initMass);
533  //what is the uncertainty on initMass?
534  Amg::MatrixX jacobianMass(3*m_nobj, 1);
535  jacobianMass.setZero();
536  for(unsigned int ii=0;ii<m_nobj;ii++) {
537  jacobianMass(3*ii + 0, 0) = (2.*(etot) * p[0][ii] / p[3][ii] - 2.* xtot);
538  jacobianMass(3*ii + 1, 0) = (2.*(etot) * p[1][ii] / p[3][ii] - 2.* ytot);
539  jacobianMass(3*ii + 2, 0) = (2.*(etot) * p[2][ii] / p[3][ii] - 2.* ztot);
540  }
541  ATH_MSG_VERBOSE( "likelihoodMass2: the jacobian mass/matrix \n" << jacobianMass << "\n" << m_covarianceInit );
542 
543 
544  Amg::MatrixX sigm = jacobianMass.transpose() * m_covarianceInit * jacobianMass;
545 
546  ATH_MSG_VERBOSE( "likelihoodMass2: sigm \n" << sigm );
547 
548  double sig = sigm(0, 0);
549  sig = sqrt(sig);
550  m_resolution = sig;
551  sig /= 2 * initMass;
552  double MassResol = sig;
553  double maxmass = initMass;
554 
555  ATH_MSG_VERBOSE( "likelihoodMass2: 2 " );
556 
557  double max = -(maxmass - initMass) * (maxmass - initMass) / 2. / MassResol / MassResol
558  - log((maxmass * maxmass - m_conMass * m_conMass) * (maxmass * maxmass - m_conMass * m_conMass) + m_conMass * m_conMass * m_conWidth * m_conWidth);
559 
560  for (int i = 1; i < 401; i++) {
561  double ytest = initMass + (m_conMass - initMass) / 400 * i;
562 
563  double val = -(ytest - initMass) * (ytest - initMass) / 2. / MassResol / MassResol
564  - log((ytest * ytest - m_conMass * m_conMass) * (ytest * ytest - m_conMass * m_conMass) + m_conMass * m_conMass * m_conWidth * m_conWidth);
565 
566  if (val > max) {
567  max = val;
568  maxmass = ytest;
569  }
570  }
571 
572  ATH_MSG_VERBOSE( "likelihoodMass2: 3 " );
573 
574  for( int i = 0 ; i < rows ; i++ )
575  delete [] p[i];
576  delete [] p;
577 
578  ATH_MSG_VERBOSE( "likelihoodMass2: maxmass " << maxmass );
579 
580  return maxmass;
581  }
582 
583  void
584  ConstraintFit::massFitRun(ConstraintFitOutput& output, double zresol)
585  {
586  ATH_MSG_VERBOSE( "massFitRun: entering ");
587 
588  if(!m_ignoreInputChecks && !m_theInput.isOK()) {
589  ATH_MSG_WARNING("ConstraintFit::massFitRun() the input is not OK.");
590  ATH_MSG_WARNING("ConstraintFit::massFitRun()(cont) Following H->4l group decision (Feb 1, 2013), this event will not be fit");
591  ATH_MSG_WARNING("ConstraintFit::massFitRun()(cont) To bypass this behaviour set the ingore flag");
592 
593  std::vector<TLorentzVector> particleList;
594  for(unsigned int i=0;i<m_nobj;i++)
595  {
596  TLorentzVector particle;
597  particle.SetXYZM(m_parametersInit(0 + i*m_parameters,0),
598  m_parametersInit(1 + i*m_parameters,0),
599  m_parametersInit(2 + i*m_parameters,0),
600  m_objmass[i]);
601  particleList.push_back(particle);
602 
603  }
604  Amg::MatrixX covOut(5*m_nobj,5*m_nobj);
605  covOut.setZero();
606  for(unsigned int iobj=0;iobj<m_nobj;iobj++)
607  for(int ii=0;ii<5;ii++)
608  for(int jj=0;jj<5;jj++)
609  covOut(ii+5*iobj,jj+5*iobj) = m_theInput.getCovarianceCartesian(iobj)(ii,jj);
610 
611  ATH_MSG_VERBOSE( "massFitRun: cov 1 ");
612 
613  // convert covariance matrix from cartesian to d0z0PhiThetaP
614  Amg::MatrixX covard0z0PhiThetaP;
615  covard0z0PhiThetaP.setZero();
616  convertCovXYZTod0z0PhiThetaP(particleList, covOut, covard0z0PhiThetaP);
617  // save output
618  output.setFitOutput(particleList, covOut, covard0z0PhiThetaP);
619 
620  // ATH_MSG_VERBOSE( "massFitRun: cov 1 " << covOut);
621 
622  return;
623  }
624 
625  ATH_MSG_VERBOSE( "massFitRun: 1 ");
626 
627  Amg::MatrixX parametersFit = m_parametersInit;
628  double mass = -1.;
629  if (!m_conHasWidth)
630  mass = m_conMass;
631  else if(zresol<0.)
632  mass = likelihoodMass2();
633  else
634  mass = likelihoodMass(zresol);
635 
636  ATH_MSG_VERBOSE( "massFitRun: resolution/mass " << m_resolution << "/" << mass);
637 
638  // do not use chi2:
639  // double chi2 = massFit(&m_parametersInit, &m_covarianceInit, mass, &parametersFit, &m_covarianceFinal);
640  massFit(m_parametersInit, m_covarianceInit, mass, parametersFit, m_covarianceFinal);
641 
642  ATH_MSG_VERBOSE( "massFitRun: 2");
643 
644  std::vector<TLorentzVector> particleList;
645  for(unsigned int i=0;i<m_nobj;i++) {
646  TLorentzVector particle;
647  particle.SetXYZM(m_parametersFinal(0 + i*m_parameters,0),
648  m_parametersFinal(1 + i*m_parameters,0),
649  m_parametersFinal(2 + i*m_parameters,0),
650  m_objmass[i]);
651  particleList.push_back(particle);
652  }
653 
654  ATH_MSG_VERBOSE( "massFitRun: 3 ");
655 
656  Amg::MatrixX covOut(5*m_nobj,5*m_nobj);
657  covOut.setZero();
658  for(unsigned int iobj=0;iobj<m_nobj;iobj++)
659  for(int ii=0;ii<2;ii++)
660  for(int jj=0;jj<2;jj++)
661  covOut(ii+5*iobj,jj+5*iobj) = m_theInput.getCovarianceCartesian(iobj)(ii,jj);
662  for(unsigned int iobj=0;iobj<m_nobj;iobj++)
663  for(int ii=2;ii<5;ii++)
664  for(int jj=2;jj<5;jj++)
665  covOut(ii+5*iobj,jj+5*iobj) = m_covarianceFinal(ii-2+iobj*3,jj-2+iobj*3);
666 
667  ATH_MSG_VERBOSE( "massFitRun: 4 ");
668 
669  // empirical corrections for constraint width
670  if(!m_conHasWidth) {
671  // correlation terms
672  for(unsigned int iobj1=0;iobj1<m_nobj;iobj1++) {
673  for(unsigned int iobj2=0;iobj2<m_nobj;iobj2++) {
674  if(iobj1 == iobj2) continue;
675  for(int ii=2;ii<5;ii++)
676  for(int jj=2;jj<5;jj++)
677  covOut(ii + 5*iobj1, jj+5*iobj2) = m_covarianceFinal(ii-2+iobj1*3,jj-2+iobj2*3);
678  }
679  }
680  }
681 
682  // ATH_MSG_DEBUG( "massFitRun: cov 2 ");
683 
684  // convert covariance matrix from cartesian to d0z0PhiThetaP
685  Amg::MatrixX covard0z0PhiThetaP;
686  covard0z0PhiThetaP.setZero();
687  convertCovXYZTod0z0PhiThetaP(particleList, covOut, covard0z0PhiThetaP);
688  // save output
689  output.setFitOutput(particleList, covOut, covard0z0PhiThetaP);
690 
691  // ATH_MSG_DEBUG( "massFitRun: cov 2 " << covOut);
692 
693 
694  ATH_MSG_VERBOSE( "massFitRun: 5 conHasWidth " << m_conHasWidth);
695 
696  // ATH_MSG_DEBUG( "massFitRun: 5.5 conHasWidth " << m_conHasWidth);
697 
698  if(m_conHasWidth) {
699  TLorentzVector aa;
700  for(unsigned int iobj = 0; iobj < m_nobj; iobj++) {
701  aa += m_theInput.getConstituentFourVector(iobj);
702  }
703  double massInit = aa.M();
704 
705  // ATH_MSG_DEBUG( "massFitRun: 5.6 massInit " << massInit);
706 
707  if(massInit<83.e3) {
708 
709  ATH_MSG_VERBOSE( "massFitRun: 6 \n" << covOut);
710 
711  for(unsigned int iobj = 0; iobj < m_nobj; iobj++) {
712  for(int ii = 0; ii < 5; ii++) {
713  for(int jj = 0; jj < 5; jj++) {
714  covOut(ii+5*iobj, jj+5*iobj) = m_theInput.getCovarianceCartesian(iobj)(ii, jj);
715  }
716  }
717  }
718 
719  ATH_MSG_VERBOSE( "massFitRun: 6.1 \n" << covOut);
720  }
721  else {
722 
723  // ATH_MSG_DEBUG( "massFitRun: 7 ");
724 
725  std::vector<double> resolInit;
726  std::vector<double> resolFinal;
727  std::vector<double> resolFinalnew;
728  resolInit.reserve(m_nobj);
729  resolFinal.reserve(m_nobj);
730  resolFinalnew.reserve(m_nobj);
731  double sumResolInit=0.;
732  const double constraintWidthSquare = /*m_resolution*m_resolution + */m_conWidth*m_conWidth/2.35/2.35;
733  for(unsigned int iobj = 0; iobj < m_nobj; iobj++) {
734 
735  // ATH_MSG_DEBUG( "massFitRun: 7 iobj " << iobj << " " << m_theInput.getNConstituents());
736 
737  resolInit.push_back((m_theInput.getConstituentCovariance(iobj))(4,4));
738 
739  // ATH_MSG_DEBUG( "massFitRun: 7 1 ");
740 
741  AmgMatrix(5,5) tmpMtx;
742  tmpMtx.setZero();
743  output.getConstituentCovarianced0z0PhiThetaP(iobj, tmpMtx);
744  resolFinal.push_back(tmpMtx(4,4));
745  // resolFinal.push_back((output.getConstituentCovarianced0z0PhiThetaP(iobj))(4,4));
746 
747  // ATH_MSG_DEBUG( "massFitRun: 7 2 ");
748 
749  sumResolInit += resolInit[iobj];
750  //std::cout << " ### "<< iobj <<"\t"<< resolInit[iobj] <<"\t"<<resolFinal[iobj]<<std::endl;
751  }
752 
753 
754  ATH_MSG_VERBOSE( "massFitRun: 8 \n" << covOut);
755 
756  for(unsigned int iobj=0;iobj<m_nobj;iobj++) {
757  resolFinalnew.push_back( resolFinal[iobj] + constraintWidthSquare *
758  resolInit[iobj]*resolInit[iobj]
759  /sumResolInit/sumResolInit );
760  covOut(4+5*iobj,4+5*iobj) = covOut(4+5*iobj,4+5*iobj)*resolFinalnew[iobj]/resolFinal[iobj];
761  covOut(3+5*iobj,3+5*iobj) = covOut(3+5*iobj,3+5*iobj)*resolFinalnew[iobj]/resolFinal[iobj];
762  covOut(2+5*iobj,2+5*iobj) = covOut(2+5*iobj,2+5*iobj)*resolFinalnew[iobj]/resolFinal[iobj];
763  //std::cout << " ### "<< iobj <<"\t"<< resolInit[iobj] <<"\t"<<resolFinal[iobj]<<"\t"<< resolFinalnew[iobj]<<std::endl;
764  }
765  ATH_MSG_VERBOSE( "massFitRun: 8.1 \n" << covOut);
766  }
767 
768  ATH_MSG_VERBOSE( "massFitRun: 9 ");
769 
770  // ATH_MSG_DEBUG( "massFitRun: cov 3 ");
771 
772  // convert covariance matrix from cartesian to d0z0PhiThetaP
773  Amg::MatrixX covard0z0PhiThetaP;
774  convertCovXYZTod0z0PhiThetaP(particleList, covOut, covard0z0PhiThetaP);
775  // save output
776  output.setFitOutput(particleList, covOut, covard0z0PhiThetaP);
777 
778  ATH_MSG_VERBOSE( "massFitRun: cov 3 \n" << covOut);
779 
780  }
781 
782  ATH_MSG_VERBOSE( "massFitRun: end ");
783  }
784 
785  double
786  ConstraintFit::massFitCalculation(const Amg::MatrixX& var, double mass, Amg::MatrixX& p0)
787  {
788  const int rows =4;
789  const int columns = m_nobj;
790 
791  double **p = 0;
792  p = new double *[rows] ;
793  for( int i = 0 ; i < rows ; i++ )
794  p[i] = new double[columns];
795 
796  for(unsigned int iobj=0; iobj<m_nobj; iobj++) {
797  p[3][iobj] = m_objmass[iobj]*m_objmass[iobj];
798  for(int i=0;i<3;i++)
799  {
800  p[i][iobj] = p0(i+3*iobj, 0);
801  p[3][iobj] += p[i][iobj]*p[i][iobj];
802  }
803  p[3][iobj] = sqrt(p[3][iobj]);
804  }
805 
806  double Etot = 0.;
807  double Xtot = 0.;
808  double Ytot = 0.;
809  double Ztot = 0.;
810  for(unsigned int ij=0;ij<m_nobj;ij++)
811  {
812  Etot += p[3][ij];
813  Xtot += p[0][ij];
814  Ytot += p[1][ij];
815  Ztot += p[2][ij];
816  }
817  double XYZtot[3] = {Xtot, Ytot, Ztot};
818  std::vector<double> constraintD;
819  for(unsigned int ij=0;ij<m_nobj;ij++)
820  for(int ik=0;ik<3;ik++)
821  constraintD.push_back((2.*(Etot) * p[ik][ij] / p[3][ij] - 2.* XYZtot[ik]));
822 
823  double constraintd = Etot * Etot - Xtot * Xtot - Ytot * Ytot - Ztot * Ztot;
824  constraintd = constraintd - mass * mass;
825  //std::cout << "constraint = " <<constraintd << " mass value " <<mass<<std::endl;
826  Amg::MatrixX D(1,3*m_nobj);
827  for(unsigned int ij=0;ij<3*m_nobj;ij++)
828  D(0, ij) = constraintD.at(ij);
829  AmgMatrix(1,1) d;
830  d(0, 0) = constraintd;
831 
832  Amg::MatrixX DVD(D * var * D.transpose());
833 
834  Amg::MatrixX DVDinverse(1,1);
835  DVDinverse=DVD.inverse();
836 
837  Amg::MatrixX VD(DVDinverse);
838  Amg::MatrixX lambda(VD * d);
839  Amg::MatrixX test(p0 - var * D.transpose() * lambda);
840  Amg::MatrixX Vp(var - var * D.transpose() * VD * D * var);
841  for(unsigned int ij=0;ij<3*m_nobj;ij++) {
842  p0(ij, 0) = test(ij, 0);
843  m_parametersFinal(ij,0) = p0(ij,0);
844  for(unsigned int jk = 0; jk < 3 * m_nobj; jk++)
845  m_covarianceFinal(ij,jk) = Vp(ij,jk);
846  }
847  double constraintValue = d(0, 0);
848 
849  for( int i = 0 ; i < rows ; i++ )
850  delete [] p[i];
851  delete [] p;
852  return constraintValue;
853  }
854 
855  double
856  ConstraintFit::massFit(const Amg::MatrixX& /*p0*/, const Amg::MatrixX& var, double mass, Amg::MatrixX& pOut,
857  Amg::MatrixX& /*varout*/)
858  {
859  Amg::MatrixX ivar = var.inverse();
860 
861  bool doIter = true;
862  int maxIterations = 20;
863  int iIter = 0;
864  double constraintValue = -1.e10;
865  ATH_MSG_VERBOSE( "massFit: transpose matrix \n" << pOut.transpose()
866  << " Iterations " << iIter
867  << " parameters \n" << m_parametersFinal
868  << " cov matrix \n" << m_covarianceFinal
869  << " Iteration " << iIter);
870  while (doIter) {
871  Amg::MatrixX p_old(pOut);
872  constraintValue = massFitCalculation(var, mass, pOut); //call the fit
873 
874  // Convergence criteria
875  // 1. the parameters should not change too much (<1.e-6 relative change)
876  // 2. the constraint should be satisfied very well (<1.e-6 absolute)
877  double maxDiff = 5.e15;
878  Amg::MatrixX diff = (pOut - p_old);
879  for (unsigned int i = 0; i < m_parameters * m_nobj; i++) {
880  diff(i, 0) = diff(i, 0) / p_old(i, 0);
881  if (maxDiff > diff(i, 0))
882  maxDiff = diff(i, 0);
883  }
884  if ((maxDiff < 1.e-4 && fabs(constraintValue) < 5.e-5) || maxIterations <= iIter)
885  doIter = false;
886  ATH_MSG_VERBOSE( "massFit: transpose matrix \n" << pOut.transpose()
887  << " Iterations " << iIter
888  << " parameters \n" << m_parametersFinal
889  << " cov matrix \n" << m_covarianceFinal
890  << " \n Iteration " << iIter << " doIter " << doIter);
891  iIter++;
892 
893  }
894  ATH_MSG_VERBOSE( "massFit: transpose matrix \n" << pOut.transpose()
895  << " Iterations " << iIter
896  << " parameters \n" << m_parametersFinal
897  << " cov matrix \n" << m_covarianceFinal
898  << " Iteration " << iIter);
899  // double chi2 = 0.;//calculateChi2(p0, var, mass);
900  return 0;
901  }
902 
904  double
905  ConstraintFit::getMassError (const ConstraintFitInput& firstInput,
906  const ConstraintFitInput& secondInput)
907  {
908  unsigned int nPartFirst = firstInput.getNConstituents();
909  unsigned int nPartSecond = secondInput.getNConstituents();
910  unsigned int nPart = nPartFirst + nPartSecond;
911  std::vector<TLorentzVector> particles(nPart);
912  std::vector<AmgMatrix(3,3)> covariances(nPart);
913  for ( unsigned int iobj = 0; iobj < firstInput.getNConstituents(); ++iobj) {
914  covariances[iobj].setZero();
915  particles[iobj] = firstInput.getConstituentFourVector(iobj);
916  firstInput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj]);
917  }
918  for ( unsigned int iobj = 0; iobj < secondInput.getNConstituents(); ++iobj) {
919  covariances[iobj + nPartFirst].setZero();
920  particles[iobj + nPartFirst] = secondInput.getConstituentFourVector(iobj);
921  secondInput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj + nPartFirst]);
922  }
923  return getMassError( particles, covariances );
924  }
925 
928  double
929  ConstraintFit::getMassError (const ConstraintFitOutput& fitOutput,
930  const ConstraintFitInput& extraInput)
931  {
932  unsigned int nPartFirst = fitOutput.getNConstituents();
933  unsigned int nPartSecond = extraInput.getNConstituents();
934  unsigned int nPart = nPartFirst + nPartSecond;
935  std::vector<TLorentzVector> particles(nPart);
936  std::vector<AmgMatrix(3,3)> covariances(nPart);
937  for ( unsigned int iobj = 0; iobj < fitOutput.getNConstituents(); ++iobj) {
938  covariances[iobj].setZero();
939  particles[iobj] = fitOutput.getConstituentFourVector(iobj);
940  fitOutput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj]);
941  }
942  for ( unsigned int iobj = 0; iobj < extraInput.getNConstituents(); ++iobj) {
943  covariances[iobj + nPartFirst].setZero();
944  particles[iobj + nPartFirst] = extraInput.getConstituentFourVector(iobj);
945  extraInput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj + nPartFirst]);
946  }
947  return getMassError( particles, covariances );
948  }
949 
950 
951 
953  double
954  ConstraintFit::getMassError (const ConstraintFitOutput& fitOutput,
955  const ConstraintFitOutput& secondFitOutput)
956  {
957  unsigned int nPartFirst = fitOutput.getNConstituents();
958  unsigned int nPartSecond = secondFitOutput.getNConstituents();
959  unsigned int nPart = nPartFirst + nPartSecond;
960  std::vector<TLorentzVector> particles(nPart);
961  std::vector<AmgMatrix(3,3)> covariances(nPart);
962  for ( unsigned int iobj = 0; iobj < fitOutput.getNConstituents(); ++iobj) {
963  covariances[iobj].setZero();
964  particles[iobj] = fitOutput.getConstituentFourVector(iobj);
965  fitOutput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj]);
966  }
967  for ( unsigned int iobj = 0; iobj < secondFitOutput.getNConstituents(); ++iobj) {
968  covariances[iobj + nPartFirst].setZero();
969  particles[iobj + nPartFirst] = secondFitOutput.getConstituentFourVector(iobj);
970  secondFitOutput.getConstituentCovariancePhiThetaP(iobj, covariances[iobj + nPartFirst]);
971  }
972  return getMassError( particles, covariances );
973  }
974 
975  double
976  ConstraintFit::getMassError(const std::vector<TLorentzVector>& particles,
977  const std::vector<AmgMatrix(3,3)>& covariances)
978  {
979  double masserror=0;
980 
981  if (particles.size() != covariances.size()) {
982  ATH_MSG_ERROR( "getMassError: Number of particles is not equal to the number of covariance matrices");
983  return 1E11;
984  }
985 
986  ATH_MSG_DEBUG ("getMassError: 1 " << particles.size() << "/" << covariances.size());
987 
988  //composite lorentz vector and its invariant mass
989  TLorentzVector combv;
990  for ( auto lv : particles ) combv += lv;
991  double invmass = combv.M();
992 
993  ATH_MSG_VERBOSE ("getMassError: 1.1 invmass " << invmass);
994 
995  //Calculation of 1x4 jacobian - used in derivation of var(M)
996  Amg::MatrixX jacobianNPtoM(1,4);
997  jacobianNPtoM.setZero();
998  jacobianNPtoM(0,0) = -1.*combv.Px()/invmass; // dM/dpxN
999  jacobianNPtoM(0,1) = -1.*combv.Py()/invmass; // dM/dpyN
1000  jacobianNPtoM(0,2) = -1.*combv.Pz()/invmass; // dM/dpzN
1001  jacobianNPtoM(0,3) = combv.E()/invmass; // dM/dEN
1002 
1003  ATH_MSG_VERBOSE ("getMassError: 2\n" << jacobianNPtoM);
1004 
1005  int iobj = -1;
1006  for ( auto lv : particles ) {
1007  ++iobj;
1008 
1009  ATH_MSG_VERBOSE ("getMassError: 1.2 iobj " << iobj);
1010 
1011  // Constituent 3*3 covariance matrix (Phi, Theta, P)
1012  const AmgMatrix(3,3)& covPhiThetaP = covariances[iobj];
1013 
1014  ATH_MSG_VERBOSE ("getMassError: 1.3 covPhiThetaP\n " << covPhiThetaP);
1015 
1016  double theta = lv.Theta();
1017  double phi = lv.Phi();
1018  double p = lv.P();
1019  double e = lv.E();
1020  double m = lv.M();
1021 
1022  ATH_MSG_VERBOSE ("getMassError: 2.1 iobj " << iobj << "covPhiThetaP\n" << covPhiThetaP);
1023 
1024  // Convert (Phi, Theta, P) to (Phi, eta, P), as covPhiEtaP
1025  AmgMatrix(3,3) jacobian1;
1026  jacobian1.setZero();
1027  jacobian1(0,0) = 1.; //
1028  jacobian1(1,1) = -1*(1./sin(theta)); // deta/dtheta
1029  jacobian1(2,2) = 1.; //
1030  AmgMatrix(3,3) covPhiEtaP = jacobian1 * covPhiThetaP * jacobian1.transpose();
1031 
1032  ATH_MSG_VERBOSE ("getMassError: 2.2 jacobian1\n" << jacobian1);
1033  ATH_MSG_VERBOSE ("getMassError: 2.4 covPhiEtaP\n" << covPhiEtaP);
1034 
1035 
1036  // Rearrange things to get (E, eta, phi, M), as covEEtaPhiM
1037  AmgMatrix(4,4) covEEtaPhiM;
1038  covEEtaPhiM.setZero();
1039  for(int i = 0; i < 3; i++) {
1040  for(int j = 0; j < 3; j++) {
1041  covEEtaPhiM(i,j) = covEEtaPhiM(j,i) = covPhiEtaP(2-i,2-j);
1042  }
1043  }
1044 
1045  ATH_MSG_VERBOSE ("getMassError: 3 covEEtaPhiM\n" << covEEtaPhiM);
1046 
1047  //convert (E, eta, phi, M) to (px,py,pz,E) representation, as covPxPyPzE
1048  AmgMatrix(4,4) jacobian2;
1049  jacobian2.setZero();
1050  jacobian2(0,0) = e/p * sin(theta)*cos(phi); // dpx/dE
1051  jacobian2(1,0) = e/p * sin(theta)*sin(phi); // dpy/dE
1052  jacobian2(2,0) = e/p * cos(theta); // dpy/dE
1053  jacobian2(3,0) = 1; // dE/dE
1054 
1055  jacobian2(0,1) = -p * cos(phi) * cos(theta)*sin(theta); // dpx/deta
1056  jacobian2(1,1) = -p * sin(phi) * cos(theta)*sin(theta); // dpy/deta
1057  jacobian2(2,1) = p * sin(theta)*sin(theta); // dpz/deta;
1058 
1059  jacobian2(0,2) = -p * sin(theta) * sin(phi); // dpx/dphi
1060  jacobian2(1,2) = p * sin(theta) * cos(phi); // dpy/dphi
1061 
1062  jacobian2(0,3) = -m/p * sin(theta)*cos(phi); // dpx/dM
1063  jacobian2(1,3) = -m/p * sin(theta)*sin(phi); // dpy/dM
1064  jacobian2(2,3) = -m/p * cos(theta); // dpz/dM
1065 
1066  AmgMatrix(4,4) covPxPyPzE = jacobian2 * covEEtaPhiM * jacobian2.transpose();
1067 
1068  ATH_MSG_VERBOSE ("getMassError: 4 jacobian2\n" << jacobian2);
1069  ATH_MSG_VERBOSE ("getMassError: 5 covPxPyPzE\n" << covPxPyPzE);
1070 
1071  //Get mass variance
1072  // (Similarity transform: C = A*B*A^T (from egammaFourMomentumError/GeneralUtils.cxx))
1073  Amg::MatrixX em = jacobianNPtoM * covPxPyPzE * jacobianNPtoM.transpose();
1074  masserror += em(0,0);
1075 
1076  ATH_MSG_VERBOSE ("getMassError: 6 em\n" << em);
1077 
1078  }
1079 
1080  //Return square root of variance.
1081  if (masserror < 0.) {
1082  ATH_MSG_WARNING( "getMassError: Mass covariance element less than zero! Returning 1E11 ...");
1083  masserror = 1E11;
1084  }
1085  return sqrt(masserror);
1086 
1087  }
1088 
1089  bool
1090  ConstraintFit::doSanityChecksOnCovariance(const TLorentzVector& vector,
1091  const AmgMatrix(5,5)& covar) const
1092  {
1093  bool isOK = true;
1094  if ( sqrt(covar(qP,qP))/vector.P()>1.0 ) {
1095  isOK = false;
1096  ATH_MSG_WARNING("ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Fractional uncertainty on P = "
1097  << sqrt(covar(qP,qP))/vector.P() << " is > 1 ...this is not ok...");
1098 
1099  }
1100 
1101  if( covar(theta,theta) > 1.e-2) {
1102  isOK = false;
1103  ATH_MSG_WARNING("ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Uncertainty on Theta = "
1104  << covar(theta,theta) << " is > 10^-2 ...this is not ok");
1105  }
1106  if( covar(phi0,phi0) > 1.e-2) {
1107  isOK = false;
1108  ATH_MSG_WARNING("ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Uncertainty on Phi = "
1109  << covar(phi0,phi0) << " is > 10^-2 ...this is not ok");
1110  }
1111 
1112  for(int i=0;i<5;i++) {
1113  for(int j=0;j<5;j++) {
1114  if(i==j) continue;
1115  if(fabs(covar(i,j))/sqrt(covar(i,i))/sqrt(covar(j,j)) > 1.) {
1116  ATH_MSG_WARNING( "ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance:: Off-diagonal term "
1117  << i << " " << j << " is > 1 - doesn't look ok");
1118  isOK = false;
1119  }
1120  }
1121  }
1122 
1123  return isOK;
1124  }
1125 
1126  void
1127  ConstraintFit::convertCovd0z0PhiThetaPToXYZ(const TLorentzVector& fourVec,
1128  const AmgMatrix(5,5)& covard0z0PhiThetaP,
1129  AmgMatrix(5,5)& covarXYZ) const
1130  {
1131  //going from d0,z0,phi,theta,P --> d0,z0,px,py,pz
1132  double phi = fourVec.Phi();
1133  double theta = fourVec.Theta();
1134  double P = fourVec.P();
1135  //std::cout << ":::::L \t"<<P<<"\t"<< phi<<"\t"<<theta<<std::endl;
1136  AmgMatrix(5,5) jacobian;
1137  jacobian.setZero();
1138  jacobian(0,0) = 1.;
1139  jacobian(1,1) = 1.;
1140  jacobian(2,2) = -P * TMath::Sin(theta) * TMath::Sin(phi);
1141  jacobian(2,3) = P * TMath::Sin(theta) * TMath::Cos(phi);
1142  jacobian(3,2) = P * TMath::Cos(theta) * TMath::Cos(phi);
1143  jacobian(3,3) = P * TMath::Cos(theta) * TMath::Sin(phi);
1144  jacobian(3,4) = -P * TMath::Sin(theta);
1145  jacobian(4,2) = TMath::Sin(theta) * TMath::Cos(phi);
1146  jacobian(4,3) = TMath::Sin(theta) * TMath::Sin(phi);
1147  jacobian(4,4) = TMath::Cos(theta);
1148 
1149  covarXYZ = jacobian.transpose() * covard0z0PhiThetaP * jacobian;
1150 
1151  // std::cout << "initial\n"<< covar << "final \n" << newcovariance << "jac\n"
1152  // << jacobian << std::endl;
1153 
1154  }
1155 
1156 
1157  void
1158  ConstraintFit::convertCovXYZTod0z0PhiThetaP(const std::vector<TLorentzVector>& particleList,
1159  const Amg::MatrixX& covarXYZ,
1160  Amg::MatrixX& covard0z0PhiThetaP) const
1161  {
1162  unsigned int matrixSize = 5 * particleList.size();
1163  Amg::MatrixX Jacobian(matrixSize, matrixSize);
1164  Jacobian.setZero();
1166  for(unsigned int ii = 0; ii < particleList.size(); ii++) {
1167  double phi = particleList.at(ii).Phi();
1168  double theta = particleList.at(ii).Theta();
1169  double P = particleList.at(ii).P();
1170  Jacobian( 5*ii, 5*ii) = 1.;
1171  Jacobian(1 + 5*ii, 1 + 5*ii) = 1.;
1172 
1173  Jacobian(2 + 5*ii, 2 + 5*ii) = -P * TMath::Sin(theta) * TMath::Sin(phi);
1174  Jacobian(2 + 5*ii, 3 + 5*ii) = P * TMath::Sin(theta) * TMath::Cos(phi);
1175  Jacobian(3 + 5*ii, 2 + 5*ii) = P * TMath::Cos(theta) * TMath::Cos(phi);
1176  Jacobian(3 + 5*ii, 3 + 5*ii) = P * TMath::Cos(theta) * TMath::Sin(phi);
1177  Jacobian(3 + 5*ii, 4 + 5*ii) = -P * TMath::Sin(theta);
1178  Jacobian(4 + 5*ii, 2 + 5*ii) = TMath::Sin(theta) * TMath::Cos(phi);
1179  Jacobian(4 + 5*ii, 3 + 5*ii) = TMath::Sin(theta) * TMath::Sin(phi);
1180  Jacobian(4 + 5*ii, 4 + 5*ii) = TMath::Cos(theta);
1181  }
1182 
1183  Amg::MatrixX Jacobianinverse(matrixSize, matrixSize);
1184  Jacobianinverse=Jacobian.inverse();
1185  covard0z0PhiThetaP = Jacobianinverse.transpose() * covarXYZ * Jacobianinverse;
1186  }
1187 
1188  float
1189  ConstraintFit::retrieve_eta_calo(const xAOD::IParticle& part) const
1190  {
1191  double eta_calo = 0;
1192  const xAOD::Egamma* eg = dynamic_cast<const xAOD::Egamma*> (&part);
1193  if (!eg) {
1194  ATH_MSG_ERROR("retrieve_eta_calo - unable to cast to Egamma");
1195  return eta_calo;
1196  }
1197  const xAOD::CaloCluster& cluster = *eg->caloCluster();
1198  if (cluster.retrieveMoment(xAOD::CaloCluster::ETACALOFRAME,
1199  eta_calo)) { }
1200  else if (cluster.isAvailable<float>("etaCalo")) {
1201  eta_calo = cluster.auxdata<float>("etaCalo");
1202  }
1203  else {
1204  ATH_MSG_ERROR("retrieve_eta_calo - etaCalo not available as auxilliary variable");
1205  ATH_MSG_WARNING("retrieve_eta_calo - using eta as etaCalo");
1206  eta_calo = cluster.eta();
1207  }
1208  return eta_calo;
1209  }
1210 }
1211 
1212 
LArG4FSStartPointFilter.part
part
Definition: LArG4FSStartPointFilter.py:21
ZMassConstraint::ConstraintFit::theta
@ theta
Definition: ConstraintFit.h:86
ZMassConstraint::ConstraintFit::m_conMass
double m_conMass
Definition: ConstraintFit.h:110
beamspotnt.var
var
Definition: bin/beamspotnt.py:1394
ZMassConstraint::ConstraintFit::m_parametersInit
Amg::MatrixX m_parametersInit
Definition: ConstraintFit.h:123
PATCore::ParticleType::UnconvertedPhoton
@ UnconvertedPhoton
Definition: PATCoreEnums.h:38
ZMassConstraint::ConstraintFitInput::getNConstituents
unsigned int getNConstituents() const
Number of particles.
Definition: ConstraintFitInput.h:22
Trk::ParticleSwitcher::particle
constexpr ParticleHypothesis particle[PARTICLEHYPOTHESES]
the array of masses
Definition: ParticleHypothesis.h:76
Amg::MatrixX
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixX
Dynamic Matrix - dynamic allocation.
Definition: EventPrimitives.h:29
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
max
#define max(a, b)
Definition: cfImp.cxx:41
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
ZMassConstraint::ConstraintFit::qP
@ qP
Definition: ConstraintFit.h:88
PATCore::ParticleType::Type
Type
Definition: PATCoreEnums.h:35
ParticleTest.eg
eg
Definition: ParticleTest.py:29
ZMassConstraint::ConstraintFit::massFitCalculation
double massFitCalculation(const Amg::MatrixX &var, double mass, Amg::MatrixX &p0)
Definition: ConstraintFit.cxx:786
ZMassConstraint::ConstraintFit::m_objmass
std::vector< double > m_objmass
Definition: ConstraintFit.h:122
ZMassConstraint::ConstraintFit::ConstraintFit
ConstraintFit(const std::string &name)
Create a proper constructor for Athena.
Definition: ConstraintFit.cxx:13
AthCommonDataStore< AthCommonMsg< AlgTool > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
hist_file_dump.d
d
Definition: hist_file_dump.py:137
ZMassConstraint::isID_MCMT
@ isID_MCMT
Definition: IConstraintFit.h:37
ZMassConstraint::ConstraintFitOutput
Definition: ConstraintFitOutput.h:13
asg
Definition: DataHandleTestTool.h:28
ZMassConstraint::isMS_MCMT
@ isMS_MCMT
Definition: IConstraintFit.h:38
mc.diff
diff
Definition: mc.SFGenPy8_MuMu_DD.py:14
ZMassConstraint::ConstraintFit::z0
@ z0
Definition: ConstraintFit.h:84
xAOD::Egamma_v1
Definition: Egamma_v1.h:56
ZMassConstraint::ConstraintFitOutput::getConstituentCovariancePhiThetaP
Amg::MatrixX getConstituentCovariancePhiThetaP(int index) const
Access to individual covariance PhiThetaP (3,3)
Definition: ConstraintFitOutput.h:31
ZMassConstraint::ConstraintFit::m_parametersFinal
Amg::MatrixX m_parametersFinal
Definition: ConstraintFit.h:125
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition: ConstAccessor.h:54
ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance
bool doSanityChecksOnCovariance(const TLorentzVector &vector, const AmgMatrix(5, 5)&covar) const
Definition: ConstraintFit.cxx:1090
TrigInDetValidation_Base.test
test
Definition: TrigInDetValidation_Base.py:144
ConstraintFitInput.h
ZMassConstraint::ConstraintFitInput::getConstituentFourVector
const TLorentzVector & getConstituentFourVector(int i) const
Access to individual particle 4-vec.
Definition: ConstraintFitInput.h:25
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
xAOD::IParticle
Class providing the definition of the 4-vector interface.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:40
ZMassConstraint::MassConstraintMuonType
MassConstraintMuonType
Definition: IConstraintFit.h:35
ZMassConstraint::ConstraintFitInput::getConstituentCovariancePhiThetaP
void getConstituentCovariancePhiThetaP(int i, AmgMatrix(3, 3)&outMatrix) const
Access to individual covariance PhiThetaP (3,3)
Definition: ConstraintFitInput.h:31
ConstraintFit.h
xAOD::Muon_v1
Class describing a Muon.
Definition: Muon_v1.h:38
AmgMatrix
#define AmgMatrix(rows, cols)
Definition: EventPrimitives.h:51
dqt_zlumi_pandas.mass
mass
Definition: dqt_zlumi_pandas.py:170
ZMassConstraint::ConstraintFit::massFitRun
void massFitRun(ConstraintFitOutput &output, double zresol=-1.)
Definition: ConstraintFit.cxx:584
ZMassConstraint::ConstraintFit::massFitInterface
int massFitInterface(const ConstraintFitInput &theInput)
Definition: ConstraintFit.cxx:384
PATCore::ParticleType::ConvertedPhoton
@ ConvertedPhoton
Definition: PATCoreEnums.h:39
ZMassConstraint::ConstraintFit::convertCovXYZTod0z0PhiThetaP
void convertCovXYZTod0z0PhiThetaP(const std::vector< TLorentzVector > &particleList, const Amg::MatrixX &covarXYZ, Amg::MatrixX &covard0z0PhiThetaP) const
Definition: ConstraintFit.cxx:1158
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:59
ZMassConstraint::ConstraintFit::m_nobj
unsigned int m_nobj
Definition: ConstraintFit.h:121
xAOD::EgammaHelpers::isConvertedPhoton
bool isConvertedPhoton(const xAOD::Egamma *eg, bool excludeTRT=false)
is the object a converted photon
Definition: EgammaxAODHelpers.cxx:26
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
CheckAppliedSFs.e3
e3
Definition: CheckAppliedSFs.py:264
ZMassConstraint::ConstraintFit::getMassError
double getMassError(const ConstraintFitInput &firstInput, const ConstraintFitInput &secondInput=ConstraintFitInput())
Calculate the mass error without fit - use just the inputs.
Definition: ConstraintFit.cxx:905
ZMassConstraint::ConstraintFit::m_theInput
ConstraintFitInput m_theInput
Definition: ConstraintFit.h:119
ZMassConstraint::ConstraintFit::m_covarianceInit
Amg::MatrixX m_covarianceInit
Definition: ConstraintFit.h:124
xAOD::CaloCluster_v1::eta
virtual double eta() const
The pseudorapidity ( ) of the particle.
Definition: CaloCluster_v1.cxx:251
lumiFormat.i
int i
Definition: lumiFormat.py:92
ZMassConstraint::ConstraintFit::convertCovd0z0PhiThetaPToXYZ
void convertCovd0z0PhiThetaPToXYZ(const TLorentzVector &fourVec, const AmgMatrix(5, 5)&covard0z0PhiThetaP, AmgMatrix(5, 5)&covarXYZ) const
Definition: ConstraintFit.cxx:1127
ZMassConstraint::ConstraintFit::d0
@ d0
Definition: ConstraintFit.h:83
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
vector
Definition: MultiHisto.h:13
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
ZMassConstraint::ConstraintFit::likelihoodMass
double likelihoodMass(double)
Definition: ConstraintFit.cxx:434
PlotPulseshapeFromCool.input
input
Definition: PlotPulseshapeFromCool.py:106
plotBeamSpotCompare.ivar
int ivar
Definition: plotBeamSpotCompare.py:383
ZMassConstraint::ConstraintFitInput
Definition: ConstraintFitInput.h:12
ZMassConstraint::ConstraintFit::retrieve_eta_calo
float retrieve_eta_calo(const xAOD::IParticle &part) const
Definition: ConstraintFit.cxx:1189
python.BuildSignatureFlags.sig
sig
Definition: BuildSignatureFlags.py:215
PhotonxAODHelpers.h
ZMassConstraint
Definition: ConstraintFit.cxx:12
plotIsoValidation.el
el
Definition: plotIsoValidation.py:197
ZMassConstraint::ConstraintFit::m_conHasWidth
bool m_conHasWidth
Definition: ConstraintFit.h:109
xAOD::CaloCluster_v1::retrieveMoment
bool retrieveMoment(MomentType type, double &value) const
Retrieve individual moment.
Definition: CaloCluster_v1.cxx:738
beamspotnt.rows
list rows
Definition: bin/beamspotnt.py:1112
ZMassConstraint::ConstraintFitOutput::getConstituentFourVector
const TLorentzVector & getConstituentFourVector(int index) const
Access to individual particle 4-vec.
Definition: ConstraintFitOutput.h:24
ZMassConstraint::ConstraintFit::massFit
double massFit(const Amg::MatrixX &, const Amg::MatrixX &var, double mass, Amg::MatrixX &pOut, Amg::MatrixX &)
Definition: ConstraintFit.cxx:856
Trk::Combined
@ Combined
Definition: TrackSummaryTool.h:32
merge.output
output
Definition: merge.py:17
ZMassConstraint::ConstraintFit::addFSRParticle
void addFSRParticle(const xAOD::IParticle &part, const TLorentzVector &fsr4Vec, ConstraintFitInput &input)
Add in FSR photon to input, (energy resolution is obtain in method)
Definition: ConstraintFit.cxx:83
ConstraintFitOutput.h
ZMassConstraint::ConstraintFit::m_conWidth
double m_conWidth
Definition: ConstraintFit.h:111
ZMassConstraint::ConstraintFit::m_parameters
unsigned int m_parameters
Definition: ConstraintFit.h:120
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
ZMassConstraint::ConstraintFit::phi0
@ phi0
Definition: ConstraintFit.h:85
xAOD::Electron_v1
Definition: Electron_v1.h:34
ZMassConstraint::ConstraintFit::m_ignoreInputChecks
bool m_ignoreInputChecks
Definition: ConstraintFit.h:113
ZMassConstraint::ConstraintFitInput::isOK
bool isOK()
Check that covariance matrix is OK.
Definition: ConstraintFitInput.h:46
xAOD::IParticle::isAvailable
bool isAvailable(const std::string &name, const std::string &clsname="") const
Check if a user property is available for reading or not.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:131
xAOD::Photon
Photon_v1 Photon
Definition of the current "egamma version".
Definition: Event/xAOD/xAODEgamma/xAODEgamma/Photon.h:17
TrackParticle.h
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
ZMassConstraint::ConstraintFit::m_resolution
double m_resolution
Definition: ConstraintFit.h:112
ZMassConstraint::ConstraintFitOutput::getNConstituents
unsigned int getNConstituents() const
Number of particles.
Definition: ConstraintFitOutput.h:21
xAOD::Photon_v1
Definition: Photon_v1.h:37
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
ZMassConstraint::ConstraintFit::doMassFit
StatusCode doMassFit(const ConstraintFitInput &input, ConstraintFitOutput &output)
Perform the constrained mass fit.
Definition: ConstraintFit.cxx:54
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
LArG4FSStartPointFilter.particles
list particles
Definition: LArG4FSStartPointFilter.py:84
xAOD::IParticle::auxdata
T & auxdata(const std::string &name, const std::string &clsname="")
Fetch an aux data variable, as a non-const reference.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:96
ZMassConstraint::ConstraintFit::initialize
StatusCode initialize()
Initialize constraint fit.
Definition: ConstraintFit.cxx:34
ZMassConstraint::ConstraintFit::m_mu_resolSFTool
ToolHandle< CP::IMuonCalibrationAndSmearingTool > m_mu_resolSFTool
Definition: ConstraintFit.h:117
ZMassConstraint::ConstraintFit::m_energyRescaler
ToolHandle< CP::IEgammaCalibrationAndSmearingTool > m_energyRescaler
Definition: ConstraintFit.h:116
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
ZMassConstraint::ConstraintFit::addParticle
void addParticle(const xAOD::Muon &part, ConstraintFitInput &input, MassConstraintMuonType muonType=isCombMCMT)
Add muon to input, must provide the resolution Scale Factor.
Definition: ConstraintFit.cxx:145
PATCore::ParticleType::Electron
@ Electron
Definition: PATCoreEnums.h:40
xAOD::track
@ track
Definition: TrackingPrimitives.h:512
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
ZMassConstraint::ConstraintFit::likelihoodMass2
double likelihoodMass2(void)
Definition: ConstraintFit.cxx:502
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
xAOD::CaloCluster_v1::ETACALOFRAME
@ ETACALOFRAME
Eta in the calo frame (for egamma)
Definition: CaloCluster_v1.h:187
ZMassConstraint::ConstraintFit::~ConstraintFit
~ConstraintFit(void)
Definition: ConstraintFit.cxx:30
CaloNoise_fillDB.mu
mu
Definition: CaloNoise_fillDB.py:53
ZMassConstraint::ConstraintFit::m_covarianceFinal
Amg::MatrixX m_covarianceFinal
Definition: ConstraintFit.h:126
test_AnalysisBaseEventLoopJob.aa
aa
Definition: test_AnalysisBaseEventLoopJob.py:37
ZMassConstraint::ConstraintFit::qOverP
@ qOverP
Definition: ConstraintFit.h:87
Generated on Sun Jun 30 2024 21:12:56 for ATLAS Offline Software by   doxygen 1.8.18