ATLAS Offline Software
Loading...
Searching...
No Matches
ConstraintFit.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
5#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 {
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
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 static const SG::ConstAccessor<float> etaCaloAcc("etaCalo");
1199 if (cluster.retrieveMoment(xAOD::CaloCluster::ETACALOFRAME,
1200 eta_calo)) { }
1201 else if (etaCaloAcc.isAvailable(cluster)) {
1202 eta_calo = etaCaloAcc(cluster);
1203 }
1204 else {
1205 ATH_MSG_ERROR("retrieve_eta_calo - etaCalo not available as auxilliary variable");
1206 ATH_MSG_WARNING("retrieve_eta_calo - using eta as etaCalo");
1207 eta_calo = cluster.eta();
1208 }
1209 return eta_calo;
1210 }
1211}
1212
1213
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::AuxElement::ConstAccessor
SG::ConstAccessor< T, ALLOC > ConstAccessor
Definition AuxElement.h:569
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:54
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
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::massFitInterface
int massFitInterface(const ConstraintFitInput &theInput)
Definition ConstraintFit.cxx:384
ZMassConstraint::ConstraintFit::massFitRun
void massFitRun(ConstraintFitOutput &output, double zresol=-1.)
Definition ConstraintFit.cxx:584
ZMassConstraint::ConstraintFit::convertCovd0z0PhiThetaPToXYZ
void convertCovd0z0PhiThetaPToXYZ(const TLorentzVector &fourVec, const AmgMatrix(5, 5)&covard0z0PhiThetaP, AmgMatrix(5, 5)&covarXYZ) const
Definition ConstraintFit.cxx:1127
ZMassConstraint::ConstraintFit::retrieve_eta_calo
float retrieve_eta_calo(const xAOD::IParticle &part) const
Definition ConstraintFit.cxx:1189
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:34
ZMassConstraint::ConstraintFit::doSanityChecksOnCovariance
bool doSanityChecksOnCovariance(const TLorentzVector &vector, const AmgMatrix(5, 5)&covar) const
Definition ConstraintFit.cxx:1090
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
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:786
ZMassConstraint::ConstraintFit::convertCovXYZTod0z0PhiThetaP
void convertCovXYZTod0z0PhiThetaP(const std::vector< TLorentzVector > &particleList, const Amg::MatrixX &covarXYZ, Amg::MatrixX &covard0z0PhiThetaP) const
Definition ConstraintFit.cxx:1158
ZMassConstraint::ConstraintFit::massFit
double massFit(const Amg::MatrixX &, const Amg::MatrixX &var, double mass, Amg::MatrixX &pOut, Amg::MatrixX &)
Definition ConstraintFit.cxx:856
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::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.14.0