Belle II Software  light-2205-abys
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 
45 using namespace std;
46 using namespace Belle2;
47 
49 static const double realNaN = std::numeric_limits<double>::quiet_NaN();
51 static const B2Vector3D vecNaN(realNaN, realNaN, realNaN);
53 static const TMatrixDSym matNaN(3, (double [])
54 {
55  realNaN, realNaN, realNaN,
56  realNaN, realNaN, realNaN,
57  realNaN, realNaN, realNaN
58 });
59 
60 // import tools from RotationTools.h
61 using RotationTools::rotateTensor;
62 using RotationTools::rotateTensorInv;
63 using RotationTools::toSymMatrix;
64 using RotationTools::toVec;
65 using RotationTools::getUnitOrthogonal;
66 
67 //-----------------------------------------------------------------
68 // Register the Module
69 //-----------------------------------------------------------------
71 
72 //-----------------------------------------------------------------
73 // Implementation
74 //-----------------------------------------------------------------
75 
76 TagVertexModule::TagVertexModule() : Module(),
77  m_Bfield(0), m_fitTruthStatus(0), m_rollbackStatus(0), m_fitPval(0), m_mcTagLifeTime(-1), m_mcPDG(0), m_mcLifeTimeReco(-1),
78  m_deltaT(0), m_deltaTErr(0), m_mcDeltaTau(0), m_mcDeltaT(0),
79  m_shiftZ(0), m_FitType(0), m_tagVl(0),
80  m_truthTagVl(0), m_tagVlErr(0), m_tagVol(0), m_truthTagVol(0), m_tagVolErr(0), m_tagVNDF(0), m_tagVChi2(0), m_tagVChi2IP(0),
81  m_verbose(true)
82 {
83  // Set module properties
84  setDescription("Tag side Vertex Fitter for modular analysis");
85 
86  // Parameter definitions
87  addParam("listName", m_listName, "name of particle list", string(""));
88  addParam("confidenceLevel", m_confidenceLevel,
89  "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",
90  0.001);
91  addParam("MCAssociation", m_useMCassociation,
92  "'': no MC association. breco: use standard Breco MC association. internal: use internal MC association", string("breco"));
93  addParam("constraintType", m_constraintType,
94  "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)",
95  string("IP"));
96  addParam("trackFindingType", m_trackFindingType,
97  "Choose how to reconstruct the tracks on the tag side: standard, standard_PXD",
98  string("standard_PXD"));
99  addParam("maskName", m_roeMaskName,
100  "Choose ROE mask to get particles from ", string(RestOfEvent::c_defaultMaskName));
101  addParam("askMCInformation", m_mcInfo,
102  "TRUE when requesting MC Information from the tracks performing the vertex fit", false);
103  addParam("reqPXDHits", m_reqPXDHits,
104  "Minimum number of PXD hits for a track to be used in the vertex fit", 0);
105  addParam("fitAlgorithm", m_fitAlgo,
106  "Fitter used for the tag vertex fit: Rave or KFit", string("Rave"));
107  addParam("useTruthInFit", m_useTruthInFit,
108  "Use the true track parameters in the vertex fit", false);
109  addParam("useRollBack", m_useRollBack,
110  "Use rolled back non-primary tracks", false);
111 }
112 
114 {
115  // magnetic field
117  // RAVE setup
119  B2INFO("TagVertexModule : magnetic field = " << m_Bfield);
120  // truth fit status will be set to 2 only if the MC info cannot be recovered
122  // roll back status will be set to 2 only if the MC info cannot be recovered
124 
125  //input
127  m_plist.isRequired(m_listName);
128  // output
129  m_verArray.registerInDataStore();
131  //check if the fitting algorithm name is set correctly
132  if (m_fitAlgo != "Rave" && m_fitAlgo != "KFit")
133  B2FATAL("TagVertexModule: invalid fitting algorithm (must be set to either Rave or KFit).");
135  B2FATAL("TagVertexModule: invalid fitting option (useRollBack and useTruthInFit cannot be simultaneously set to true).");
136  //temporary while the one track fit is broken
137  if (m_trackFindingType == "singleTrack" || m_trackFindingType == "singleTrack_PXD")
138  B2FATAL("TagVertexModule : the singleTrack option is temporarily broken.");
139 }
140 
142 {
143  //TODO: set magnetic field for each run
144  //m_Bfield = BFieldMap::Instance().getBField(m_BeamSpotCenter).Z();
145 }
146 
148 {
149  if (!m_plist) {
150  B2ERROR("TagVertexModule: ParticleList " << m_listName << " not found");
151  return;
152  }
153 
154  // output
156 
157  std::vector<unsigned int> toRemove;
158 
159  for (unsigned i = 0; i < m_plist->getListSize(); ++i) {
161 
162  Particle* particle = m_plist->getParticle(i);
163  if (m_useMCassociation == "breco" || m_useMCassociation == "internal") BtagMCVertex(particle);
164  bool ok = doVertexFit(particle);
165  if (ok) deltaT(particle);
166 
169  toRemove.push_back(particle->getArrayIndex());
170  } else {
171  // save information in the Vertex StoreArray
172  TagVertex* ver = m_verArray.appendNew();
173  // create relation: Particle <-> Vertex
174  particle->addRelationTo(ver);
175  // fill Vertex with content
176  if (ok) {
177  ver->setTagVertex(m_tagV);
180  ver->setDeltaT(m_deltaT);
182  ver->setMCTagVertex(m_mcTagV);
183  ver->setMCTagBFlavor(m_mcPDG);
185  ver->setMCDeltaT(m_mcDeltaT);
186  ver->setFitType(m_FitType);
187  ver->setNTracks(m_tagParticles.size());
188  ver->setTagVl(m_tagVl);
190  ver->setTagVlErr(m_tagVlErr);
191  ver->setTagVol(m_tagVol);
194  ver->setTagVNDF(m_tagVNDF);
195  ver->setTagVChi2(m_tagVChi2);
205  } else {
206  ver->setTagVertex(m_tagV);
207  ver->setTagVertexPval(-1.);
208  ver->setDeltaT(m_deltaT);
210  ver->setMCTagVertex(m_mcTagV);
211  ver->setMCTagBFlavor(0.);
213  ver->setMCDeltaT(m_mcDeltaT);
214  ver->setFitType(m_FitType);
215  ver->setNTracks(m_tagParticles.size());
216  ver->setTagVl(m_tagVl);
218  ver->setTagVlErr(m_tagVlErr);
219  ver->setTagVol(m_tagVol);
222  ver->setTagVNDF(-1111.);
223  ver->setTagVChi2(-1111.);
224  ver->setTagVChi2IP(-1111.);
233  }
234  }
235  }
236  m_plist->removeParticles(toRemove);
237 
238  //free memory allocated by rave. initialize() would be enough, except that we must clean things up before program end...
239  //
241 }
242 
244 {
245  //reset the fit truth status in case it was set to 2 in a previous fit
246 
248 
249  //reset the roll back status in case it was set to 2 in a previous fit
250 
252 
253  //set constraint type, reset pVal and B field
254 
255  m_fitPval = 1;
256 
257  if (!(Breco->getRelatedTo<RestOfEvent>())) {
258  m_FitType = -1;
259  return false;
260  }
261 
262  if (m_Bfield == 0) {
263  B2ERROR("TagVertex: No magnetic field");
264  return false;
265  }
266 
267  // recover beam spot info
268 
269  m_BeamSpotCenter = m_beamSpotDB->getIPPosition();
270  m_BeamSpotCov.ResizeTo(3, 3);
271  m_BeamSpotCov = m_beamSpotDB->getCovVertex();
272 
273  //make the beam spot bigger for the standard constraint
274 
275  double beta = PCmsLabTransform().getBoostVector().Mag();
276  double bg = beta / sqrt(1 - beta * beta);
277 
278  //TODO: What's the origin of these numbers?
279  double cut = 8.717575e-02 * bg;
280  m_shiftZ = 4.184436e+02 * bg * 0.0001;
281 
282  //tube length here set to 20 * 2 * c tau beta gamma ~= 0.5 cm, should be enough to not bias the decay
283  //time but should still help getting rid of some pions from kshorts
284  m_constraintCov.ResizeTo(3, 3);
286  else if (m_constraintType == "tube") tie(m_constraintCenter, m_constraintCov) = findConstraintBTube(Breco, 1000 * cut);
287  else if (m_constraintType == "boost") tie(m_constraintCenter, m_constraintCov) = findConstraintBoost(cut * 200000.);
288  else if (m_constraintType == "breco") tie(m_constraintCenter, m_constraintCov) = findConstraint(Breco, cut * 2000.);
289  else if (m_constraintType == "noConstraint") m_constraintCenter = B2Vector3D(); //zero vector
290  else {
291  B2ERROR("TagVertex: Invalid constraintType selected");
292  return false;
293  }
294 
295  if (m_constraintCenter == vecNaN) {
296  B2ERROR("TagVertex: No correct fit constraint");
297  return false;
298  }
299 
300  /* 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
301  high efficiency, the next algorithm less restrictive is used. I.e, if standard_PXD does not work, the program tries with standard.
302  */
303 
304  m_FitType = 0;
305  double minPVal = (m_fitAlgo != "KFit") ? 0.001 : 0.;
306  bool ok = false;
307 
308  if (m_trackFindingType == "standard_PXD") {
310  if (m_tagParticles.size() > 0) {
311  ok = makeGeneralFit();
312  m_FitType = 3;
313  }
314  }
315 
316  if (ok == false || m_fitPval < minPVal || m_trackFindingType == "standard") {
318  ok = m_tagParticles.size() > 0;
319  if (ok) {
320  ok = makeGeneralFit();
321  m_FitType = 4;
322  }
323  }
324 
325  if ((ok == false || (m_fitPval <= 0. && m_fitAlgo == "Rave")) && m_constraintType != "noConstraint") {
327  ok = (m_constraintCenter != vecNaN);
328  if (ok) {
330  ok = (m_tagParticles.size() > 0);
331  }
332  if (ok) {
333  ok = makeGeneralFit();
334  m_FitType = 5;
335  }
336  }
337 
338  return ok;
339 }
340 
341 pair<B2Vector3D, TMatrixDSym> TagVertexModule::findConstraint(const Particle* Breco, double cut) const
342 {
343  if (Breco->getPValue() < 0.) return make_pair(vecNaN, matNaN);
344 
345  TMatrixDSym beamSpotCov(3);
346  beamSpotCov = m_beamSpotDB->getCovVertex();
347 
349 
350  double pmag = Breco->getMomentumMagnitude();
351  double xmag = (B2Vector3D(Breco->getVertex()) - m_BeamSpotCenter).Mag();
352 
353 
354  TMatrixDSym TerrMatrix = Breco->getMomentumVertexErrorMatrix();
355  TMatrixDSym PerrMatrix(7);
356 
357  for (int i = 0; i < 3; ++i) {
358  for (int j = 0; j < 3; ++j) {
359  if (i == j) {
360  PerrMatrix(i, j) = (beamSpotCov(i, j) + TerrMatrix(i, j)) * pmag / xmag;
361  } else {
362  PerrMatrix(i, j) = TerrMatrix(i, j);
363  }
364  PerrMatrix(i + 4, j + 4) = TerrMatrix(i + 4, j + 4);
365  }
366  }
367 
368  PerrMatrix(3, 3) = 0.;
369 
370  //Copy Breco, but use errors as are in PerrMatrix
371  Particle* Breco2 = ParticleCopy::copyParticle(Breco);
372  Breco2->setMomentumVertexErrorMatrix(PerrMatrix);
373 
374 
375  const Particle* BRecoRes = doVertexFitForBTube(Breco2, "kalman");
376  if (BRecoRes->getPValue() < 0) return make_pair(vecNaN, matNaN); //problems
377 
378  // Overall error matrix
379  TMatrixDSym errFinal = TMatrixDSym(Breco->getVertexErrorMatrix() + BRecoRes->getVertexErrorMatrix());
380 
381  // TODO : to be developed the extraction of the momentum from the rave fitted track
382 
383  // Get expected pBtag 4-momentum using transverse-momentum conservation
384  ROOT::Math::XYZVector BvertDiff = pmag * (Breco->getVertex() - BRecoRes->getVertex()).Unit();
385  ROOT::Math::PxPyPzMVector pBrecEstimate(BvertDiff.x(), BvertDiff.y(), BvertDiff.z(), Breco->getPDGMass());
386  ROOT::Math::PxPyPzMVector pBtagEstimate = PCmsLabTransform::labToCms(pBrecEstimate);
387  pBtagEstimate.SetPxPyPzE(-pBtagEstimate.px(), -pBtagEstimate.py(), -pBtagEstimate.pz(), pBtagEstimate.E());
388  pBtagEstimate = PCmsLabTransform::cmsToLab(pBtagEstimate);
389 
390  // rotate err-matrix such that pBrecEstimate goes to eZ
391  TMatrixD TubeZ = rotateTensorInv(pBrecEstimate.Vect(), errFinal);
392 
393  TubeZ(2, 2) = cut * cut;
394  TubeZ(2, 0) = 0; TubeZ(0, 2) = 0;
395  TubeZ(2, 1) = 0; TubeZ(1, 2) = 0;
396 
397 
398  // rotate err-matrix such that eZ goes to pBtagEstimate
399  TMatrixD Tube = rotateTensor(B2Vector3D(pBtagEstimate.Vect()), TubeZ);
400 
401  // Standard algorithm needs no shift
402  return make_pair(m_BeamSpotCenter, toSymMatrix(Tube));
403 }
404 
405 pair<B2Vector3D, TMatrixDSym> TagVertexModule::findConstraintBTube(const Particle* Breco, double cut)
406 {
407  //Use Breco as the creator of the B tube.
408  if ((Breco->getVertexErrorMatrix()(2, 2)) == 0.0) {
409  B2WARNING("In TagVertexModule::findConstraintBTube: cannot get a proper vertex for BReco. BTube constraint replaced by Boost.");
410  return findConstraintBoost(cut);
411  }
412 
413 
414  //vertex fit will give the intersection between the beam spot and the trajectory of the B
415  //(base of the BTube, or primary vtx cov matrix)
416  const Particle* tubecreatorBCopy = doVertexFitForBTube(Breco, "avf");
417  if (tubecreatorBCopy->getPValue() < 0) return make_pair(vecNaN, matNaN); //if problems
418 
419 
420  //get direction of B tag = opposite direction of B rec in CMF
421  ROOT::Math::PxPyPzEVector pBrec = tubecreatorBCopy->get4Vector();
422 
423  //if we want the true info, replace the 4vector by the true one
424  if (m_useTruthInFit) {
425  const MCParticle* mcBr = Breco->getRelated<MCParticle>();
426  if (mcBr)
427  pBrec = mcBr->get4Vector();
428  else
429  m_fitTruthStatus = 2;
430  }
431  ROOT::Math::PxPyPzEVector pBtag = PCmsLabTransform::labToCms(pBrec);
432  pBtag.SetPxPyPzE(-pBtag.px(), -pBtag.py(), -pBtag.pz(), pBtag.E());
433  pBtag = PCmsLabTransform::cmsToLab(pBtag);
434 
435  //To create the B tube, strategy is: take the primary vtx cov matrix, and add to it a cov
436  //matrix corresponding to an very big error in the direction of the B tag
437  TMatrixDSym pv = tubecreatorBCopy->getVertexErrorMatrix();
438 
439  //print some stuff if wanted
440  if (m_verbose) {
441  B2DEBUG(10, "Brec decay vertex before fit: " << printVector(Breco->getVertex()));
442  B2DEBUG(10, "Brec decay vertex after fit: " << printVector(tubecreatorBCopy->getVertex()));
443  B2DEBUG(10, "Brec direction before fit: " << printVector(float(1. / Breco->getP()) * Breco->getMomentum()));
444  B2DEBUG(10, "Brec direction after fit: " << printVector(float(1. / tubecreatorBCopy->getP()) * tubecreatorBCopy->getMomentum()));
445  B2DEBUG(10, "IP position: " << printVector(ROOT::Math::XYZVector(m_BeamSpotCenter.x(), m_BeamSpotCenter.y(),
446  m_BeamSpotCenter.z())));
447  B2DEBUG(10, "IP covariance: " << printMatrix(m_BeamSpotCov));
448  B2DEBUG(10, "Brec primary vertex: " << printVector(tubecreatorBCopy->getVertex()));
449  B2DEBUG(10, "Brec PV covariance: " << printMatrix(pv));
450  B2DEBUG(10, "BTag direction: " << printVector((1. / pBtag.P())*pBtag.Vect()));
451  }
452 
453  //make a long error matrix along BTag direction
454  TMatrixD longerror(3, 3); longerror(2, 2) = cut * cut;
455 
456 
457  // make rotation matrix from z axis to BTag line of flight
458  TMatrixD longerrorRotated = rotateTensor(B2Vector3D(pBtag.Vect()), longerror);
459 
460  //pvNew will correspond to the covariance matrix of the B tube
461  TMatrixD pvNew = TMatrixD(pv) + longerrorRotated;
462 
463  //set the constraint
464  ROOT::Math::XYZVector constraintCenter = tubecreatorBCopy->getVertex();
465 
466  //if we want the true info, set the centre of the constraint to the primary vertex
467  if (m_useTruthInFit) {
468  const MCParticle* mcBr = Breco->getRelated<MCParticle>();
469  if (mcBr) {
470  constraintCenter = mcBr->getProductionVertex();
471  }
472  }
473 
474  if (m_verbose) {
475  B2DEBUG(10, "IPTube covariance: " << printMatrix(pvNew));
476  }
477 
478  //The following is done to do the BTube constraint with a virtual track
479  //(ie KFit way)
480 
481  m_tagMomentum = pBtag;
482 
483  m_pvCov.ResizeTo(pv);
484  m_pvCov = pv;
485 
486  return make_pair(B2Vector3D(constraintCenter), toSymMatrix(pvNew));
487 }
488 
489 pair<B2Vector3D, TMatrixDSym> TagVertexModule::findConstraintBoost(double cut, double shiftAlongBoost) const
490 {
491  //make a long error matrix along boost direction
492  TMatrixD longerror(3, 3); longerror(2, 2) = cut * cut;
494  TMatrixD longerrorRotated = rotateTensor(boostDir, longerror);
495 
496  //Extend error of BeamSpotCov matrix in the boost direction
497  TMatrixDSym beamSpotCov = m_beamSpotDB->getCovVertex();
498  TMatrixD Tube = TMatrixD(beamSpotCov) + longerrorRotated;
499 
500  // Standard algorithm needs no shift
501  B2Vector3D constraintCenter = m_BeamSpotCenter;
502 
503  // The constraint used in the Single Track Fit needs to be shifted in the boost direction.
504  if (shiftAlongBoost > -1000) {
505  constraintCenter += shiftAlongBoost * boostDir;
506  }
507 
508  return make_pair(constraintCenter, toSymMatrix(Tube));
509 }
510 
512 static double getProperLifeTime(const MCParticle* mc)
513 {
514  double beta = mc->getMomentum().Mag() / mc->getEnergy();
515  return 1e3 * mc->getLifetime() * sqrt(1 - pow(beta, 2));
516 }
517 
519 {
520  //fill vector with mcB (intended order: Reco, Tag)
521  vector<const MCParticle*> mcBs;
522  for (const MCParticle& mc : m_mcParticles) {
523  if (abs(mc.getPDG()) == abs(Breco->getPDGCode()))
524  mcBs.push_back(&mc);
525  }
526  //too few Bs
527  if (mcBs.size() < 2) return;
528 
529  if (mcBs.size() > 2) {
530  B2WARNING("TagVertexModule:: Too many Bs found in MC");
531  }
532 
533  auto isReco = [&](const MCParticle * mc) {
534  return (m_useMCassociation == "breco") ? (mc == Breco->getRelated<MCParticle>())
535  : compBrecoBgen(Breco, mc); //internal association
536  };
537 
538  //nothing matched?
539  if (!isReco(mcBs[0]) && !isReco(mcBs[1])) {
540  return;
541  }
542 
543  //first is Tag, second Reco -> swap the order
544  if (!isReco(mcBs[0]) && isReco(mcBs[1]))
545  swap(mcBs[0], mcBs[1]);
546 
547  //both matched -> use closest vertex dist as Reco
548  if (isReco(mcBs[0]) && isReco(mcBs[1])) {
549  double dist0 = (mcBs[0]->getDecayVertex() - B2Vector3D(Breco->getVertex())).Mag2();
550  double dist1 = (mcBs[1]->getDecayVertex() - B2Vector3D(Breco->getVertex())).Mag2();
551  if (dist0 > dist1)
552  swap(mcBs[0], mcBs[1]);
553  }
554 
555  m_mcVertReco = mcBs[0]->getDecayVertex();
556  m_mcLifeTimeReco = getProperLifeTime(mcBs[0]);
557  m_mcTagV = mcBs[1]->getDecayVertex();
558  m_mcTagLifeTime = getProperLifeTime(mcBs[1]);
559  m_mcPDG = mcBs[1]->getPDG();
560 }
561 
562 // static
563 bool TagVertexModule::compBrecoBgen(const Particle* Breco, const MCParticle* Bgen)
564 {
565 
566  bool isDecMode = true;
567 
568  const std::vector<Belle2::Particle*> recDau = Breco->getDaughters();
569  const std::vector<Belle2::MCParticle*> genDau = Bgen->getDaughters();
570 
571  if (recDau.size() > 0 && genDau.size() > 0) {
572  for (auto dauRec : recDau) {
573  bool isDau = false;
574  for (auto dauGen : genDau) {
575  if (dauGen->getPDG() == dauRec->getPDGCode())
576  isDau = compBrecoBgen(dauRec, dauGen) ;
577  }
578  if (!isDau) isDecMode = false;
579  }
580  } else {
581  if (recDau.size() == 0) { //&& genDau.size()==0){
582  if (Bgen->getPDG() != Breco->getPDGCode()) isDecMode = false;;
583  } else {isDecMode = false;}
584  }
585 
586  return isDecMode;
587 }
588 
589 // STANDARD FIT ALGORITHM
590 /* This algorithm basically takes all the tracks coming from the Rest Of Events and send them to perform a multi-track fit
591  The option to request PXD hits for the tracks can be chosen by the user.
592  */
593 std::vector<const Particle*> TagVertexModule::getTagTracks_standardAlgorithm(const Particle* Breco, int reqPXDHits) const
594 {
595  std::vector<const Particle*> fitParticles;
596  const RestOfEvent* roe = Breco->getRelatedTo<RestOfEvent>();
597  if (!roe) return fitParticles;
598  //load all particles from the ROE
599  std::vector<const Particle*> ROEParticles = roe->getChargedParticles(m_roeMaskName, 0, false);
600  if (ROEParticles.size() == 0) return fitParticles;
601 
602  for (auto& ROEParticle : ROEParticles) {
603  HitPatternVXD roeTrackPattern = ROEParticle->getTrackFitResult()->getHitPatternVXD();
604 
605  if (roeTrackPattern.getNPXDHits() >= reqPXDHits) {
606  fitParticles.push_back(ROEParticle);
607  }
608  }
609  return fitParticles;
610 }
611 
612 vector<ParticleAndWeight> TagVertexModule::getParticlesWithoutKS(const vector<const Particle*>& tagParticles,
613  double massWindowWidth) const
614 {
615  vector<ParticleAndWeight> particleAndWeights;
616 
617  ParticleAndWeight particleAndWeight;
618  particleAndWeight.mcParticle = 0;
619  particleAndWeight.weight = -1111.;
620 
621  // remove tracks from KS
622  for (unsigned i = 0; i < tagParticles.size(); ++i) {
623  const Particle* particle1 = tagParticles.at(i);
624  if (!particle1) continue;
625  ROOT::Math::PxPyPzEVector mom1 = particle1->get4Vector();
626  if (!isfinite(mom1.mag2())) continue;
627 
628  //is from Ks decay?
629  bool isKsDau = false;
630  for (unsigned j = 0; j < tagParticles.size(); ++j) {
631  if (i == j) continue;
632  const Particle* particle2 = tagParticles.at(j);
633  if (!particle2) continue;
634  ROOT::Math::PxPyPzEVector mom2 = particle2->get4Vector();
635  if (!isfinite(mom2.mag2())) continue;
636  double mass = (mom1 + mom2).M();
637  if (abs(mass - Const::K0Mass) < massWindowWidth) {
638  isKsDau = true;
639  break;
640  }
641  }
642  //if from Ks decay, skip
643  if (isKsDau) continue;
644 
645  particleAndWeight.particle = particle1;
646 
647  if (m_useMCassociation == "breco" || m_useMCassociation == "internal")
648  particleAndWeight.mcParticle = particle1->getRelatedTo<MCParticle>();
649 
650  particleAndWeights.push_back(particleAndWeight);
651 
652  }
653 
654  return particleAndWeights;
655 }
656 
658 {
659  if (m_fitAlgo == "Rave") return makeGeneralFitRave();
660  else if (m_fitAlgo == "KFit") return makeGeneralFitKFit();
661  return false;
662 }
663 
664 void TagVertexModule::fillParticles(vector<ParticleAndWeight>& particleAndWeights)
665 {
666  unsigned n = particleAndWeights.size();
667  sort(particleAndWeights.begin(), particleAndWeights.end(),
668  [](const ParticleAndWeight & a, const ParticleAndWeight & b) { return a.weight > b.weight; });
669 
670  m_raveParticles.resize(n);
671  m_raveWeights.resize(n);
672  m_raveMCParticles.resize(n);
673 
674  for (unsigned i = 0; i < n; ++i) {
675  m_raveParticles.at(i) = particleAndWeights.at(i).particle;
676  m_raveMCParticles.at(i) = particleAndWeights.at(i).mcParticle;
677  m_raveWeights.at(i) = particleAndWeights.at(i).weight;
678  }
679 }
680 
681 void TagVertexModule::fillTagVinfo(const B2Vector3D& tagVpos, const TMatrixDSym& tagVposErr)
682 {
683  m_tagV = tagVpos;
684 
685  if (m_constraintType != "noConstraint") {
686  TMatrixDSym tubeInv = m_constraintCov;
687  tubeInv.Invert();
688  TVectorD dV = toVec(m_tagV - m_BeamSpotCenter);
689  m_tagVChi2IP = tubeInv.Similarity(dV);
690  }
691 
692  m_tagVErrMatrix.ResizeTo(tagVposErr);
693  m_tagVErrMatrix = tagVposErr;
694 }
695 
697 {
698  // apply constraint
700  if (m_constraintType != "noConstraint")
703 
704  //feed rave with tracks without Kshorts
705  vector<ParticleAndWeight> particleAndWeights = getParticlesWithoutKS(m_tagParticles);
706 
707  for (const auto& pw : particleAndWeights) {
708  try {
709  if (m_useTruthInFit) {
710  if (pw.mcParticle) {
712  rFit.addTrack(&tfr);
713  } else
714  m_fitTruthStatus = 2;
715  } else if (m_useRollBack) {
716  if (pw.mcParticle) {
718  rFit.addTrack(&tfr);
719  } else
720  m_rollbackStatus = 2;
721  } else {
722  rFit.addTrack(pw.particle->getTrackFitResult());
723  }
724  } catch (const rave::CheckedFloatException&) {
725  B2ERROR("Exception caught in TagVertexModule::makeGeneralFitRave(): Invalid inputs (nan/inf)?");
726  }
727  }
728 
729  //perform fit
730 
731  int isGoodFit(-1);
732 
733  try {
734  isGoodFit = rFit.fit("avf");
735  // if problems
736  if (isGoodFit < 1) return false;
737  } catch (const rave::CheckedFloatException&) {
738  B2ERROR("Exception caught in TagVertexModule::makeGeneralFitRave(): Invalid inputs (nan/inf)?");
739  return false;
740  }
741 
742  //save the track info for later use
743 
744  for (unsigned int i(0); i < particleAndWeights.size() && isGoodFit >= 1; ++i)
745  particleAndWeights.at(i).weight = rFit.getWeight(i);
746 
747  //Tracks are sorted from highest rave weight to lowest
748 
749  fillParticles(particleAndWeights);
750 
751  //if the fit is good, save the infos related to the vertex
752  fillTagVinfo(rFit.getPos(0), rFit.getCov(0));
753 
754  //fill quality variables
755  m_tagVNDF = rFit.getNdf(0);
756  m_tagVChi2 = rFit.getChi2(0);
757  m_fitPval = rFit.getPValue();
758 
759  return true;
760 }
761 
763 {
764  //initialize KFit
766  kFit.setMagneticField(m_Bfield);
767 
768  // apply constraint
769  if (m_constraintType != "noConstraint") {
770  if (m_constraintType == "tube") {
771  CLHEP::HepSymMatrix err(7, 0);
772  //copy m_pvCov to the end of err matrix
773  err.sub(5, ROOTToCLHEP::getHepSymMatrix(m_pvCov));
774  kFit.setIpTubeProfile(
775  ROOTToCLHEP::getHepLorentzVector(m_tagMomentum),
776  ROOTToCLHEP::getPoint3DFromB2Vector(m_constraintCenter),
777  err,
778  0.);
779  } else {
780  kFit.setIpProfile(ROOTToCLHEP::getPoint3DFromB2Vector(m_constraintCenter),
781  ROOTToCLHEP::getHepSymMatrix(m_constraintCov));
782  }
783  }
784 
785  //feed KFit with tracks without Kshorts
786  vector<ParticleAndWeight> particleAndWeights = getParticlesWithoutKS(m_tagParticles);
787 
788  int nTracksAdded = 0;
789  for (auto& pawi : particleAndWeights) {
790  int addedOK = 1;
791  if (m_useTruthInFit) {
792  if (pawi.mcParticle) {
793  addedOK = kFit.addTrack(
794  ROOTToCLHEP::getHepLorentzVector(pawi.mcParticle->get4Vector()),
795  ROOTToCLHEP::getPoint3DFromB2Vector(getTruePoca(pawi)),
796  ROOTToCLHEP::getHepSymMatrix(pawi.particle->getMomentumVertexErrorMatrix()),
797  pawi.particle->getCharge());
798  } else {
799  m_fitTruthStatus = 2;
800  }
801  } else if (m_useRollBack) {
802  if (pawi.mcParticle) {
803  addedOK = kFit.addTrack(
804  ROOTToCLHEP::getHepLorentzVector(pawi.mcParticle->get4Vector()),
805  ROOTToCLHEP::getPoint3DFromB2Vector(getRollBackPoca(pawi)),
806  ROOTToCLHEP::getHepSymMatrix(pawi.particle->getMomentumVertexErrorMatrix()),
807  pawi.particle->getCharge());
808  } else {
809  m_rollbackStatus = 2;
810  }
811  } else {
812  addedOK = kFit.addParticle(pawi.particle);
813  }
814 
815  if (addedOK == 0) {
816  ++nTracksAdded;
817  pawi.weight = 1.;
818  } else {
819  B2WARNING("TagVertexModule::makeGeneralFitKFit: failed to add a track");
820  pawi.weight = 0.;
821  }
822  }
823 
824  //perform fit if there are enough tracks
825  if ((nTracksAdded < 2 && m_constraintType == "noConstraint") || nTracksAdded < 1)
826  return false;
827 
828  int isGoodFit = kFit.doFit();
829  if (isGoodFit != 0) return false;
830 
831  //save the track info for later use
832  //Tracks are sorted by weight, ie pushing the tracks with 0 weight (from KS) to the end of the list
833  fillParticles(particleAndWeights);
834 
835  //Save the infos related to the vertex
836  fillTagVinfo(CLHEPToROOT::getXYZVector(kFit.getVertex()),
837  CLHEPToROOT::getTMatrixDSym(kFit.getVertexError()));
838 
839  m_tagVNDF = kFit.getNDF();
840  m_tagVChi2 = kFit.getCHIsq();
841  m_fitPval = TMath::Prob(m_tagVChi2, m_tagVNDF);
842 
843  return true;
844 }
845 
847 {
848 
850  B2Vector3D boostDir = boost.Unit();
851  double bg = boost.Mag() / sqrt(1 - boost.Mag2());
852  double c = Const::speedOfLight / 1000.; // cm ps-1
853 
854  //Reconstructed DeltaL & DeltaT in the boost direction
855  B2Vector3D dVert = B2Vector3D(Breco->getVertex()) - m_tagV; //reconstructed vtxReco - vtxTag
856  double dl = dVert.Dot(boostDir);
857  m_deltaT = dl / (bg * c);
858 
859  //Truth DeltaL & approx DeltaT in the boost direction
860  B2Vector3D MCdVert = m_mcVertReco - m_mcTagV; //truth vtxReco - vtxTag
861  double MCdl = MCdVert.Dot(boostDir);
862  m_mcDeltaT = MCdl / (bg * c);
863 
864  // MCdeltaTau=tauRec-tauTag
866  if (m_mcLifeTimeReco == -1 || m_mcTagLifeTime == -1)
867  m_mcDeltaTau = realNaN;
868 
869  TVectorD bVec = toVec(boostDir);
870 
871  //TagVertex error in boost dir
872  m_tagVlErr = sqrt(m_tagVErrMatrix.Similarity(bVec));
873 
874  //bReco error in boost dir
875  double bRecoErrL = sqrt(Breco->getVertexErrorMatrix().Similarity(bVec));
876 
877  //Delta t error
878  m_deltaTErr = hypot(m_tagVlErr, bRecoErrL) / (bg * c);
879 
880  m_tagVl = m_tagV.Dot(boostDir);
881  m_truthTagVl = m_mcTagV.Dot(boostDir);
882 
883  // calculate tagV component and error in the direction orthogonal to the boost
884  B2Vector3D oboost = getUnitOrthogonal(boostDir);
885  TVectorD oVec = toVec(oboost);
886 
887  //TagVertex error in boost-orthogonal dir
888  m_tagVolErr = sqrt(m_tagVErrMatrix.Similarity(oVec));
889 
890  m_tagVol = m_tagV.Dot(oboost);
891  m_truthTagVol = m_mcTagV.Dot(oboost);
892 }
893 
894 Particle* TagVertexModule::doVertexFitForBTube(const Particle* motherIn, std::string fitType) const
895 {
896  //make a copy of motherIn to not modify the original object
897  Particle* mother = ParticleCopy::copyParticle(motherIn);
898 
899  //Here rave is used to find the upsilon(4S) vtx as the intersection
900  //between the mother B trajectory and the beam spot
902 
904  rsg.addTrack(mother);
905  int nvert = rsg.fit(fitType);
906  if (nvert != 1) {
907  mother->setPValue(-1); //error
908  return mother;
909  } else {
910  rsg.updateDaughters();
911  return mother;
912  }
913 }
914 
915 
916 
918 {
919  if (!paw.mcParticle) {
920  B2ERROR("In TagVertexModule::getTrackWithTrueCoordinate: no MC particle set");
921  return TrackFitResult();
922  }
923 
924  const TrackFitResult* tfr(paw.particle->getTrackFitResult());
925 
926  return TrackFitResult(getTruePoca(paw),
927  paw.mcParticle->getMomentum(),
928  tfr->getCovariance6(),
929  tfr->getChargeSign(),
930  tfr->getParticleType(),
931  tfr->getPValue(),
932  m_Bfield, 0, 0, tfr->getNDF());
933 }
934 
935 // static
937 {
938  if (!paw.mcParticle) {
939  B2ERROR("In TagVertexModule::getTruePoca: no MC particle set");
940  return B2Vector3D(0., 0., 0.);
941  }
942 
944  paw.mcParticle->getMomentum(),
946 }
947 
949 {
950  const TrackFitResult* tfr(paw.particle->getTrackFitResult());
951 
952  return TrackFitResult(getRollBackPoca(paw),
953  tfr->getMomentum(),
954  tfr->getCovariance6(),
955  tfr->getChargeSign(),
956  tfr->getParticleType(),
957  tfr->getPValue(),
958  m_Bfield, 0, 0, tfr->getNDF());
959 }
960 
962 {
963  if (!paw.mcParticle) {
964  B2ERROR("In TagVertexModule::getTruePoca: no MC particle set");
965  return B2Vector3D(0., 0., 0.);
966  }
967 
969 }
970 
972 {
973  m_raveParticles.resize(0);
974  m_raveMCParticles.resize(0);
975  m_tagParticles.resize(0);
976  m_raveWeights.resize(0);
977 
978  m_fitPval = realNaN;
979  m_tagV = vecNaN;
980  m_tagVErrMatrix.ResizeTo(matNaN);
981  m_tagVErrMatrix = matNaN;
982  m_mcTagV = vecNaN;
983  m_mcVertReco = vecNaN;
984  m_deltaT = realNaN;
985  m_deltaTErr = realNaN;
986  m_mcDeltaTau = realNaN;
987  m_constraintCov.ResizeTo(matNaN);
988  m_constraintCov = matNaN;
989  m_constraintCenter = vecNaN;
990  m_tagVl = realNaN;
991  m_truthTagVl = realNaN;
992  m_tagVlErr = realNaN;
993  m_tagVol = realNaN;
994  m_truthTagVol = realNaN;
995  m_tagVolErr = realNaN;
996  m_tagVNDF = realNaN;
997  m_tagVChi2 = realNaN;
998  m_tagVChi2IP = realNaN;
999  m_pvCov.ResizeTo(matNaN);
1000  m_pvCov = matNaN;
1001  m_tagMomentum = ROOT::Math::PxPyPzEVector(realNaN, realNaN, realNaN, realNaN);
1002 }
1003 
1004 //The following functions are just here to help printing stuff
1005 
1006 // static
1007 std::string TagVertexModule::printVector(const ROOT::Math::XYZVector& vec)
1008 {
1009  std::ostringstream oss;
1010  int w = 14;
1011  oss << "(" << std::setw(w) << vec.x() << ", " << std::setw(w) << vec.y() << ", " << std::setw(w) << vec.z() << ")" << std::endl;
1012  return oss.str();
1013 }
1014 
1015 // static
1016 std::string TagVertexModule::printMatrix(const TMatrixD& mat)
1017 {
1018  std::ostringstream oss;
1019  int w = 14;
1020  for (int i = 0; i < mat.GetNrows(); ++i) {
1021  for (int j = 0; j < mat.GetNcols(); ++j) {
1022  oss << std::setw(w) << mat(i, j) << " ";
1023  }
1024  oss << endl;
1025  }
1026  return oss.str();
1027 }
1028 
1029 // static
1030 std::string TagVertexModule::printMatrix(const TMatrixDSym& mat)
1031 {
1032  std::ostringstream oss;
1033  int w = 14;
1034  for (int i = 0; i < mat.GetNrows(); ++i) {
1035  for (int j = 0; j < mat.GetNcols(); ++j) {
1036  oss << std::setw(w) << mat(i, j) << " ";
1037  }
1038  oss << endl;
1039  }
1040  return oss.str();
1041 }
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:429
DataType y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:421
DataType z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:423
DataType Mag() const
The magnitude (rho in spherical coordinate system).
Definition: B2Vector3.h:153
DataType x() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:419
DataType Mag2() const
The magnitude squared (rho^2 in spherical coordinate system).
Definition: B2Vector3.h:151
B2Vector3< DataType > Unit() const
Unit vector parallel to this.
Definition: B2Vector3.h:263
DataType Dot(const B2Vector3< DataType > &p) const
Scalar product.
Definition: B2Vector3.h:284
static B2Vector3D getFieldInTesla(const B2Vector3D &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:70
static const double K0Mass
neutral kaon mass
Definition: Const.h:573
static const double speedOfLight
[cm/ns]
Definition: Const.h:575
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:34
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:191
TVector3 getMomentum() const
Return momentum.
Definition: MCParticle.h:200
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:209
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:114
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:74
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
Definition: Particle.cc:479
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:596
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:560
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:403
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:660
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:636
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:488
void setMomentumVertexErrorMatrix(const TMatrixFSym &errMatrix)
Sets 7x7 error matrix.
Definition: Particle.cc:425
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:497
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:331
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:452
const TrackFitResult * getTrackFitResult() const
Returns the pointer to the TrackFitResult that was used to create this Particle (ParticleType == c_Tr...
Definition: Particle.cc:846
double getMomentumMagnitude() const
Returns momentum magnitude.
Definition: Particle.h:506
double getP() const
Returns momentum magnitude (same as getMomentumMagnitude but with shorter name)
Definition: Particle.h:515
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.
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to 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_shiftZ
parameter for testing the systematic error from the IP measurement
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
B2Vector3D getRollBackPoca(ParticleAndWeight const &paw)
This shifts the position of tracks by the vector difference of mother B and production point of track...
bool m_useTruthInFit
Set to true if the tag fit is to be made with the TRUE tag track momentum and position.
void fillTagVinfo(const B2Vector3D &tagVpos, const TMatrixDSym &tagVposErr)
Fill tagV vertex info.
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.
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)
static B2Vector3D getTruePoca(ParticleAndWeight const &paw)
This finds the point on the true particle trajectory closest to the measured track position.
std::vector< ParticleAndWeight > getParticlesWithoutKS(const std::vector< const Particle * > &tagParticles, double massWindowWidth=0.01) const
Get a list of pions from a list of pions removing the Kshorts Warning: this assumes all the particles...
virtual void event() override
Event processor.
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::pair< B2Vector3D, TMatrixDSym > findConstraint(const Particle *Breco, double cut) const
calculate the constraint for the vertex fit on the tag side using Breco information
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
void BtagMCVertex(const Particle *Breco)
get the vertex of the MC B particle associated to Btag.
bool m_mcInfo
true if user wants to retrieve MC information out from the tracks used in the fit
B2Vector3D m_tagV
tag side fit result
std::string m_useMCassociation
No MC association or standard Breco particle or internal MCparticle association.
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
B2Vector3D 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.
virtual void beginRun() override
Called when entering a new run.
B2Vector3D m_constraintCenter
centre position of the constraint for the tag Vertex fit
double m_mcDeltaTau
generated DeltaT
std::pair< B2Vector3D, TMatrixDSym > findConstraintBoost(double cut, double shiftAlongBoost=-2000.) const
calculate the standard constraint for the vertex fit on the tag side
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.
B2Vector3D m_mcTagV
generated tag side vertex
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.
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
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...
std::pair< B2Vector3D, 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...
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)
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.
B2Vector3D m_mcVertReco
generated Breco decay vertex
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.
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.
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
int getNDF() const
Getter for number of degrees of freedom of the track fit.
TVector3 getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
Const::ParticleType getParticleType() const
Getter for ParticleType of the mass hypothesis of the track fit.
The Unit class.
Definition: Unit.h:40
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
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:502
B2Vector3D poca(B2Vector3D const &trackPos, B2Vector3D const &trackP, B2Vector3D 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:23
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