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