Belle II Software  release-05-01-25
MillepedeCollectorModule.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: tadeas *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <alignment/modules/MillepedeCollector/MillepedeCollectorModule.h>
12 
13 #include <alignment/dataobjects/MilleData.h>
14 #include <alignment/GblMultipleScatteringController.h>
15 #include <alignment/GlobalDerivatives.h>
16 #include <alignment/GlobalLabel.h>
17 #include <alignment/GlobalParam.h>
18 #include <alignment/GlobalTimeLine.h>
19 #include <alignment/Manager.h>
20 #include <alignment/reconstruction/AlignableCDCRecoHit.h>
21 #include <alignment/reconstruction/AlignablePXDRecoHit.h>
22 #include <alignment/reconstruction/AlignableSVDRecoHit.h>
23 #include <alignment/reconstruction/AlignableSVDRecoHit2D.h>
24 #include <alignment/reconstruction/AlignableBKLMRecoHit.h>
25 #include <alignment/reconstruction/AlignableEKLMRecoHit.h>
26 #include <analysis/dataobjects/ParticleList.h>
27 #include <analysis/utility/ReferenceFrame.h>
28 #include <framework/core/FileCatalog.h>
29 #include <framework/database/DBObjPtr.h>
30 #include <framework/dataobjects/EventT0.h>
31 #include <framework/dataobjects/FileMetaData.h>
32 #include <framework/datastore/StoreArray.h>
33 #include <framework/datastore/StoreObjPtr.h>
34 #include <framework/dbobjects/BeamParameters.h>
35 #include <framework/particledb/EvtGenDatabasePDG.h>
36 #include <framework/pcore/ProcHandler.h>
37 #include <mdst/dbobjects/BeamSpot.h>
38 #include <mdst/dataobjects/Track.h>
39 #include <tracking/trackFitting/fitter/base/TrackFitter.h>
40 #include <tracking/trackFitting/measurementCreator/adder/MeasurementAdder.h>
41 
42 #include <genfit/FullMeasurement.h>
43 #include <genfit/GblFitter.h>
44 #include <genfit/KalmanFitterInfo.h>
45 #include <genfit/PlanarMeasurement.h>
46 #include <genfit/Track.h>
47 
48 #include <TMath.h>
49 #include <TH1F.h>
50 #include <TTree.h>
51 #include <TDecompSVD.h>
52 
53 using namespace std;
54 using namespace Belle2;
55 using namespace alignment;
56 
57 //-----------------------------------------------------------------
58 // Register the Module
59 //-----------------------------------------------------------------
60 REG_MODULE(MillepedeCollector)
61 
62 //-----------------------------------------------------------------
63 // Implementation
64 //-----------------------------------------------------------------
65 
67 {
68  setPropertyFlags(c_ParallelProcessingCertified);
69  setDescription("Calibration data collector for Millepede Algorithm");
70 
71  // Configure input sample types
72  addParam("tracks", m_tracks, "Names of collections of RecoTracks (already fitted with DAF) for calibration", vector<string>({""}));
73  addParam("particles", m_particles, "Names of particle list of single particles", vector<string>());
74  addParam("vertices", m_vertices,
75  "Name of particle list of (mother) particles with daughters for calibration using vertex constraint", vector<string>());
76  addParam("primaryVertices", m_primaryVertices,
77  "Name of particle list of (mother) particles with daughters for calibration using vertex + IP profile constraint",
78  vector<string>());
79  addParam("twoBodyDecays", m_twoBodyDecays,
80  "Name of particle list of (mother) particles with daughters for calibration using vertex + mass constraint",
81  vector<string>());
82  addParam("primaryTwoBodyDecays", m_primaryTwoBodyDecays,
83  "Name of particle list of (mother) particles with daughters for calibration using vertex + IP profile + kinematics constraint",
84  vector<string>());
85  addParam("primaryMassTwoBodyDecays", m_primaryMassTwoBodyDecays,
86  "Name of particle list of (mother) particles with daughters for calibration using vertex + mass constraint",
87  vector<string>());
88  addParam("primaryMassVertexTwoBodyDecays", m_primaryMassVertexTwoBodyDecays,
89  "Name of particle list of (mother) particles with daughters for calibration using vertex + IP profile + mass constraint",
90  vector<string>());
91 
92  addParam("stableParticleWidth", m_stableParticleWidth,
93  "Width (in GeV/c/c) to use for invariant mass constraint for 'stable' particles (like K short). Temporary until proper solution is found.",
94  double(0.002));
95  // Configure output
96  addParam("doublePrecision", m_doublePrecision, "Use double (=true) or single/float (=false) precision for writing binary files",
97  bool(false));
98  addParam("useGblTree", m_useGblTree, "Store GBL trajectories in a tree instead of output to binary files",
99  bool(true));
100  addParam("absFilePaths", m_absFilePaths, "Use absolute paths to remember binary files. Only applies if useGblTree=False",
101  bool(false));
102 
103  // Configure global parameters
104  addParam("components", m_components,
105  "Specify which DB objects are calibrated, like ['BeamSpot', 'CDCTimeWalks'] or leave empty to use all components available.",
106  m_components);
107  addParam("calibrateVertex", m_calibrateVertex,
108  "For primary vertices / two body decays, beam spot vertex calibration derivatives are added",
109  bool(true));
110  addParam("calibrateKinematics", m_calibrateKinematics,
111  "For primary two body decays, beam spot kinematics calibration derivatives are added",
112  bool(true));
113 
114  //Configure GBL fit of individual tracks
115  addParam("externalIterations", m_externalIterations, "Number of external iterations of GBL fitter",
116  int(0));
117  addParam("internalIterations", m_internalIterations, "String defining internal GBL iterations for outlier down-weighting",
118  string(""));
119  addParam("recalcJacobians", m_recalcJacobians, "Up to which external iteration propagation Jacobians should be re-calculated",
120  int(0));
121 
122  addParam("minPValue", m_minPValue, "Minimum p-value to write out a (combined) trajectory. Set <0 to write out all.",
123  double(-1.));
124 
125  // Configure CDC specific options
126  addParam("fitTrackT0", m_fitTrackT0, "Add local parameter for track T0 fit in GBL",
127  bool(true));
128  addParam("updateCDCWeights", m_updateCDCWeights, "Update L/R weights from previous DAF fit result",
129  bool(true));
130  addParam("minCDCHitWeight", m_minCDCHitWeight, "Minimum (DAF) CDC hit weight for usage by GBL",
131  double(1.0E-6));
132  addParam("minUsedCDCHitFraction", m_minUsedCDCHitFraction, "Minimum used CDC hit fraction to write out a trajectory",
133  double(0.85));
134 
135  addParam("hierarchyType", m_hierarchyType, "Type of (VXD only now) hierarchy: 0 = None, 1 = Flat, 2 = Half-Shells, 3 = Full",
136  int(3));
137  addParam("enablePXDHierarchy", m_enablePXDHierarchy, "Enable PXD in hierarchy (flat or full)",
138  bool(true));
139  addParam("enableSVDHierarchy", m_enableSVDHierarchy, "Enable SVD in hierarchy (flat or full)",
140  bool(true));
141 
142  addParam("enableWireByWireAlignment", m_enableWireByWireAlignment, "Enable global derivatives for wire-by-wire alignment",
143  bool(false));
144  addParam("enableWireSagging", m_enableWireSagging, "Enable global derivatives for wire sagging",
145  bool(false));
146 
147  // Time dependence
148  addParam("events", m_eventNumbers,
149  "List of (event, run, exp) with event numbers at which payloads can change for timedep calibration.",
150  m_eventNumbers);
151  // Time dependence config
152  addParam("timedepConfig", m_timedepConfig,
153  "list{ {list{param1, param2, ...}, list{(ev1, run1, exp1), ...}}, ... }.",
154  m_timedepConfig);
155 
156  // Custom mass+width config
157  addParam("customMassConfig", m_customMassConfig,
158  "dict{ list_name: (mass, width), ... } with custom mass and width to use as external measurement.",
159  m_customMassConfig);
160 }
161 
162 void MillepedeCollectorModule::prepare()
163 {
165  emd.isRequired();
166 
167  StoreObjPtr<EventT0> eventT0;
168  //eventT0.isRequired();
169 
170  if (m_tracks.empty() &&
171  m_particles.empty() &&
172  m_vertices.empty() &&
173  m_primaryVertices.empty() &&
174  m_twoBodyDecays.empty() &&
175  m_primaryTwoBodyDecays.empty() &&
176  m_primaryMassTwoBodyDecays.empty() &&
177  m_primaryMassVertexTwoBodyDecays.empty())
178  B2ERROR("You have to specify either arrays of single tracks or particle lists of single single particles or mothers with vertex constrained daughters.");
179 
180  if (!m_tracks.empty()) {
181  for (auto arrayName : m_tracks)
182  continue;
183  // StoreArray<RecoTrack>::required(arrayName);
184  }
185 
186  if (!m_particles.empty() || !m_vertices.empty() || !m_primaryVertices.empty()) {
187  StoreArray<RecoTrack> recoTracks;
188  StoreArray<Track> tracks;
189  StoreArray<TrackFitResult> trackFitResults;
190 
191  //recoTracks.isRequired();
192  //tracks.isRequired();
193  //trackFitResults.isRequired();
194  }
195 
196  for (auto listName : m_particles) {
197  StoreObjPtr<ParticleList> list(listName);
198  //list.isRequired();
199  }
200 
201  for (auto listName : m_vertices) {
202  StoreObjPtr<ParticleList> list(listName);
203  //list.isRequired();
204  }
205 
206  for (auto listName : m_primaryVertices) {
207  StoreObjPtr<ParticleList> list(listName);
208  //list.isRequired();
209  }
210 
211  // Register Mille output
212  registerObject<MilleData>("mille", new MilleData(m_doublePrecision, m_absFilePaths));
213 
214  auto gblDataTree = new TTree("GblDataTree", "GblDataTree");
215  gblDataTree->Branch<std::vector<gbl::GblData>>("GblData", &m_currentGblData, 32000, 99);
216  registerObject<TTree>("GblDataTree", gblDataTree);
217 
218  registerObject<TH1I>("ndf", new TH1I("ndf", "ndf", 200, 0, 200));
219  registerObject<TH1F>("chi2_per_ndf", new TH1F("chi2_per_ndf", "chi2 divided by ndf", 200, 0., 50.));
220  registerObject<TH1F>("pval", new TH1F("pval", "pval", 100, 0., 1.));
221 
222  registerObject<TH1F>("cdc_hit_fraction", new TH1F("cdc_hit_fraction", "cdc_hit_fraction", 100, 0., 1.));
223  registerObject<TH1F>("evt0", new TH1F("evt0", "evt0", 400, -100., 100.));
224 
225  // Configure the (VXD) hierarchy before being built
226  if (m_hierarchyType == 0)
227  Belle2::alignment::VXDGlobalParamInterface::s_hierarchyType = VXDGlobalParamInterface::c_None;
228  else if (m_hierarchyType == 1)
229  Belle2::alignment::VXDGlobalParamInterface::s_hierarchyType = VXDGlobalParamInterface::c_Flat;
230  else if (m_hierarchyType == 2)
231  Belle2::alignment::VXDGlobalParamInterface::s_hierarchyType = VXDGlobalParamInterface::c_HalfShells;
232  else if (m_hierarchyType == 3)
233  Belle2::alignment::VXDGlobalParamInterface::s_hierarchyType = VXDGlobalParamInterface::c_Full;
234 
237 
238  std::vector<EventMetaData> events;
239  for (auto& ev_run_exp : m_eventNumbers) {
240  events.push_back(EventMetaData(std::get<0>(ev_run_exp), std::get<1>(ev_run_exp), std::get<2>(ev_run_exp)));
241  }
242 
243  // This will also build the hierarchy for the first time:
244  if (!m_timedepConfig.empty() && m_eventNumbers.empty()) {
245  auto autoEvents = Belle2::alignment::timeline::setupTimedepGlobalLabels(m_timedepConfig);
247  } else if (m_timedepConfig.empty() && !m_eventNumbers.empty()) {
249  } else if (m_timedepConfig.empty() && m_eventNumbers.empty()) {
251  } else {
252  B2ERROR("Cannot set both, event list and timedep config.");
253  }
254 
255 // Belle2::alignment::GlobalCalibrationManager::getInstance().writeConstraints("constraints.txt");
256 
257  AlignableCDCRecoHit::s_enableTrackT0LocalDerivative = m_fitTrackT0;
258  AlignableCDCRecoHit::s_enableWireSaggingGlobalDerivative = m_enableWireSagging;
259  AlignableCDCRecoHit::s_enableWireByWireAlignmentGlobalDerivatives = m_enableWireByWireAlignment;
260 }
261 
262 void MillepedeCollectorModule::collect()
263 {
265  alignment::GlobalCalibrationManager::getInstance().preCollect(*emd);
266  StoreObjPtr<EventT0> eventT0;
267 
268  if (!m_useGblTree) {
269  // Open new file on request (at start or after being closed)
270  auto mille = getObjectPtr<MilleData>("mille");
271  if (!mille->isOpen())
272  mille->open(getUniqueMilleName());
273  }
274 
275  std::shared_ptr<genfit::GblFitter> gbl(new genfit::GblFitter());
276  double chi2 = -1.;
277  double lostWeight = -1.;
278  int ndf = -1;
279  float evt0 = -9999.;
280 
281  for (auto arrayName : m_tracks) {
282  StoreArray<RecoTrack> recoTracks(arrayName);
283  if (!recoTracks.isValid())
284  continue;
285 
286  for (auto& recoTrack : recoTracks) {
287 
288  if (!fitRecoTrack(recoTrack))
289  continue;
290 
291  auto& track = RecoTrackGenfitAccess::getGenfitTrack(recoTrack);
292  if (!track.hasFitStatus())
293  continue;
294  genfit::GblFitStatus* fs = dynamic_cast<genfit::GblFitStatus*>(track.getFitStatus());
295  if (!fs)
296  continue;
297 
298  if (!fs->isFittedWithReferenceTrack())
299  continue;
300 
301  using namespace gbl;
302  GblTrajectory trajectory(gbl->collectGblPoints(&track, track.getCardinalRep()), fs->hasCurvature());
303 
304  trajectory.fit(chi2, ndf, lostWeight);
305  getObjectPtr<TH1I>("ndf")->Fill(ndf);
306  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
307  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
308  if (eventT0.isValid() && eventT0->hasEventT0()) {
309  evt0 = eventT0->getEventT0();
310  getObjectPtr<TH1F>("evt0")->Fill(evt0);
311  }
312 
313  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(trajectory);
314 
315  }
316 
317  }
318 
319  for (auto listName : m_particles) {
320  StoreObjPtr<ParticleList> list(listName);
321  if (!list.isValid())
322  continue;
323 
324  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
325  for (auto& track : getParticlesTracks({list->getParticle(iParticle)}, false)) {
326  auto gblfs = dynamic_cast<genfit::GblFitStatus*>(track->getFitStatus());
327 
328  gbl::GblTrajectory trajectory(gbl->collectGblPoints(track, track->getCardinalRep()), gblfs->hasCurvature());
329 
330  trajectory.fit(chi2, ndf, lostWeight);
331  getObjectPtr<TH1I>("ndf")->Fill(ndf);
332  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
333  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
334  if (eventT0.isValid() && eventT0->hasEventT0()) {
335  evt0 = eventT0->getEventT0();
336  getObjectPtr<TH1F>("evt0")->Fill(evt0);
337  }
338 
339  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(trajectory);
340 
341  }
342  }
343  }
344 
345  for (auto listName : m_vertices) {
346  StoreObjPtr<ParticleList> list(listName);
347  if (!list.isValid())
348  continue;
349 
350  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
351  auto mother = list->getParticle(iParticle);
352  std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
353 
354  for (auto& track : getParticlesTracks(mother->getDaughters()))
355  daughters.push_back({
356  gbl->collectGblPoints(track, track->getCardinalRep()),
357  getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2)
358  });
359 
360  if (daughters.size() > 1) {
361  gbl::GblTrajectory combined(daughters);
362 
363  combined.fit(chi2, ndf, lostWeight);
364  getObjectPtr<TH1I>("ndf")->Fill(ndf);
365  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
366  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
367  if (eventT0.isValid() && eventT0->hasEventT0()) {
368  evt0 = eventT0->getEventT0();
369  getObjectPtr<TH1F>("evt0")->Fill(evt0);
370  }
371 
372 
373  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
374 
375  B2RESULT("Vertex-constrained fit NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
376 
377  }
378  }
379  }
380 
381  for (auto listName : m_primaryVertices) {
382  StoreObjPtr<ParticleList> list(listName);
383  if (!list.isValid())
384  continue;
385 
386  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
387  auto mother = list->getParticle(iParticle);
388  std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
389 
390  TMatrixD extProjection(5, 3);
391  TMatrixD locProjection(3, 5);
392 
393  bool first(true);
394  for (auto& track : getParticlesTracks(mother->getDaughters())) {
395  if (first) {
396  // For first trajectory only
397  extProjection = getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2);
398  locProjection = getLocalToGlobalTransform(track->getFittedState()).GetSub(0, 2, 0, 4);
399  first = false;
400  }
401  daughters.push_back({
402  gbl->collectGblPoints(track, track->getCardinalRep()),
403  getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2)
404  });
405  }
406 
407  if (daughters.size() > 1) {
408  auto beam = getPrimaryVertexAndCov();
409 
410  TMatrixDSym vertexCov(get<TMatrixDSym>(beam));
411  TMatrixDSym vertexPrec(get<TMatrixDSym>(beam).Invert());
412  B2Vector3D vertexResidual = - (mother->getVertex() - get<B2Vector3D>(beam));
413 
414  TVectorD extMeasurements(3);
415  extMeasurements[0] = vertexResidual[0];
416  extMeasurements[1] = vertexResidual[1];
417  extMeasurements[2] = vertexResidual[2];
418 
419  TMatrixD extDeriv(3, 3);
420  extDeriv.Zero();
421  // beam vertex constraint
422  extDeriv(0, 0) = 1.;
423  extDeriv(1, 1) = 1.;
424  extDeriv(2, 2) = 1.;
425 
426  if (m_calibrateVertex) {
427  TMatrixD derivatives(3, 3);
428  derivatives.Zero();
429  derivatives(0, 0) = 1.;
430  derivatives(1, 1) = 1.;
431  derivatives(2, 2) = 1.;
432 
433  std::vector<int> labels;
434  GlobalLabel label = GlobalLabel::construct<BeamSpot>(0, 0);
435  labels.push_back(label.setParameterId(1));
436  labels.push_back(label.setParameterId(2));
437  labels.push_back(label.setParameterId(3));
438 
439  // Allow to disable BeamSpot externally
440  alignment::GlobalDerivatives globals(labels, derivatives);
441  // Add derivatives for vertex calibration to first point of first trajectory
442  // NOTE: use GlobalDerivatives operators vector<int> and TMatrixD which filter
443  // the derivatives to not pass those with zero labels (usefull to get rid of some params)
444  std::vector<int> lab(globals); TMatrixD der(globals);
445 
446  // Transformation from local system at (vertex) point to global (vx,vy,vz)
447  // of the (decay) vertex
448  //
449  // d(q/p,u',v',u,v)/d(vy,vy,vz) = dLocal_dExt
450  //
451  //
452  // Note its transpose is its "inverse" in the sense that
453  //
454  // dloc/dext * (dloc/dext)^T = diag(0, 0, 0, 0, 1, 1)
455  //
456  //
457  // N.B. typical dLocal_dExt matrix (5x3):
458  //
459  // | 0 | 1 | 2 |
460  // --------------------------------------------
461  // 0 | 0 0 0
462  // 1 | 0 0 0
463  // 2 | 0 0 0
464  // 3 | -0.02614 -0.9997 0
465  // 4 | 0 0 1
466  //
467  // Therefore one can simplify things by only taking the last two rows/columns in vectors/matrices
468  // and vertex measurement can be expressed as standard 2D measurement in GBL.
469  //
470  TMatrixD dLocal_dExt = extProjection;
471  TMatrixD dExt_dLocal = locProjection;
472 
473  TVectorD locRes = dLocal_dExt * extMeasurements;
474  // Do not use inverted covariance - seems to have issues with numeric precision
475  TMatrixD locCov = dLocal_dExt * vertexCov * dExt_dLocal;
476  // Invert here only the 2D sub-matrix (rest is zero due to the foÅ•m of dLocal_dExt)
477  TMatrixD locPrec = locCov.GetSub(3, 4, 3, 4).Invert();
478  TMatrixDSym locPrec2D(2); locPrec2D.Zero();
479  for (int i = 0; i < 2; ++i)
480  for (int j = 0; j < 2; ++j)
481  locPrec2D(i, j) = locPrec(i, j);
482 
483  // Take the 2 last components also for residuals and global derivatives
484  // (in local system of vertex point - defined during fitRecoTrack(..., particle) and using
485  // the (hopefully) updated momentum and position seed after vertex fit by modularAnalysis
486  TVectorD locRes2D = locRes.GetSub(3, 4);
487  TMatrixD locDerivs2D = (extProjection * der).GetSub(3, 4, 0, 2);
488 
489  // Attach the primary beamspot vertex position as a measurement at 1st point
490  // of first trajectory (and optionaly also the global derivatives for beamspot alignment
491  daughters[0].first[0].addMeasurement(locRes2D, locPrec2D);
492  if (!lab.empty()) {
493  daughters[0].first[0].addGlobals(lab, locDerivs2D);
494  }
495 
496  gbl::GblTrajectory combined(daughters);
497  //combined.printTrajectory(100);
498  //combined.printPoints(100);
499 
500  combined.fit(chi2, ndf, lostWeight);
501  getObjectPtr<TH1I>("ndf")->Fill(ndf);
502  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
503  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
504  if (eventT0.isValid() && eventT0->hasEventT0()) {
505  evt0 = eventT0->getEventT0();
506  getObjectPtr<TH1F>("evt0")->Fill(evt0);
507  }
508 
509  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
510  B2RESULT("Beam vertex constrained fit results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
511 
512  } else {
513 
514  gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, vertexPrec);
515 
516  combined.fit(chi2, ndf, lostWeight);
517  getObjectPtr<TH1I>("ndf")->Fill(ndf);
518  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
519  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
520  if (eventT0.isValid() && eventT0->hasEventT0()) {
521  evt0 = eventT0->getEventT0();
522  getObjectPtr<TH1F>("evt0")->Fill(evt0);
523  }
524 
525  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
526 
527  B2RESULT("Beam vertex constrained fit results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
528 
529  }
530  }
531  }
532  }
533 
534  for (auto listName : m_twoBodyDecays) {
535  StoreObjPtr<ParticleList> list(listName);
536  if (!list.isValid())
537  continue;
538 
539  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
540 
541  auto mother = list->getParticle(iParticle);
542  auto track12 = getParticlesTracks(mother->getDaughters());
543  if (track12.size() != 2) {
544  B2ERROR("Did not get 2 fitted tracks. Skipping this mother.");
545  continue;
546  }
547 
548  auto pdgdb = EvtGenDatabasePDG::Instance();
549  double motherMass = mother->getPDGMass();
550  double motherWidth = pdgdb->GetParticle(mother->getPDGCode())->Width();
551 
552  updateMassWidthIfSet(listName, motherMass, motherWidth);
553 
554  //TODO: what to take as width for "real" particles? -> make a param for default detector mass resolution??
555  if (motherWidth == 0.) {
556  motherWidth = m_stableParticleWidth * Unit::GeV;
557  B2WARNING("Using artificial width for " << pdgdb->GetParticle(mother->getPDGCode())->GetName() << " : " << motherWidth << " GeV");
558  }
559 
560  auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
561  std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
562 
563  daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
564  daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
565 
566  TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
567  TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
568 
569  TVectorD extMeasurements(1);
570  extMeasurements[0] = massResidual[0];
571 
572  TMatrixD extDeriv(1, 9);
573  extDeriv.Zero();
574  extDeriv(0, 8) = 1.;
575 
576  gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, massPrec);
577 
578  combined.fit(chi2, ndf, lostWeight);
579  //combined.printTrajectory(1000);
580  //combined.printPoints(1000);
581  getObjectPtr<TH1I>("ndf")->Fill(ndf);
582  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
583  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
584  if (eventT0.isValid() && eventT0->hasEventT0()) {
585  evt0 = eventT0->getEventT0();
586  getObjectPtr<TH1F>("evt0")->Fill(evt0);
587  }
588 
589 
590  B2RESULT("Mass(PDG) + vertex constrained fit results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
591 
592  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
593 
594  }
595  }
596 
597  for (auto listName : m_primaryMassTwoBodyDecays) {
598  StoreObjPtr<ParticleList> list(listName);
599  if (!list.isValid())
600  continue;
601 
603 
604  double motherMass = beam->getMass();
605  double motherWidth = sqrt((beam->getCovHER() + beam->getCovLER())(0, 0));
606 
607  updateMassWidthIfSet(listName, motherMass, motherWidth);
608 
609  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
610 
611  auto mother = list->getParticle(iParticle);
612  auto track12 = getParticlesTracks(mother->getDaughters());
613  if (track12.size() != 2) {
614  B2ERROR("Did not get 2 fitted tracks. Skipping this mother.");
615  continue;
616  }
617 
618  auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
619  std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
620 
621  daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
622  daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
623 
624  TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
625  TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
626 
627  TVectorD extMeasurements(1);
628  extMeasurements[0] = massResidual[0];
629 
630  TMatrixD extDeriv(1, 9);
631  extDeriv.Zero();
632  extDeriv(0, 8) = 1.;
633 
634  gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, massPrec);
635 
636  combined.fit(chi2, ndf, lostWeight);
637  getObjectPtr<TH1I>("ndf")->Fill(ndf);
638  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
639  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
640  if (eventT0.isValid() && eventT0->hasEventT0()) {
641  evt0 = eventT0->getEventT0();
642  getObjectPtr<TH1F>("evt0")->Fill(evt0);
643  }
644 
645 
646  B2RESULT("Mass constrained fit results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
647 
648  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
649 
650  }
651  }
652 
653  for (auto listName : m_primaryMassVertexTwoBodyDecays) {
654  StoreObjPtr<ParticleList> list(listName);
655  if (!list.isValid())
656  continue;
657 
659 
660  double motherMass = beam->getMass();
661  double motherWidth = sqrt((beam->getCovHER() + beam->getCovLER())(0, 0));
662 
663  updateMassWidthIfSet(listName, motherMass, motherWidth);
664 
665  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
666 
667  auto mother = list->getParticle(iParticle);
668  auto track12 = getParticlesTracks(mother->getDaughters());
669  if (track12.size() != 2) {
670  B2ERROR("Did not get 2 fitted tracks. Skipping this mother.");
671  continue;
672  }
673 
674  auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
675  std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
676 
677  daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
678  daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
679 
680  TMatrixDSym vertexPrec(get<TMatrixDSym>(getPrimaryVertexAndCov()).Invert());
681  B2Vector3D vertexResidual = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()));
682 
683  TMatrixDSym massPrec(1); massPrec(0, 0) = 1. / motherWidth / motherWidth;
684  TVectorD massResidual(1); massResidual = - (mother->getMass() - motherMass);
685 
686  TMatrixDSym extPrec(4); extPrec.Zero();
687  extPrec.SetSub(0, 0, vertexPrec);
688  extPrec(3, 3) = massPrec(0, 0);
689 
690  TVectorD extMeasurements(4);
691  extMeasurements[0] = vertexResidual[0];
692  extMeasurements[1] = vertexResidual[1];
693  extMeasurements[2] = vertexResidual[2];
694  extMeasurements[3] = massResidual[0];
695 
696  TMatrixD extDeriv(4, 9);
697  extDeriv.Zero();
698  extDeriv(0, 0) = 1.;
699  extDeriv(1, 1) = 1.;
700  extDeriv(2, 2) = 1.;
701  extDeriv(3, 8) = 1.;
702 
703  gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, extPrec);
704 
705  combined.fit(chi2, ndf, lostWeight);
706  getObjectPtr<TH1I>("ndf")->Fill(ndf);
707  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
708  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
709  if (eventT0.isValid() && eventT0->hasEventT0()) {
710  evt0 = eventT0->getEventT0();
711  getObjectPtr<TH1F>("evt0")->Fill(evt0);
712  }
713 
714 
715  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
716 
717  B2RESULT("Mass + vertex constrained fit results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
718 
719  }
720  }
721 
722  for (auto listName : m_primaryTwoBodyDecays) {
723  B2WARNING("This should NOT be used for production of calibration constants for the real detector (yet)!");
724 
725  StoreObjPtr<ParticleList> list(listName);
726  if (!list.isValid())
727  continue;
728 
730 
731  // For the error of invariant mass M = 2 * sqrt(E_HER * E_LER) (for m_e ~ 0)
732  double M = beam->getMass();
733  double E_HER = beam->getHER().E();
734  double E_LER = beam->getLER().E();
735 
736  double pz = (beam->getHER().Vect() + beam->getLER().Vect())[2];
737  double E = (beam->getHER() + beam->getLER()).E();
738 
739  double motherMass = beam->getMass();
740  double motherWidth = sqrt((E_HER / M) * (E_HER / M) * beam->getCovLER()(0, 0) + (E_LER / M) * (E_LER / M) * beam->getCovHER()(0,
741  0));
742 
743  updateMassWidthIfSet(listName, motherMass, motherWidth);
744 
745  for (unsigned int iParticle = 0; iParticle < list->getListSize(); ++iParticle) {
746 
747  B2WARNING("Two body decays with full kinematic constraint not yet correct - need to resolve strange covariance provided by BeamParameters!");
748 
749  auto mother = list->getParticle(iParticle);
750 
751  auto track12 = getParticlesTracks(mother->getDaughters());
752  if (track12.size() != 2) {
753  B2ERROR("Did not get exactly 2 fitted tracks. Skipping this mother in list " << listName);
754  continue;
755  }
756 
757  auto dfdextPlusMinus = getTwoBodyToLocalTransform(*mother, motherMass);
758  std::vector<std::pair<std::vector<gbl::GblPoint>, TMatrixD> > daughters;
759 
760  daughters.push_back({gbl->collectGblPoints(track12[0], track12[0]->getCardinalRep()), dfdextPlusMinus.first});
761  daughters.push_back({gbl->collectGblPoints(track12[1], track12[1]->getCardinalRep()), dfdextPlusMinus.second});
762 
763  TMatrixDSym extCov(7); extCov.Zero();
764 
765  // 3x3 IP vertex covariance
766  extCov.SetSub(0, 0, get<TMatrixDSym>(getPrimaryVertexAndCov()));
767 
768  // 3x3 boost vector covariance
769  //NOTE: BeamSpot return covarince in variables (E, theta_x, theta_y)
770  // We need to transform it to our variables (px, py, pz)
771 
772  TMatrixD dBoost_dVect(3, 3);
773  dBoost_dVect(0, 0) = 0.; dBoost_dVect(0, 1) = 1. / pz; dBoost_dVect(0, 2) = 0.;
774  dBoost_dVect(1, 0) = 0.; dBoost_dVect(1, 1) = 0.; dBoost_dVect(1, 2) = 1. / pz;
775  dBoost_dVect(2, 0) = pz / E; dBoost_dVect(2, 1) = 0.; dBoost_dVect(2, 2) = 0.;
776 
777  TMatrixD dVect_dBoost(3, 3);
778  dVect_dBoost(0, 0) = 0.; dVect_dBoost(0, 1) = 0.; dVect_dBoost(0, 2) = E / pz;
779  dVect_dBoost(1, 0) = pz; dVect_dBoost(1, 1) = 0.; dVect_dBoost(1, 2) = 0.;
780  dVect_dBoost(2, 0) = 0.; dVect_dBoost(2, 1) = pz; dVect_dBoost(2, 2) = 0.;
781 
782  TMatrixD covBoost(3, 3);
783  for (int i = 0; i < 3; ++i) {
784  for (int j = i; j < 3; ++j) {
785  covBoost(j, i) = covBoost(i, j) = (beam->getCovHER() + beam->getCovLER())(i, j);
786  }
787  }
788  //TODO: Temporary fix: if theta_x, theta_y covariance is zero, use arbitrary 10mrad^2
789 // if (covBoost(1, 1) == 0.) covBoost(1, 1) = 1.;
790 // if (covBoost(2, 2) == 0.) covBoost(2, 2) = 1.;
791  if (covBoost(1, 1) == 0.) covBoost(1, 1) = 1.e-4;
792  if (covBoost(2, 2) == 0.) covBoost(2, 2) = 1.e-4;
793 
794  TMatrixD covVect = dBoost_dVect * covBoost * dVect_dBoost;
795 
796  extCov.SetSub(3, 3, covVect);
797 
798  extCov(6, 6) = motherWidth * motherWidth;
799  auto extPrec = extCov; extPrec.Invert();
800 
801  TVectorD extMeasurements(7);
802  extMeasurements[0] = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()))[0];
803  extMeasurements[1] = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()))[1];
804  extMeasurements[2] = - (mother->getVertex() - get<B2Vector3D>(getPrimaryVertexAndCov()))[2];
805  extMeasurements[3] = - (mother->getMomentum() - (beam->getHER().Vect() + beam->getLER().Vect()))[0];
806  extMeasurements[4] = - (mother->getMomentum() - (beam->getHER().Vect() + beam->getLER().Vect()))[1];
807  extMeasurements[5] = - (mother->getMomentum() - (beam->getHER().Vect() + beam->getLER().Vect()))[2];
808  extMeasurements[6] = - (mother->getMass() - motherMass);
809 
810  B2INFO("mother mass = " << mother->getMass() << " and beam mass = " << beam->getMass());
811 
812  TMatrixD extDeriv(7, 9);
813  extDeriv.Zero();
814  // beam vertex constraint
815  extDeriv(0, 0) = 1.;
816  extDeriv(1, 1) = 1.;
817  extDeriv(2, 2) = 1.;
818  // beam kinematics constraint
819  extDeriv(3, 3) = 1.;
820  extDeriv(4, 4) = 1.;
821  extDeriv(5, 5) = 1.;
822  // beam inv. mass constraint
823  extDeriv(6, 8) = 1;
824 
825  if (m_calibrateVertex || m_calibrateKinematics) {
826  B2WARNING("Primary vertex+kinematics calibration not (yet?) fully implemented!");
827  B2WARNING("This code is highly experimental and has (un)known issues!");
828 
829  // up to d(x,y,z,px,py,pz,theta,phi,M)/d(vx,vy,vz,theta_x,theta_y,E)
830  TMatrixD derivatives(9, 6);
831  std::vector<int> labels;
832  derivatives.Zero();
833 
834  if (m_calibrateVertex) {
835  derivatives(0, 0) = 1.;
836  derivatives(1, 1) = 1.;
837  derivatives(2, 2) = 1.;
838  GlobalLabel label = GlobalLabel::construct<BeamSpot>(0, 0);
839  labels.push_back(label.setParameterId(1));
840  labels.push_back(label.setParameterId(2));
841  labels.push_back(label.setParameterId(3));
842  } else {
843  labels.push_back(0);
844  labels.push_back(0);
845  labels.push_back(0);
846  }
847 
848  if (m_calibrateKinematics) {
849  derivatives(3, 3) = mother->getMomentumMagnitude();
850  derivatives(4, 4) = mother->getMomentumMagnitude();
851  derivatives(8, 5) = (beam->getLER().E() + beam->getHER().E()) / beam->getMass();
852 
853  GlobalLabel label = GlobalLabel::construct<BeamSpot>(0, 0);
854  labels.push_back(label.setParameterId(4)); //theta_x
855  labels.push_back(label.setParameterId(5)); //theta_y
856  labels.push_back(label.setParameterId(6)); //E
857 
858  } else {
859  labels.push_back(0);
860  labels.push_back(0);
861  labels.push_back(0);
862  }
863 
864  // Allow to disable BeamSpot externally
865  alignment::GlobalDerivatives globals(labels, derivatives);
866 
867  // Add derivatives for vertex calibration to first point of first trajectory
868  // NOTE: use GlobalDerivatives operators vector<int> and TMatrixD which filter
869  // the derivatives to not pass those with zero labels (usefull to get rid of some params)
870  std::vector<int> lab(globals); TMatrixD der(globals);
871 
872  // I want: dlocal/dext = dlocal/dtwobody * dtwobody/dext = dfdextPlusMinus * dtwobody/dext
873  TMatrixD dTwoBody_dExt(9, 7);
874  dTwoBody_dExt.Zero();
875  // beam vertex constraint
876  dTwoBody_dExt(0, 0) = 1.;
877  dTwoBody_dExt(1, 1) = 1.;
878  dTwoBody_dExt(2, 2) = 1.;
879  // beam kinematics constraint
880  dTwoBody_dExt(3, 3) = 1.;
881  dTwoBody_dExt(4, 4) = 1.;
882  dTwoBody_dExt(5, 5) = 1.;
883  // beam inv. mass constraint
884  dTwoBody_dExt(8, 6) = 1.;
885 
886  const TMatrixD dLocal_dExt = dfdextPlusMinus.first * dTwoBody_dExt;
887  TMatrixD dLocal_dExt_T = dLocal_dExt; dLocal_dExt_T.T();
888 
889  // The 5x7 transformation matrix d(q/p,u',v',u,v)/d(vx,vy,vz,px,py,pz,M) needs to be "inverted"
890  // to transform the covariance of the beamspot and boost vector of SuperKEKB into the local system
891  // of one GBL point - such that Millepede can align the beamspot (or even beam kinematics) if requested.
892  //
893  // I tested also other methods, but only the Singular Value Decomposition gives nice-enough results,
894  // with almost no code:
895  //
896  TDecompSVD svd(dLocal_dExt_T);
897  TMatrixD dExt_dLocal = svd.Invert().T();
898  //
899  // (dLocal_dExt * dExt_dLocal).Print(); // Check how close we are to unit matrix
900  //
901  // 5x5 matrix is as follows
902  //
903  // | 0 | 1 | 2 | 3 | 4 |
904  // ----------------------------------------------------------------------
905  // 0 | 1 -2.58e-17 6.939e-18 1.571e-17 -1.649e-19
906  // 1 | 1.787e-14 1 5.135e-16 -3.689e-16 -2.316e-18
907  // 2 | -1.776e-15 -7.806e-17 1 5.636e-17 6.193e-18
908  // 3 | -2.453e-15 7.26e-18 2.009e-16 1 -1.14e-16
909  // 4 | -1.689e-14 -9.593e-17 -2.317e-15 -3.396e-17 1
910  //
911  // It took me half a day to find out how to do this with 2 lines of code (3 with the include).
912  // Source: ROOT macro example - actually found at:
913  // <https://root.cern.ch/root/html/tutorials/matrix/solveLinear.C.html>
914  for (int i = 0; i < 7; ++i) {
915  for (int j = 0; j < 5; ++j) {
916  if (fabs(dExt_dLocal(i, j)) < 1.e-6)
917  dExt_dLocal(i, j) = 0.;
918  }
919  }
920  const TVectorD locRes = dLocal_dExt * extMeasurements;
921  const TMatrixD locPrec = dLocal_dExt * extPrec * dExt_dLocal;
922 
923  TMatrixDSym locPrecSym(5); locPrecSym.Zero();
924  for (int i = 0; i < 5; ++i) {
925  for (int j = i; j < 5; ++j) {
926  //locPrecSym(j, i) = locPrecSym(i, j) = locPrec(i, j);
927  locPrecSym(j, i) = locPrecSym(i, j) = (fabs(locPrec(i, j)) > 1.e-6) ? locPrec(i, j) : 0.;
928  }
929  }
930 
931  daughters[0].first[0].addMeasurement(locRes, locPrecSym);
932  if (!lab.empty())
933  daughters[0].first[0].addGlobals(lab, dfdextPlusMinus.first * der);
934 
935  //TODO: Understand this: either find a bug somewhere or improve the parametrization or .... ?
936  // This should be enough, but the parametrization seems to fail for nearly horizontal pairs...
937  //gbl::GblTrajectory combined(daughters);
938  // This should not be needed, it actually seems to make worse Chi2/NDF, but GBL does not fail.
939  // The measurement added just to be able to add the global derivatives (done just above) is redundant
940  // to the external measurement added here:
941  gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, extPrec);
942  //combined.printTrajectory(1000);
943  //combined.printPoints(1000);
944 
945  combined.fit(chi2, ndf, lostWeight);
946  getObjectPtr<TH1I>("ndf")->Fill(ndf);
947  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
948  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
949  if (eventT0.isValid() && eventT0->hasEventT0()) {
950  evt0 = eventT0->getEventT0();
951  getObjectPtr<TH1F>("evt0")->Fill(evt0);
952  }
953 
954 
955  B2RESULT("Full kinematic-constrained fit (calibration version) results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
956 
957  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
958 
959  } else {
960 
961  gbl::GblTrajectory combined(daughters, extDeriv, extMeasurements, extPrec);
962  //combined.printTrajectory(1000);
963  //combined.printPoints(1000);
964 
965  combined.fit(chi2, ndf, lostWeight);
966  getObjectPtr<TH1I>("ndf")->Fill(ndf);
967  getObjectPtr<TH1F>("chi2_per_ndf")->Fill(chi2 / double(ndf));
968  getObjectPtr<TH1F>("pval")->Fill(TMath::Prob(chi2, ndf));
969  if (eventT0.isValid() && eventT0->hasEventT0()) {
970  evt0 = eventT0->getEventT0();
971  getObjectPtr<TH1F>("evt0")->Fill(evt0);
972  }
973 
974 
975  B2RESULT("Full kinematic-constrained fit results NDF = " << ndf << " Chi2/NDF = " << chi2 / double(ndf));
976 
977  if (TMath::Prob(chi2, ndf) > m_minPValue) storeTrajectory(combined);
978  }
979  }
980  }
981 }
982 
983 void MillepedeCollectorModule::closeRun()
984 {
985  // We close the file at end of run, producing
986  // one file per run (and process id) which is more
987  // convenient than one large binary block.
988  auto mille = getObjectPtr<MilleData>("mille");
989  if (mille->isOpen())
990  mille->close();
991 }
992 
993 void MillepedeCollectorModule::finish()
994 {
996 
997  StoreObjPtr<FileMetaData> fileMetaData("", DataStore::c_Persistent);
998  if (!fileMetaData.isValid()) {
999  B2ERROR("Cannot register binaries in FileCatalog.");
1000  return;
1001  }
1002 
1003 
1004  const std::vector<string> parents = {fileMetaData->getLfn()};
1005  for (auto binary : getObjectPtr<MilleData>("mille")->getFiles()) {
1006  FileMetaData milleMetaData(*fileMetaData);
1007  // We reset filename to be set directly by the registerFile procedure
1008  milleMetaData.setLfn("");
1009  milleMetaData.setParents(parents);
1010  FileCatalog::Instance().registerFile(binary, milleMetaData);
1011  }
1012 
1013 }
1014 
1015 void MillepedeCollectorModule::storeTrajectory(gbl::GblTrajectory& trajectory)
1016 {
1017  if (m_useGblTree) {
1018  if (trajectory.isValid())
1019  m_currentGblData = trajectory.getData();
1020  else
1021  m_currentGblData.clear();
1022 
1023  if (!m_currentGblData.empty())
1024  getObjectPtr<TTree>("GblDataTree")->Fill();
1025  } else {
1026  getObjectPtr<MilleData>("mille")->fill(trajectory);
1027  }
1028 }
1029 
1030 std::string MillepedeCollectorModule::getUniqueMilleName()
1031 {
1033  string name = getName();
1034 
1035  name += "-e" + to_string(emd->getExperiment());
1036  name += "-r" + to_string(emd->getRun());
1037  name += "-ev" + to_string(emd->getEvent());
1038 
1039  if (ProcHandler::parallelProcessingUsed())
1040  name += "-pid" + to_string(ProcHandler::EvtProcID());
1041 
1042  name += ".mille";
1043 
1044  return name;
1045 }
1046 
1047 bool MillepedeCollectorModule::fitRecoTrack(RecoTrack& recoTrack, Particle* particle)
1048 {
1049  try {
1050  // For already fitted tracks, try to get fitted (DAF) weights for CDC
1051  if (m_updateCDCWeights && recoTrack.getNumberOfCDCHits() && recoTrack.getTrackFitStatus()
1052  && recoTrack.getTrackFitStatus()->isFitted()) {
1053  double sumCDCWeights = recoTrack.getNumberOfCDCHits(); // start with full weights
1054  // Do the hits synchronisation
1055  auto relatedRecoHitInformation =
1057 
1058  for (RecoHitInformation& recoHitInformation : relatedRecoHitInformation) {
1059 
1060  if (recoHitInformation.getFlag() == RecoHitInformation::c_pruned) {
1061  B2FATAL("Found pruned point in RecoTrack. Pruned tracks cannot be used in MillepedeCollector.");
1062  }
1063 
1064  if (recoHitInformation.getTrackingDetector() != RecoHitInformation::c_CDC) continue;
1065 
1066  const genfit::TrackPoint* trackPoint = recoTrack.getCreatedTrackPoint(&recoHitInformation);
1067  if (trackPoint) {
1068  if (not trackPoint->hasFitterInfo(recoTrack.getCardinalRepresentation()))
1069  continue;
1070  auto kalmanFitterInfo = dynamic_cast<genfit::KalmanFitterInfo*>(trackPoint->getFitterInfo());
1071  if (not kalmanFitterInfo) {
1072  continue;
1073  } else {
1074  std::vector<double> weights = kalmanFitterInfo->getWeights();
1075  if (weights.size() == 2) {
1076  if (weights.at(0) > weights.at(1))
1077  recoHitInformation.setRightLeftInformation(RecoHitInformation::c_left);
1078  else if (weights.at(0) < weights.at(1))
1079  recoHitInformation.setRightLeftInformation(RecoHitInformation::c_right);
1080 
1081  double weightLR = weights.at(0) + weights.at(1);
1082  if (weightLR < m_minCDCHitWeight) recoHitInformation.setUseInFit(false);
1083  sumCDCWeights += weightLR - 1.; // reduce weight sum if weightLR<1
1084  }
1085  }
1086  }
1087  }
1088 
1089  double usedCDCHitFraction = sumCDCWeights / double(recoTrack.getNumberOfCDCHits());
1090  getObjectPtr<TH1F>("cdc_hit_fraction")->Fill(usedCDCHitFraction);
1091  if (usedCDCHitFraction < m_minUsedCDCHitFraction)
1092  return false;
1093  }
1094  } catch (...) {
1095  B2ERROR("Error in checking DAF weights from previous fit to resolve hit ambiguity. Why? Failed fit points in DAF? Skip track to be sure.");
1096  return false;
1097  }
1098 
1099  std::shared_ptr<genfit::GblFitter> gbl(new genfit::GblFitter());
1100  gbl->setOptions(m_internalIterations, true, true, m_externalIterations, m_recalcJacobians);
1101  gbl->setTrackSegmentController(new GblMultipleScatteringController);
1102 
1103  MeasurementAdder factory("", "", "", "", "");
1104 
1105  // We need the store arrays
1111 
1112  // Create the genfit::MeasurementFactory
1113  genfit::MeasurementFactory<genfit::AbsMeasurement> genfitMeasurementFactory;
1114 
1115  // Add producer for alignable RecoHits to factory
1116  if (pxdHits.isOptional()) {
1119  genfitMeasurementFactory.addProducer(Const::PXD, PXDProducer);
1120  }
1121 
1122  if (svdHits.isOptional()) {
1125  genfitMeasurementFactory.addProducer(Const::SVD, SVDProducer);
1126  }
1127 
1128  if (cdcHits.isOptional()) {
1131  genfitMeasurementFactory.addProducer(Const::CDC, CDCProducer);
1132  }
1133 
1134  if (bklmHits.isOptional()) {
1137  genfitMeasurementFactory.addProducer(Const::BKLM, BKLMProducer);
1138  }
1139 
1140  if (eklmHits.isOptional()) {
1143  genfitMeasurementFactory.addProducer(Const::EKLM, EKLMProducer);
1144  }
1145 
1146 
1147  // Create the measurement creators
1148  std::vector<std::shared_ptr<PXDBaseMeasurementCreator>> pxdMeasurementCreators = { std::shared_ptr<PXDBaseMeasurementCreator>(new PXDCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1149  std::vector<std::shared_ptr<SVDBaseMeasurementCreator>> svdMeasurementCreators = { std::shared_ptr<SVDBaseMeasurementCreator>(new SVDCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1150  // TODO: Create a new MeasurementCreator based on SVDBaseMeasurementCreator (or on SVDCoordinateMeasurementCreator), which does the combination on the fly.
1151 
1152  std::vector<std::shared_ptr<CDCBaseMeasurementCreator>> cdcMeasurementCreators = { std::shared_ptr<CDCBaseMeasurementCreator>(new CDCCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1153  std::vector<std::shared_ptr<BKLMBaseMeasurementCreator>> bklmMeasurementCreators = { std::shared_ptr<BKLMBaseMeasurementCreator>(new BKLMCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1154  std::vector<std::shared_ptr<EKLMBaseMeasurementCreator>> eklmMeasurementCreators = { std::shared_ptr<EKLMBaseMeasurementCreator>(new EKLMCoordinateMeasurementCreator(genfitMeasurementFactory)) };
1155 
1156  // TODO: Or put it in here and leave the svdMeasurementCreators empty.
1157  std::vector<std::shared_ptr<BaseMeasurementCreator>> additionalMeasurementCreators = {};
1158  factory.resetMeasurementCreators(pxdMeasurementCreators, svdMeasurementCreators, cdcMeasurementCreators, bklmMeasurementCreators,
1159  eklmMeasurementCreators, additionalMeasurementCreators);
1160  factory.addMeasurements(recoTrack);
1161 
1162  auto& gfTrack = RecoTrackGenfitAccess::getGenfitTrack(recoTrack);
1163 
1164  int currentPdgCode = TrackFitter::createCorrectPDGCodeForChargedStable(Const::muon, recoTrack);
1165  if (particle)
1166  currentPdgCode = particle->getPDGCode();
1167 
1168  genfit::AbsTrackRep* trackRep = RecoTrackGenfitAccess::createOrReturnRKTrackRep(recoTrack, currentPdgCode);
1169  gfTrack.setCardinalRep(gfTrack.getIdForRep(trackRep));
1170 
1171  if (particle) {
1172  B2Vector3D vertexPos = particle->getVertex();
1173  B2Vector3D vertexMom = particle->getMomentum();
1174  gfTrack.setStateSeed(vertexPos, vertexMom);
1175 
1176  genfit::StateOnPlane vertexSOP(gfTrack.getCardinalRep());
1177  B2Vector3D vertexRPhiDir(vertexPos[0], vertexPos[1], 0);
1178  B2Vector3D vertexZDir(0, 0, vertexPos[2]);
1179  //FIXME: This causes problem to current GBL version in genfit -> needs update of GBL to re-enable
1180  // genfit::SharedPlanePtr vertexPlane(new genfit::DetPlane(vertexPos, vertexRPhiDir, vertexZDir));
1181  //This works instead fine:
1182  genfit::SharedPlanePtr vertexPlane(new genfit::DetPlane(vertexPos, vertexMom));
1183 
1184  vertexSOP.setPlane(vertexPlane);
1185  vertexSOP.setPosMom(vertexPos, vertexMom);
1186  TMatrixDSym vertexCov(5);
1187  vertexCov.UnitMatrix();
1188  // By using negative covariance no measurement is added to GBL. But this first point
1189  // is then used as additional point in trajectory at the assumed point of its fitted vertex
1190  vertexCov *= -1.;
1191  genfit::MeasuredStateOnPlane mop(vertexSOP, vertexCov);
1192  genfit::FullMeasurement* vertex = new genfit::FullMeasurement(mop, Const::IR);
1193  gfTrack.insertMeasurement(vertex, 0);
1194  }
1195 
1196  try {
1197  for (unsigned int i = 0; i < gfTrack.getNumPoints() - 1; ++i) {
1198  //if (gfTrack.getPointWithMeasurement(i)->getNumRawMeasurements() != 1)
1199  // continue;
1200  genfit::PlanarMeasurement* planarMeas1 = dynamic_cast<genfit::PlanarMeasurement*>(gfTrack.getPointWithMeasurement(
1201  i)->getRawMeasurement(0));
1202  genfit::PlanarMeasurement* planarMeas2 = dynamic_cast<genfit::PlanarMeasurement*>(gfTrack.getPointWithMeasurement(
1203  i + 1)->getRawMeasurement(0));
1204 
1205  if (planarMeas1 != NULL && planarMeas2 != NULL &&
1206  planarMeas1->getDetId() == planarMeas2->getDetId() &&
1207  planarMeas1->getPlaneId() != -1 && // -1 is default plane id
1208  planarMeas1->getPlaneId() == planarMeas2->getPlaneId()) {
1209  Belle2::AlignableSVDRecoHit* hit1 = dynamic_cast<Belle2::AlignableSVDRecoHit*>(planarMeas1);
1210  Belle2::AlignableSVDRecoHit* hit2 = dynamic_cast<Belle2::AlignableSVDRecoHit*>(planarMeas2);
1211  if (hit1 && hit2) {
1212  Belle2::AlignableSVDRecoHit* hitU(NULL);
1213  Belle2::AlignableSVDRecoHit* hitV(NULL);
1214  // We have to decide U/V now (else AlignableSVDRecoHit2D could throw FATAL)
1215  if (hit1->isU() && !hit2->isU()) {
1216  hitU = hit1;
1217  hitV = hit2;
1218  } else if (!hit1->isU() && hit2->isU()) {
1219  hitU = hit2;
1220  hitV = hit1;
1221  } else {
1222  continue;
1223  }
1225  // insert measurement before point i (increases number of currect point to i+1)
1226  gfTrack.insertMeasurement(hit, i);
1227  // now delete current point (at its original place, we have the new 2D recohit)
1228  gfTrack.deletePoint(i + 1);
1229  gfTrack.deletePoint(i + 1);
1230  }
1231  }
1232  }
1233  } catch (std::exception& e) {
1234  B2ERROR(e.what());
1235  B2ERROR("SVD Cluster combination failed. This is symptomatic of pruned tracks. MillepedeCollector cannot process pruned tracks.");
1236  return false;
1237  }
1238 
1239  try {
1240  gbl->processTrackWithRep(&gfTrack, gfTrack.getCardinalRep(), true);
1241  } catch (genfit::Exception& e) {
1242  B2ERROR(e.what());
1243  return false;
1244  } catch (...) {
1245  B2ERROR("GBL fit failed.");
1246  return false;
1247  }
1248 
1249  return true;
1250 }
1251 
1252 std::vector< genfit::Track* > MillepedeCollectorModule::getParticlesTracks(std::vector<Particle*> particles, bool addVertexPoint)
1253 {
1254  std::vector< genfit::Track* > tracks;
1255  for (auto particle : particles) {
1256  auto belle2Track = particle->getTrack();
1257  if (!belle2Track) {
1258  B2WARNING("No Belle2::Track for particle (particle->X");
1259  continue;
1260  }
1261 // auto trackFitResult = belle2Track->getTrackFitResult(Const::chargedStableSet.find(abs(particle->getPDGCode())));
1262 // if (!trackFitResult) {
1263 // B2INFO("No track fit result for track");
1264 // continue;
1265 // }
1266 // auto recoTrack = trackFitResult->getRelatedFrom<RecoTrack>();
1267  auto recoTrack = belle2Track->getRelatedTo<RecoTrack>();
1268 
1269  if (!recoTrack) {
1270  B2WARNING("No related RecoTrack for Belle2::Track (particle->Track->X)");
1271  continue;
1272  }
1273 
1274  // If any track fails, fail completely
1275  if (!fitRecoTrack(*recoTrack, (addVertexPoint) ? particle : nullptr))
1276  return {};
1277 
1278  auto& track = RecoTrackGenfitAccess::getGenfitTrack(*recoTrack);
1279 
1280  if (!track.hasFitStatus()) {
1281  B2WARNING("Track has no fit status");
1282  continue;
1283  }
1284  genfit::GblFitStatus* fs = dynamic_cast<genfit::GblFitStatus*>(track.getFitStatus());
1285  if (!fs) {
1286  B2WARNING("Track FitStatus is not GblFitStatus.");
1287  continue;
1288  }
1289  if (!fs->isFittedWithReferenceTrack()) {
1290  B2WARNING("Track is not fitted with reference track.");
1291  continue;
1292  }
1293 
1294  tracks.push_back(&track);
1295  }
1296 
1297  return tracks;
1298 }
1299 
1300 std::pair<TMatrixD, TMatrixD> MillepedeCollectorModule::getTwoBodyToLocalTransform(Particle& mother,
1301  double motherMass)
1302 {
1303  std::vector<TMatrixD> result;
1304 
1305  double px = mother.getMomentum()[0];
1306  double py = mother.getMomentum()[1];
1307  double pz = mother.getMomentum()[2];
1308  double pt = sqrt(px * px + py * py);
1309  double p = mother.getMomentumMagnitude();
1310  double M = motherMass;
1311  double m = mother.getDaughter(0)->getPDGMass();
1312 
1313  if (mother.getNDaughters() != 2
1314  || m != mother.getDaughter(1)->getPDGMass()) B2FATAL("Only two same-mass daughters (V0->f+f- decays) allowed.");
1315 
1316  // Rotation matrix from mother reference system to lab system
1317  TMatrixD mother2lab(3, 3);
1318  mother2lab(0, 0) = px * pz / pt / p; mother2lab(0, 1) = - py / pt; mother2lab(0, 2) = px / p;
1319  mother2lab(1, 0) = py * pz / pt / p; mother2lab(1, 1) = px / pt; mother2lab(1, 2) = py / p;
1320  mother2lab(2, 0) = - pt / p; mother2lab(2, 1) = 0; mother2lab(2, 2) = pz / p;
1321  auto lab2mother = mother2lab; lab2mother.Invert();
1322 
1323  // Need to rotate and boost daughters' momenta to know which goes forward (+sign in decay model)
1324  // and to get the angles theta, phi of the decaying daughter system in mothers' reference frame
1325  RestFrame boostedFrame(&mother);
1326  TLorentzVector fourVector1(lab2mother * mother.getDaughter(0)->get4Vector().Vect(), mother.getDaughter(0)->get4Vector().T());
1327  TLorentzVector fourVector2(lab2mother * mother.getDaughter(1)->get4Vector().Vect(), mother.getDaughter(1)->get4Vector().T());
1328 
1329  auto mom1 = boostedFrame.getMomentum(fourVector1).Vect();
1330  auto mom2 = boostedFrame.getMomentum(fourVector2).Vect();
1331  // One momentum has opposite direction (otherwise should be same in CMS of mother), but which?
1332  double sign = 1.;
1333  auto avgMom = 0.5 * (mom1 - mom2);
1334  if (avgMom[2] < 0.) {
1335  avgMom *= -1.;
1336  // switch meaning of plus/minus trajectories
1337  sign = -1.;
1338  }
1339 
1340  double theta = atan2(avgMom.Perp(), avgMom[2]);
1341  double phi = atan2(avgMom[1], avgMom[0]);
1342  if (phi < 0.) phi += 2. * TMath::Pi();
1343 
1344  double alpha = M / 2. / m;
1345  double c1 = m * sqrt(alpha * alpha - 1.);
1346  double c2 = 0.5 * sqrt((alpha * alpha - 1.) / alpha / alpha * (p * p + M * M));
1347 
1348  double p3 = p * p * p;
1349  double pt3 = pt * pt * pt;
1350 
1351 
1352  for (auto& track : getParticlesTracks(mother.getDaughters())) {
1353 
1354 
1355  TMatrixD R = mother2lab;
1356  B2Vector3D P(sign * c1 * sin(theta) * cos(phi),
1357  sign * c1 * sin(theta) * sin(phi),
1358  p / 2. + sign * c2 * cos(theta));
1359 
1360  TMatrixD dRdpx(3, 3);
1361  dRdpx(0, 0) = - pz * (pow(px, 4.) - pow(py, 4.) - py * py * pz * pz) / pt3 / p3;
1362  dRdpx(0, 1) = px * py / pt3;
1363  dRdpx(0, 2) = (py * py + pz * pz) / p3;
1364 
1365  dRdpx(1, 0) = - px * py * pz * (2. * px * px + 2. * py * py + pz * pz) / pt3 / p3;
1366  dRdpx(1, 1) = - py * py / pt3;
1367  dRdpx(1, 2) = px * py / p3;
1368 
1369  dRdpx(2, 0) = - px * pz * pz / pt / p3;
1370  dRdpx(2, 1) = 0.;
1371  dRdpx(2, 2) = - px * pz / p3;
1372 
1373  TMatrixD dRdpy(3, 3);
1374  dRdpy(0, 0) = - px * py * pz * (2. * px * px + 2. * py * py + pz * pz) / pt3 / p3;
1375  dRdpy(0, 1) = - px * px / pt3;
1376  dRdpy(0, 2) = px * pz / p3;
1377 
1378  dRdpy(1, 0) = - pz * (- pow(px, 4.) - px * px * pz * pz + pow(py, 4.)) / pt3 / p3;
1379  dRdpy(1, 1) = px * py / pt3;
1380  dRdpy(1, 2) = (px * px + pz * pz) / p3;
1381 
1382  dRdpy(2, 0) = - py * pz * pz / pt / p3;
1383  dRdpy(2, 1) = 0.;
1384  dRdpy(2, 2) = - py * pz / p3;
1385 
1386  TMatrixD dRdpz(3, 3);
1387  dRdpz(0, 0) = px * pt / p3;
1388  dRdpz(0, 1) = 0.;
1389  dRdpz(0, 2) = - px * pz / p3;
1390 
1391  dRdpz(1, 0) = py * pt / p3;
1392  dRdpz(1, 1) = 0.;
1393  dRdpz(1, 2) = py * pz / p3;
1394 
1395  dRdpz(2, 0) = pz * pt / p3;
1396  dRdpz(2, 1) = 0.;
1397  dRdpz(2, 2) = (px * px + py * py) / p3;
1398 
1399  B2Vector3D dpdpx = dRdpx * P;
1400  B2Vector3D dpdpy = dRdpy * P;
1401  B2Vector3D dpdpz = dRdpz * P;
1402 
1403  B2Vector3D dpdtheta = R * B2Vector3D(sign * c1 * cos(theta) * cos(phi),
1404  sign * c1 * cos(theta) * sin(phi),
1405  p / 2. + sign * c2 * (- sin(phi)));
1406 
1407 
1408  B2Vector3D dpdphi = R * B2Vector3D(sign * c1 * sin(theta) * (- sin(phi)),
1409  sign * c1 * sin(theta) * cos(phi),
1410  0.);
1411 
1412  double dc1dM = m * M / (2. * sqrt(M * M - 4. * m * m));
1413  double dc2dM = M * (4. * m * m * p * p + pow(M, 4)) / (2 * M * M * M * sqrt((M * M - 4. * m * m) * (p * p + M * M)));
1414 
1415  B2Vector3D dpdM = R * B2Vector3D(sign * sin(theta) * cos(phi) * dc1dM,
1416  sign * sin(theta) * sin(phi) * dc1dM,
1417  sign * cos(theta) * dc2dM);
1418 
1419  TMatrixD dpdz(3, 6);
1420  dpdz(0, 0) = dpdpx(0); dpdz(0, 1) = dpdpy(0); dpdz(0, 2) = dpdpz(0); dpdz(0, 3) = dpdtheta(0); dpdz(0, 4) = dpdphi(0);
1421  dpdz(0, 5) = dpdM(0);
1422  dpdz(1, 0) = dpdpx(1); dpdz(1, 1) = dpdpy(1); dpdz(1, 2) = dpdpz(1); dpdz(1, 3) = dpdtheta(1); dpdz(1, 4) = dpdphi(1);
1423  dpdz(1, 5) = dpdM(1);
1424  dpdz(2, 0) = dpdpx(2); dpdz(2, 1) = dpdpy(2); dpdz(2, 2) = dpdpz(2); dpdz(2, 3) = dpdtheta(2); dpdz(2, 4) = dpdphi(2);
1425  dpdz(2, 5) = dpdM(2);
1426 
1427  TMatrixD dqdv = getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 0, 2);
1428  TMatrixD dqdp = getGlobalToLocalTransform(track->getFittedState()).GetSub(0, 4, 3, 5);
1429  TMatrixD dfdvz(5, 9);
1430  dfdvz.SetSub(0, 0, dqdv);
1431  dfdvz.SetSub(0, 3, dqdp * dpdz);
1432 
1433  result.push_back(dfdvz);
1434 
1435  // switch sign for second trajectory
1436  sign *= -1.;
1437  }
1438 
1439  return {result[0], result[1]};
1440 }
1441 
1442 TMatrixD MillepedeCollectorModule::getGlobalToLocalTransform(const genfit::MeasuredStateOnPlane& msop)
1443 {
1444  auto state = msop;
1445  const B2Vector3D& U(state.getPlane()->getU());
1446  const B2Vector3D& V(state.getPlane()->getV());
1447  const B2Vector3D& O(state.getPlane()->getO());
1448  const B2Vector3D& W(state.getPlane()->getNormal());
1449 
1450  const double* state5 = state.getState().GetMatrixArray();
1451 
1452  double spu = 1.;
1453 
1454  const TVectorD& auxInfo = state.getAuxInfo();
1455  if (auxInfo.GetNrows() == 2
1456  || auxInfo.GetNrows() == 1) // backwards compatibility with old RKTrackRep
1457  spu = state.getAuxInfo()(0);
1458 
1459  TVectorD state7(7);
1460 
1461  state7[0] = O.X() + state5[3] * U.X() + state5[4] * V.X(); // x
1462  state7[1] = O.Y() + state5[3] * U.Y() + state5[4] * V.Y(); // y
1463  state7[2] = O.Z() + state5[3] * U.Z() + state5[4] * V.Z(); // z
1464 
1465  state7[3] = spu * (W.X() + state5[1] * U.X() + state5[2] * V.X()); // a_x
1466  state7[4] = spu * (W.Y() + state5[1] * U.Y() + state5[2] * V.Y()); // a_y
1467  state7[5] = spu * (W.Z() + state5[1] * U.Z() + state5[2] * V.Z()); // a_z
1468 
1469  // normalize dir
1470  double norm = 1. / sqrt(state7[3] * state7[3] + state7[4] * state7[4] + state7[5] * state7[5]);
1471  for (unsigned int i = 3; i < 6; ++i) state7[i] *= norm;
1472 
1473  state7[6] = state5[0]; // q/p
1474 
1475  const double AtU = state7[3] * U.X() + state7[4] * U.Y() + state7[5] * U.Z();
1476  const double AtV = state7[3] * V.X() + state7[4] * V.Y() + state7[5] * V.Z();
1477  const double AtW = state7[3] * W.X() + state7[4] * W.Y() + state7[5] * W.Z();
1478 
1479  // J_Mp matrix is d(q/p,u',v',u,v) / d(x,y,z,px,py,pz) (in is 6x6)
1480 
1481  const double qop = state7[6];
1482  const double p = state.getCharge() / qop; // momentum
1483 
1484  TMatrixD J_Mp_6x5(6, 5);
1485  J_Mp_6x5.Zero();
1486 
1487  //d(u)/d(x,y,z)
1488  J_Mp_6x5(0, 3) = U.X(); // [0][3]
1489  J_Mp_6x5(1, 3) = U.Y(); // [1][3]
1490  J_Mp_6x5(2, 3) = U.Z(); // [2][3]
1491  //d(v)/d(x,y,z)
1492  J_Mp_6x5(0, 4) = V.X(); // [0][4]
1493  J_Mp_6x5(1, 4) = V.Y(); // [1][4]
1494  J_Mp_6x5(2, 4) = V.Z(); // [2][4]
1495 
1496  // d(q/p)/d(px,py,pz)
1497  double fact = (-1.) * qop / p;
1498  J_Mp_6x5(3, 0) = fact * state7[3]; // [3][0]
1499  J_Mp_6x5(4, 0) = fact * state7[4]; // [4][0]
1500  J_Mp_6x5(5, 0) = fact * state7[5]; // [5][0]
1501  // d(u')/d(px,py,pz)
1502  fact = 1. / (p * AtW * AtW);
1503  J_Mp_6x5(3, 1) = fact * (U.X() * AtW - W.X() * AtU); // [3][1]
1504  J_Mp_6x5(4, 1) = fact * (U.Y() * AtW - W.Y() * AtU); // [4][1]
1505  J_Mp_6x5(5, 1) = fact * (U.Z() * AtW - W.Z() * AtU); // [5][1]
1506  // d(v')/d(px,py,pz)
1507  J_Mp_6x5(3, 2) = fact * (V.X() * AtW - W.X() * AtV); // [3][2]
1508  J_Mp_6x5(4, 2) = fact * (V.Y() * AtW - W.Y() * AtV); // [4][2]
1509  J_Mp_6x5(5, 2) = fact * (V.Z() * AtW - W.Z() * AtV); // [5][2]
1510 
1511  return J_Mp_6x5.T();
1512 }
1513 
1514 TMatrixD MillepedeCollectorModule::getLocalToGlobalTransform(const genfit::MeasuredStateOnPlane& msop)
1515 {
1516  auto state = msop;
1517  // get vectors and aux variables
1518  const B2Vector3D& U(state.getPlane()->getU());
1519  const B2Vector3D& V(state.getPlane()->getV());
1520  const B2Vector3D& W(state.getPlane()->getNormal());
1521 
1522  const TVectorD& state5(state.getState());
1523  double spu = 1.;
1524 
1525  const TVectorD& auxInfo = state.getAuxInfo();
1526  if (auxInfo.GetNrows() == 2
1527  || auxInfo.GetNrows() == 1) // backwards compatibility with old RKTrackRep
1528  spu = state.getAuxInfo()(0);
1529 
1530  TVectorD pTilde(3);
1531  pTilde[0] = spu * (W.X() + state5(1) * U.X() + state5(2) * V.X()); // a_x
1532  pTilde[1] = spu * (W.Y() + state5(1) * U.Y() + state5(2) * V.Y()); // a_y
1533  pTilde[2] = spu * (W.Z() + state5(1) * U.Z() + state5(2) * V.Z()); // a_z
1534 
1535  const double pTildeMag = sqrt(pTilde[0] * pTilde[0] + pTilde[1] * pTilde[1] + pTilde[2] * pTilde[2]);
1536  const double pTildeMag2 = pTildeMag * pTildeMag;
1537 
1538  const double utpTildeOverpTildeMag2 = (U.X() * pTilde[0] + U.Y() * pTilde[1] + U.Z() * pTilde[2]) / pTildeMag2;
1539  const double vtpTildeOverpTildeMag2 = (V.X() * pTilde[0] + V.Y() * pTilde[1] + V.Z() * pTilde[2]) / pTildeMag2;
1540 
1541  //J_pM matrix is d(x,y,z,px,py,pz) / d(q/p,u',v',u,v) (out is 6x6)
1542 
1543  const double qop = state5(0);
1544  const double p = state.getCharge() / qop; // momentum
1545 
1546  TMatrixD J_pM_5x6(5, 6);
1547  J_pM_5x6.Zero();
1548 
1549  // d(px,py,pz)/d(q/p)
1550  double fact = -1. * p / (pTildeMag * qop);
1551  J_pM_5x6(0, 3) = fact * pTilde[0]; // [0][3]
1552  J_pM_5x6(0, 4) = fact * pTilde[1]; // [0][4]
1553  J_pM_5x6(0, 5) = fact * pTilde[2]; // [0][5]
1554  // d(px,py,pz)/d(u')
1555  fact = p * spu / pTildeMag;
1556  J_pM_5x6(1, 3) = fact * (U.X() - pTilde[0] * utpTildeOverpTildeMag2); // [1][3]
1557  J_pM_5x6(1, 4) = fact * (U.Y() - pTilde[1] * utpTildeOverpTildeMag2); // [1][4]
1558  J_pM_5x6(1, 5) = fact * (U.Z() - pTilde[2] * utpTildeOverpTildeMag2); // [1][5]
1559  // d(px,py,pz)/d(v')
1560  J_pM_5x6(2, 3) = fact * (V.X() - pTilde[0] * vtpTildeOverpTildeMag2); // [2][3]
1561  J_pM_5x6(2, 4) = fact * (V.Y() - pTilde[1] * vtpTildeOverpTildeMag2); // [2][4]
1562  J_pM_5x6(2, 5) = fact * (V.Z() - pTilde[2] * vtpTildeOverpTildeMag2); // [2][5]
1563  // d(x,y,z)/d(u)
1564  J_pM_5x6(3, 0) = U.X(); // [3][0]
1565  J_pM_5x6(3, 1) = U.Y(); // [3][1]
1566  J_pM_5x6(3, 2) = U.Z(); // [3][2]
1567  // d(x,y,z)/d(v)
1568  J_pM_5x6(4, 0) = V.X(); // [4][0]
1569  J_pM_5x6(4, 1) = V.Y(); // [4][1]
1570  J_pM_5x6(4, 2) = V.Z(); // [4][2]
1571 
1572  return J_pM_5x6.T();
1573 
1574 }
1575 
1576 tuple<B2Vector3D, TMatrixDSym> MillepedeCollectorModule::getPrimaryVertexAndCov() const
1577 {
1578  DBObjPtr<BeamSpot> beam;
1579  return {beam->getIPPosition(), beam->getSizeCovMatrix()};
1580 }
1581 
1582 void MillepedeCollectorModule::updateMassWidthIfSet(string listName, double& mass, double& width)
1583 {
1584  if (m_customMassConfig.find(listName) != m_customMassConfig.end()) {
1585  auto massWidth = m_customMassConfig.at(listName);
1586  mass = std::get<0>(massWidth);
1587  width = std::get<1>(massWidth);
1588  }
1589 }
1590 
Belle2::CalibrationCollectorModule
Calibration collector module base class.
Definition: CalibrationCollectorModule.h:44
genfit::Exception
Exception class for error handling in GENFIT (provides storage for diagnostic information)
Definition: Exception.h:48
genfit::TrackPoint
Object containing AbsMeasurement and AbsFitterInfo objects.
Definition: TrackPoint.h:46
Belle2::alignment::GlobalCalibrationManager::initialize
void initialize(const std::vector< std::string > &components={}, const std::vector< EventMetaData > &timeSlices={})
Initialize the manager with given configuration (from MillepedeCollector)
Definition: Manager.cc:45
Belle2::RecoTrack::getTrackFitStatus
const genfit::FitStatus * getTrackFitStatus(const genfit::AbsTrackRep *representation=nullptr) const
Return the track fit status for the given representation or for the cardinal one. You are not allowed...
Definition: RecoTrack.h:537
genfit::SharedPlanePtr
std::shared_ptr< genfit::DetPlane > SharedPlanePtr
Shared Pointer to a DetPlane.
Definition: SharedPlanePtr.h:40
Belle2::alignment::VXDGlobalParamInterface::s_hierarchyType
static E_VXDHierarchyType s_hierarchyType
What type of hierarchy to use for VXD?
Definition: GlobalParam.h:92
Belle2::MeasurementAdder::addMeasurements
bool addMeasurements(RecoTrack &recoTrack) const
After you have filled the internal storage with measurement creators (either by providing your own or...
Definition: MeasurementAdder.cc:162
Belle2::RecoTrack::getCardinalRepresentation
genfit::AbsTrackRep * getCardinalRepresentation() const
Get a pointer to the cardinal track representation. You are not allowed to modify or delete it!
Definition: RecoTrack.h:547
genfit::MeasuredStateOnPlane
#StateOnPlane with additional covariance matrix.
Definition: MeasuredStateOnPlane.h:39
genfit::GblFitStatus
FitStatus for use with GblFitter.
Definition: GblFitStatus.h:39
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::CDCHit
Class containing the result of the unpacker in raw data and the result of the digitizer in simulation...
Definition: CDCHit.h:51
Belle2::CDCCoordinateMeasurementCreator
CoordinateMeasurementCreator< RecoHitInformation::UsedCDCHit, Const::CDC > CDCCoordinateMeasurementCreator
Needed for templating.
Definition: CoordinateMeasurementCreator.h:46
genfit::StateOnPlane
A state with arbitrary dimension defined in a DetPlane.
Definition: StateOnPlane.h:47
genfit::AbsTrackRep
Abstract base class for a track representation.
Definition: AbsTrackRep.h:66
Belle2::B2Vector3::Z
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:434
genfit::MeasurementProducer
Template class for a measurement producer module.
Definition: MeasurementProducer.h:76
Belle2::EKLMAlignmentHit
This dataobject is used only for EKLM alignment.
Definition: EKLMAlignmentHit.h:39
Belle2::FileMetaData::setParents
void setParents(const std::vector< std::string > &parents)
Parents setter.
Definition: FileMetaData.h:173
gbl::GblTrajectory
GBL trajectory.
Definition: GblTrajectory.h:48
genfit::MeasurementFactory< genfit::AbsMeasurement >
Belle2::GlobalLabel
Class to convert to/from global labels for Millepede II to/from detector & parameter identificators.
Definition: GlobalLabel.h:51
Belle2::RelationsInterface::getRelationsTo
RelationVector< TO > getRelationsTo(const std::string &name="", const std::string &namedRelation="") const
Get the relations that point from this object to another store array.
Definition: RelationsObject.h:199
genfit::PlanarMeasurement
Measurement class implementing a planar hit geometry (1 or 2D).
Definition: PlanarMeasurement.h:44
Belle2::FileMetaData
Metadata information about a file.
Definition: FileMetaData.h:39
Belle2::MeasurementAdder::resetMeasurementCreators
void resetMeasurementCreators(const std::vector< std::shared_ptr< PXDBaseMeasurementCreator >> &pxdMeasurementCreators, const std::vector< std::shared_ptr< SVDBaseMeasurementCreator >> &svdMeasurementCreators, const std::vector< std::shared_ptr< CDCBaseMeasurementCreator >> &cdcMeasurementCreators, const std::vector< std::shared_ptr< BKLMBaseMeasurementCreator >> &bklmMeasurementCreators, const std::vector< std::shared_ptr< EKLMBaseMeasurementCreator >> &eklmMeasurementCreators, const std::vector< std::shared_ptr< BaseMeasurementCreator >> &additionalMeasurementCreators)
If you want to use non-default settings for the store arrays, you can create your own instances of th...
Definition: MeasurementAdder.cc:105
gbl
Namespace for the general broken lines package.
Definition: BorderedBandMatrix.h:44
Belle2::BKLMCoordinateMeasurementCreator
CoordinateMeasurementCreator< RecoHitInformation::UsedBKLMHit, Const::BKLM > BKLMCoordinateMeasurementCreator
Hit to reco hit measurement creator for the BKLM.
Definition: CoordinateMeasurementCreator.h:52
Belle2::DBObjPtr
Class for accessing objects in the database.
Definition: DBObjPtr.h:31
Belle2::B2Vector3< double >
Belle2::RecoTrack::getCreatedTrackPoint
const genfit::TrackPoint * getCreatedTrackPoint(const RecoHitInformation *recoHitInformation) const
Get a pointer to the TrackPoint that was created from this hit.
Definition: RecoTrack.cc:220
Belle2::MilleData
Mergeable class holding list of so far opened mille binaries and providing the binaries.
Definition: MilleData.h:17
Belle2::AlignableEKLMRecoHit
Alignable EKLM hit.
Definition: AlignableEKLMRecoHit.h:40
genfit::TrackPoint::getFitterInfo
AbsFitterInfo * getFitterInfo(const AbsTrackRep *rep=nullptr) const
Get fitterInfo for rep. Per default, use cardinal rep.
Definition: TrackPoint.cc:170
Belle2::FileMetaData::setLfn
void setLfn(const std::string &lfn)
Setter for LFN.
Definition: FileMetaData.h:145
Belle2::RecoTrack
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:78
Belle2::SVDCoordinateMeasurementCreator
CoordinateMeasurementCreator< RecoHitInformation::UsedSVDHit, Const::SVD > SVDCoordinateMeasurementCreator
Hit to reco hit measurement creator for the SVD.
Definition: CoordinateMeasurementCreator.h:48
genfit::KalmanFitterInfo
Collects information needed and produced by a AbsKalmanFitter implementations and is specific to one ...
Definition: KalmanFitterInfo.h:44
Belle2::B2Vector3D
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:507
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::AlignableSVDRecoHit2D
This class is used to transfer SVD information to the track fit.
Definition: AlignableSVDRecoHit2D.h:30
Belle2::SVDRecoHit::isU
bool isU() const
Is the coordinate u or v?
Definition: SVDRecoHit.h:100
Belle2::GblMultipleScatteringController
TrackSegmentController for use with GblFitter in Belle2.
Definition: GblMultipleScatteringController.h:36
Belle2::RestFrame
Rest frame of a particle.
Definition: ReferenceFrame.h:141
Belle2::alignment::GlobalCalibrationManager::writeConstraints
void writeConstraints(std::string txtFilename)
Write-out complete hierarchy to a text file.
Definition: Manager.cc:162
genfit::MeasurementFactory::addProducer
void addProducer(int detID, AbsMeasurementProducer< measurement_T > *hitProd)
Register a producer module to the factory.
Definition: MeasurementFactory.h:97
Belle2::alignment::VXDGlobalParamInterface::s_enablePXD
static bool s_enablePXD
Enable PXD in hierarchy?
Definition: GlobalParam.h:94
Belle2::RecoHitInformation
This class stores additional information to every CDC/SVD/PXD hit stored in a RecoTrack.
Definition: RecoHitInformation.h:48
Belle2::alignment::GlobalDerivatives
Class for easier manipulation with global derivatives (and their labels)
Definition: GlobalDerivatives.h:37
Belle2::PXDCluster
The PXD Cluster class This class stores all information about reconstructed PXD clusters The position...
Definition: PXDCluster.h:41
Belle2::alignment::GlobalCalibrationManager::getInstance
static GlobalCalibrationManager & getInstance()
Get instance of the Manager auto& gcm = GlobalCalibrationManager::getInstance();.
Definition: Manager.cc:20
Belle2::StoreArray::isValid
bool isValid() const
Check wether the array was registered.
Definition: StoreArray.h:298
Belle2::AlignableBKLMRecoHit
Alignable BKLM hit.
Definition: AlignableBKLMRecoHit.h:40
Belle2::EKLMCoordinateMeasurementCreator
CoordinateMeasurementCreator< RecoHitInformation::UsedEKLMHit, Const::EKLM > EKLMCoordinateMeasurementCreator
Hit to reco hit measurement creator for the EKLM.
Definition: CoordinateMeasurementCreator.h:54
Belle2::MeasurementAdder
Algorithm class to translate the added detector hits (e.g.
Definition: MeasurementAdder.h:78
Belle2::MillepedeCollectorModule
Calibration data collector for Millepede Algorithm.
Definition: MillepedeCollectorModule.h:42
Belle2::StoreArray::getPtr
TClonesArray * getPtr() const
Raw access to the underlying TClonesArray.
Definition: StoreArray.h:321
Belle2::RecoTrack::getNumberOfCDCHits
unsigned int getNumberOfCDCHits() const
Return the number of cdc hits.
Definition: RecoTrack.h:417
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::SVDCluster
The SVD Cluster class This class stores all information about reconstructed SVD clusters.
Definition: SVDCluster.h:38
Belle2::AlignablePXDRecoHit
This class is used to transfer PXD information to the track fit.
Definition: AlignablePXDRecoHit.h:30
Belle2::RestFrame::getMomentum
virtual TLorentzVector getMomentum(const TLorentzVector &vector) const override
Get Lorentz vector in rest frame System.
Definition: ReferenceFrame.cc:49
Belle2::EventMetaData
Store event, run, and experiment numbers.
Definition: EventMetaData.h:43
Belle2::alignment::VXDGlobalParamInterface::s_enableSVD
static bool s_enableSVD
Enable SVD in hierarchy?
Definition: GlobalParam.h:96
genfit::FullMeasurement
Measurement class implementing a measurement of all track parameters.
Definition: FullMeasurement.h:42
genfit::DetPlane
Detector plane.
Definition: DetPlane.h:59
Belle2::AlignableSVDRecoHit
This class is used to transfer SVD information to the track fit.
Definition: AlignableSVDRecoHit.h:30
gbl::GblTrajectory::isValid
bool isValid() const
Retrieve validity of trajectory.
Definition: GblTrajectory.cc:289
Belle2::AlignableCDCRecoHit
This class is used to transfer CDC information to the track fit and Millepede.
Definition: AlignableCDCRecoHit.h:32
Belle2::StoreArray
Accessor to arrays stored in the data store.
Definition: ECLMatchingPerformanceExpertModule.h:33
Belle2::RecoTrack::getStoreArrayNameOfRecoHitInformation
const std::string & getStoreArrayNameOfRecoHitInformation() const
Name of the store array of the reco hit informations.
Definition: RecoTrack.h:646
Belle2::B2Vector3::X
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:430
genfit::FitStatus::isFitted
bool isFitted() const
Has the track been fitted?
Definition: FitStatus.h:94
Belle2::PXDCoordinateMeasurementCreator
CoordinateMeasurementCreator< RecoHitInformation::UsedPXDHit, Const::PXD > PXDCoordinateMeasurementCreator
Hit to reco hit measurement creator for the PXD.
Definition: CoordinateMeasurementCreator.h:50
genfit::GblFitter
Generic GBL implementation.
Definition: GblFitter.h:53
Belle2::BKLMHit2d
Store one BKLM strip hit as a ROOT object.
Definition: BKLMHit2d.h:42
Belle2::B2Vector3::Y
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:432
gbl::GblTrajectory::fit
unsigned int fit(double &Chi2, int &Ndf, double &lostWeight, std::string optionList="")
Perform fit of (valid) trajectory.
Definition: GblTrajectory.cc:1043
Belle2::StoreObjPtr::isValid
bool isValid() const
Check whether the object was created.
Definition: StoreObjPtr.h:120