Belle II Software development
ParticleKinematicFitterModule.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8#include <sstream>
9#define SSTR( x ) dynamic_cast< std::ostringstream && >(( std::ostringstream() << std::dec << x ) ).str()
10
11// kinfitter
12#include <analysis/modules/ParticleKinematicFitter/ParticleKinematicFitterModule.h>
13#include <analysis/OrcaKinFit/BaseFitObject.h>
14#include <analysis/OrcaKinFit/BaseFitter.h>
15#include <analysis/OrcaKinFit/OPALFitterGSL.h>
16#include <analysis/OrcaKinFit/JetFitObject.h>
17#include <analysis/OrcaKinFit/NewtonFitterGSL.h>
18#include <analysis/OrcaKinFit/NewFitterGSL.h>
19#include <analysis/OrcaKinFit/ParticleFitObject.h>
20#include <analysis/OrcaKinFit/PxPyPzMFitObject.h>
21#include <analysis/OrcaKinFit/TextTracer.h>
22
23#include <mdst/dataobjects/ECLCluster.h>
24
25#include <TDatabasePDG.h>
26
27// framework utilities
28#include <framework/gearbox/Const.h>
29#include <framework/logging/Logger.h>
30
31// analysis dataobjects
32#include <analysis/dataobjects/Particle.h>
33
34// analysis utilities
35#include <analysis/utility/PCmsLabTransform.h>
36#include <analysis/utility/ParticleCopy.h>
37#include <analysis/ClusterUtility/ClusterUtils.h>
38
39#include <cmath>
40
41using namespace CLHEP;
42using namespace std;
43using namespace Belle2;
44using namespace Belle2::OrcaKinFit;
45
46//-----------------------------------------------------------------
47// Register module
48//-----------------------------------------------------------------
49
50REG_MODULE(ParticleKinematicFitter);
51
52//-----------------------------------------------------------------
53// Implementation
54//-----------------------------------------------------------------
55
57 DataStore::c_Event)
58{
59 setDescription("Kinematic fitter for modular analysis");
61
62 // Add parameters
63 addParam("listName", m_listName, "Name of particle list.", string(""));
64 addParam("kinematicFitter", m_kinematicFitter, "Available: OrcaKinFit.", string("OrcaKinFit"));
65 addParam("orcaFitterEngine", m_orcaFitterEngine, "OrcaKinFit engine: NewFitterGSL, NewtonFitterGSL, OPALFitterGSL.",
66 string("OPALFitterGSL"));
67 addParam("orcaTracer", m_orcaTracer, "OrcaKinFit tracer: None, Text.", string("None"));
68 addParam("orcaConstraint", m_orcaConstraint, "OrcaKinFit constraint: HardBeam, RecoilMass.", string("HardBeam"));
69 addParam("debugFitter", m_debugFitter, "Switch on/off internal debugging output if available.", false);
70 addParam("debugFitterLevel", m_debugFitterLevel, "Internal debugging output level if available.", 10);
71 addParam("addUnmeasuredPhoton", m_addUnmeasuredPhoton, "Add one unmeasured photon (-3C).", false);
72 addParam("fixUnmeasuredToHER", m_fixUnmeasuredPhotonToHER, "fix the momentum of the unmeasured photon to HER (+2C).", false);
73 addParam("fixUnmeasuredToLER", m_fixUnmeasuredPhotonToLER, "fix the momentum of the unmeasured photon to LER (+2C).", false);
74 addParam("add3CPhoton", m_add3CPhoton, "Add one photon with unmeasured energy (-1C).", false);
75 addParam("liftPhotonTheta", m_liftPhotonTheta, "Lift theta constraint of 3CPhoton. Valid when add3CPhoton is true.", false);
76 addParam("decayString", m_decayString, "Specifies which daughter particles are included in the kinematic fit.", string(""));
77 addParam("updateMother", m_updateMother, "Update the mother kinematics.", true);
78 addParam("updateDaughters", m_updateDaughters, "Update the daughter kinematics.", false);
79 addParam("recoilMass", m_recoilMass, "Recoil mass in GeV. RecoilMass constraint only.", 0.0);
80 addParam("invMass", m_invMass, "Invariant mass in GeV. Mass constraint only.", 0.0);
81 addParam("variablePrefix", m_prefix, "Prefix attached to extra info variables.", string(""));
82 addParam("decayStringForDirectionOnlyParticles", m_decayStringForDirectionOnlyParticles,
83 "DecayString specifying the particles to use only direction information in the fit", std::string(""));
84 addParam("decayStringForAlternateMassParticles", m_decayStringForAlternateMassParticles,
85 "DecayString specifying the particles where an alternate mass hypothesis is used", std::string(""));
86 addParam("decayStringForNeutronVsAntiNeutron", m_decayStringForNeutronVsAntiNeutron,
87 "DecayString specifying the charged particle used to tag whether n or nbar. If tag particle has negative charge, PDG sign of n/nbar is flipped from default given in alternateMassHypos",
88 std::string(""));
89 addParam("alternateMassHypos", m_listAlternateMassHypo,
90 "integer list of pdg values for particles where different mass hypothesis is used in the fit");
91
92}
93
95{
96 m_eventextrainfo.isRequired();
97
98 if (m_decayString != "") {
100 B2INFO("ParticleKinematicFitter: Using specified decay string: " << m_decayString);
101 }
102
105 }
106
109 if (m_listAlternateMassHypo.size() == 0)
110 B2FATAL("decayStringForAlternateMassParticles specified but alternateMassHypos not provided.");
111 } else if (m_listAlternateMassHypo.size() > 0) {
112 B2FATAL("alternateMassHypos specified but decayStringForAlternateMassParticles not provided.");
113 }
114
117 }
118
119 m_plist.isRequired(m_listName);
120}
121
123{
124 B2INFO("ParticleKinematicFitterModule::terminate");
125}
126
128{
129 B2DEBUG(17, "ParticleKinematicFitterModule::event");
130
131 unsigned int n = m_plist->getListSize();
132
133 for (unsigned i = 0; i < n; i++) {
134 Particle* particle = m_plist->getParticle(i);
135
136 if (m_updateDaughters == true) {
137 if (m_decayString.empty()) ParticleCopy::copyDaughters(particle);
138 else B2ERROR("Daughters update works only when all daughters are selected. Daughters will not be updated");
139 }
140
141 bool ok = doKinematicFit(particle);
142
143 if (!ok) particle->setPValue(-1.);
144 }
145}
146
148{
149 B2DEBUG(17, "ParticleKinematicFitterModule::doKinematicFit");
150
151 bool ok = false;
152
153 // fitting with OrcaKinFit
154 if (m_kinematicFitter == "OrcaKinFit") {
155
156 // select subset of particles for the fit
157 if (m_decayString != "") {
158 B2FATAL("ParticleKinematicFitterModule: OrcaKinFit does not support yet selection of daughters via decay string!") ;
159 }
160
161 // check requested fit engine
162 if (!(m_orcaFitterEngine == "OPALFitterGSL" or
163 m_orcaFitterEngine == "NewtonFitterGSL" or
164 m_orcaFitterEngine == "NewFitterGSL")) {
165 B2FATAL("ParticleKinematicFitterModule: " << m_orcaFitterEngine << " is an invalid OrcaKinFit fitter engine!");
166 }
167
168 // check requested constraint
169 if (!(m_orcaConstraint == "HardBeam" or
170 m_orcaConstraint == "RecoilMass" or
171 m_orcaConstraint == "Mass")) {
172 B2FATAL("ParticleKinematicFitterModule: " << m_orcaConstraint << " is an invalid OrcaKinFit constraint!");
173 }
174
175 // basic check is good, go to fitting routine
176 ok = doOrcaKinFitFit(mother);
177 } else { // invalid fitter
178 B2FATAL("ParticleKinematicFitter: " << m_kinematicFitter << " is an invalid kinematic fitter!");
179 }
180
181 if (!ok) return false;
182
183 return true;
184
185}
186
188{
189 if (mother->getNDaughters() <= 1) {
190 B2WARNING("ParticleKinematicFitterModule: Cannot fit with " << mother->getNDaughters() << " daughters.");
191 return false;
192 }
193
194 //use getMdstSource as unique identifier
195 std::vector<int> directionOnlyParticlesIdentifiers;
197 std::vector<const Particle*> selparticles = m_decaydescriptorForDirectionOnlyParticles.getSelectionParticles(mother);
198 for (unsigned int i = 0; i < selparticles.size(); i++) {
199 directionOnlyParticlesIdentifiers.push_back(selparticles[i]->getMdstSource());
200 }
201 }
202
203 //use getMdstSource as unique identifier
204 std::vector<int> alternateMassHypoIdentifiers;
206 std::vector<const Particle*> selparticles = m_decaydescriptorForAlternateMassParticles.getSelectionParticles(mother);
207 for (unsigned int i = 0; i < selparticles.size(); i++) {
208 alternateMassHypoIdentifiers.push_back(selparticles[i]->getMdstSource());
209 }
210 if (alternateMassHypoIdentifiers.size() != m_listAlternateMassHypo.size()) {
211 B2FATAL("alternateMassHypoIdentifiers size must match listAlternateMassHypo");
212 }
213 }
214
215 //Determine if n or nbar based on charge of tag particle indicated
216 //If tag particle has negative charge, sign of n and nbar are flipped.
217 //Required to use optimal position resolution
218 bool flipNeutronPDGsign = 0;
220 std::vector<const Particle*> selparticles = m_decaydescriptorForNeutronVsAntiNeutron.getSelectionParticles(mother);
221 if (selparticles.size() != 1) {
222 B2FATAL("Select only one particle to tag neutron vs. antineutron");
223 }
224 if (selparticles[0]->getCharge() < 0) flipNeutronPDGsign = 1;
225 }
226
227
228 // fill particles
229 std::vector<Particle*> particleChildren;
230 bool validChildren = fillFitParticles(mother, particleChildren);
231
232 if (!validChildren) {
233 B2WARNING("ParticleKinematicFitterModule: Cannot find valid children for the fit.");
234 return false;
235 }
236
237 // set fit engine
238 BaseFitter* pfitter;
239
240 // internal debugger
241 int debugfitter = 0;
242 if (m_debugFitter) debugfitter = m_debugFitterLevel;
243
244 // choose minimization
245 if (m_orcaFitterEngine == "OPALFitterGSL") {
246 pfitter = new OPALFitterGSL(); // OPAL fitter has no debugger
247 } else if (m_orcaFitterEngine == "NewtonFitterGSL") {
248 pfitter = new NewtonFitterGSL();
249 (dynamic_cast<NewtonFitterGSL*>(pfitter))->setDebug(debugfitter);
250 } else if (m_orcaFitterEngine == "NewFitterGSL") {
251 pfitter = new NewFitterGSL();
252 (dynamic_cast<NewFitterGSL*>(pfitter))->setDebug(debugfitter);
253 } else {
254 B2FATAL("ParticleKinematicFitterModule: " << m_orcaFitterEngine << " is an invalid OrcaKinFit fitter engine!");
255 return false;
256 }
257
258 if (!pfitter) return false;
259 BaseFitter& fitter(*pfitter);
260
261 // reset fitter
262 resetFitter(fitter);
263
264 // set constraints (not connected to a fitter or particles at this point!)
266
267
268 // add fit particles from particle list to the fitter and to all constraints
269 for (unsigned iChild = 0; iChild < particleChildren.size(); iChild++) {
270
271 bool useDirectionOnly = false;
272 if (directionOnlyParticlesIdentifiers.size() > 0) {
273 if (std::find(directionOnlyParticlesIdentifiers.begin(), directionOnlyParticlesIdentifiers.end(),
274 particleChildren[iChild]->getMdstSource()) != directionOnlyParticlesIdentifiers.end()) useDirectionOnly = true;
275 }
276
277 int massHypo = 0;
278 if (alternateMassHypoIdentifiers.size() > 0) {
279 for (unsigned int i = 0; i < alternateMassHypoIdentifiers.size(); i++) {
280 if (alternateMassHypoIdentifiers[i] == particleChildren[iChild]->getMdstSource()) {
281 massHypo = m_listAlternateMassHypo[i];
282 break;
283 }
284 }
285 //Always use direction only for neutrons
286 if (abs(massHypo) == Const::neutron.getPDGCode()) {
287 if (useDirectionOnly == false)
288 B2WARNING("Neutron mass hypothesis assigned to fit particle but directionOnly flag not specified for same particle. Setting candidate to useDirectionOnly.");
289 useDirectionOnly = true;
290 if (flipNeutronPDGsign) massHypo = -massHypo;
291 }
292 }
293 addParticleToOrcaKinFit(fitter, particleChildren[iChild], iChild, useDirectionOnly, massHypo);
294 }
295
296 // add unmeasured photon to the fitter and to all constraints
298
299 // add constraints to the fitter
301
302 // add tracers to the fitter
303 addTracerToFitter(fitter);
304
305 //store information before the fit
306 storeOrcaKinFitParticles("Measured", fitter, particleChildren, mother);
307
308 double prob = fitter.fit();
309 double chi2 = fitter.getChi2();
310 int niter = fitter.getIterations();
311 int ndof = fitter.getDoF();
312 int errorcode = fitter.getError();
313
314 B2DEBUG(17, "ParticleKinematicFitterModule: -------------------------------------------");
315 B2DEBUG(17, "ParticleKinematicFitterModule: Fit result of OrcaKinFit using " << m_orcaFitterEngine);
316 B2DEBUG(17, "ParticleKinematicFitterModule: prob " << prob);
317 B2DEBUG(17, "ParticleKinematicFitterModule: chi2 " << chi2);
318 B2DEBUG(17, "ParticleKinematicFitterModule: iterations " << niter);
319 B2DEBUG(17, "ParticleKinematicFitterModule: ndf " << ndof);
320 B2DEBUG(17, "ParticleKinematicFitterModule: errorcode " << errorcode);
321 B2DEBUG(17, "ParticleKinematicFitterModule: -------------------------------------------");
322
323 // default update mother information
324 if (m_updateMother) updateOrcaKinFitMother(fitter, particleChildren, mother);
325
326 // update daughter information if that is requested
327 if (m_updateDaughters) updateOrcaKinFitDaughters(fitter, mother);
328
329 // store information after the fit
330 storeOrcaKinFitParticles("Fitted", fitter, particleChildren, mother);
331
332 //store general fit results
333 mother->addExtraInfo(m_prefix + "OrcaKinFitProb", prob);
334 mother->setPValue(prob);
335 mother->addExtraInfo(m_prefix + "OrcaKinFitChi2", chi2);
336 mother->addExtraInfo(m_prefix + "OrcaKinFitErrorCode", errorcode);
337
338 // if we added an unmeasured photon, add the kinematics to the mother - at some point we may want to create a particle list from this?
339 std::vector <BaseFitObject*>* fitObjectContainer = fitter.getFitObjects();
340 for (auto fo : *fitObjectContainer) {
342 const std::string name = fo->getName();
343 if (name.find("Unmeasured") != std::string::npos) {
344 auto* fitobject = static_cast<ParticleFitObject*>(fo);
345 ROOT::Math::PxPyPzEVector tlv = getLorentzVector(fitobject);
346 mother->addExtraInfo(m_prefix + "OrcaKinFit" + name + "Theta", tlv.Theta());
347 mother->addExtraInfo(m_prefix + "OrcaKinFit" + name + "Phi", tlv.Phi());
348 mother->addExtraInfo(m_prefix + "OrcaKinFit" + name + "E", tlv.E());
349
350 // Uncertainty
351 // const double err0 = getFitObjectError(fitobject, 0);
352 mother->addExtraInfo(m_prefix + "OrcaKinFit" + name + "ErrorTheta", getFitObjectError(fitobject, 1));
353 mother->addExtraInfo(m_prefix + "OrcaKinFit" + name + "ErrorPhi", getFitObjectError(fitobject, 2));
354 mother->addExtraInfo(m_prefix + "OrcaKinFit" + name + "ErrorE", getFitObjectError(fitobject, 0));
355 }
356 }
357 delete fo;
358 }
359
360 delete pfitter;
361 delete m_textTracer;
362 return true;
363}
364
365bool ParticleKinematicFitterModule::fillFitParticles(Particle* mother, std::vector<Particle*>& particleChildren)
366{
367 for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
368 auto* child = const_cast<Particle*>(mother->getDaughter(ichild));
369
370 if (child->getNDaughters() > 0) {
371 bool err = fillFitParticles(child, particleChildren);
372 if (!err) {
373 B2WARNING("ParticleKinematicFitterModule: Cannot find valid children for the fit.");
374 return false;
375 }
376 } else if (child->getPValue() > 0) {
377 particleChildren.push_back(child);
378 } else {
379 B2ERROR("Daughter with PDG code " << child->getPDGCode() << " does not have a valid p-value: p=" << child->getPValue() << ", E=" <<
380 child->getEnergy() << " GeV");
381 return false; // error matrix not valid
382 }
383 }
384 return true;
385}
386
388{
389 TMatrixFSym MomentumVertexErrorMatrix(7);
390 for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
391 auto* child = const_cast<Particle*>(mother->getDaughter(ichild));
392
393 if (child->getPValue() > 0) {
394 MomentumVertexErrorMatrix += child->getMomentumVertexErrorMatrix();
395 } else if (child->getNDaughters() > 0) {
396 AddFour(child);
397 MomentumVertexErrorMatrix += child->getMomentumVertexErrorMatrix();
398 } else {
399 B2ERROR("Daughter with PDG code " << child->getPDGCode() << " does not have a valid p-value: p=" << child->getPValue() << ", E=" <<
400 child->getEnergy() << " GeV");
401 return false; // error matrix not valid
402 }
403 }
404 mother->setMomentumVertexErrorMatrix(MomentumVertexErrorMatrix);
405 return true;
406}
407
409 bool useOnlyDirection, int massHypoPDG)
410{
411 B2DEBUG(17, "ParticleKinematicFitterModule: adding a particle to the fitter!");
412
413 if (m_add3CPhoton && index == 0) {
414 if (particle -> getPDGCode() != Const::photon.getPDGCode()) {
415 B2ERROR("In 3C Kinematic fit, the first daughter should be the Unmeasured Photon!");
416 }
417
418 const ECLCluster* cluster = particle->getECLCluster();
419 double startingE = cluster->getEnergy(particle->getECLClusterEHypothesisBit());
420 double startingPhi = cluster->getPhi();
421 double startingTheta = cluster->getTheta();
422
423 ClusterUtils clutls;
424 const auto EPhiThetaCov = clutls.GetCovarianceMatrix3x3FromCluster(cluster);
425 double startingEError = sqrt(fabs(EPhiThetaCov[0][0]));
426 double startingPhiError = sqrt(fabs(EPhiThetaCov[1][1]));
427 double startingThetaError = sqrt(fabs(EPhiThetaCov[2][2]));
428
429 B2DEBUG(17, startingPhi << " " << startingTheta << " " << startingPhiError << " " << startingThetaError);
430 // create a fit object
431 ParticleFitObject* pfitobject;
432 // memory allocated: it will be deallocated via "delete fo" in doOrcaKinFitFit
433 pfitobject = new JetFitObject(startingE, startingTheta, startingPhi, startingEError, startingThetaError, startingPhiError, 0.);
434 pfitobject->setParam(0, startingE, false, false);
436 pfitobject->setParam(1, startingTheta, false, false);
437 else
438 pfitobject->setParam(1, startingTheta, true, false);
439 pfitobject->setParam(2, startingPhi, true, false);
440
441 std::string fitObjectName = "Unmeasured3C";
442 pfitobject->setName(fitObjectName.c_str());
443 ParticleFitObject& fitobject = *pfitobject;
444
445 // add this fit object (=particle) to the constraints
446 addFitObjectToConstraints(fitobject);
447
448 // add fit particle to the fitter
449 fitter.addFitObject(fitobject);
450
451 } else if (useOnlyDirection || massHypoPDG) {
452
453 ParticleFitObject* pfitobject;
454
455 const ECLCluster* cluster = particle->getECLCluster();
456
457 if (particle->getCharge() != 0 or !cluster) {
458 B2FATAL("ParticleKinematicFitterModule: Direction only and alternate mass options only implemented for neutral particles with ECL cluster");
459 }
460
461 double mass = particle->getPDGMass();
462 if (massHypoPDG != 0) {
463 if (TDatabasePDG::Instance()->GetParticle(massHypoPDG) == nullptr) {
464 B2FATAL("ParticleKinematicFitterModule: " << massHypoPDG << " is an unknown PDG code!");
465 }
466 mass = TDatabasePDG::Instance()->GetParticle(massHypoPDG)->Mass();
467 }
468 double clusterE = cluster->getEnergy(particle->getECLClusterEHypothesisBit());
469 double startingE = sqrt(clusterE * clusterE + mass * mass);
470 double startingPhi = cluster->getPhi();
471 double startingTheta = cluster->getTheta();
472
473 ClusterUtils clutls;
474 const auto EPhiThetaCov = clutls.GetCovarianceMatrix3x3FromCluster(cluster, massHypoPDG);
475 double startingEError = sqrt(fabs(EPhiThetaCov[0][0]));
476 double startingPhiError = sqrt(fabs(EPhiThetaCov[1][1]));
477 double startingThetaError = sqrt(fabs(EPhiThetaCov[2][2]));
478
479 // memory allocated: it will be deallocated via "delete fo" in doOrcaKinFitFit
480 pfitobject = new JetFitObject(startingE, startingTheta, startingPhi, startingEError, startingThetaError, startingPhiError, mass);
481 pfitobject->setParam(0, startingE, true, false);
482 if (useOnlyDirection) pfitobject->setParam(0, startingE, false, false);
483 pfitobject->setParam(1, startingTheta, true, false);
484 pfitobject->setParam(2, startingPhi, true, false);
485
486 std::string fitObjectName = "particle_" + SSTR(index);
487 pfitobject->setName(fitObjectName.c_str());
488 ParticleFitObject& fitobject = *pfitobject;
489
490 // add this fit object (=particle) to the constraints
491 addFitObjectToConstraints(fitobject);
492
493 // add fit particle to the fitter
494 fitter.addFitObject(fitobject);
495
496 } else {
497 // four vector
498 CLHEP::HepLorentzVector clheplorentzvector = getCLHEPLorentzVector(particle);
499
500 // error matrix
501 CLHEP::HepSymMatrix clhepmomentumerrormatrix = getCLHEPMomentumErrorMatrix(particle);
502
503 // create the fit object (ParticleFitObject is the base class)
504 ParticleFitObject* pfitobject;
505 // memory allocated: it will be deallocated via "delete fo" in doOrcaKinFitFit
506 pfitobject = new PxPyPzMFitObject(clheplorentzvector, clhepmomentumerrormatrix);
507 std::string fitObjectName = "particle_" + SSTR(index);
508 pfitobject->setName(fitObjectName.c_str());
509 ParticleFitObject& fitobject = *pfitobject;
510
511 // add this fit object (=particle) to the constraints
512 addFitObjectToConstraints(fitobject);
513
514 // add fit particle to the fitter
515 fitter.addFitObject(fitobject);
516 }
517
518 return;
519}
520
522{
523 CLHEP::HepSymMatrix covMatrix(4);
524 TMatrixFSym errMatrix = particle->getMomentumErrorMatrix();
525
526 for (int i = 0; i < 4; i++) {
527 for (int j = 0; j < 4; j++) {
528 covMatrix[i][j] = errMatrix[i][j];
529 }
530 }
531
532 return covMatrix;
533}
534
536{
537 CLHEP::HepSymMatrix covMatrix(7);
538 TMatrixFSym errMatrix = particle->getMomentumVertexErrorMatrix();
539
540 for (int i = 0; i < 7; i++) {
541 for (int j = 0; j < 7; j++) {
542 covMatrix[i][j] = errMatrix[i][j];
543 }
544 }
545
546 return covMatrix;
547}
548
550{
551 CLHEP::HepLorentzVector mom(particle->getPx(), particle->getPy(), particle->getPz(), particle->get4Vector().E());
552 return mom;
553}
554
556{
557 ROOT::Math::PxPyPzEVector mom(fitobject->getPx(), fitobject->getPy(), fitobject->getPz(), fitobject->getE());
558 return mom;
559}
560
561// TMatrixFSym ParticleKinematicFitterModule::getTMatrixFSymMomentumErrorMatrix(ParticleFitObject* fitobject) //fix the warning
563{
564 TMatrixFSym errMatrix(4);
565
566 for (int i = 0; i < 4; i++) {
567 for (int j = i; j < 4; j++) {
568 errMatrix[i][j] = 0.0;
569 }
570 }
571
572 return errMatrix;
573}
574
575// TMatrixFSym ParticleKinematicFitterModule::getTMatrixFSymMomentumVertexErrorMatrix(ParticleFitObject* fitobject) //fix the warning
577{
578 TMatrixFSym errMatrix(7);
579
580 for (int i = 0; i < 7; i++) {
581 for (int j = i; j < 7; j++) {
582 errMatrix[i][j] = 0.0;
583 }
584 }
585
586 return errMatrix;
587}
588
590{
591 if (m_orcaConstraint == "HardBeam") {
592 ROOT::Math::PxPyPzEVector constraints4vector(m_hardConstraintPx.getValue(),
593 m_hardConstraintPy.getValue(),
594 m_hardConstraintPz.getValue(),
595 m_hardConstraintE.getValue());
596 return constraints4vector;
597 } else {
598 B2FATAL("ParticleKinematicFitterModule: " << m_orcaConstraint << " is an invalid OrcaKinFit constraint!");
599 }
600
601 // should not reach this point...
602 return ROOT::Math::PxPyPzEVector(0., 0., 0., 0.);
603}
604
606{
607 if (m_orcaConstraint == "HardBeam") {
609 const ROOT::Math::PxPyPzEVector boost = T.getBeamFourMomentum();
610
611 m_hardConstraintPx = MomentumConstraint(0, 1, 0, 0, boost.Px());
612 m_hardConstraintPy = MomentumConstraint(0, 0, 1, 0, boost.Py());
613 m_hardConstraintPz = MomentumConstraint(0, 0, 0, 1, boost.Pz());
614 m_hardConstraintE = MomentumConstraint(1, 0, 0, 0, boost.E());
615
616 m_hardConstraintPx.resetFOList();
617 m_hardConstraintPy.resetFOList();
618 m_hardConstraintPz.resetFOList();
619 m_hardConstraintE.resetFOList();
620
621 m_hardConstraintPx.setName("Sum(p_x) [hard]");
622 m_hardConstraintPy.setName("Sum(p_y) [hard]");
623 m_hardConstraintPz.setName("Sum(p_z) [hard]");
624 m_hardConstraintE.setName("Sum(E) [hard]");
625
626 } else if (m_orcaConstraint == "RecoilMass") {
628 const ROOT::Math::PxPyPzEVector boost = T.getBeamFourMomentum();
629
630 m_hardConstraintRecoilMass = RecoilMassConstraint(m_recoilMass, boost.Px(), boost.Py(), boost.Pz(), boost.E());
631
632 m_hardConstraintRecoilMass.resetFOList();
633 m_hardConstraintRecoilMass.setName("Recoil Mass [hard]");
634
635 } else if (m_orcaConstraint == "Mass") {
637
638 m_hardConstraintMass.resetFOList();
639 m_hardConstraintMass.setName("Mass [hard]");
640 } else {
641 B2FATAL("ParticleKinematicFitterModule: " << m_orcaConstraint << " is an invalid OrcaKinFit constraint!");
642 }
643}
644
646{
647 B2DEBUG(17, "ParticleKinematicFitterModule: Resetting the fitter");
648 fitter.reset();
649}
650
652{
653 if (m_orcaConstraint == "HardBeam") {
654 m_hardConstraintPx.addToFOList(fitobject);
655 m_hardConstraintPy.addToFOList(fitobject);
656 m_hardConstraintPz.addToFOList(fitobject);
657 m_hardConstraintE.addToFOList(fitobject);
658 } else if (m_orcaConstraint == "RecoilMass") {
659 m_hardConstraintRecoilMass.addToFOList(fitobject);
660 } else if (m_orcaConstraint == "Mass") {
661 m_hardConstraintMass.addToFOList(fitobject);
662 } else {
663 B2FATAL("ParticleKinematicFitterModule: " << m_orcaConstraint << " is an invalid OrcaKinFit constraint!");
664 }
665}
666
668{
669 if (m_orcaConstraint == "HardBeam") {
670 fitter.addConstraint(m_hardConstraintPx);
671 fitter.addConstraint(m_hardConstraintPy);
672 fitter.addConstraint(m_hardConstraintPz);
673 fitter.addConstraint(m_hardConstraintE);
674 } else if (m_orcaConstraint == "RecoilMass") {
675 fitter.addConstraint(m_hardConstraintRecoilMass);
676 } else if (m_orcaConstraint == "Mass") {
677 fitter.addConstraint(m_hardConstraintMass);
678 }
679
680 else {
681 B2FATAL("ParticleKinematicFitterModule: " << m_orcaConstraint << " is an invalid OrcaKinFit constraint!");
682 }
683}
684
686{
687 if (m_orcaTracer == "Text") {
688 m_textTracer = new TextTracer(std::cout);
689 fitter.setTracer(m_textTracer);
690 } else if (m_orcaTracer != "None") {
691 B2FATAL("ParticleKinematicFitterModule: " << m_orcaTracer << " is an invalid OrcaKinFit tracer!");
692 }
693}
694
696{
697 B2DEBUG(17, "ParticleKinematicFitterModule::addUnmeasuredGammaToOrcaKinFit: adding an unmeasured photon to the fitter!");
698 // Initialize photon using the existing constraints
699 ROOT::Math::PxPyPzEVector tlv = getLorentzVectorConstraints();
700 double startingE = tlv.E();
701 double startingPhi = tlv.Phi();
702 double startingTheta = tlv.Theta();
703 bool paramFlag = false;
704
705 // create a fit object
706 ParticleFitObject* pfitobject;
707
708 std::string fitObjectName = "UnmeasuredAlongBeam";
709
711 startingTheta = 41.5e-3; // TODO: Read beam crossing angle from db if it's available
712 startingPhi = 0.0;
713 paramFlag = true;
714 } else if (m_fixUnmeasuredPhotonToLER) {
715 startingTheta = TMath::Pi() - 41.5e-3;
716 startingPhi = 0.0;
717 paramFlag = true;
718 } else {
719 fitObjectName = "Unmeasured";
720 }
721
722 // memory allocated: it will be deallocated via "delete fo" in doOrcaKinFitFit
723 pfitobject = new JetFitObject(startingE, startingTheta, startingPhi, 0.0, 0.0, 0.0, 0.);
724 pfitobject->setParam(0, startingE, false, false);
725 pfitobject->setParam(1, startingTheta, paramFlag, paramFlag);
726 pfitobject->setParam(2, startingPhi, paramFlag, paramFlag);
727
728 pfitobject->setName(fitObjectName.c_str());
729 ParticleFitObject& fitobject = *pfitobject;
730
731 // add this fit object (=particle) to the constraints
732 addFitObjectToConstraints(fitobject);
733
734 // add fit particle to the fitter
735 fitter.addFitObject(fitobject);
736}
737
739{
740 std::vector <Particle*> bDau = mother->getDaughters();
741 std::vector <BaseFitObject*>* fitObjectContainer = fitter.getFitObjects();
742
743 const unsigned nd = bDau.size();
744 unsigned l = 0;
745 std::vector<std::vector<unsigned>> pars;
746 std::vector<Particle*> allparticles;
747 for (unsigned ichild = 0; ichild < nd; ichild++) {
748 const Particle* daughter = mother->getDaughter(ichild);
749 std::vector<unsigned> pard;
750 if (daughter->getNDaughters() > 0) {
751 updateMapOfTrackAndDaughter(l, pars, pard, allparticles, daughter);
752 pars.push_back(pard);
753 allparticles.push_back(bDau[ichild]);
754 } else {
755 pard.push_back(l);
756 pars.push_back(pard);
757 allparticles.push_back(bDau[ichild]);
758 l++;
759 }
760 }
761
762 if (l == fitObjectContainer->size() - m_addUnmeasuredPhoton) {
763
764 if (fitter.getError() == 0) {
765 for (unsigned iDaug = 0; iDaug < allparticles.size(); iDaug++) {
766 ROOT::Math::PxPyPzEVector tlv ;
767 TMatrixFSym errMatrixU(7);
768 if (pars[iDaug].size() > 0) {
769 for (unsigned int iChild : pars[iDaug]) {
770 BaseFitObject* fo = fitObjectContainer->at(iChild);
771 auto* fitobject = static_cast<ParticleFitObject*>(fo);
772 ROOT::Math::PxPyPzEVector tlv_sub = getLorentzVector(fitobject);
773 TMatrixFSym errMatrixU_sub = getCovMat7(fitobject);
774 tlv = tlv + tlv_sub;
775 errMatrixU = errMatrixU + errMatrixU_sub;
776 }
777 } else {
778 B2FATAL("ParticleKinematicFitterModule: no fitObject could be used to update the daughter!");
779 }
780 ROOT::Math::XYZVector pos = allparticles[iDaug]->getVertex(); // we don't update the vertex yet
781 TMatrixFSym errMatrix = allparticles[iDaug]->getMomentumVertexErrorMatrix();
782 TMatrixFSym errMatrixMom = allparticles[iDaug]->getMomentumErrorMatrix();
783 TMatrixFSym errMatrixVer = allparticles[iDaug]->getVertexErrorMatrix();
784
785 for (int i = 0; i < 3; i++) {
786 for (int j = i; j < 3; j++) {
787 errMatrixU[i + 4][j + 4] = errMatrixVer[i][j];
788 }
789 }
790
791 allparticles[iDaug]->set4Vector(tlv);
792 allparticles[iDaug]->setVertex(pos);
793 allparticles[iDaug]->setMomentumVertexErrorMatrix(errMatrixU);
794 }
795 }
796
797 return true;
798 } else {
799 B2ERROR("updateOrcaKinFitDaughters: Cannot update daughters, mismatch between number of daughters and number of fitobjects!");
800 return false;
801 }
802}
803
804void ParticleKinematicFitterModule::updateMapOfTrackAndDaughter(unsigned& l, std::vector<std::vector<unsigned>>& pars,
805 std::vector<unsigned>& parm, std::vector<Particle*>& allparticles, const Particle* daughter)
806{
807 std::vector <Particle*> dDau = daughter->getDaughters();
808 for (unsigned ichild = 0; ichild < daughter->getNDaughters(); ichild++) {
809 const Particle* child = daughter->getDaughter(ichild);
810 std::vector<unsigned> pard;
811 if (child->getNDaughters() > 0) {
812 updateMapOfTrackAndDaughter(l, pars, pard, allparticles, child);
813 parm.insert(parm.end(), pard.begin(), pard.end());
814 pars.push_back(pard);
815 allparticles.push_back(dDau[ichild]);
816 } else {
817 pard.push_back(l);
818 parm.push_back(l);
819 pars.push_back(pard);
820 allparticles.push_back(dDau[ichild]);
821 l++;
822 }
823 }
824}
825
826void ParticleKinematicFitterModule::updateOrcaKinFitMother(BaseFitter& fitter, std::vector<Particle*>& particleChildren,
827 Particle* mother)
828{
829 // get old values
830 ROOT::Math::XYZVector pos = mother->getVertex();
831 TMatrixFSym errMatrix = mother->getMomentumVertexErrorMatrix();
832 float pvalue = mother->getPValue();
833
834 // update momentum vector
835 ROOT::Math::PxPyPzEVector momnew(0., 0., 0., 0.);
836
837 std::vector <BaseFitObject*>* fitObjectContainer = fitter.getFitObjects();
838 for (unsigned iChild = 0; iChild < particleChildren.size(); iChild++) {
839 BaseFitObject* fo = fitObjectContainer->at(iChild);
840 auto* fitobject = static_cast<ParticleFitObject*>(fo);
841 ROOT::Math::PxPyPzEVector tlv = getLorentzVector(fitobject);
842 momnew += tlv;
843 }
844
845 // update
846 // TODO: use pvalue of the fit or the old one of the mother? use fit covariance matrix?
847 // Maybe here should use the pvalue and errmatrix of the fit ----Yu Hu
848 mother->updateMomentum(momnew, pos, errMatrix, pvalue);
849}
850
852 std::vector<Particle*>& particleChildren, Particle* mother)
853{
854 bool updated = false;
855 std::vector <BaseFitObject*>* fitObjectContainer = fitter.getFitObjects();
856
857 for (unsigned iChild = 0; iChild < particleChildren.size(); iChild++) {
858 BaseFitObject* fo = fitObjectContainer->at(iChild);
859 auto* fitobject = static_cast<ParticleFitObject*>(fo);
860 ROOT::Math::PxPyPzEVector tlv = getLorentzVector(fitobject);
861
862 // name of extra variables
863 std::string extraVariableParticlePx = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_Px";
864 std::string extraVariableParticlePy = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_Py";
865 std::string extraVariableParticlePz = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_Pz";
866 std::string extraVariableParticleE = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_E";
867 std::string extraVariableParticlePxErr = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_PxErr";
868 std::string extraVariableParticlePyErr = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_PyErr";
869 std::string extraVariableParticlePzErr = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_PzErr";
870 std::string extraVariableParticleEErr = m_prefix + "OrcaKinFit" + fitSuffix + "_" + SSTR(iChild) + "_EErr";
871
872 mother->addExtraInfo(extraVariableParticlePx, tlv.Px());
873 mother->addExtraInfo(extraVariableParticlePy, tlv.Py());
874 mother->addExtraInfo(extraVariableParticlePz, tlv.Pz());
875 mother->addExtraInfo(extraVariableParticleE, tlv.E());
876 mother->addExtraInfo(extraVariableParticlePxErr, getFitObjectError(fitobject, 0));
877 mother->addExtraInfo(extraVariableParticlePyErr, getFitObjectError(fitobject, 1));
878 mother->addExtraInfo(extraVariableParticlePzErr, getFitObjectError(fitobject, 2));
879 mother->addExtraInfo(extraVariableParticleEErr, -1.0);
880
881 }
882
883 return updated;
884}
885
887{
888 //check if it is a PxPyPzMFitObject
889 auto* pxpypzmfitobject = static_cast<PxPyPzMFitObject*>(fitobject);
890 if (pxpypzmfitobject) {
891 return fitobject->getError(ilocal);
892 } else {
893 B2FATAL("ParticleKinematicFitterModule: not implemented yet");
894 }
895}
896
897// this returns the local covariance matrix
899{
900
901 //check if it is a PxPyPzMFitObject
902 auto* pxpypzmfitobject = static_cast<PxPyPzMFitObject*>(fitobject);
903 if (pxpypzmfitobject) {
904
905 TMatrixFSym errMatrix(3);
906
907 //loop over the i-j local variables.
908 for (int i = 0; i < 3; i++) {
909 for (int j = 0; j < 3; j++) {
910 errMatrix[i][j] = pxpypzmfitobject->getCov(i, j);
911 }
912 }
913
914 return errMatrix;
915 } else {
916 B2FATAL("ParticleKinematicFitterModule: not implemented yet");
917 }
918}
919
921{
922 TMatrixFSym fitCovMatrix(3);
923
924 if (strcmp(fitobject->getParamName(0), "E") == 0) {
925 //check if it is a JetFitObject
926 auto* jetfitObject = static_cast<JetFitObject*>(fitobject);
927 if (jetfitObject) {
928
929 fitCovMatrix = getFitObjectCovMat(fitobject);
930 ROOT::Math::PxPyPzEVector tlv = getLorentzVector(fitobject);
931
932 const double energy = tlv.E();
933 const double theta = tlv.Theta();
934 const double phi = tlv.Phi();
935
936 const double st = sin(theta);
937 const double ct = cos(theta);
938 const double sp = sin(phi);
939 const double cp = cos(phi);
940
941 // updated covariance matrix is: A * cov * A^T where A is the Jacobi matrix (use Similarity)
942 TMatrixF A(7, 3);
943 A(0, 0) = cp * st ; // dpx/dE
944 A(0, 1) = energy * cp * ct ; // dpx/dtheta
945 A(0, 2) = -energy * sp * st ; // dpx/dphi
946 A(1, 0) = sp * st ; // dpy/dE
947 A(1, 1) = energy * sp * ct ; // dpz/dtheta
948 A(1, 2) = energy * cp * st ; // dpy/dphi
949 A(2, 0) = ct ; // dpz/dE
950 A(2, 1) = -energy * st ; // dpz/dtheta
951 A(2, 2) = 0 ; // dpz/dphi
952 A(3, 0) = 1.0; // dE/dE
953 A(3, 1) = 0.0; // dE/dphi
954 A(3, 2) = 0.0; // dE/dtheta
955
956 TMatrixFSym D = fitCovMatrix.Similarity(A);
957 return D;
958
959 } else {
960 B2FATAL("ParticleKinematicFitterModule: not implemented yet");
961 }
962 } else {
963 //check if it is a PxPyPzMFitObject
964 auto* pxpypzmfitobject = static_cast<PxPyPzMFitObject*>(fitobject);
965 if (pxpypzmfitobject) {
966
967 fitCovMatrix = getFitObjectCovMat(fitobject);
968
969 // updated covariance matrix is: A * cov * A^T where A is the Jacobi matrix (use Similarity)
970 ROOT::Math::PxPyPzEVector tlv = getLorentzVector(fitobject);
971 TMatrixF A(7, 3);
972 A[0][0] = 1.; // px/dpx
973 A[0][1] = 0.; // px/dpy
974 A[0][2] = 0.; // px/dpz
975 A[1][0] = 0.; // py/dpx
976 A[1][1] = 1.; // py/dpy
977 A[1][2] = 0.; // py/dpz
978 A[2][0] = 0.; // pz/dpx
979 A[2][1] = 0.; // pz/dpy
980 A[2][2] = 1.; // pz/dpz
981 if (tlv.E() > 0.0) {
982 A[3][0] = tlv.Px() / tlv.E(); // E/dpx, E=sqrt(px^2 + py^2 + pz^2 + m^2)
983 A[3][1] = tlv.Py() / tlv.E(); // E/dpy
984 A[3][2] = tlv.Pz() / tlv.E(); // E/dpz
985 }
986
987 TMatrixFSym D = fitCovMatrix.Similarity(A);
988
989 return D;
990 } else {
991 B2FATAL("ParticleKinematicFitterModule: not implemented yet");
992 }
993 }
994}
Class to provide momentum-related information from ECLClusters.
const TMatrixDSym GetCovarianceMatrix3x3FromCluster(const ECLCluster *cluster, int particleHypo=Const::photon.getPDGCode())
Returns 3x3 covariance matrix (E, theta, phi)
static const ParticleType neutron
neutron particle
Definition Const.h:675
static const ParticleType photon
photon particle
Definition Const.h:673
In the store you can park objects that have to be accessed by various modules.
Definition DataStore.h:51
ECL cluster data.
Definition ECLCluster.h:27
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
virtual const char * getParamName(int ilocal) const =0
Get name of parameter ilocal.
virtual double getError(int ilocal) const
Get error of parameter ilocal.
virtual bool setParam(int ilocal, double par_, bool measured_, bool fixed_=false)
Set value and measured flag of parameter i; return: significant change.
void setName(const char *name_)
Set object's name.
Abstract base class for fitting engines of kinematic fits.
Definition BaseFitter.h:47
Class for jets with (E, eta, phi) in kinematic fits.
Implements constraint 0 = mass1 - mass2 - m.
Implements a constraint of the form efact*sum(E)+pxfact*sum(px)+pyfact*sum(py)+pzfact*sum(pz)=value.
A kinematic fitter using the Newton-Raphson method to solve the equations.
Description of the fit algorithm and interface:
virtual double getPx() const
Return px.
virtual double getPz() const
Return pz.
virtual double getPy() const
Return py.
virtual double getE() const
Return E.
bool fillFitParticles(Particle *mother, std::vector< Particle * > &particleChildren)
Fills valid particle's children (with valid error matrix) in the vector of Particles that will enter ...
void addFitObjectToConstraints(ParticleFitObject &fitobject)
Adds Orca fit object to the constraints.
void updateMapOfTrackAndDaughter(unsigned &l, std::vector< std::vector< unsigned > > &pars, std::vector< unsigned > &pard, std::vector< Particle * > &allparticles, const Particle *daughter)
update the map of daughter and tracks, find out which tracks belong to each daughter.
bool m_liftPhotonTheta
lift theta constraint of the 3CPhoton.
bool AddFour(Particle *mother)
Added four vectors and calculated a covariance matrix for a combined particles.
std::string m_orcaConstraint
Constraint (softBeam, hardBeam (default))
std::string m_decayStringForNeutronVsAntiNeutron
n or nbar particle tag selection
int m_debugFitterLevel
internal debugging level (for New and Newton fitter only)
virtual void initialize() override
Initialize the Module.
CLHEP::HepSymMatrix getCLHEPMomentumVertexErrorMatrix(Particle *particle)
Returns particle's 7x7 momentum-vertex-error matrix as a HepSymMatrix.
MomentumConstraint m_hardConstraintE
hard beam constraint E
std::string m_prefix
prefix attached to extrainfo names
DecayDescriptor m_decaydescriptorForDirectionOnlyParticles
Decay descriptor of direction only particles selection.
bool m_debugFitter
activate internal debugging (for New and Newton fitter only)
void addParticleToOrcaKinFit(BaseFitter &fitter, Particle *particle, const int index, bool useOnlyDirection, int massHypoPDG)
Adds given particle to the OrcaKinFit.
void addUnmeasuredGammaToOrcaKinFit(BaseFitter &fitter)
Adds an unmeasured gamma (E, phi, theta) to the fit (-3C) stored as EventExtraInfo TODO.
std::string m_decayStringForAlternateMassParticles
alternate mass particles selection
MomentumConstraint m_hardConstraintPz
hard beam constraint pz
float getFitObjectError(ParticleFitObject *fitobject, int ilocal)
Returns fit object error on the parameter ilocal.
ROOT::Math::PxPyPzEVector getLorentzVectorConstraints()
Get constraints (at whatever stage before/after fitting)
std::vector< int > m_listAlternateMassHypo
index of particles where only direction is used
bool doOrcaKinFitFit(Particle *p)
Kinematic fit using OrcaKinFit.
bool m_add3CPhoton
add one photon with unmeasured energy to the fit (costs 1 constraints)
TMatrixFSym getTMatrixFSymMomentumErrorMatrix()
Returns particle's 7x7 momentum-error matrix as a TMatrixFSym.
DecayDescriptor m_decaydescriptorForAlternateMassParticles
Decay descriptor of alternate particles selection.
void addConstraintsToFitter(BaseFitter &fitter)
Adds Orca fit object to the constraints.
void addTracerToFitter(BaseFitter &fitter)
Adds tracer to the fitter.
RecoilMassConstraint m_hardConstraintRecoilMass
hard recoil mass constraint
TMatrixFSym getTMatrixFSymMomentumVertexErrorMatrix()
Returns particle's 7x7 momentum-vertex-error matrix as a TMatrixFSym.
StoreObjPtr< EventExtraInfo > m_eventextrainfo
StoreObjPtr for the EventExtraInfo in this mode.
MomentumConstraint m_hardConstraintPy
hard beam constraint py
double m_recoilMass
Recoil mass for RecoilMass constraint.
CLHEP::HepLorentzVector getCLHEPLorentzVector(Particle *particle)
Returns particle's 4-momentum as a HepLorentzVector.
TextTracer * m_textTracer
internal text output variable
std::string m_decayStringForDirectionOnlyParticles
direction only particles selection
DecayDescriptor m_decaydescriptor
Decay descriptor of decays to look for.
TMatrixFSym getCovMat7(ParticleFitObject *fitobject)
Returns covariance matrix.
ROOT::Math::PxPyPzEVector getLorentzVector(ParticleFitObject *fitobject)
bool updateOrcaKinFitDaughters(BaseFitter &fitter, Particle *mother)
Update the daughters: momentum is sum of daughters TODO update covariance matrix.
void setConstraints()
Sets constraints, this is not connect to particles or a fitter at this stage.
bool m_addUnmeasuredPhoton
add one unmeasured photon to the fit (costs 3 constraints)
TMatrixFSym getFitObjectCovMat(ParticleFitObject *fitobject)
Returns covariance matrix.
bool storeOrcaKinFitParticles(const std::string &fitSuffix, BaseFitter &fitter, std::vector< Particle * > &particleChildren, Particle *mother)
store fit object information as ExtraInfo
CLHEP::HepSymMatrix getCLHEPMomentumErrorMatrix(Particle *particle)
Returns particle's 4x4 momentum-error matrix as a HepSymMatrix.
bool m_fixUnmeasuredPhotonToHER
fix the momentum of the unmeasured photon to HER
DecayDescriptor m_decaydescriptorForNeutronVsAntiNeutron
Decay descriptor of n or nbar particle tag selection.
StoreObjPtr< ParticleList > m_plist
StoreObjPtr for the particle list.
void updateOrcaKinFitMother(BaseFitter &fitter, std::vector< Particle * > &particleChildren, Particle *mother)
Update the mother: momentum is sum of daughters TODO update covariance matrix.
void resetFitter(BaseFitter &fitter)
Resets all objects associated with the OrcaKinFit fitter.
bool m_fixUnmeasuredPhotonToLER
fix the momentum of the unmeasured photon to LER
bool doKinematicFit(Particle *p)
Main steering routine for any kinematic fitter.
MomentumConstraint m_hardConstraintPx
hard beam constraint px
Class to hold Lorentz transformations from/to CMS and boost vector.
ROOT::Math::PxPyPzEVector getBeamFourMomentum() const
Returns LAB four-momentum of e+e-, i.e.
Class to store reconstructed particles.
Definition Particle.h:76
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition Particle.h:687
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition Particle.h:651
unsigned getNDaughters(void) const
Returns number of daughter particles.
Definition Particle.h:747
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition Particle.cc:1421
std::vector< Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition Particle.cc:668
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition Particle.cc:424
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition Particle.h:377
void updateMomentum(const ROOT::Math::PxPyPzEVector &p4, const ROOT::Math::XYZVector &vertex, const TMatrixFSym &errMatrix, double pValue)
Sets Lorentz vector, position, 7x7 error matrix and p-value.
Definition Particle.h:397
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition Particle.cc:451
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition Particle.cc:662
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
void copyDaughters(Particle *mother)
Function copies all (grand-)^n-daughter particles of the argument mother Particle.
Abstract base class for different kinds of events.
STL namespace.