Belle II Software light-2406-ragdoll
TagVertexModule.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
9#include <analysis/modules/TagVertex/TagVertexModule.h>
10
11//to help printing out stuff
12#include<sstream>
13
14// framework aux
15#include <framework/gearbox/Unit.h>
16#include <framework/gearbox/Const.h>
17#include <framework/logging/Logger.h>
18
19// dataobjects
20#include <analysis/dataobjects/RestOfEvent.h>
21#include <analysis/dataobjects/FlavorTaggerInfo.h>
22
23// utilities
24#include <analysis/utility/PCmsLabTransform.h>
25#include <analysis/variables/TrackVariables.h>
26#include <analysis/utility/ParticleCopy.h>
27#include <analysis/utility/CLHEPToROOT.h>
28#include <analysis/utility/ROOTToCLHEP.h>
29#include <analysis/utility/DistanceTools.h>
30#include <analysis/utility/RotationTools.h>
31
32// vertex fitting
33#include <analysis/VertexFitting/KFit/VertexFitKFit.h>
34
35// mdst dataobject
36#include <mdst/dataobjects/HitPatternVXD.h>
37
38// Magnetic field
39#include <framework/geometry/BFieldManager.h>
40
41#include <TVector.h>
42#include <TRotation.h>
43#include <Math/Vector4D.h>
44
45using namespace std;
46using namespace Belle2;
47
49static const double realNaN = std::numeric_limits<double>::quiet_NaN();
50
52static const ROOT::Math::XYZVector vecNaN(realNaN, realNaN, realNaN);
53
55static const double arrayNaN[] = {
56 realNaN, realNaN, realNaN,
57 realNaN, realNaN, realNaN,
58 realNaN, realNaN, realNaN,
59};
60
62static const TMatrixDSym matNaN(3, arrayNaN);
63
64// import tools from RotationTools.h
65using RotationTools::rotateTensor;
66using RotationTools::rotateTensorInv;
67using RotationTools::toSymMatrix;
68using RotationTools::toVec;
69using RotationTools::getUnitOrthogonal;
70
71//-----------------------------------------------------------------
72// Register the Module
73//-----------------------------------------------------------------
75
76//-----------------------------------------------------------------
77// Implementation
78//-----------------------------------------------------------------
79
81 m_Bfield(0), m_fitTruthStatus(0), m_rollbackStatus(0), m_fitPval(0), m_mcTagLifeTime(-1), m_mcPDG(0), m_mcLifeTimeReco(-1),
82 m_deltaT(0), m_deltaTErr(0), m_mcDeltaTau(0), m_mcDeltaT(0),
83 m_FitType(0), m_tagVl(0),
84 m_truthTagVl(0), m_tagVlErr(0), m_tagVol(0), m_truthTagVol(0), m_tagVolErr(0), m_tagVNDF(0), m_tagVChi2(0), m_tagVChi2IP(0),
85 m_kFitReqReducedChi2(5),
86 m_verbose(true)
87{
88 // Set module properties
89 setDescription("Tag side Vertex Fitter for modular analysis");
90
91 // Parameter definitions
92 addParam("listName", m_listName, "name of particle list", string(""));
93 addParam("confidenceLevel", m_confidenceLevel,
94 "required confidence level of fit to keep particles in the list. Note that even with confidenceLevel == 0.0, errors during the fit might discard Particles in the list. confidenceLevel = -1 if an error occurs during the fit",
95 0.001);
96 addParam("MCAssociation", m_useMCassociation,
97 "'': no MC association. breco: use standard Breco MC association. internal: use internal MC association", string("breco"));
98 addParam("constraintType", m_constraintType,
99 "Choose the type of the constraint: noConstraint, IP (tag tracks constrained to be within the beam spot), tube (long tube along the BTag line of flight, only for fully reconstruced B rec), boost (long tube along the Upsilon(4S) boost direction), (breco)",
100 string("tube"));
101 addParam("trackFindingType", m_trackFindingType,
102 "Choose how to reconstruct the tracks on the tag side: standard, standard_PXD",
103 string("standard_PXD"));
104 addParam("maskName", m_roeMaskName,
105 "Choose ROE mask to get particles from ", string(RestOfEvent::c_defaultMaskName));
106 addParam("askMCInformation", m_mcInfo,
107 "TRUE when requesting MC Information from the tracks performing the vertex fit", false);
108 addParam("reqPXDHits", m_reqPXDHits,
109 "Minimum number of PXD hits for a track to be used in the vertex fit", 0);
110 addParam("fitAlgorithm", m_fitAlgo,
111 "Fitter used for the tag vertex fit: Rave or KFit", string("KFit"));
112 addParam("kFitReqReducedChi2", m_kFitReqReducedChi2,
113 "The required chi2/ndf to accept the kFit result, if it is higher, iteration procedure is applied", 5.0);
114 addParam("useTruthInFit", m_useTruthInFit,
115 "Use the true track parameters in the vertex fit", false);
116 addParam("useRollBack", m_useRollBack,
117 "Use rolled back non-primary tracks", false);
118}
119
121{
122 // magnetic field
123 m_Bfield = BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(m_BeamSpotCenter)).Z();
124 // RAVE setup
126 B2INFO("TagVertexModule : magnetic field = " << m_Bfield);
127 // truth fit status will be set to 2 only if the MC info cannot be recovered
129 // roll back status will be set to 2 only if the MC info cannot be recovered
131
132 //input
134 m_plist.isRequired(m_listName);
135 // output
136 m_verArray.registerInDataStore();
138 //check if the fitting algorithm name is set correctly
139 if (m_fitAlgo != "Rave" && m_fitAlgo != "KFit")
140 B2FATAL("TagVertexModule: invalid fitting algorithm (must be set to either Rave or KFit).");
142 B2FATAL("TagVertexModule: invalid fitting option (useRollBack and useTruthInFit cannot be simultaneously set to true).");
143 //temporary while the one track fit is broken
144 if (m_trackFindingType == "singleTrack" || m_trackFindingType == "singleTrack_PXD")
145 B2FATAL("TagVertexModule : the singleTrack option is temporarily broken.");
146}
147
149{
150 //TODO: set magnetic field for each run
151 //m_Bfield = BFieldMap::Instance().getBField(m_BeamSpotCenter).Z();
152}
153
155{
156 if (!m_plist) {
157 B2ERROR("TagVertexModule: ParticleList " << m_listName << " not found");
158 return;
159 }
160
161 // output
163
164 std::vector<unsigned int> toRemove;
165
166 for (unsigned i = 0; i < m_plist->getListSize(); ++i) {
168
169 Particle* particle = m_plist->getParticle(i);
170 if (m_useMCassociation == "breco" || m_useMCassociation == "internal") BtagMCVertex(particle);
171 bool ok = doVertexFit(particle);
172 if (ok) deltaT(particle);
173
176 toRemove.push_back(particle->getArrayIndex());
177 } else {
178 // save information in the Vertex StoreArray
179 TagVertex* ver = m_verArray.appendNew();
180 // create relation: Particle <-> Vertex
181 particle->addRelationTo(ver);
182 // fill Vertex with content
183 if (ok) {
184 ver->setTagVertex(m_tagV);
187 ver->setDeltaT(m_deltaT);
193 ver->setFitType(m_FitType);
194 ver->setNTracks(m_tagParticles.size());
195 ver->setTagVl(m_tagVl);
198 ver->setTagVol(m_tagVol);
201 ver->setTagVNDF(m_tagVNDF);
212 } else {
213 ver->setTagVertex(m_tagV);
214 ver->setTagVertexPval(-1.);
215 ver->setDeltaT(m_deltaT);
218 ver->setMCTagBFlavor(0.);
221 ver->setFitType(m_FitType);
222 ver->setNTracks(m_tagParticles.size());
223 ver->setTagVl(m_tagVl);
226 ver->setTagVol(m_tagVol);
229 ver->setTagVNDF(-1111.);
230 ver->setTagVChi2(-1111.);
231 ver->setTagVChi2IP(-1111.);
240 }
241 }
242 }
243 m_plist->removeParticles(toRemove);
244
245 //free memory allocated by rave. initialize() would be enough, except that we must clean things up before program end...
246 //
248}
249
251{
252 //reset the fit truth status in case it was set to 2 in a previous fit
253
255
256 //reset the roll back status in case it was set to 2 in a previous fit
257
259
260 //set constraint type, reset pVal and B field
261
262 m_fitPval = 1;
263
264 if (!(Breco->getRelatedTo<RestOfEvent>())) {
265 m_FitType = -1;
266 return false;
267 }
268
269 if (m_Bfield == 0) {
270 B2ERROR("TagVertex: No magnetic field");
271 return false;
272 }
273
274 // recover beam spot info
275
276 m_BeamSpotCenter = m_beamSpotDB->getIPPosition();
277 m_BeamSpotCov.ResizeTo(3, 3);
278 m_BeamSpotCov = m_beamSpotDB->getCovVertex();
279
280 //make the beam spot bigger for the standard constraint
281
282 double beta = PCmsLabTransform().getBoostVector().Mag();
283 double bg = beta / sqrt(1 - beta * beta);
284
285 //TODO: What's the origin of these numbers?
286 double tauB = 1.519; //B0 lifetime in ps
287 double c = Const::speedOfLight / 1000.; // cm ps-1
288 double lB0 = tauB * bg * c;
289
290 //tube length here set to 20 * 2 * c tau beta gamma ~= 0.5 cm, should be enough to not bias the decay
291 //time but should still help getting rid of some pions from kshorts
292 m_constraintCov.ResizeTo(3, 3);
294 else if (m_constraintType == "tube") tie(m_constraintCenter, m_constraintCov) = findConstraintBTube(Breco, 200 * lB0);
295 else if (m_constraintType == "boost") tie(m_constraintCenter, m_constraintCov) = findConstraintBoost(200 * lB0);
296 else if (m_constraintType == "breco") tie(m_constraintCenter, m_constraintCov) = findConstraint(Breco, 200 * lB0);
297 else if (m_constraintType == "noConstraint") m_constraintCenter = ROOT::Math::XYZVector(); //zero vector
298 else {
299 B2ERROR("TagVertex: Invalid constraintType selected");
300 return false;
301 }
302
303 if (m_constraintCenter == vecNaN) {
304 B2ERROR("TagVertex: No correct fit constraint");
305 return false;
306 }
307
308 /* Depending on the user's choice, one of the possible algorithms is chosen for the fit. In case the algorithm does not converge, in order to assure
309 high efficiency, the next algorithm less restrictive is used. I.e, if standard_PXD does not work, the program tries with standard.
310 */
311
312 m_FitType = 0;
313 double minPVal = (m_fitAlgo != "KFit") ? 0.001 : 0.;
314 bool ok = false;
315
316 if (m_trackFindingType == "standard_PXD") {
318 if (m_tagParticles.size() > 0) {
319 ok = makeGeneralFit();
320 m_FitType = 3;
321 }
322 }
323
324 if (ok == false || m_fitPval < minPVal || m_trackFindingType == "standard") {
326 ok = m_tagParticles.size() > 0;
327 if (ok) {
328 ok = makeGeneralFit();
329 m_FitType = 4;
330 }
331 }
332
333 if ((ok == false || (m_fitPval <= 0. && m_fitAlgo == "Rave")) && m_constraintType != "noConstraint") {
335 ok = (m_constraintCenter != vecNaN);
336 if (ok) {
338 ok = (m_tagParticles.size() > 0);
339 }
340 if (ok) {
341 ok = makeGeneralFit();
342 m_FitType = 5;
343 }
344 }
345
346 return ok;
347}
348
349pair<ROOT::Math::XYZVector, TMatrixDSym> TagVertexModule::findConstraint(const Particle* Breco, double cut) const
350{
351 if (Breco->getPValue() < 0.) return make_pair(vecNaN, matNaN);
352
353 TMatrixDSym beamSpotCov(3);
354 beamSpotCov = m_beamSpotDB->getCovVertex();
355
357
358 double pmag = Breco->getMomentumMagnitude();
359 double xmag = (Breco->getVertex() - m_BeamSpotCenter).R();
360
361
362 TMatrixDSym TerrMatrix = Breco->getMomentumVertexErrorMatrix();
363 TMatrixDSym PerrMatrix(7);
364
365 for (int i = 0; i < 3; ++i) {
366 for (int j = 0; j < 3; ++j) {
367 if (i == j) {
368 PerrMatrix(i, j) = (beamSpotCov(i, j) + TerrMatrix(i, j)) * pmag / xmag;
369 } else {
370 PerrMatrix(i, j) = TerrMatrix(i, j);
371 }
372 PerrMatrix(i + 4, j + 4) = TerrMatrix(i + 4, j + 4);
373 }
374 }
375
376 PerrMatrix(3, 3) = 0.;
377
378 //Copy Breco, but use errors as are in PerrMatrix
379 Particle* Breco2 = ParticleCopy::copyParticle(Breco);
380 Breco2->setMomentumVertexErrorMatrix(PerrMatrix);
381
382
383 const Particle* BRecoRes = doVertexFitForBTube(Breco2, "kalman");
384 if (BRecoRes->getPValue() < 0) return make_pair(vecNaN, matNaN); //problems
385
386 // Overall error matrix
387 TMatrixDSym errFinal = TMatrixDSym(Breco->getVertexErrorMatrix() + BRecoRes->getVertexErrorMatrix());
388
389 // TODO : to be developed the extraction of the momentum from the rave fitted track
390
391 // Get expected pBtag 4-momentum using transverse-momentum conservation
392 ROOT::Math::XYZVector BvertDiff = pmag * (Breco->getVertex() - BRecoRes->getVertex()).Unit();
393 ROOT::Math::PxPyPzMVector pBrecEstimate(BvertDiff.X(), BvertDiff.Y(), BvertDiff.Z(), Breco->getPDGMass());
394 ROOT::Math::PxPyPzMVector pBtagEstimate = PCmsLabTransform::labToCms(pBrecEstimate);
395 pBtagEstimate.SetPxPyPzE(-pBtagEstimate.px(), -pBtagEstimate.py(), -pBtagEstimate.pz(), pBtagEstimate.E());
396 pBtagEstimate = PCmsLabTransform::cmsToLab(pBtagEstimate);
397
398 // rotate err-matrix such that pBrecEstimate goes to eZ
399 TMatrixD TubeZ = rotateTensorInv(pBrecEstimate.Vect(), errFinal);
400
401 TubeZ(2, 2) = cut * cut;
402 TubeZ(2, 0) = 0; TubeZ(0, 2) = 0;
403 TubeZ(2, 1) = 0; TubeZ(1, 2) = 0;
404
405
406 // rotate err-matrix such that eZ goes to pBtagEstimate
407 TMatrixD Tube = rotateTensor(B2Vector3D(pBtagEstimate.Vect()), TubeZ);
408
409 // Standard algorithm needs no shift
410 return make_pair(m_BeamSpotCenter, toSymMatrix(Tube));
411}
412
413pair<ROOT::Math::XYZVector, TMatrixDSym> TagVertexModule::findConstraintBTube(const Particle* Breco, double cut)
414{
415 //Use Breco as the creator of the B tube.
416 if ((Breco->getVertexErrorMatrix()(2, 2)) == 0.0) {
417 B2WARNING("In TagVertexModule::findConstraintBTube: cannot get a proper vertex for BReco. BTube constraint replaced by Boost.");
418 return findConstraintBoost(cut);
419 }
420
421
422 //vertex fit will give the intersection between the beam spot and the trajectory of the B
423 //(base of the BTube, or primary vtx cov matrix)
424 const Particle* tubecreatorBCopy = doVertexFitForBTube(Breco, "avf");
425 if (tubecreatorBCopy->getPValue() < 0) return make_pair(vecNaN, matNaN); //if problems
426
427
428 //get direction of B tag = opposite direction of B rec in CMF
429 ROOT::Math::PxPyPzEVector pBrec = tubecreatorBCopy->get4Vector();
430
431 //if we want the true info, replace the 4vector by the true one
432 if (m_useTruthInFit) {
433 const MCParticle* mcBr = Breco->getRelated<MCParticle>();
434 if (mcBr)
435 pBrec = mcBr->get4Vector();
436 else
438 }
439 ROOT::Math::PxPyPzEVector pBtag = PCmsLabTransform::labToCms(pBrec);
440 pBtag.SetPxPyPzE(-pBtag.px(), -pBtag.py(), -pBtag.pz(), pBtag.E());
441 pBtag = PCmsLabTransform::cmsToLab(pBtag);
442
443 //To create the B tube, strategy is: take the primary vtx cov matrix, and add to it a cov
444 //matrix corresponding to an very big error in the direction of the B tag
445 TMatrixDSym pv = tubecreatorBCopy->getVertexErrorMatrix();
446
447 //print some stuff if wanted
448 if (m_verbose) {
449 B2DEBUG(10, "Brec decay vertex before fit: " << printVector(Breco->getVertex()));
450 B2DEBUG(10, "Brec decay vertex after fit: " << printVector(tubecreatorBCopy->getVertex()));
451 B2DEBUG(10, "Brec direction before fit: " << printVector(float(1. / Breco->getP()) * Breco->getMomentum()));
452 B2DEBUG(10, "Brec direction after fit: " << printVector(float(1. / tubecreatorBCopy->getP()) * tubecreatorBCopy->getMomentum()));
453 B2DEBUG(10, "IP position: " << printVector(ROOT::Math::XYZVector(m_BeamSpotCenter.X(), m_BeamSpotCenter.Y(),
454 m_BeamSpotCenter.Z())));
455 B2DEBUG(10, "IP covariance: " << printMatrix(m_BeamSpotCov));
456 B2DEBUG(10, "Brec primary vertex: " << printVector(tubecreatorBCopy->getVertex()));
457 B2DEBUG(10, "Brec PV covariance: " << printMatrix(pv));
458 B2DEBUG(10, "BTag direction: " << printVector((1. / pBtag.P())*pBtag.Vect()));
459 }
460
461 //make a long error matrix along BTag direction
462 TMatrixD longerror(3, 3); longerror(2, 2) = cut * cut;
463
464
465 // make rotation matrix from z axis to BTag line of flight
466 TMatrixD longerrorRotated = rotateTensor(B2Vector3D(pBtag.Vect()), longerror);
467
468 //pvNew will correspond to the covariance matrix of the B tube
469 TMatrixD pvNew = TMatrixD(pv) + longerrorRotated;
470
471 //set the constraint
472 ROOT::Math::XYZVector constraintCenter = tubecreatorBCopy->getVertex();
473
474 //if we want the true info, set the centre of the constraint to the primary vertex
475 if (m_useTruthInFit) {
476 const MCParticle* mcBr = Breco->getRelated<MCParticle>();
477 if (mcBr) {
478 constraintCenter = mcBr->getProductionVertex();
479 }
480 }
481
482 if (m_verbose) {
483 B2DEBUG(10, "IPTube covariance: " << printMatrix(pvNew));
484 }
485
486 //The following is done to do the BTube constraint with a virtual track
487 //(ie KFit way)
488
489 m_tagMomentum = pBtag;
490
491 m_pvCov.ResizeTo(pv);
492 m_pvCov = pv;
493
494 return make_pair(constraintCenter, toSymMatrix(pvNew));
495}
496
497pair<ROOT::Math::XYZVector, TMatrixDSym> TagVertexModule::findConstraintBoost(double cut) const
498{
499 double d = 20e-4; //average transverse distance flown by B0
500
501 //make a long error matrix along boost direction
502 TMatrixD longerror(3, 3); longerror(2, 2) = cut * cut;
503 longerror(0, 0) = longerror(1, 1) = d * d;
504
506
507 TMatrixD longerrorRotated = rotateTensor(boostDir, longerror);
508
509 //Extend error of BeamSpotCov matrix in the boost direction
510 TMatrixDSym beamSpotCov = m_beamSpotDB->getCovVertex();
511 TMatrixD Tube = TMatrixD(beamSpotCov) + longerrorRotated;
512
513 // Standard algorithm needs no shift
514 ROOT::Math::XYZVector constraintCenter = m_BeamSpotCenter;
515
516 return make_pair(constraintCenter, toSymMatrix(Tube));
517}
518
520static double getProperLifeTime(const MCParticle* mc)
521{
522 double beta = mc->getMomentum().R() / mc->getEnergy();
523 return 1e3 * mc->getLifetime() * sqrt(1 - pow(beta, 2));
524}
525
527{
528 //fill vector with mcB (intended order: Reco, Tag)
529 vector<const MCParticle*> mcBs;
530 for (const MCParticle& mc : m_mcParticles) {
531 if (abs(mc.getPDG()) == abs(Breco->getPDGCode()))
532 mcBs.push_back(&mc);
533 }
534 //too few Bs
535 if (mcBs.size() < 2) return;
536
537 if (mcBs.size() > 2) {
538 B2WARNING("TagVertexModule:: Too many Bs found in MC");
539 }
540
541 auto isReco = [&](const MCParticle * mc) {
542 return (m_useMCassociation == "breco") ? (mc == Breco->getRelated<MCParticle>())
543 : compBrecoBgen(Breco, mc); //internal association
544 };
545
546 //nothing matched?
547 if (!isReco(mcBs[0]) && !isReco(mcBs[1])) {
548 return;
549 }
550
551 //first is Tag, second Reco -> swap the order
552 if (!isReco(mcBs[0]) && isReco(mcBs[1]))
553 swap(mcBs[0], mcBs[1]);
554
555 //both matched -> use closest vertex dist as Reco
556 if (isReco(mcBs[0]) && isReco(mcBs[1])) {
557 double dist0 = (mcBs[0]->getDecayVertex() - Breco->getVertex()).Mag2();
558 double dist1 = (mcBs[1]->getDecayVertex() - Breco->getVertex()).Mag2();
559 if (dist0 > dist1)
560 swap(mcBs[0], mcBs[1]);
561 }
562
563 m_mcVertReco = mcBs[0]->getDecayVertex();
564 m_mcLifeTimeReco = getProperLifeTime(mcBs[0]);
565 m_mcTagV = mcBs[1]->getDecayVertex();
566 m_mcTagLifeTime = getProperLifeTime(mcBs[1]);
567 m_mcPDG = mcBs[1]->getPDG();
568}
569
570// static
572{
573
574 bool isDecMode = true;
575
576 const std::vector<Belle2::Particle*> recDau = Breco->getDaughters();
577 const std::vector<Belle2::MCParticle*> genDau = Bgen->getDaughters();
578
579 if (recDau.size() > 0 && genDau.size() > 0) {
580 for (auto dauRec : recDau) {
581 bool isDau = false;
582 for (auto dauGen : genDau) {
583 if (dauGen->getPDG() == dauRec->getPDGCode())
584 isDau = compBrecoBgen(dauRec, dauGen) ;
585 }
586 if (!isDau) isDecMode = false;
587 }
588 } else {
589 if (recDau.size() == 0) { //&& genDau.size()==0){
590 if (Bgen->getPDG() != Breco->getPDGCode()) isDecMode = false;;
591 } else {isDecMode = false;}
592 }
593
594 return isDecMode;
595}
596
597// STANDARD FIT ALGORITHM
598/* This algorithm basically takes all the tracks coming from the Rest Of Events and send them to perform a multi-track fit
599 The option to request PXD hits for the tracks can be chosen by the user.
600 */
601std::vector<const Particle*> TagVertexModule::getTagTracks_standardAlgorithm(const Particle* Breco, int reqPXDHits) const
602{
603 std::vector<const Particle*> fitParticles;
604 const RestOfEvent* roe = Breco->getRelatedTo<RestOfEvent>();
605 if (!roe) return fitParticles;
606 //load all particles from the ROE
607 std::vector<const Particle*> ROEParticles = roe->getChargedParticles(m_roeMaskName, 0, false);
608 if (ROEParticles.size() == 0) return fitParticles;
609
610 for (auto& ROEParticle : ROEParticles) {
611 HitPatternVXD roeTrackPattern = ROEParticle->getTrackFitResult()->getHitPatternVXD();
612
613 if (roeTrackPattern.getNPXDHits() >= reqPXDHits) {
614 fitParticles.push_back(ROEParticle);
615 }
616 }
617 return fitParticles;
618}
619
620
621vector<ParticleAndWeight> TagVertexModule::getParticlesAndWeights(const vector<const Particle*>& tagParticles) const
622{
623 vector<ParticleAndWeight> particleAndWeights;
624
625 for (const Particle* particle : tagParticles) {
626 ROOT::Math::PxPyPzEVector mom = particle->get4Vector();
627 if (!isfinite(mom.mag2())) continue;
628
629 ParticleAndWeight particleAndWeight;
630 particleAndWeight.mcParticle = 0;
631 particleAndWeight.weight = -1111.;
632 particleAndWeight.particle = particle;
633
634 if (m_useMCassociation == "breco" || m_useMCassociation == "internal")
635 particleAndWeight.mcParticle = particle->getRelatedTo<MCParticle>();
636
637 particleAndWeights.push_back(particleAndWeight);
638 }
639
640 return particleAndWeights;
641}
642
644{
645 if (m_fitAlgo == "Rave") return makeGeneralFitRave();
646 else if (m_fitAlgo == "KFit") return makeGeneralFitKFit();
647 return false;
648}
649
650void TagVertexModule::fillParticles(vector<ParticleAndWeight>& particleAndWeights)
651{
652 unsigned n = particleAndWeights.size();
653 sort(particleAndWeights.begin(), particleAndWeights.end(),
654 [](const ParticleAndWeight & a, const ParticleAndWeight & b) { return a.weight > b.weight; });
655
656 m_raveParticles.resize(n);
657 m_raveWeights.resize(n);
658 m_raveMCParticles.resize(n);
659
660 for (unsigned i = 0; i < n; ++i) {
661 m_raveParticles.at(i) = particleAndWeights.at(i).particle;
662 m_raveMCParticles.at(i) = particleAndWeights.at(i).mcParticle;
663 m_raveWeights.at(i) = particleAndWeights.at(i).weight;
664 }
665}
666
667void TagVertexModule::fillTagVinfo(const ROOT::Math::XYZVector& tagVpos, const TMatrixDSym& tagVposErr)
668{
669 m_tagV = tagVpos;
670
671 if (m_constraintType != "noConstraint") {
672 TMatrixDSym tubeInv = m_constraintCov;
673 tubeInv.Invert();
674 TVectorD dV = toVec(m_tagV - m_BeamSpotCenter);
675 m_tagVChi2IP = tubeInv.Similarity(dV);
676 }
677
678 m_tagVErrMatrix.ResizeTo(tagVposErr);
679 m_tagVErrMatrix = tagVposErr;
680}
681
683{
684 // apply constraint
686 if (m_constraintType != "noConstraint")
689
690 //feed rave with tracks without Kshorts
691 vector<ParticleAndWeight> particleAndWeights = getParticlesAndWeights(m_tagParticles);
692
693 for (const auto& pw : particleAndWeights) {
694 try {
695 if (m_useTruthInFit) {
696 if (pw.mcParticle) {
698 rFit.addTrack(&tfr);
699 } else
701 } else if (m_useRollBack) {
702 if (pw.mcParticle) {
704 rFit.addTrack(&tfr);
705 } else
707 } else {
708 rFit.addTrack(pw.particle->getTrackFitResult());
709 }
710 } catch (const rave::CheckedFloatException&) {
711 B2ERROR("Exception caught in TagVertexModule::makeGeneralFitRave(): Invalid inputs (nan/inf)?");
712 }
713 }
714
715 //perform fit
716
717 int isGoodFit(-1);
718
719 try {
720 isGoodFit = rFit.fit("avf");
721 // if problems
722 if (isGoodFit < 1) return false;
723 } catch (const rave::CheckedFloatException&) {
724 B2ERROR("Exception caught in TagVertexModule::makeGeneralFitRave(): Invalid inputs (nan/inf)?");
725 return false;
726 }
727
728 //save the track info for later use
729
730 for (unsigned int i(0); i < particleAndWeights.size() && isGoodFit >= 1; ++i)
731 particleAndWeights.at(i).weight = rFit.getWeight(i);
732
733 //Tracks are sorted from highest rave weight to lowest
734
735 fillParticles(particleAndWeights);
736
737 //if the fit is good, save the infos related to the vertex
738 fillTagVinfo(ROOT::Math::XYZVector(rFit.getPos(0)), rFit.getCov(0));
739
740 //fill quality variables
741 m_tagVNDF = rFit.getNdf(0);
742 m_tagVChi2 = rFit.getChi2(0);
743 m_fitPval = rFit.getPValue();
744
745 return true;
746}
747
748
749analysis::VertexFitKFit TagVertexModule::doSingleKfit(vector<ParticleAndWeight>& particleAndWeights)
750{
751 //initialize KFit
753 kFit.setMagneticField(m_Bfield);
754
755 // apply constraint
756 if (m_constraintType != "noConstraint") {
757 if (m_constraintType == "tube") {
758 CLHEP::HepSymMatrix err(7, 0);
759 //copy m_pvCov to the end of err matrix
760 err.sub(5, ROOTToCLHEP::getHepSymMatrix(m_pvCov));
761 kFit.setIpTubeProfile(
762 ROOTToCLHEP::getHepLorentzVector(m_tagMomentum),
763 ROOTToCLHEP::getPoint3DFromB2Vector(m_constraintCenter),
764 err,
765 0.);
766 } else {
767 kFit.setIpProfile(ROOTToCLHEP::getPoint3DFromB2Vector(m_constraintCenter),
768 ROOTToCLHEP::getHepSymMatrix(m_constraintCov));
769 }
770 }
771
772
773 for (auto& pawi : particleAndWeights) {
774 int addedOK = 1;
775 if (m_useTruthInFit) {
776 if (pawi.mcParticle) {
777 addedOK = kFit.addTrack(
778 ROOTToCLHEP::getHepLorentzVector(pawi.mcParticle->get4Vector()),
779 ROOTToCLHEP::getPoint3DFromB2Vector(getTruePoca(pawi)),
780 ROOTToCLHEP::getHepSymMatrix(pawi.particle->getMomentumVertexErrorMatrix()),
781 pawi.particle->getCharge());
782 } else {
784 }
785 } else if (m_useRollBack) {
786 if (pawi.mcParticle) {
787 addedOK = kFit.addTrack(
788 ROOTToCLHEP::getHepLorentzVector(pawi.mcParticle->get4Vector()),
789 ROOTToCLHEP::getPoint3DFromB2Vector(getRollBackPoca(pawi)),
790 ROOTToCLHEP::getHepSymMatrix(pawi.particle->getMomentumVertexErrorMatrix()),
791 pawi.particle->getCharge());
792 } else {
794 }
795 } else {
796 addedOK = kFit.addParticle(pawi.particle);
797 }
798
799 if (addedOK == 0) {
800 pawi.weight = 1.;
801 } else {
802 B2WARNING("TagVertexModule::makeGeneralFitKFit: failed to add a track");
803 pawi.weight = 0.;
804 }
805 }
806
807
808 int nTracksAdded = kFit.getTrackCount();
809
810 //perform fit if there are enough tracks
811 if ((nTracksAdded < 2 && m_constraintType == "noConstraint") || nTracksAdded < 1)
813
814 int isGoodFit = kFit.doFit();
815 if (isGoodFit != 0) return analysis::VertexFitKFit();
816
817 return kFit;
818}
819
820
821static int getLargestChi2ID(const analysis::VertexFitKFit& kFit)
822{
823 int largest_chi2_trackid = -1;
824 double largest_track_chi2 = -1;
825 for (int i = 0; i < kFit.getTrackCount(); ++i) {
826 double track_chi2 = kFit.getTrackCHIsq(i);
827 if (track_chi2 > largest_track_chi2) {
828 largest_track_chi2 = track_chi2;
829 largest_chi2_trackid = i;
830 }
831 }
832 return largest_chi2_trackid;
833}
834
835
836
837//uses m_tagMomentum, m_constraintCenter, m_constraintCov, m_tagParticles
839{
840 //feed KFit with tracks without Kshorts
841 vector<ParticleAndWeight> particleAndWeights = getParticlesAndWeights(m_tagParticles);
842
844
845
846 // iterative procedure which removes tracks with high chi2
847 for (int iteration_counter = 0; iteration_counter < 100; ++iteration_counter) {
848 analysis::VertexFitKFit kFitTemp = doSingleKfit(particleAndWeights);
849 if (!kFitTemp.isFitted() || isnan(kFitTemp.getCHIsq()))
850 return false;
851
852 double reduced_chi2 = kFitTemp.getCHIsq() / kFitTemp.getNDF();
853 int nTracks = kFitTemp.getTrackCount();
854
855 if (nTracks != int(particleAndWeights.size()))
856 B2ERROR("TagVertexModule: Different number of tracks in kFit and particles");
857
858 if (reduced_chi2 <= m_kFitReqReducedChi2 || nTracks <= 1 || (nTracks <= 2 && m_constraintType == "noConstraint")) {
859 kFit = kFitTemp;
860 break;
861 } else { // remove particle with highest chi2/ndf and continue
862 int badTrackID = getLargestChi2ID(kFitTemp);
863 if (0 <= badTrackID && badTrackID < int(particleAndWeights.size()))
864 particleAndWeights.erase(particleAndWeights.begin() + badTrackID);
865 else
866 B2ERROR("TagVertexModule: Obtained badTrackID is not within limits");
867 }
868
869 }
870
871 //save the track info for later use
872 //Tracks are sorted by weight, i.e. pushing the tracks with 0 weight (from KS) to the end of the list
873 fillParticles(particleAndWeights);
874
875 //Save the infos related to the vertex
876 fillTagVinfo(CLHEPToROOT::getXYZVector(kFit.getVertex()),
877 CLHEPToROOT::getTMatrixDSym(kFit.getVertexError()));
878
879 m_tagVNDF = kFit.getNDF();
880 m_tagVChi2 = kFit.getCHIsq();
881 m_fitPval = TMath::Prob(m_tagVChi2, m_tagVNDF);
882
883 return true;
884}
885
887{
888
890 ROOT::Math::XYZVector boostDir = ROOT::Math::XYZVector(boost.Unit());
891 double bg = boost.Mag() / sqrt(1 - boost.Mag2());
892 double c = Const::speedOfLight / 1000.; // cm ps-1
893
894 //Reconstructed DeltaL & DeltaT in the boost direction
895 ROOT::Math::XYZVector dVert = Breco->getVertex() - m_tagV; //reconstructed vtxReco - vtxTag
896 double dl = dVert.Dot(boostDir);
897 m_deltaT = dl / (bg * c);
898
899 //Truth DeltaL & approx DeltaT in the boost direction
900 ROOT::Math::XYZVector MCdVert = m_mcVertReco - m_mcTagV; //truth vtxReco - vtxTag
901 double MCdl = MCdVert.Dot(boostDir);
902 m_mcDeltaT = MCdl / (bg * c);
903
904 // MCdeltaTau=tauRec-tauTag
906 if (m_mcLifeTimeReco == -1 || m_mcTagLifeTime == -1)
907 m_mcDeltaTau = realNaN;
908
909 TVectorD bVec = toVec(B2Vector3D(boostDir));
910
911 //TagVertex error in boost dir
912 m_tagVlErr = sqrt(m_tagVErrMatrix.Similarity(bVec));
913
914 //bReco error in boost dir
915 double bRecoErrL = sqrt(Breco->getVertexErrorMatrix().Similarity(bVec));
916
917 //Delta t error
918 m_deltaTErr = hypot(m_tagVlErr, bRecoErrL) / (bg * c);
919
920 m_tagVl = m_tagV.Dot(boostDir);
921 m_truthTagVl = m_mcTagV.Dot(boostDir);
922
923 // calculate tagV component and error in the direction orthogonal to the boost
924 B2Vector3D oboost = getUnitOrthogonal(B2Vector3D(boostDir));
925 TVectorD oVec = toVec(oboost);
926
927 //TagVertex error in boost-orthogonal dir
928 m_tagVolErr = sqrt(m_tagVErrMatrix.Similarity(oVec));
929
930 m_tagVol = m_tagV.Dot(oboost);
931 m_truthTagVol = m_mcTagV.Dot(oboost);
932}
933
934Particle* TagVertexModule::doVertexFitForBTube(const Particle* motherIn, std::string fitType) const
935{
936 //make a copy of motherIn to not modify the original object
937 Particle* mother = ParticleCopy::copyParticle(motherIn);
938
939 //Here rave is used to find the upsilon(4S) vtx as the intersection
940 //between the mother B trajectory and the beam spot
942
944 rsg.addTrack(mother);
945 int nvert = rsg.fit(fitType);
946 if (nvert != 1) {
947 mother->setPValue(-1); //error
948 return mother;
949 } else {
950 rsg.updateDaughters();
951 return mother;
952 }
953}
954
955
956
958{
959 if (!paw.mcParticle) {
960 B2ERROR("In TagVertexModule::getTrackWithTrueCoordinate: no MC particle set");
961 return TrackFitResult();
962 }
963
964 const TrackFitResult* tfr(paw.particle->getTrackFitResult());
965
966 return TrackFitResult(getTruePoca(paw),
967 paw.mcParticle->getMomentum(),
968 tfr->getCovariance6(),
969 tfr->getChargeSign(),
970 tfr->getParticleType(),
971 tfr->getPValue(),
972 m_Bfield, 0, 0, tfr->getNDF());
973}
974
975// static
976ROOT::Math::XYZVector TagVertexModule::getTruePoca(ParticleAndWeight const& paw)
977{
978 if (!paw.mcParticle) {
979 B2ERROR("In TagVertexModule::getTruePoca: no MC particle set");
980 return ROOT::Math::XYZVector(0., 0., 0.);
981 }
982
984 paw.mcParticle->getMomentum(),
986}
987
989{
990 const TrackFitResult* tfr(paw.particle->getTrackFitResult());
991
993 tfr->getMomentum(),
994 tfr->getCovariance6(),
995 tfr->getChargeSign(),
996 tfr->getParticleType(),
997 tfr->getPValue(),
998 m_Bfield, 0, 0, tfr->getNDF());
999}
1000
1002{
1003 if (!paw.mcParticle) {
1004 B2ERROR("In TagVertexModule::getTruePoca: no MC particle set");
1005 return ROOT::Math::XYZVector(0., 0., 0.);
1006 }
1007
1008 return paw.particle->getTrackFitResult()->getPosition() - paw.mcParticle->getProductionVertex() + ROOT::Math::XYZVector(m_mcTagV);
1009}
1010
1012{
1013 m_raveParticles.resize(0);
1014 m_raveMCParticles.resize(0);
1015 m_tagParticles.resize(0);
1016 m_raveWeights.resize(0);
1017
1018 m_fitPval = realNaN;
1019 m_tagV = vecNaN;
1020 m_tagVErrMatrix.ResizeTo(matNaN);
1021 m_tagVErrMatrix = matNaN;
1022 m_mcTagV = vecNaN;
1023 m_mcVertReco = vecNaN;
1024 m_deltaT = realNaN;
1025 m_deltaTErr = realNaN;
1026 m_mcDeltaTau = realNaN;
1027 m_constraintCov.ResizeTo(matNaN);
1028 m_constraintCov = matNaN;
1029 m_constraintCenter = vecNaN;
1030 m_tagVl = realNaN;
1031 m_truthTagVl = realNaN;
1032 m_tagVlErr = realNaN;
1033 m_tagVol = realNaN;
1034 m_truthTagVol = realNaN;
1035 m_tagVolErr = realNaN;
1036 m_tagVNDF = realNaN;
1037 m_tagVChi2 = realNaN;
1038 m_tagVChi2IP = realNaN;
1039 m_pvCov.ResizeTo(matNaN);
1040 m_pvCov = matNaN;
1041 m_tagMomentum = ROOT::Math::PxPyPzEVector(realNaN, realNaN, realNaN, realNaN);
1042}
1043
1044//The following functions are just here to help printing stuff
1045
1046// static
1047std::string TagVertexModule::printVector(const ROOT::Math::XYZVector& vec)
1048{
1049 std::ostringstream oss;
1050 int w = 14;
1051 oss << "(" << std::setw(w) << vec.X() << ", " << std::setw(w) << vec.Y() << ", " << std::setw(w) << vec.Z() << ")" << std::endl;
1052 return oss.str();
1053}
1054
1055// static
1056std::string TagVertexModule::printMatrix(const TMatrixD& mat)
1057{
1058 std::ostringstream oss;
1059 int w = 14;
1060 for (int i = 0; i < mat.GetNrows(); ++i) {
1061 for (int j = 0; j < mat.GetNcols(); ++j) {
1062 oss << std::setw(w) << mat(i, j) << " ";
1063 }
1064 oss << endl;
1065 }
1066 return oss.str();
1067}
1068
1069// static
1070std::string TagVertexModule::printMatrix(const TMatrixDSym& mat)
1071{
1072 std::ostringstream oss;
1073 int w = 14;
1074 for (int i = 0; i < mat.GetNrows(); ++i) {
1075 for (int j = 0; j < mat.GetNcols(); ++j) {
1076 oss << std::setw(w) << mat(i, j) << " ";
1077 }
1078 oss << endl;
1079 }
1080 return oss.str();
1081}
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:159
DataType Mag2() const
The magnitude squared (rho^2 in spherical coordinate system).
Definition: B2Vector3.h:157
B2Vector3< DataType > Unit() const
Unit vector parallel to this.
Definition: B2Vector3.h:269
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:61
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
Hit pattern of the VXD within a track.
Definition: HitPatternVXD.h:37
unsigned short getNPXDHits() const
Get total number of hits in the PXD.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition: MCParticle.cc:52
ROOT::Math::XYZVector getProductionVertex() const
Return production vertex position.
Definition: MCParticle.h:189
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:207
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition: MCParticle.h:198
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
Class to hold Lorentz transformations from/to CMS and boost vector.
static ROOT::Math::PxPyPzMVector labToCms(const ROOT::Math::PxPyPzMVector &vec)
Transforms Lorentz vector into CM System.
static ROOT::Math::PxPyPzMVector cmsToLab(const ROOT::Math::PxPyPzMVector &vec)
Transforms Lorentz vector into Laboratory System.
B2Vector3D getBoostVector() const
Returns boost vector (beta=p/E)
Class to store reconstructed particles.
Definition: Particle.h:75
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
Definition: Particle.cc:447
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:667
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:631
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:637
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:604
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:547
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:393
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:366
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:420
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
Definition: Particle.cc:854
double getMomentumMagnitude() const
Returns momentum magnitude.
Definition: Particle.h:569
double getP() const
Returns momentum magnitude (same as getMomentumMagnitude but with shorter name)
Definition: Particle.h:578
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
This is a general purpose class for collecting reconstructed MDST data objects that are not used in r...
Definition: RestOfEvent.h:57
std::vector< const Particle * > getChargedParticles(const std::string &maskName=c_defaultMaskName, unsigned int pdg=0, bool unpackComposite=true) const
Get charged particles from ROE mask.
Definition: RestOfEvent.cc:108
static constexpr const char * c_defaultMaskName
Default mask name.
Definition: RestOfEvent.h:60
bool isRequired(const std::string &name="")
Ensure this array/object has been registered previously.
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:140
int m_fitTruthStatus
Store info about whether the fit was performed with the truth info 0 fit performed with measured para...
TMatrixDSym m_constraintCov
constraint to be used in the tag vertex fit
double m_truthTagVol
MC tagV component in the direction orthogonal to the boost.
std::vector< const Particle * > m_tagParticles
tracks of the rest of the event
double m_tagVl
tagV component in the boost direction
bool doVertexFit(const Particle *Breco)
central method for the tag side vertex fit
std::vector< const Particle * > getTagTracks_standardAlgorithm(const Particle *Breco, int nPXDHits) const
performs the fit using the standard algorithm - using all tracks in RoE The user can specify a reques...
double m_truthTagVl
MC tagV component in the boost direction
std::pair< ROOT::Math::XYZVector, TMatrixDSym > findConstraintBoost(double cut) const
calculate the standard constraint for the vertex fit on the tag side
bool m_useTruthInFit
Set to true if the tag fit is to be made with the TRUE tag track momentum and position.
std::vector< double > m_raveWeights
Store the weights used by Rave in the vtx fit so that they can be accessed later.
TrackFitResult getTrackWithRollBackCoordinates(ParticleAndWeight const &paw)
If the fit has to be done with the rolled back tracks, Rave or KFit is fed with a track where the pos...
virtual void initialize() override
Initialize the Module.
void fillTagVinfo(const ROOT::Math::XYZVector &tagVpos, const TMatrixDSym &tagVposErr)
Fill tagV vertex info.
ROOT::Math::XYZVector m_mcVertReco
generated Breco decay vertex
TMatrixDSym m_pvCov
covariance matrix of the PV (useful with tube and KFit)
double m_mcDeltaT
generated DeltaT with boost-direction approximation
static std::string printVector(const ROOT::Math::XYZVector &vec)
Print a XYZVector (useful for debugging)
ROOT::Math::XYZVector m_tagV
tag side fit result
virtual void event() override
Event processor.
std::pair< ROOT::Math::XYZVector, TMatrixDSym > findConstraintBTube(const Particle *Breco, double cut)
calculate constraint for the vertex fit on the tag side using the B tube (cylinder along the expected...
Particle * doVertexFitForBTube(const Particle *mother, std::string fitType) const
it returns an intersection between B rec and beam spot (= origin of BTube)
std::string m_listName
Breco particle list name.
std::vector< ParticleAndWeight > getParticlesAndWeights(const std::vector< const Particle * > &tagParticles) const
Get a list of particles with attached weight and associated MC particle.
bool m_useRollBack
Set to true if the tag fit is to be made with the tag track position rolled back to mother B.
double m_tagVlErr
Error of the tagV component in the boost direction
std::string m_roeMaskName
ROE particles from this mask will be used for vertex fitting.
double m_tagVChi2
chi^2 value of the tag vertex fit result
static ROOT::Math::XYZVector getTruePoca(ParticleAndWeight const &paw)
This finds the point on the true particle trajectory closest to the measured track position.
void BtagMCVertex(const Particle *Breco)
get the vertex of the MC B particle associated to Btag.
std::pair< ROOT::Math::XYZVector, TMatrixDSym > findConstraint(const Particle *Breco, double cut) const
calculate the constraint for the vertex fit on the tag side using Breco information
bool m_mcInfo
true if user wants to retrieve MC information out from the tracks used in the fit
double m_kFitReqReducedChi2
The required chi2/ndf to accept the kFit result, if it is higher, iteration procedure is applied.
std::string m_useMCassociation
No MC association or standard Breco particle or internal MCparticle association.
ROOT::Math::XYZVector m_constraintCenter
centre position of the constraint for the tag Vertex fit
double m_tagVolErr
Error of the tagV component in the direction orthogonal to the boost.
double m_mcTagLifeTime
generated tag side life time of B-decay
ROOT::Math::XYZVector m_BeamSpotCenter
Beam spot position.
double m_tagVNDF
Number of degrees of freedom in the tag vertex fit.
double m_deltaTErr
reconstructed DeltaT error
std::string m_fitAlgo
Algorithm used for the tag fit (Rave or KFit)
bool makeGeneralFit()
TO DO: tag side vertex fit in the case of semileptonic tag side decay.
ROOT::Math::XYZVector getRollBackPoca(ParticleAndWeight const &paw)
This shifts the position of tracks by the vector difference of mother B and production point of track...
virtual void beginRun() override
Called when entering a new run.
double m_mcDeltaTau
generated DeltaT
double m_fitPval
P value of the tag side fit result.
StoreArray< TagVertex > m_verArray
StoreArray of TagVertexes.
DBObjPtr< BeamSpot > m_beamSpotDB
Beam spot database object.
void deltaT(const Particle *Breco)
calculate DeltaT and MC-DeltaT (rec - tag) in ps from Breco and Btag vertices DT = Dl / gamma beta c ...
std::vector< const Particle * > m_raveParticles
tracks given to rave for the track fit (after removing Kshorts
void fillParticles(std::vector< ParticleAndWeight > &particleAndWeights)
Fill sorted list of particles into external variable.
TagVertexModule()
Constructor.
int m_reqPXDHits
N of PXD hits for a track to be used.
double m_confidenceLevel
required fit confidence level
double m_tagVChi2IP
IP component of the chi^2 of the tag vertex fit result.
double m_tagVol
tagV component in the direction orthogonal to the boost
analysis::VertexFitKFit doSingleKfit(std::vector< ParticleAndWeight > &particleAndWeights)
performs single KFit on particles stored in particleAndWeights this function can be iterated several ...
std::string m_constraintType
Choose constraint: noConstraint, IP, tube, boost, (breco)
void resetReturnParams()
Reset all parameters that are computed in each event and then used to compute tuple variables.
ROOT::Math::PxPyPzEVector m_tagMomentum
B tag momentum computed from fully reconstructed B sig.
bool makeGeneralFitRave()
make the vertex fit on the tag side: RAVE AVF tracks coming from Ks removed all other tracks used
TrackFitResult getTrackWithTrueCoordinates(ParticleAndWeight const &paw) const
If the fit has to be done with the truth info, Rave is fed with a track where the momentum is replace...
StoreArray< MCParticle > m_mcParticles
StoreArray of MCParticles.
double m_mcLifeTimeReco
generated Breco life time
static bool compBrecoBgen(const Particle *Breco, const MCParticle *Bgen)
compare Breco with the two MC B particles
double m_Bfield
magnetic field from data base
TMatrixDSym m_tagVErrMatrix
Error matrix of the tag side fit result.
int m_mcPDG
generated tag side B flavor
StoreObjPtr< ParticleList > m_plist
input particle list
static std::string printMatrix(const TMatrixD &mat)
Print a TMatrix (useful for debugging)
ROOT::Math::XYZVector m_mcTagV
generated tag side vertex
int m_rollbackStatus
Store info about whether the fit was performed with the rolled back tracks 0 fit performed with measu...
std::string m_trackFindingType
Choose how to find the tag tracks: standard, standard_PXD.
bool m_verbose
choose if you want to print extra infos
int m_FitType
fit algo used
bool makeGeneralFitKFit()
make the vertex fit on the tag side: KFit tracks coming from Ks removed all other tracks used
std::vector< const MCParticle * > m_raveMCParticles
Store the MC particles corresponding to each track used by Rave in the vtx fit.
TMatrixDSym m_BeamSpotCov
size of the beam spot == covariance matrix on the beam spot position
double m_deltaT
reconstructed DeltaT
TagVertex data object: contains Btag Vertex and DeltaT.
Definition: TagVertex.h:34
void setConstraintType(const std::string &constraintType)
Set the type of the constraint for the tag fit.
Definition: TagVertex.cc:312
void setTagVlErr(float TagVlErr)
Set the error of the tagV component in the boost direction.
Definition: TagVertex.cc:250
void setTruthTagVl(float TruthTagVl)
Set the MC tagV component in the boost direction.
Definition: TagVertex.cc:245
void setTagVertex(const B2Vector3D &TagVertex)
Set BTag Vertex.
Definition: TagVertex.cc:185
void setConstraintCenter(const B2Vector3D &constraintCenter)
Set the centre of the constraint for the tag fit.
Definition: TagVertex.cc:301
void setTruthTagVol(float TruthTagVol)
Set the tagV component in the direction orthogonal to the boost.
Definition: TagVertex.cc:260
void setMCTagBFlavor(int mcTagBFlavor)
Set generated Btag PDG code.
Definition: TagVertex.cc:215
void setTagVolErr(float TagVolErr)
Set the error of the tagV component in the direction orthogonal to the boost.
Definition: TagVertex.cc:265
void setMCTagVertex(const B2Vector3D &mcTagVertex)
Set generated BTag Vertex.
Definition: TagVertex.cc:210
void setTagVNDF(float TagVNDF)
Set the number of degrees of freedom in the tag vertex fit.
Definition: TagVertex.cc:270
void setDeltaTErr(float DeltaTErr)
Set DeltaTErr.
Definition: TagVertex.cc:205
void setNTracks(int nTracks)
Set number of tracks used in the fit.
Definition: TagVertex.cc:235
void setTagVChi2(float TagVChi2)
Set the chi^2 value of the tag vertex fit result.
Definition: TagVertex.cc:275
void setMCDeltaT(float mcDeltaT)
Set generated DeltaT (in kin.
Definition: TagVertex.cc:225
void setRollBackStatus(int backStatus)
Set the status of the fit performed with the rolled back tracks.
Definition: TagVertex.cc:336
void setVertexFitMCParticles(const std::vector< const MCParticle * > &vtxFitMCParticles)
Set a vector of pointers to the MC p'cles corresponding to the tracks in the tag vtx fit.
Definition: TagVertex.cc:291
void setTagVol(float TagVol)
Set the tagV component in the direction orthogonal to the boost.
Definition: TagVertex.cc:255
void setDeltaT(float DeltaT)
Set DeltaT.
Definition: TagVertex.cc:200
void setRaveWeights(const std::vector< double > &raveWeights)
Set the weights used by Rave in the tag vtx fit.
Definition: TagVertex.cc:296
void setTagVertexPval(float TagVertexPval)
Set BTag Vertex P value.
Definition: TagVertex.cc:195
void setMCDeltaTau(float mcDeltaTau)
Set generated DeltaT.
Definition: TagVertex.cc:220
void setTagVl(float TagVl)
Set the tagV component in the boost direction.
Definition: TagVertex.cc:240
void setTagVertexErrMatrix(const TMatrixDSym &TagVertexErrMatrix)
Set BTag Vertex (3x3) error matrix.
Definition: TagVertex.cc:190
void setConstraintCov(const TMatrixDSym &constraintCov)
Set the covariance matrix of the constraint for the tag fit.
Definition: TagVertex.cc:306
void setFitType(float FitType)
Set fit algo type.
Definition: TagVertex.cc:230
void setVertexFitParticles(const std::vector< const Particle * > &vtxFitParticles)
Set a vector of pointers to the tracks used in the tag vtx fit.
Definition: TagVertex.cc:285
void setTagVChi2IP(float TagVChi2IP)
Set the IP component of the chi^2 value of the tag vertex fit result.
Definition: TagVertex.cc:280
void setFitTruthStatus(int truthStatus)
Set the status of the fit performed with the truth info of the tracks.
Definition: TagVertex.cc:331
Values of the result of a track fit with a given particle hypothesis.
float getNDF() const
Getter for number of degrees of freedom of the track fit.
short getChargeSign() const
Return track charge (1 or -1).
double getPValue() const
Getter for Chi2 Probability of the track fit.
TMatrixDSym getCovariance6() const
Position and Momentum Covariance Matrix.
Const::ParticleType getParticleType() const
Getter for ParticleType of the mass hypothesis of the track fit.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
ROOT::Math::XYZVector getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
The Unit class.
Definition: Unit.h:40
virtual int getNDF(void) const
Get an NDF of the fit.
Definition: KFitBase.cc:114
bool isFitted(void) const
Return false if fit is not performed yet or performed fit is failed; otherwise true.
Definition: KFitBase.cc:728
int getTrackCount(void) const
Get the number of added tracks.
Definition: KFitBase.cc:107
void unsetBeamSpot()
unset beam spot constraint
Definition: RaveSetup.cc:80
static void initialize(int verbosity=1, double MagneticField=1.5)
Set everything up so everything needed for vertex fitting is there.
Definition: RaveSetup.cc:33
static RaveSetup * getInstance()
get the pointer to the instance to get/set any of options stored in RaveSetup
Definition: RaveSetup.h:43
void setBeamSpot(const B2Vector3D &beamSpot, const TMatrixDSym &beamSpotCov)
The beam spot position and covariance is known you can set it here so that and a vertex in the beam s...
Definition: RaveSetup.cc:72
void reset()
frees memory allocated by initialize().
Definition: RaveSetup.cc:58
The RaveVertexFitter class is part of the RaveInterface together with RaveSetup.
TMatrixDSym getCov(VecSize vertexId=0) const
get the covariance matrix (3x3) of the of the fitted vertex position.
int fit(std::string options="default")
do the vertex fit with all tracks previously added with the addTrack or addMother function.
void addTrack(const Particle *const aParticlePtr)
add a track (in the format of a Belle2::Particle) to set of tracks that should be fitted to a vertex
double getNdf(VecSize vertexId=0) const
get the number of degrees of freedom (NDF) of the fitted vertex.
double getChi2(VecSize vertexId=0) const
get the χ² of the fitted vertex.
B2Vector3D getPos(VecSize vertexId=0) const
get the position of the fitted vertex.
void updateDaughters()
update the Daughters particles
double getWeight(int trackId, VecSize vertexId=0) const
get the weight Rave assigned to a specific input track.
double getPValue(VecSize vertexId=0) const
get the p value of the fitted vertex.
VertexFitKFit is a derived class from KFitBase to perform vertex-constraint kinematical fit.
Definition: VertexFitKFit.h:34
double getCHIsq(void) const override
Get a chi-square of the fit.
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:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
ROOT::Math::XYZVector poca(ROOT::Math::XYZVector const &trackPos, ROOT::Math::XYZVector const &trackP, ROOT::Math::XYZVector const &vtxPos)
Returns the Point Of Closest Approach of a track to a vertex.
Particle * copyParticle(const Particle *original)
Function takes argument Particle and creates a copy of it and copies of all its (grand-)^n-daughters.
Definition: ParticleCopy.cc:18
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
this struct is used to store and sort the tag tracks
const Particle * particle
tag track fit result with pion mass hypo, for sorting purposes
const MCParticle * mcParticle
mc particle matched to the tag track, for sorting purposes
double weight
rave weight associated to the track, for sorting purposes