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