Belle II Software  light-2212-foldex
B2BIIConvertMdstModule.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 <b2bii/modules/B2BIIMdstInput/B2BIIConvertMdstModule.h>
10 
11 #include <framework/datastore/StoreObjPtr.h>
12 #include <framework/datastore/RelationArray.h>
13 #include <framework/database/Database.h>
14 #include <framework/pcore/ProcHandler.h>
15 
16 #include <mdst/dataobjects/HitPatternVXD.h>
17 #include <mdst/dataobjects/HitPatternCDC.h>
18 
19 // Belle II utilities
20 #include <framework/gearbox/Unit.h>
21 #include <framework/gearbox/Const.h>
22 #include <framework/geometry/B2Vector3.h>
23 #include <analysis/dataobjects/ParticleExtraInfoMap.h>
24 
25 // Belle II dataobjects
26 #include <framework/dataobjects/EventMetaData.h>
27 #include <framework/dataobjects/Helix.h>
28 #include <framework/dataobjects/UncertainHelix.h>
29 
30 // Belle utilities
31 #include <b2bii/utility/BelleMdstToGenHepevt.h>
32 
33 // ROOT
34 #include <Math/RotationY.h>
35 #include <Math/Vector3D.h>
36 #include <Math/Vector4D.h>
37 
38 #include <limits>
39 #include <algorithm>
40 #include <queue>
41 #include <utility>
42 
43 #ifdef HAVE_EID
44 #include "belle_legacy/eid/eid.h"
45 #endif
46 
47 #ifdef HAVE_KID_ACC
48 #include "belle_legacy/kid/kid_acc.h"
49 #include "belle_legacy/kid/kid_cdc.h"
50 #endif
51 
52 #ifdef HAVE_FINDKS
53 #include "belle_legacy/findKs/findKs.h"
54 #endif
55 
56 #ifdef HAVE_NISKSFINDER
57 #include "belle_legacy/nisKsFinder/nisKsFinder.h"
58 #endif
59 
60 #ifdef HAVE_GOODLAMBDA
61 #include "belle_legacy/findLambda/findLambda.h"
62 #endif
63 
64 #include "belle_legacy/benergy/BeamEnergy.h"
65 #include "belle_legacy/ip/IpProfile.h"
66 #include "belle_legacy/tables/evtcls.h"
67 #include "belle_legacy/tables/trg.h"
68 
69 
70 #include <cmath>
71 #include <cfloat>
72 #include <bitset>
73 using namespace Belle2;
74 
76 
77 bool approximatelyEqual(float a, float b, float epsilon)
78 {
79  return fabs(a - b) <= ((fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon);
80 }
81 
82 double adjustAngleRange(double phi)
83 {
84  phi = phi - int(phi / TMath::TwoPi()) * TMath::TwoPi();
85  return phi - int(phi / TMath::Pi()) * TMath::TwoPi();
86 }
87 
88 void fill7x7ErrorMatrix(const TrackFitResult* tfr, TMatrixDSym& error7x7, const double mass, const double bField)
89 {
90  short charge = tfr->getChargeSign();
91 
92  double d0 = tfr->getD0();
93  double phi0 = tfr->getPhi0();
94  double omega = tfr->getOmega();
95  //double z0 = tfr->getZ0();
96  double tanl = tfr->getTanLambda();
97 
98  double alpha = tfr->getHelix().getAlpha(bField);
99 
100  double cosPhi0 = TMath::Cos(phi0);
101  double sinPhi0 = TMath::Sin(phi0);
102 
103  double rho;
104  if (omega != 0)
105  rho = 1.0 / alpha / omega;
106  else
107  rho = (DBL_MAX);
108 
109  double energy = TMath::Sqrt(mass * mass + (1.0 + tanl * tanl) * rho * rho);
110 
111  const int iPx = 0;
112  const int iPy = 1;
113  const int iPz = 2;
114  const int iE = 3;
115  const int iX = 4;
116  const int iY = 5;
117  const int iZ = 6;
118 
119  const int iD0 = 0;
120  const int iPhi0 = 1;
121  const int iOmega = 2;
122  const int iZ0 = 3;
123  const int iTanl = 4;
124 
125  TMatrixD jacobian(7, 5);
126  jacobian.Zero();
127 
128  jacobian(iPx, iPhi0) = - fabs(rho) * sinPhi0;
129  jacobian(iPx, iOmega) = - charge * rho * rho * cosPhi0 * alpha;
130  jacobian(iPy, iPhi0) = fabs(rho) * cosPhi0;
131  jacobian(iPy, iOmega) = - charge * rho * rho * sinPhi0 * alpha;
132  jacobian(iPz, iOmega) = - charge * rho * rho * tanl * alpha;
133  jacobian(iPz, iTanl) = fabs(rho);
134  if (omega != 0 && energy != 0) {
135  jacobian(iE, iOmega) = - (1.0 + tanl * tanl) * rho * rho / omega / energy;
136  jacobian(iE, iTanl) = tanl * rho * rho / energy;
137  } else {
138  jacobian(iE, iOmega) = (DBL_MAX);
139  jacobian(iE, iTanl) = (DBL_MAX);
140  }
141  jacobian(iX, iD0) = sinPhi0;
142  jacobian(iX, iPhi0) = d0 * cosPhi0;
143  jacobian(iY, iD0) = - cosPhi0;
144  jacobian(iY, iPhi0) = d0 * sinPhi0;
145  jacobian(iZ, iZ0) = 1.0;
146 
147  TMatrixDSym error5x5 = tfr->getCovariance5();
148 
149  error7x7 = error5x5.Similarity(jacobian);
150 }
151 //-----------------------------------------------------------------
152 // Register the Module
153 //-----------------------------------------------------------------
154 REG_MODULE(B2BIIConvertMdst);
155 
156 //-----------------------------------------------------------------
157 // Implementation
158 //-----------------------------------------------------------------
159 
161  m_mcMatchingMode(c_Direct)
162 {
163  //Set module properties
164  setDescription("Converts Belle mDST objects (Panther tables and records) to Belle II mDST objects.");
165 
166  addParam("convertBeamParameters", m_convertBeamParameters,
167  "Convert beam parameters or use information stored in "
168  "Belle II database.", true);
169  addParam("use6x6CovarianceMatrix4Tracks", m_use6x6CovarianceMatrix4Tracks,
170  "Use 6x6 (position, momentum) covariance matrix for charged tracks instead of 5x5 (helix parameters) covariance matrix", false);
171  addParam("mcMatchingMode", m_mcMatchingModeString,
172  "MC matching mode: 'Direct', or 'GeneratorLevel'",
173  std::string("Direct"));
174  addParam("matchType2E9oE25Threshold", m_matchType2E9oE25Threshold,
175  "clusters with a E9/E25 value above this threshold are classified as neutral even if tracks are matched to their connected region (matchType == 2)",
176  -1.1);
177 
178  addParam("convertEvtcls", m_convertEvtcls, "Flag to switch on conversion of Mdst_evtcls", true);
179  addParam("nisKsInfo", m_nisEnable, "Flag to switch on conversion of nisKsFinder info", true);
180  addParam("RecTrg", m_convertRecTrg, "Flag to switch on conversion of rectrg_summary3", false);
181  addParam("TrkExtra", m_convertTrkExtra, " Flag to switch on conversion of first_x,y,z and last_x,y,z from Mdst_trk_fit", true);
182 
183  m_realData = false;
184 
185  B2DEBUG(1, "B2BIIConvertMdst: Constructor done.");
186 }
187 
188 
190 {
191 }
192 
194 {
195  // Initialize Belle II DataStore
197  if (m_mcMatchingModeString == "Direct")
199  else if (m_mcMatchingModeString == "GeneratorLevel")
201  else
202  B2FATAL("Unknown MC matching mode: " << m_mcMatchingModeString);
203  B2INFO("B2BIIConvertMdst: initialized.");
204  if (!m_nisEnable)
205  B2WARNING("nisKsFinder output has been disabled. ksnbVLike, ksnbNoLam, ksnbStandard will not be converted.");
206 }
207 
209 {
210  B2DEBUG(99, "[B2BIIConvertMdstModule::initializeDataStore] initialization of DataStore started");
211 
212  // list here all converted Belle2 objects
213  m_eclClusters.registerInDataStore();
214  m_klmClusters.registerInDataStore();
215  m_tracks.registerInDataStore();
216  m_trackFitResults.registerInDataStore();
217  m_v0s.registerInDataStore();
218  m_particles.registerInDataStore();
219 
221  extraInfoMap.registerInDataStore();
222 
223  if (m_convertEvtcls || m_convertRecTrg) m_evtInfo.registerInDataStore();
224 
225  if (m_convertTrkExtra) m_belleTrkExtra.registerInDataStore();
226 
227  StoreObjPtr<ParticleList> gammaParticleList("gamma:mdst");
228  gammaParticleList.registerInDataStore();
229  StoreObjPtr<ParticleList> pi0ParticleList("pi0:mdst");
230  pi0ParticleList.registerInDataStore();
231  StoreObjPtr<ParticleList> kShortParticleList("K_S0:mdst");
232  kShortParticleList.registerInDataStore();
233  StoreObjPtr<ParticleList> kLongParticleList("K_L0:mdst");
234  kLongParticleList.registerInDataStore();
235  StoreObjPtr<ParticleList> lambdaParticleList("Lambda0:mdst");
236  lambdaParticleList.registerInDataStore();
237  StoreObjPtr<ParticleList> antiLambdaParticleList("anti-Lambda0:mdst");
238  antiLambdaParticleList.registerInDataStore();
239  StoreObjPtr<ParticleList> gammaConversionsParticleList("gamma:v0mdst");
240  gammaConversionsParticleList.registerInDataStore();
241 
242  m_pidLikelihoods.registerInDataStore();
243 
244  // needs to be registered, even if running over data, since this information is available only at the begin_run function
245  m_mcParticles.registerInDataStore();
246 
247  //list here all Relations between Belle2 objects
248  m_tracks.registerRelationTo(m_mcParticles);
249  m_tracks.registerRelationTo(m_pidLikelihoods);
250  if (m_convertTrkExtra) m_tracks.registerRelationTo(m_belleTrkExtra);
251  m_eclClusters.registerRelationTo(m_mcParticles);
252  m_tracks.registerRelationTo(m_eclClusters);
253  m_klmClusters.registerRelationTo(m_tracks);
254  m_klmClusters.registerRelationTo(m_eclClusters);
255  m_particles.registerRelationTo(m_mcParticles);
256  m_particles.registerRelationTo(m_pidLikelihoods);
257 
258  B2DEBUG(99, "[B2BIIConvertMdstModule::initializeDataStore] initialization of DataStore ended");
259 }
260 
261 
263 {
264  B2DEBUG(99, "B2BIIConvertMdst: beginRun called.");
265 
267  //BeamEnergy class updated by fixmdst module in beginRun()
268  Belle::BeamEnergy::begin_run();
270  Belle::BeamEnergy::dump();
271 
272  // load IP data from DB server
273  Belle::IpProfile::begin_run();
274  convertIPProfile(true);
275  Belle::IpProfile::dump();
276  bool usableIP = Belle::IpProfile::usable();
277  B2DEBUG(99, "B2BIIConvertMdst: IpProfile is usable = " << usableIP);
278  }
279 
280  //init eID
281 #ifdef HAVE_EID
282  Belle::eid::init_data();
283  Belle::eid::show_use("ALL");
284 #endif
285 }
286 
287 
289 {
291  // Are we running on MC or DATA?
292  Belle::Belle_event_Manager& evman = Belle::Belle_event_Manager::get_manager();
293  Belle::Belle_event& evt = evman[0];
294 
295  if (evt.ExpMC() == 2)
296  m_realData = false; // <- this is MC sample
297  else
298  m_realData = true; // <- this is real data sample
299 
300  // 0. Convert IPProfile to BeamSpot
302 
303  // Make sure beam parameters are correct: if they are not found in the
304  // database or different from the ones in the database we need to override them
305  if (!m_beamSpotDB || !(m_beamSpot == *m_beamSpotDB) ||
308  B2INFO("No database entry for this run yet, create one");
310  IntervalOfValidity iov(event->getExperiment(), event->getRun(), event->getExperiment(), event->getRun());
311  Database::Instance().storeData("CollisionBoostVector", &m_collisionBoostVector, iov);
312  Database::Instance().storeData("CollisionInvariantMass", &m_collisionInvM, iov);
313  Database::Instance().storeData("BeamSpot", &m_beamSpot, iov);
314  }
315  if (m_realData) {
316  B2ERROR("BeamParameters from condition database are different from converted "
317  "ones, overriding database. Did you make sure the globaltag B2BII is used?");
318  } else {
319  B2INFO("BeamSpot, BoostVector, and InvariantMass from condition database are different from converted "
320  "ones, overriding database");
321  }
323  B2FATAL("Cannot reliably override the Database content in parallel processing "
324  "mode, please run the conversion in single processing mode");
325  }
327  DBStore::Instance().addConstantOverride("CollisionInvariantMass", new CollisionInvariantMass(m_collisionInvM), true);
328  DBStore::Instance().addConstantOverride("BeamSpot", new BeamSpot(m_beamSpot), true);
329  }
330  }
331 
332  // 1. Convert MC information
334 
335  // 2. Convert ECL information
337 
338  // 3. Convert KLM information
340 
341  // 4. Convert Tracking information
343 
344  // 5. Set Track -> ECLCluster relations
346 
347  // 6. Set KLMCluster -> Track, ECLCluster relations
349 
350  // 7. Convert Gamma information
352 
353  // 8. Convert Pi0 information
355 
356  // 9. Convert V0s
358 
359  // 10. Convert KLong information
361 
362  // 11. Convert Evtcls panther table information
364 
365  // 12. Convert trigger information from rectrg_summary3
367 
368 }
369 
370 
371 //-----------------------------------------------------------------------------
372 // CONVERT TABLES
373 //-----------------------------------------------------------------------------
375 {
376  const double Eher = Belle::BeamEnergy::E_HER();
377  const double Eler = Belle::BeamEnergy::E_LER();
378  const double crossingAngle = Belle::BeamEnergy::Cross_angle();
379  const double angleLer = M_PI; //parallel to negative z axis (different from Belle II!)
380  const double angleHer = crossingAngle; //in positive z and x direction, verified to be consistent with Upsilon(4S) momentum
381  const double mass_e = Const::electronMass; //mass of electron: 0.0 in basf, 0.000510998902 in basf2;
382  TMatrixDSym covariance(0); //0 entries = no error
383  HepLorentzVector p_beam = Belle::BeamEnergy::p_beam(); // Testing only
384 
385  // Get four momentum of LER and HER
386  ROOT::Math::PxPyPzEVector P_her(0.0, 0.0, TMath::Sqrt(Eher * Eher - mass_e * mass_e), Eher);
387  ROOT::Math::RotationY rotateAroundYAxis(angleHer);
388  P_her = rotateAroundYAxis(P_her);
389  ROOT::Math::PxPyPzEVector P_ler(0.0, 0.0, TMath::Sqrt(Eler * Eler - mass_e * mass_e), Eler);
390  rotateAroundYAxis.SetAngle(angleLer);
391  P_ler = rotateAroundYAxis(P_ler);
392 
393  // Get four momentum of beam
394  ROOT::Math::PxPyPzEVector P_beam = P_her + P_ler;
395 
396  m_collisionBoostVector.setBoost(B2Vector3D(-P_beam.BoostToCM()), covariance);
397  m_collisionInvM.setMass(P_beam.M(), 0.0, 0.0);
398 
399  B2DEBUG(99, "Beam Energy: E_HER = " << Eher << "; E_LER = " << Eler << "; angle = " << crossingAngle);
400  B2DEBUG(99, "Beam Momentum (pre-convert) : P_X = " << p_beam.px() << "; P_Y = " << p_beam.py() << "; P_Z = " << p_beam.pz());
401  B2DEBUG(99, "Beam Momentum (post-convert) : P_X = " << P_beam.Px() << "; P_Y = " << P_beam.Py() << "; P_Z = " << P_beam.Pz());
402 }
403 
405 {
406  if (!Belle::IpProfile::usable()) {
407  // No IPProfile for this run ...
408  if (beginRun) {
409  // no IPProfile, set vertex to NaN without errors for the full run
411  TVector3(std::numeric_limits<double>::quiet_NaN(),
412  std::numeric_limits<double>::quiet_NaN(),
413  std::numeric_limits<double>::quiet_NaN()
414  ), TMatrixTSym<double>()
415  );
416  }
417  return;
418  }
419  HepPoint3D ip;
420  CLHEP::HepSymMatrix ipErr;
421  if (beginRun) {
422  // use event independent average in begin run
423  ip = Belle::IpProfile::position();
424  ipErr = Belle::IpProfile::position_err();
425  } else {
426  // update evtbin
427  Belle::IpProfile::set_evtbin_number();
428  // check if it changed, if not there's nothing to do
429  if (Belle::IpProfile::EvtBinNo() == m_lastIPProfileBin) return;
430  // get event dependent position and error
431  ip = Belle::IpProfile::e_position();
432  ipErr = Belle::IpProfile::e_position_err();
433  }
434  // reset last ipprofile bin
435  m_lastIPProfileBin = Belle::IpProfile::EvtBinNo();
436 
437  TMatrixDSym cov(ipErr.num_col());
438  for (int i = 0; i < ipErr.num_row(); ++i) {
439  for (int j = 0; j < ipErr.num_col(); ++j) {
440  cov(i, j) = ipErr(i + 1, j + 1);
441  }
442  }
443  m_beamSpot.setIP(TVector3(ip.x(), ip.y(), ip.z()), cov);
444 }
445 
447 {
448  // Relations
449  RelationArray tracksToMCParticles(m_tracks, m_mcParticles);
450 
451  // Loop over all Belle charged tracks
452  Belle::Mdst_charged_Manager& m = Belle::Mdst_charged_Manager::get_manager();
453  for (Belle::Mdst_charged_Manager::iterator chargedIterator = m.begin(); chargedIterator != m.end(); ++chargedIterator) {
454  Belle::Mdst_charged belleTrack = *chargedIterator;
455 
456  auto track = m_tracks.appendNew();
457 
458  // convert MDST_Charged -> Track
459  convertMdstChargedObject(belleTrack, track);
460 
461  convertPIDData(belleTrack, track);
462 
463  if (m_realData)
464  continue;
465 
466  // create Track -> MCParticle relation
467  // step 1: MDSTCharged -> Gen_hepevt
468  const Belle::Gen_hepevt& hep0 = get_hepevt(belleTrack);
469  if (hep0 == 0)
470  continue;
471  const Belle::Gen_hepevt* hep = nullptr;
472  switch (m_mcMatchingMode) {
473  case c_Direct:
474  hep = &hep0;
475  break;
476  case c_GeneratorLevel:
477  hep = &gen_level(hep0);
478  break;
479  }
480  // step 2: Gen_hepevt -> MCParticle
481  if (genHepevtToMCParticle.count(hep->get_ID()) > 0) {
482  int matchedMCParticle = genHepevtToMCParticle[hep->get_ID()];
483 
484  // step 3: set the relation
485  tracksToMCParticles.add(track->getArrayIndex(), matchedMCParticle);
486 
487  testMCRelation(*hep, m_mcParticles[matchedMCParticle], "Track");
488  } else {
489  B2DEBUG(99, "Can not find MCParticle corresponding to this gen_hepevt (Panther ID = " << hep->get_ID() << ")");
490  B2DEBUG(99, "Gen_hepevt: Panther ID = " << hep->get_ID() << "; idhep = " << hep->idhep() << "; isthep = " << hep->isthep());
491  }
492  }
493 }
494 
496 {
497  // Create and initialize K_S0 particle list
498  StoreObjPtr<ParticleList> ksPList("K_S0:mdst");
499  ksPList.create();
500  ksPList->initialize(310, ksPList.getName());
501 
502  // Create and initialize Lambda0 and anti-Lamda0 particle list
503  StoreObjPtr<ParticleList> lambda0PList("Lambda0:mdst");
504  lambda0PList.create();
505  lambda0PList->initialize(3122, lambda0PList.getName());
506 
507  StoreObjPtr<ParticleList> antiLambda0PList("anti-Lambda0:mdst");
508  antiLambda0PList.create();
509  antiLambda0PList->initialize(-3122, antiLambda0PList.getName());
510 
511  antiLambda0PList->bindAntiParticleList(*lambda0PList);
512 
513  // Create and initialize converted gamma particle list
514  StoreObjPtr<ParticleList> convGammaPList("gamma:v0mdst");
515  convGammaPList.create();
516  convGammaPList->initialize(22, convGammaPList.getName());
517 
518  // Loop over all Belle Vee2 candidates
519  Belle::Mdst_vee2_Manager& m = Belle::Mdst_vee2_Manager::get_manager();
520  for (Belle::Mdst_vee2_Manager::iterator vee2Iterator = m.begin(); vee2Iterator != m.end(); ++vee2Iterator) {
521  Belle::Mdst_vee2 belleV0 = *vee2Iterator;
522 
523  // +ve track
524  Belle::Mdst_charged belleTrackP = belleV0.chgd(0);
525  // -ve track
526  Belle::Mdst_charged belleTrackM = belleV0.chgd(1);
527 
528  // type of V0
531  int belleHypP = -1;
532  int belleHypM = -1;
533 
534  switch (belleV0.kind()) {
535  case 1 : // K0s -> pi+ pi-
536  pTypeP = Const::pion;
537  pTypeM = Const::pion;
538  belleHypP = 2;
539  belleHypM = 2;
540  break;
541  case 2 : // Lambda -> p+ pi-
542  pTypeP = Const::proton;
543  pTypeM = Const::pion;
544  belleHypP = 4;
545  belleHypM = 2;
546  break;
547  case 3 : // anti-Lambda -> pi+ anti-p-
548  pTypeP = Const::pion;
549  pTypeM = Const::proton;
550  belleHypP = 2;
551  belleHypM = 4;
552  break;
553  case 4 : // gamma -> e+ e-
554  pTypeP = Const::electron;
555  pTypeM = Const::electron;
556  belleHypP = 0;
557  belleHypM = 0;
558  break;
559  default :
560  B2WARNING("Conversion of vee2 candidate of unknown kind! kind = " << belleV0.kind());
561  }
562 
563  // This part is copied from Relation.cc in BASF
564  int trackID[2] = {0, 0};
565  unsigned nTrack = 0;
566  Belle::Mdst_charged_Manager& charged_mag = Belle::Mdst_charged_Manager::get_manager();
567  for (std::vector<Belle::Mdst_charged>::iterator chgIterator = charged_mag.begin(); chgIterator != charged_mag.end();
568  ++chgIterator) {
569  if (belleV0.chgd(0).get_ID() >= 1 && trackID[0] == 0 && belleV0.chgd(0).get_ID() == chgIterator->get_ID()) {
570  trackID[0] = (int)(chgIterator->get_ID()); //+ve trac
571  ++nTrack;
572  }
573  if (belleV0.chgd(1).get_ID() >= 1 && trackID[1] == 0 && belleV0.chgd(1).get_ID() == chgIterator->get_ID()) {
574  trackID[1] = (int)(chgIterator->get_ID()); //-ve trac
575  ++nTrack;
576  }
577  if (nTrack == 2)
578  break;
579  }
580 
581  HepPoint3D dauPivot(belleV0.vx(), belleV0.vy(), belleV0.vz());
582  int trackFitPIndex = -1;
583  int trackFitMIndex = -1;
584  Particle daughterP, daughterM;
585  CLHEP::HepLorentzVector momentumP;
586  CLHEP::HepSymMatrix error7x7P(7, 0);
587  HepPoint3D positionP;
588  TMatrixFSym errMatrixP(7);
589  CLHEP::HepLorentzVector momentumM;
590  CLHEP::HepSymMatrix error7x7M(7, 0);
591  HepPoint3D positionM;
592  TMatrixFSym errMatrixM(7);
593  CLHEP::HepSymMatrix error5x5(5, 0);
594  if (trackID[0] >= 1) {
595  if (belleV0.daut()) {
596  std::vector<float> helixParam(5);
597  std::vector<float> helixError(15);
598  belleVeeDaughterHelix(belleV0, 1, helixParam, helixError);
599 
600  auto trackFitP = m_trackFitResults.appendNew(helixParam, helixError, pTypeP, 0.5, -1, -1, 0);
601  trackFitPIndex = trackFitP->getArrayIndex();
602 
603  belleVeeDaughterToCartesian(belleV0, 1, pTypeP, momentumP, positionP, error7x7P);
604  TrackFitResult* tmpTFR = new TrackFitResult(createTrackFitResult(momentumP, positionP, error7x7P, 1, pTypeP, 0.5, -1, -1, 0));
605  // TrackFitResult internaly stores helix parameters at pivot = (0,0,0) so the momentum of the Particle will be wrong again.
606  // Overwrite it.
607 
608  for (unsigned i = 0; i < 7; i++)
609  for (unsigned j = 0; j < 7; j++)
610  errMatrixP(i, j) = error7x7P[i][j];
611 
612  daughterP = Particle(trackID[0] - 1, tmpTFR, pTypeP);
613  daughterP.updateMomentum(ROOT::Math::PxPyPzEVector(momentumP.px(), momentumP.py(), momentumP.pz(), momentumP.e()),
614  ROOT::Math::XYZVector(positionP.x(), positionP.y(), positionP.z()),
615  errMatrixP, 0.5);
616  delete tmpTFR;
617  } else {
618  Belle::Mdst_trk_fit& trk_fit = charged_mag[trackID[0] - 1].trk().mhyp(belleHypP);
619  double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
620 
621  std::vector<float> helixParam(5);
622  std::vector<float> helixError(15);
623  convertHelix(trk_fit, HepPoint3D(0., 0., 0.), helixParam, helixError);
624 
625  // Checking for invalid helix curvature with parameter 2 equal to 0:
626  if (helixParam[2] == 0) {
627  B2WARNING("Helix parameter for curvature == 0. Skipping Track! The parameter is: " << helixParam[2] << "...");
628  continue;
629  }
630 
631  auto trackFitP = m_trackFitResults.appendNew(helixParam, helixError, pTypeP, pValue, -1, -1, 0);
632 
633  trackFitPIndex = trackFitP->getArrayIndex();
634 
635  daughterP = Particle(trackID[0] - 1, trackFitP, pTypeP);
636  // set momentum/positions at pivot = V0 decay vertex
637  getHelixParameters(trk_fit, pTypeP.getMass(), dauPivot,
638  helixParam, error5x5,
639  momentumP, positionP, error7x7P);
640 
641  for (unsigned i = 0; i < 7; i++)
642  for (unsigned j = 0; j < 7; j++)
643  errMatrixP(i, j) = error7x7P[i][j];
644 
645  daughterP.updateMomentum(ROOT::Math::PxPyPzEVector(momentumP.px(), momentumP.py(), momentumP.pz(), momentumP.e()),
646  ROOT::Math::XYZVector(positionP.x(), positionP.y(), positionP.z()),
647  errMatrixP, pValue);
648  }
649  }
650  if (trackID[1] >= 1) {
651  if (belleV0.daut()) {
652  std::vector<float> helixParam(5);
653  std::vector<float> helixError(15);
654  belleVeeDaughterHelix(belleV0, -1, helixParam, helixError);
655 
656  auto trackFitM = m_trackFitResults.appendNew(helixParam, helixError, pTypeM, 0.5, -1, -1, 0);
657  trackFitMIndex = trackFitM->getArrayIndex();
658 
659  belleVeeDaughterToCartesian(belleV0, -1, pTypeM, momentumM, positionM, error7x7M);
660  TrackFitResult* tmpTFR = new TrackFitResult(createTrackFitResult(momentumM, positionM, error7x7M, -1, pTypeM, 0.5, -1, -1, 0));
661  // TrackFitResult internaly stores helix parameters at pivot = (0,0,0) so the momentum of the Particle will be wrong again.
662  // Overwrite it.
663  for (unsigned i = 0; i < 7; i++)
664  for (unsigned j = 0; j < 7; j++)
665  errMatrixM(i, j) = error7x7M[i][j];
666 
667  daughterM = Particle(trackID[1] - 1, tmpTFR, pTypeM);
668  daughterM.updateMomentum(ROOT::Math::PxPyPzEVector(momentumM.px(), momentumM.py(), momentumM.pz(), momentumM.e()),
669  ROOT::Math::XYZVector(positionM.x(), positionM.y(), positionM.z()),
670  errMatrixM, 0.5);
671  delete tmpTFR;
672  } else {
673  Belle::Mdst_trk_fit& trk_fit = charged_mag[trackID[1] - 1].trk().mhyp(belleHypM);
674  double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
675 
676  std::vector<float> helixParam(5);
677  std::vector<float> helixError(15);
678  convertHelix(trk_fit, HepPoint3D(0., 0., 0.), helixParam, helixError);
679 
680  // Checking for invalid helix curvature with parameter 2 equal to 0:
681  if (helixParam[2] == 0) {
682  B2WARNING("Helix parameter for curvature == 0. Skipping Track! The parameter is: " << helixParam[2] << "...");
683  continue;
684  }
685 
686  auto trackFitM = m_trackFitResults.appendNew(helixParam, helixError, pTypeM, pValue, -1, -1, 0);
687 
688  trackFitMIndex = trackFitM->getArrayIndex();
689 
690  daughterM = Particle(trackID[1] - 1, trackFitM, pTypeM);
691  // set momentum/positions at pivot = V0 decay vertex
692  getHelixParameters(trk_fit, pTypeM.getMass(), dauPivot,
693  helixParam, error5x5,
694  momentumM, positionM, error7x7M);
695 
696  for (unsigned i = 0; i < 7; i++)
697  for (unsigned j = 0; j < 7; j++)
698  errMatrixM(i, j) = error7x7M[i][j];
699 
700  daughterM.updateMomentum(ROOT::Math::PxPyPzEVector(momentumM.px(), momentumM.py(), momentumM.pz(), momentumM.e()),
701  ROOT::Math::XYZVector(positionM.x(), positionM.y(), positionM.z()),
702  errMatrixM, pValue);
703  }
704  }
705 
706  Track* trackP = m_tracks[trackID[0] - 1];
707  Track* trackM = m_tracks[trackID[1] - 1];
708 
709  TrackFitResult* trackFitP = m_trackFitResults[trackFitPIndex];
710  TrackFitResult* trackFitM = m_trackFitResults[trackFitMIndex];
711 
712  m_v0s.appendNew(std::make_pair(trackP, trackFitP), std::make_pair(trackM, trackFitM));
713 
714  // create Ks Particle and add it to the 'K_S0:mdst' ParticleList
715  const PIDLikelihood* pidP = trackP->getRelated<PIDLikelihood>();
716  const PIDLikelihood* pidM = trackM->getRelated<PIDLikelihood>();
717  const MCParticle* mcParticleP = trackP->getRelated<MCParticle>();
718  const MCParticle* mcParticleM = trackM->getRelated<MCParticle>();
719 
720  Particle* newDaugP = m_particles.appendNew(daughterP);
721  if (pidP)
722  newDaugP->addRelationTo(pidP);
723  if (mcParticleP)
724  newDaugP->addRelationTo(mcParticleP);
725  Particle* newDaugM = m_particles.appendNew(daughterM);
726  if (pidM)
727  newDaugM->addRelationTo(pidM);
728  if (mcParticleM)
729  newDaugM->addRelationTo(mcParticleM);
730 
731  ROOT::Math::PxPyPzEVector v0Momentum(belleV0.px(), belleV0.py(), belleV0.pz(), belleV0.energy());
732  ROOT::Math::XYZVector v0Vertex(belleV0.vx(), belleV0.vy(), belleV0.vz());
733 
734  /*
735  * Documentation of Mdst_vee2 vertex fit:
736  * /sw/belle/belle/b20090127_0910/share/tables/mdst.tdf (L96-L125)
737  */
738  auto appendVertexFitInfo = [](Belle::Mdst_vee2 & _belle_V0, Particle & _belle2_V0) {
739  // Add chisq of vertex fit. chiSquared=10^10 means the fit fails.
740  _belle2_V0.addExtraInfo("chiSquared", _belle_V0.chisq());
741  // Ndf of the vertex kinematic fit is 1
742  _belle2_V0.addExtraInfo("ndf", 1);
743  // Add p-value to extra Info
744  double prob = TMath::Prob(_belle_V0.chisq(), 1);
745  _belle2_V0.setPValue(prob);
746  };
747 
748  Particle* newV0 = nullptr;
749  if (belleV0.kind() == 1) { // K0s -> pi+ pi-
750  Particle KS(v0Momentum, 310);
751  KS.appendDaughter(newDaugP);
752  KS.appendDaughter(newDaugM);
753  KS.setVertex(v0Vertex);
754  appendVertexFitInfo(belleV0, KS);
755  newV0 = m_particles.appendNew(KS);
756  ksPList->addParticle(newV0);
757 
758  // append extra info: goodKs flag
759  Belle::FindKs belleKSFinder;
760  belleKSFinder.candidates(belleV0, Belle::IpProfile::position(1));
761  newV0->addExtraInfo("goodKs", belleKSFinder.goodKs());
762 
763  /*
764  std::cout << " ---- B1 Ks ---- " << std::endl;
765  std::cout << " momentum = " << std::endl;
766  v0Momentum.Print();
767  std::cout << " position = " << std::endl;
768  v0Vertex.Print();
769  std::cout << " ---- B2 Ks ---- " << std::endl;
770  std::cout << " momentum = " << std::endl;
771  newKS->get4Vector().Print();
772  std::cout << " position = " << std::endl;
773  newKS->getVertex().Print();
774  std::cout << " ---- B1 Ks.child(0) ---- " << std::endl;
775  std::cout << " momentum = " << momentumP << std::endl;
776  std::cout << " position = " << positionP << std::endl;
777  std::cout << " error7x7 = " << error7x7P << std::endl;
778  std::cout << " ---- B2 Ks.child(0) ---- " << std::endl;
779  std::cout << " momentum = " << std::endl;
780  newKS->getDaughter(0)->get4Vector().Print();
781  std::cout << " position = " << std::endl;
782  newKS->getDaughter(0)->getVertex().Print();
783  std::cout << " error7x7 = " << std::endl;
784  newKS->getDaughter(0)->getMomentumVertexErrorMatrix().Print();
785  std::cout << " ---- B1 Ks.child(1) ---- " << std::endl;
786  std::cout << " momentum = " << momentumM << std::endl;
787  std::cout << " position = " << positionM << std::endl;
788  std::cout << " error7x7 = " << error7x7M << std::endl;
789  std::cout << " ---- B2 Ks.child(1) ---- " << std::endl;
790  std::cout << " momentum = " << std::endl;
791  newKS->getDaughter(1)->get4Vector().Print();
792  std::cout << " position = " << std::endl;
793  newKS->getDaughter(1)->getVertex().Print();
794  std::cout << " error7x7 = " << std::endl;
795  newKS->getDaughter(1)->getMomentumVertexErrorMatrix().Print();
796  */
797  } else if (belleV0.kind() == 2) { // Lambda -> p+ pi-
798  Particle Lambda0(v0Momentum, 3122);
799  Lambda0.appendDaughter(newDaugP);
800  Lambda0.appendDaughter(newDaugM);
801  Lambda0.setVertex(v0Vertex);
802  appendVertexFitInfo(belleV0, Lambda0);
803  newV0 = m_particles.appendNew(Lambda0);
804  lambda0PList->addParticle(newV0);
805 
806  // GoodLambda flag as extra info
807  Belle::FindLambda lambdaFinder;
808  lambdaFinder.candidates(belleV0, Belle::IpProfile::position(1));
809  newV0->addExtraInfo("goodLambda", lambdaFinder.goodLambda());
810  } else if (belleV0.kind() == 3) { // anti-Lambda -> pi+ anti-p
811  Particle antiLambda0(v0Momentum, -3122);
812  antiLambda0.appendDaughter(newDaugM);
813  antiLambda0.appendDaughter(newDaugP);
814  antiLambda0.setVertex(v0Vertex);
815  appendVertexFitInfo(belleV0, antiLambda0);
816  newV0 = m_particles.appendNew(antiLambda0);
817  antiLambda0PList->addParticle(newV0);
818 
819  // GoodLambda flag as extra info
820  Belle::FindLambda lambdaFinder;
821  lambdaFinder.candidates(belleV0, Belle::IpProfile::position(1));
822  newV0->addExtraInfo("goodLambda", lambdaFinder.goodLambda());
823  } else if (belleV0.kind() == 4) { // gamma -> e+ e-
824  Particle gamma(v0Momentum, 22);
825  gamma.appendDaughter(newDaugP);
826  gamma.appendDaughter(newDaugM);
827  gamma.setVertex(v0Vertex);
828  appendVertexFitInfo(belleV0, gamma);
829  newV0 = m_particles.appendNew(gamma);
830  convGammaPList->addParticle(newV0);
831  }
832  // append extra info: nisKsFinder quality indicators
833  if (m_nisEnable) {
834  if (belleV0.kind() <= 3) { // K_S0, Lambda, anti-Lambda
835  Belle::nisKsFinder ksnb;
836  double protIDP = atcPID(pidP, 2, 4);
837  double protIDM = atcPID(pidM, 2, 4);
838  ksnb.candidates(belleV0, Belle::IpProfile::position(1), momentumP, protIDP, protIDM);
839  // K_S0 and Lambda (inverse cut on ksnbNoLam for Lambda selection).
840  newV0->addExtraInfo("ksnbVLike", ksnb.nb_vlike());
841  newV0->addExtraInfo("ksnbNoLam", ksnb.nb_nolam());
842  // K_S0 only
843  if (belleV0.kind() == 1)
844  newV0->addExtraInfo("ksnbStandard", ksnb.standard());
845  }
846  }
847  }
848 }
849 
851 {
852  if (m_realData)
853  return;
854 
855  // clear the Gen_hepevt_ID <-> MCParticleGraphPosition map
856  genHepevtToMCParticle.clear();
857 
858  // check if the Gen_hepevt table has any entries
859  Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
860  if (genMgr.count() == 0)
861  return;
862 
863  typedef std::pair<MCParticleGraph::GraphParticle*, Belle::Gen_hepevt> halfFamily;
864  halfFamily currFamily;
865  halfFamily family;
866  std::queue < halfFamily > heritancesQueue;
867 
868  // Add motherless particles. The root particle is often the first one,
869  // but this is not correct in general case; thus, all particles are added,
870  // including the beam-background ones.
872  for (Belle::Gen_hepevt_Manager::iterator genIterator = genMgr.begin();
873  genIterator != genMgr.end(); ++genIterator) {
874  Belle::Gen_hepevt hep = *genIterator;
875  // Select particles without mother.
876  if (!(hep.moFirst() == 0 && hep.moLast() == 0))
877  continue;
878 
879  int position = m_particleGraph.size();
881  genHepevtToMCParticle[hep.get_ID()] = position;
882  MCParticleGraph::GraphParticle* graphParticle = &m_particleGraph[position];
883  convertGenHepevtObject(hep, graphParticle);
884  for (int iDaughter = hep.daFirst(); iDaughter <= hep.daLast();
885  ++iDaughter) {
886  if (iDaughter == 0) {
887  B2DEBUG(95, "Trying to access generated daughter with Panther ID == 0");
888  continue;
889  }
890  currFamily.first = graphParticle;
891  currFamily.second = genMgr(Belle::Panther_ID(iDaughter));
892  heritancesQueue.push(currFamily);
893  }
894  }
895 
896  //now we can go through the queue:
897  while (!heritancesQueue.empty()) {
898  currFamily = heritancesQueue.front(); //get the first entry from the queue
899  heritancesQueue.pop(); //remove the entry.
900 
901  MCParticleGraph::GraphParticle* currMother = currFamily.first;
902  Belle::Gen_hepevt& currDaughter = currFamily.second;
903 
904  // Skip particles with idhep = 0.
905  if (currDaughter.idhep() == 0)
906  continue;
907 
908  //putting the daughter in the graph:
909  int position = m_particleGraph.size();
911  genHepevtToMCParticle[currDaughter.get_ID()] = position;
912 
913  MCParticleGraph::GraphParticle* graphDaughter = &m_particleGraph[position];
914  convertGenHepevtObject(currDaughter, graphDaughter);
915 
916  //add relation between mother and daughter to graph:
917  currMother->decaysInto((*graphDaughter));
918 
919  int nGrandChildren = currDaughter.daLast() - currDaughter.daFirst() + 1;
920 
921  if (nGrandChildren > 0 && currDaughter.daFirst() != 0) {
922  for (int igrandchild = currDaughter.daFirst(); igrandchild <= currDaughter.daLast(); ++igrandchild) {
923  if (igrandchild == 0) {
924  B2DEBUG(95, "Trying to access generated daughter with Panther ID == 0");
925  continue;
926  }
927 
928  family.first = graphDaughter;
929  family.second = genMgr(Belle::Panther_ID(igrandchild));
930  heritancesQueue.push(family);
931  }
932  }
933  }
934 
936 }
937 
939 {
940  // Relations
941  RelationArray eclClustersToMCParticles(m_eclClusters, m_mcParticles);
942 
943  // Clear the mdstEcl <-> ECLCluster map
944  mdstEclToECLCluster.clear();
945 
946  // Loop over all Belle Mdst_ecl
947  Belle::Mdst_ecl_Manager& ecl_manager = Belle::Mdst_ecl_Manager::get_manager();
948  Belle::Mdst_ecl_aux_Manager& ecl_aux_manager = Belle::Mdst_ecl_aux_Manager::get_manager();
949 
950  for (Belle::Mdst_ecl_Manager::iterator eclIterator = ecl_manager.begin(); eclIterator != ecl_manager.end(); ++eclIterator) {
951 
952  // Pull Mdst_ecl from manager
953  Belle::Mdst_ecl mdstEcl = *eclIterator;
954  Belle::Mdst_ecl_aux mdstEclAux(ecl_aux_manager(mdstEcl.get_ID()));
955 
956  // Create Belle II ECLCluster
957  auto B2EclCluster = m_eclClusters.appendNew();
958 
959  // Convert Mdst_ecl -> ECLCluster and create map of indices
960  convertMdstECLObject(mdstEcl, mdstEclAux, B2EclCluster);
961  mdstEclToECLCluster[mdstEcl.get_ID()] = B2EclCluster->getArrayIndex();
962 
963  // set ConnectedRegionID and ClusterID to
964  // cluster's array index + 1 and 1, respectively
965  B2EclCluster->setConnectedRegionId(B2EclCluster->getArrayIndex() + 1);
966  B2EclCluster->setClusterId(1);
967 
968  if (m_realData)
969  continue;
970 
971  // Create ECLCluster -> MCParticle relation
972  // Step 1: MDST_ECL -> Gen_hepevt
973  const Belle::Gen_hepevt& hep0 = get_hepevt(mdstEcl);
974  if (hep0 == 0)
975  continue;
976  const Belle::Gen_hepevt* hep = nullptr;
977  switch (m_mcMatchingMode) {
978  case c_Direct:
979  hep = &hep0;
980  break;
981  case c_GeneratorLevel:
982  hep = &gen_level(hep0);
983  break;
984  }
985  // Step 2: Gen_hepevt -> MCParticle
986  if (genHepevtToMCParticle.count(hep->get_ID()) > 0) {
987  int matchedMCParticleID = genHepevtToMCParticle[hep->get_ID()];
988  // Step 3: set the relation
989  eclClustersToMCParticles.add(B2EclCluster->getArrayIndex(), matchedMCParticleID);
990  testMCRelation(*hep, m_mcParticles[matchedMCParticleID], "ECLCluster");
991  } else {
992  B2DEBUG(79, "Cannot find MCParticle corresponding to this gen_hepevt (Panther ID = " << hep->get_ID() << ")");
993  B2DEBUG(79, "Gen_hepevt: Panther ID = " << hep->get_ID() << "; idhep = " << hep->idhep() << "; isthep = " << hep->isthep());
994  }
995  }
996 }
997 
999 {
1000  // There was no MC matching in Belle for KLM Clusters
1001 
1002  // Clear the mdstKlm <-> KLMCluster map
1003  mdstKlmToKLMCluster.clear();
1004 
1005  // Loop over all Belle Mdst_klm_cluster
1006  Belle::Mdst_klm_cluster_Manager& klm_cluster_manager = Belle::Mdst_klm_cluster_Manager::get_manager();
1007 
1008  for (Belle::Mdst_klm_cluster_Manager::iterator klmC_Ite = klm_cluster_manager.begin(); klmC_Ite != klm_cluster_manager.end();
1009  ++klmC_Ite) {
1010 
1011  // Pull Mdst_ecl from manager
1012  Belle::Mdst_klm_cluster mdstKlm_cluster = *klmC_Ite;
1013 
1014  // Create Belle II ECLCluster
1015  auto B2KlmCluster = m_klmClusters.appendNew();
1016 
1017  // Convert Mdst_klm_cluster -> KLMCluster and create map of indices
1018  convertMdstKLMObject(mdstKlm_cluster, B2KlmCluster);
1019  mdstKlmToKLMCluster[mdstKlm_cluster.get_ID()] = B2KlmCluster->getArrayIndex();
1020 
1021  }
1022 }
1023 
1025 {
1026  // Relations
1027  RelationArray particlesToMCParticles(m_particles, m_mcParticles);
1028 
1029  // Clear the mdstGamma <-> Particle map
1030  mdstGammaToParticle.clear();
1031 
1032  // Create and initialize particle list
1033  StoreObjPtr<ParticleList> plist("gamma:mdst");
1034  plist.create();
1035  plist->initialize(22, "gamma:mdst");
1036 
1037  // Loop over all Belle Mdst_gamma
1038  Belle::Mdst_gamma_Manager& gamma_manager = Belle::Mdst_gamma_Manager::get_manager();
1039 
1040  for (Belle::Mdst_gamma_Manager::iterator gammaIterator = gamma_manager.begin(); gammaIterator != gamma_manager.end();
1041  ++gammaIterator) {
1042 
1043  // Pull Mdst_gamma from manager and Mdst_ecl from pointer to Mdst_ecl
1044  Belle::Mdst_gamma mdstGamma = *gammaIterator;
1045  Belle::Mdst_ecl mdstEcl = mdstGamma.ecl();
1046  if (!mdstEcl)
1047  continue;
1048 
1049  // Get ECLCluster from map
1050  ECLCluster* B2EclCluster = m_eclClusters[mdstEclToECLCluster[mdstEcl.get_ID()]];
1051  if (!B2EclCluster)
1052  continue;
1053 
1054  // Create Particle from ECLCluster, add to StoreArray, create gamma map entry
1055  Particle* B2Gamma = m_particles.appendNew(B2EclCluster);
1056  mdstGammaToParticle[mdstGamma.get_ID()] = B2Gamma->getArrayIndex();
1057 
1058  // Add particle to particle list
1059  plist->addParticle(B2Gamma);
1060 
1061  if (m_realData)
1062  continue;
1063 
1064  // Relation to MCParticle
1065  MCParticle* matchedMCParticle = B2EclCluster->getRelated<MCParticle>();
1066  if (matchedMCParticle)
1067  B2Gamma->addRelationTo(matchedMCParticle);
1068  }
1069 }
1070 
1072 {
1073  // Create and initialize particle list
1074  StoreObjPtr<ParticleList> plist("pi0:mdst");
1075  plist.create();
1076  plist->initialize(111, "pi0:mdst");
1077 
1078  // Loop over all Mdst_pi0
1079  Belle::Mdst_pi0_Manager& pi0_manager = Belle::Mdst_pi0_Manager::get_manager();
1080  for (Belle::Mdst_pi0_Manager::iterator pi0Iterator = pi0_manager.begin(); pi0Iterator != pi0_manager.end(); ++pi0Iterator) {
1081 
1082  // Pull Mdst_pi0 from manager and Mdst_gammas from pointers to Mdst_gammas
1083  Belle::Mdst_pi0 mdstPi0 = *pi0Iterator;
1084  Belle::Mdst_gamma mdstGamma1 = mdstPi0.gamma(0);
1085  Belle::Mdst_gamma mdstGamma2 = mdstPi0.gamma(1);
1086  if (!mdstGamma1 || !mdstGamma2)
1087  continue;
1088 
1089  ROOT::Math::PxPyPzEVector p4(mdstPi0.px(), mdstPi0.py(), mdstPi0.pz(), mdstPi0.energy());
1090 
1091  // Create Particle from LorentzVector and PDG code, add to StoreArray
1092  Particle* B2Pi0 = m_particles.appendNew(p4, 111);
1093 
1094  // Get Belle II photons from map
1095  Particle* B2Gamma1 = m_particles[mdstGammaToParticle[mdstGamma1.get_ID()]];
1096  Particle* B2Gamma2 = m_particles[mdstGammaToParticle[mdstGamma2.get_ID()]];
1097  if (!B2Gamma1 || !B2Gamma2)
1098  continue;
1099 
1100  // Append photons as pi0 daughters
1101  B2Pi0->appendDaughter(B2Gamma1);
1102  B2Pi0->appendDaughter(B2Gamma2);
1103 
1104  // Add chisq of mass-constrained Kfit
1105  B2Pi0->addExtraInfo("chiSquared", mdstPi0.chisq());
1106 
1107  // Ndf of a pi0 mass-constrained kinematic fit is 1
1108  B2Pi0->addExtraInfo("ndf", 1);
1109 
1110  // Add p-value to extra Info
1111  double prob = TMath::Prob(mdstPi0.chisq(), 1);
1112  B2Pi0->setPValue(prob);
1113 
1114  // Add particle to particle list
1115  plist->addParticle(B2Pi0);
1116  }
1117 }
1118 
1120 {
1121  // Relations
1122  RelationArray particlesToMCParticles(m_particles, m_mcParticles);
1123 
1124 
1125  // Create and initialize particle list
1126  StoreObjPtr<ParticleList> plist("K_L0:mdst");
1127  plist.create();
1128  plist->initialize(Const::Klong.getPDGCode(), "K_L0:mdst");
1129 
1130  Belle::Mdst_klong_Manager& klong_manager = Belle::Mdst_klong_Manager::get_manager();
1131  for (Belle::Mdst_klong_Manager::iterator klong_Ite = klong_manager.begin(); klong_Ite != klong_manager.end(); ++klong_Ite) {
1132 
1133  // Pull Mdst_klong from manager and Mdst_klm from pointer to Mdst_klm
1134  Belle::Mdst_klong mdstKlong = *klong_Ite;
1135  Belle::Mdst_klm_cluster mdstKlm = mdstKlong.klmc();
1136 
1137  if (!mdstKlm)
1138  continue;
1139 
1140 
1141  // Get KLMCluster from map
1142  KLMCluster* B2KlmCluster = m_klmClusters[mdstKlmToKLMCluster[mdstKlm.get_ID()]];
1143  if (!B2KlmCluster)
1144  continue;
1145 
1146  // Extract cluster position from Klong and save it in KLMCluster
1147  B2KlmCluster->setClusterPosition(mdstKlong.cos_x(), mdstKlong.cos_y(), mdstKlong.cos_z());
1148 
1149  // Create Particle from KLMCluster, add to StoreArray, create Klong map entry
1150  Particle* B2Klong = m_particles.appendNew(B2KlmCluster);
1151  mdstKlongToParticle[mdstKlong.get_ID()] = B2Klong->getArrayIndex();
1152 
1153  // Add particle to particle list
1154  plist->addParticle(B2Klong);
1155  }
1156 
1157  // (Vague) MC Matching
1158  // There was no MC matching for KLongs in Belle , but a hack:
1159  // Check if MC KLong and reconstructed KLong (only without ecl) are within 15 degree for phi and theta, we set a relation
1160  // for the best reconstructed KLong to the MC KLong.
1161  // Taken and adapted from http://belle.kek.jp/secured/wiki/doku.php?id=physics:ckm:kleff
1162 
1163  if (!m_realData) {
1164 
1165  Belle::Gen_hepevt_Manager& GenMgr = Belle::Gen_hepevt_Manager::get_manager();
1166  const double dang(15. / 180.*M_PI); // check reconstructed candidates within 15 degrees
1167 
1168  for (Belle::Gen_hepevt_Manager::iterator klong_hep_it = GenMgr.begin(); klong_hep_it != GenMgr.end(); ++klong_hep_it) {
1169 
1170  if (abs((*klong_hep_it).idhep()) == Const::Klong.getPDGCode() && klong_hep_it->isthep() > 0) {
1171 
1172  CLHEP::HepLorentzVector gp4(klong_hep_it->PX(), klong_hep_it->PY(), klong_hep_it->PZ(), klong_hep_it->E());
1173  double sum(0.0);
1174  int bestRecKlongID(0);
1175 
1176  for (Belle::Mdst_klong_Manager::iterator klong_rec_it = klong_manager.begin(); klong_rec_it != klong_manager.end();
1177  ++klong_rec_it) {
1178 
1179  // if((*klong_rec_it).klmc().ecl())continue; // check only klm cand.
1180  if ((*klong_rec_it).ecl())
1181  continue; // check only klm cand.
1182  CLHEP::Hep3Vector klp3(klong_rec_it->cos_x(), klong_rec_it->cos_y(), klong_rec_it->cos_z());
1183 
1184  if (cos(gp4.theta() - klp3.theta()) > cos(dang) && cos(gp4.phi() - klp3.phi()) > cos(dang)) {
1185 
1186  double tmp_sum = cos(gp4.theta() - klp3.theta()) + cos(gp4.phi() - klp3.phi());
1187  if (tmp_sum > sum) {
1188  bestRecKlongID = mdstKlongToParticle[(*klong_rec_it).get_ID()];
1189  sum = tmp_sum;
1190  }
1191  }
1192 
1193  }
1194  if (sum > 0.0) {
1195  int matchedMCParticleID = genHepevtToMCParticle[(*klong_hep_it).get_ID()];
1196  particlesToMCParticles.add(bestRecKlongID, matchedMCParticleID);
1197  testMCRelation((*klong_hep_it), m_mcParticles[matchedMCParticleID], "m_particles");
1198  }
1199  }
1200  }
1201  }
1202 }
1203 
1205 {
1206  // Create StoreObj if it is not valid
1207  if (not m_evtInfo.isValid()) {
1208  m_evtInfo.create();
1209  }
1210  // Pull Evtcls_flag(2) from manager
1211  Belle::Evtcls_flag_Manager& EvtFlagMgr = Belle::Evtcls_flag_Manager::get_manager();
1212  Belle::Evtcls_flag2_Manager& EvtFlag2Mgr = Belle::Evtcls_flag2_Manager::get_manager();
1213 
1214  // Pull Evtcls_hadronic_flag from manager
1215  Belle::Evtcls_hadronic_flag_Manager& EvtHadFlagMgr = Belle::Evtcls_hadronic_flag_Manager::get_manager();
1216 
1217  std::string name = "evtcls_flag";
1218  std::string name_had = "evtcls_hadronic_flag";
1219  // Only one entry in each event
1220  std::vector<Belle::Evtcls_flag>::iterator eflagIterator = EvtFlagMgr.begin();
1221  std::vector<Belle::Evtcls_flag2>::iterator eflag2Iterator = EvtFlag2Mgr.begin();
1222  std::vector<Belle::Evtcls_hadronic_flag>::iterator ehadflagIterator = EvtHadFlagMgr.begin();
1223 
1224  // Converting evtcls_flag(2)
1225  std::vector<int> flag(20);
1226  for (int index = 0; index < 20; ++index) {
1227  // flag(14, 16): not filled
1228  if (index == 14 || index == 16) continue;
1229  std::string iVar = name + std::to_string(index);
1230  // 0-9 corresponding to evtcls_flag.flag(0-9)
1231  if (index < 10) {
1232  m_evtInfo->addExtraInfo(iVar, (*eflagIterator).flag(index));
1233  } else {
1234  // 10-19 corresponding to evtcls_flag2.flag(0-9)
1235  m_evtInfo->addExtraInfo(iVar, (*eflag2Iterator).flag(index - 10));
1236  }
1237  B2DEBUG(99, "evtcls_flag(" << index << ") = " << m_evtInfo->getExtraInfo(iVar));
1238  }
1239 
1240  // Converting evtcls_hadronic_flag
1241  for (int index = 0; index < 6; ++index) {
1242  std::string iVar = name_had + std::to_string(index);
1243  m_evtInfo->addExtraInfo(iVar, (*ehadflagIterator).hadronic_flag(index));
1244  B2DEBUG(99, "evtcls_hadronic_flag(" << index << ") = " << m_evtInfo->getExtraInfo(iVar));
1245  }
1246 
1247 }
1248 
1250 {
1251 
1252  // Create StoreObj if it is not valid
1253  if (not m_evtInfo.isValid()) {
1254  m_evtInfo.create();
1255  }
1256  // Load event info
1258 
1259  if (event->getExperiment() <= 27) { // Check if it's SVD 1
1260  // Pull rectrg_summary from namager
1261  Belle::Rectrg_summary_Manager& RecTrgSummaryMgr = Belle::Rectrg_summary_Manager::get_manager();
1262  std::vector<Belle::Rectrg_summary>::iterator eflagIterator = RecTrgSummaryMgr.begin();
1263  std::string name_summary = "rectrg_summary_m_final";
1264 
1265  // Converting m_final(2) from rectrg_summary
1266  for (int index = 0; index < 2; ++index) {
1267  std::string iVar = name_summary + std::to_string(index);
1268  m_evtInfo->addExtraInfo(iVar, (*eflagIterator).final(index));
1269  B2DEBUG(99, "m_final(" << index << ") = " << m_evtInfo->getExtraInfo(iVar));
1270  }
1271  } else { // For SVD2
1272  // Pull rectrg_summary3 from manager
1273  Belle::Rectrg_summary3_Manager& RecTrgSummary3Mgr = Belle::Rectrg_summary3_Manager::get_manager();
1274 
1275  std::string name_summary3 = "rectrg_summary3_m_final";
1276  // Only one entry in each event
1277  std::vector<Belle::Rectrg_summary3>::iterator eflagIterator3 = RecTrgSummary3Mgr.begin();
1278 
1279  // Converting m_final(3)
1280  for (int index = 0; index < 3; ++index) {
1281  std::string iVar = name_summary3 + std::to_string(index);
1282  m_evtInfo->addExtraInfo(iVar, (*eflagIterator3).final(index));
1283  B2DEBUG(99, "m_final(" << index << ") = " << m_evtInfo->getExtraInfo(iVar));
1284  }
1285  }
1286 
1287 }
1288 
1289 
1290 //-----------------------------------------------------------------------------
1291 // CONVERT OBJECTS
1292 //-----------------------------------------------------------------------------
1293 
1294 #ifdef HAVE_KID_ACC
1295 double B2BIIConvertMdstModule::acc_pid(const Belle::Mdst_charged& chg, int idp)
1296 {
1297  static Belle::kid_acc acc_pdf(0);
1298  //static kid_acc acc_pdf(1);
1299 
1300  const double pmass[5] = { 0.00051099907, 0.105658389, 0.13956995, 0.493677, 0.93827231 };
1301 
1302  CLHEP::Hep3Vector mom(chg.px(), chg.py(), chg.pz());
1303  double cos_theta = mom.cosTheta();
1304  double pval = mom.mag();
1305 
1306  double npe = chg.acc().photo_electron();
1307  double beta = pval / sqrt(pval * pval + pmass[idp] * pmass[idp]);
1308  double pdfval = acc_pdf.npe2pdf(cos_theta, beta, npe);
1309 
1310  return pdfval;
1311 }
1312 
1313 // this is CDC_prob5
1314 double B2BIIConvertMdstModule::cdc_pid(const Belle::Mdst_charged& chg, int idp)
1315 {
1316  CLHEP::Hep3Vector mom(chg.px(), chg.py(), chg.pz());
1317  double pval = mom.mag();
1318 
1319  Belle::kid_cdc kidCdc(5);
1320  float factor0 = kidCdc.factor0();
1321  float factor1 = kidCdc.factor1(idp, pval);
1322 
1323  if (factor0 == 1.0 && factor1 == 1.0) return chg.trk().pid(idp);
1324  //
1325  double m = chg.trk().dEdx() / factor0;
1326  double e = chg.trk().dEdx_exp(idp) * factor1;
1327  double s = chg.trk().sigma_dEdx(idp);
1328  double val = 1. / sqrt(2.*M_PI) / s * exp(-0.5 * (m - e) * (m - e) / s / s);
1329 
1330  return val;
1331 }
1332 #endif
1333 
1334 void B2BIIConvertMdstModule::setLikelihoods(PIDLikelihood* pid, Const::EDetector det, double likelihoods[c_nHyp],
1335  bool discard_allzero)
1336 {
1337  if (discard_allzero) {
1338  const double max_l = *std::max_element(likelihoods, likelihoods + c_nHyp);
1339  if (max_l <= 0.0) {
1340  return; //likelihoods broken, ignore
1341  }
1342  }
1343 
1344  for (int i = 0; i < c_nHyp; i++) {
1345  float logl = log(likelihoods[i]);
1346  pid->setLogLikelihood(det, c_belleHyp_to_chargedStable[i], logl);
1347  }
1348  //copy proton likelihood to deuterons
1349  pid->setLogLikelihood(det, Const::deuteron, pid->getLogL(Const::proton, det));
1350 }
1351 
1352 void B2BIIConvertMdstModule::convertPIDData(const Belle::Mdst_charged& belleTrack, const Track* track)
1353 {
1354  PIDLikelihood* pid = m_pidLikelihoods.appendNew();
1355  track->addRelationTo(pid);
1356 
1357  //convert data handled by atc_pid: dE/dx (-> CDC), TOF (-> TOP), ACC ( -> ARICH)
1358  //this should result in the same likelihoods used when creating atc_pid(3, 1, 5, ..., ...)
1359  //and calling prob(const Mdst_charged & chg).
1360 
1361  double likelihoods[c_nHyp];
1362  double accL[c_nHyp];
1363  double tofL[c_nHyp];
1364  double cdcL[c_nHyp];
1365  for (int i = 0; i < c_nHyp; i++) {
1366  accL[i] = tofL[i] = cdcL[i] = 1.0; // cppcheck-suppress unreadVariable
1367  }
1368 #ifdef HAVE_KID_ACC
1369  //accq0 = 3, as implemented in acc_prob3()
1370  const auto& acc = belleTrack.acc();
1371  if (acc and acc.quality() == 0) {
1372  for (int i = 0; i < c_nHyp; i++)
1373  accL[i] = likelihoods[i] = acc_pid(belleTrack, i); // cppcheck-suppress unreadVariable
1374  setLikelihoods(pid, Const::ARICH, likelihoods, true);
1375  }
1376 #endif
1377 
1378  //tofq0 = 1, as implemented in tof_prob1()
1379  //uses p1 / (p1 + p2) to create probability, so this should map directly to likelihoods
1380  const Belle::Mdst_tof& tof = belleTrack.tof();
1381  if (tof and tof.quality() == 0) {
1382  for (int i = 0; i < c_nHyp; i++)
1383  tofL[i] = likelihoods[i] = tof.pid(i); // cppcheck-suppress unreadVariable
1384  setLikelihoods(pid, Const::TOP, likelihoods, true);
1385  }
1386 
1387  // cdcq0 = 5, as implemented in cdc_prob0() (which is used for all values of cdcq0!)
1388  //uses p1 / (p1 + p2) to create probability, so this should map directly to likelihoods
1389  // eID actually uses cdc_pid (cdc_prob5)
1390  const Belle::Mdst_trk& trk = belleTrack.trk();
1391  if (trk.dEdx() > 0) {
1392  for (int i = 0; i < c_nHyp; i++) {
1393  likelihoods[i] = trk.pid(i);
1394  cdcL[i] = cdc_pid(belleTrack, i); // cppcheck-suppress unreadVariable
1395  }
1396  setLikelihoods(pid, Const::CDC, likelihoods, true);
1397  }
1398 
1399 
1400  // eid
1401  // eid is combination of atc_pid and ecl related information
1402  // since atc_pid part is already converted above only the ECL part
1403  // is converted
1404  // ECL pdfs are available only for electrons and hadrons (assumed to be pions)
1405  // likelihoods for others are set to 0
1406 
1407 #ifdef HAVE_EID
1408  Belle::eid electronID(belleTrack);
1409  float eclID_e_pdf = electronID.pdf_e_ecl();
1410  float eclID_h_pdf = electronID.pdf_h_ecl();
1411  float atcID_e_pdf = electronID.atc_pid_pdf(true, accL, tofL, cdcL);
1412  float atcID_h_pdf = electronID.atc_pid_pdf(false, accL, tofL, cdcL);
1413 
1414  // eID
1415  float eclProb = eclID_e_pdf / (eclID_e_pdf + eclID_h_pdf);
1416  float atcProb = atcID_e_pdf / (atcID_e_pdf + atcID_h_pdf);
1417 
1418  if (atcProb > 0.999999) atcProb = 0.999999;
1419  // combine the two probabilities.
1420  double eidCombinedSig = eclProb * atcProb;
1421  double eidCombinedBkg = (1. - eclProb) * (1. - atcProb);
1422 
1423  likelihoods[0] = eidCombinedSig;
1424  likelihoods[1] = 0; // no muons
1425  likelihoods[2] = eidCombinedBkg;
1426  likelihoods[3] = 0; // no kaons
1427  likelihoods[4] = 0; // no protons
1428 
1429  setLikelihoods(pid, Const::ECL, likelihoods, true);
1430 
1431  //Hep3Vector mom(belleTrack.px(), belleTrack.py(), belleTrack.pz());
1432  //B2INFO(" p = " << mom.mag() << " le_ecl = " << electronID.le_ecl());
1433 #endif
1434 
1435  //muid
1436  //Note that though it says "_likelihood()" on the label, those are
1437  //actually likelihood ratios of the type L(hyp) / (L(mu) + L(pi) + L(K)),
1438  //which are set in the FixMdst module.
1439  int muid_trackid = belleTrack.muid_ID();
1440  if (muid_trackid) {
1441  //Using approach 2. from http://belle.kek.jp/secured/muid/usage_muid.html since
1442  //it's much simpler than what Muid_mdst does.
1443  Belle::Mdst_klm_mu_ex_Manager& ex_mgr = Belle::Mdst_klm_mu_ex_Manager::get_manager();
1444  Belle::Mdst_klm_mu_ex& ex = ex_mgr(Belle::Panther_ID(muid_trackid));
1445 
1446  //filter out tracks with insufficient #hits (equal to cut on Muid_mdst::Chi_2())
1447  if (ex.Chi_2() > 0) {
1448  likelihoods[0] = 0; //no electrons
1449  likelihoods[1] = ex.Muon_likelihood();
1450  likelihoods[2] = ex.Pion_likelihood();
1451  likelihoods[3] = ex.Kaon_likelihood();
1452  likelihoods[4] = 0; //no protons
1453  //Miss_likelihood should only be != 0 for tracks that do not pass the Chi_2 cut.
1454 
1455  // in some cases the ex.XYZ_likelihood() < 0; Set it to 0 in these cases.
1456  for (int i = 0; i < 5; i++)
1457  if (likelihoods[i] < 0)
1458  likelihoods[i] = 0;
1459 
1460  //note: discard_allzero = false since all likelihoods = 0 usually means that Junk_likelihood is 1
1461  // PIDLikelihood::getProbability(hyp) will correctly return 0 then.
1462  setLikelihoods(pid, Const::KLM, likelihoods);
1463 
1464  /*
1465  const double tolerance = 1e-7;
1466  if (fabs(pid->getProbability(Const::muon, nullptr, Const::KLM) - ex.Muon_likelihood()) > tolerance ||
1467  fabs(pid->getProbability(Const::pion, nullptr, Const::KLM) - ex.Pion_likelihood()) > tolerance ||
1468  fabs(pid->getProbability(Const::kaon, nullptr, Const::KLM) - ex.Kaon_likelihood()) > tolerance) {
1469 
1470  B2INFO("muons: " << pid->getProbability(Const::muon, nullptr, Const::KLM) << " " << ex.Muon_likelihood());
1471  B2INFO("pion: " << pid->getProbability(Const::pion, nullptr, Const::KLM) << " " << ex.Pion_likelihood());
1472  B2INFO("kaon: " << pid->getProbability(Const::kaon, nullptr, Const::KLM) << " " << ex.Kaon_likelihood());
1473  B2INFO("miss/junk: " << ex.Miss_likelihood() << " " << ex.Junk_likelihood());
1474  }
1475  */
1476  }
1477  }
1478 }
1479 
1480 int B2BIIConvertMdstModule::getHelixParameters(const Belle::Mdst_trk_fit& trk_fit,
1481  const double mass,
1482  const HepPoint3D& newPivot,
1483  std::vector<float>& helixParams,
1484  CLHEP::HepSymMatrix& error5x5,
1485  CLHEP::HepLorentzVector& momentum,
1486  HepPoint3D& position,
1487  CLHEP::HepSymMatrix& error7x7, const double dPhi)
1488 {
1489  const HepPoint3D pivot(trk_fit.pivot_x(),
1490  trk_fit.pivot_y(),
1491  trk_fit.pivot_z());
1492 
1493  CLHEP::HepVector a(5);
1494  a[0] = trk_fit.helix(0);
1495  a[1] = trk_fit.helix(1);
1496  a[2] = trk_fit.helix(2);
1497  a[3] = trk_fit.helix(3);
1498  a[4] = trk_fit.helix(4);
1499  CLHEP::HepSymMatrix Ea(5, 0);
1500  Ea[0][0] = trk_fit.error(0);
1501  Ea[1][0] = trk_fit.error(1);
1502  Ea[1][1] = trk_fit.error(2);
1503  Ea[2][0] = trk_fit.error(3);
1504  Ea[2][1] = trk_fit.error(4);
1505  Ea[2][2] = trk_fit.error(5);
1506  Ea[3][0] = trk_fit.error(6);
1507  Ea[3][1] = trk_fit.error(7);
1508  Ea[3][2] = trk_fit.error(8);
1509  Ea[3][3] = trk_fit.error(9);
1510  Ea[4][0] = trk_fit.error(10);
1511  Ea[4][1] = trk_fit.error(11);
1512  Ea[4][2] = trk_fit.error(12);
1513  Ea[4][3] = trk_fit.error(13);
1514  Ea[4][4] = trk_fit.error(14);
1515 
1516  Belle::Helix helix(pivot, a, Ea);
1517 
1518  int charge = 0;
1519  if (helix.kappa() > 0)
1520  charge = 1;
1521  else
1522  charge = -1;
1523 
1524  if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
1525  helix.pivot(newPivot);
1526  momentum = helix.momentum(dPhi, mass, position, error7x7);
1527  } else {
1528  if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
1529  helix.pivot(HepPoint3D(0., 0., 0.));
1530  momentum = helix.momentum(dPhi, mass, position, error7x7);
1531  } else {
1532  momentum = helix.momentum(dPhi, mass, position, error7x7);
1533  }
1534  }
1535 
1536  convertHelix(helix, helixParams, error5x5);
1537 
1538  return charge;
1539 }
1540 
1541 void B2BIIConvertMdstModule::convertHelix(const Belle::Mdst_trk_fit& trk_fit,
1542  const HepPoint3D& newPivot,
1543  std::vector<float>& helixParams, std::vector<float>& helixError)
1544 {
1545  const HepPoint3D pivot(trk_fit.pivot_x(),
1546  trk_fit.pivot_y(),
1547  trk_fit.pivot_z());
1548 
1549  CLHEP::HepVector a(5);
1550  a[0] = trk_fit.helix(0);
1551  a[1] = trk_fit.helix(1);
1552  a[2] = trk_fit.helix(2);
1553  a[3] = trk_fit.helix(3);
1554  a[4] = trk_fit.helix(4);
1555  CLHEP::HepSymMatrix Ea(5, 0);
1556  Ea[0][0] = trk_fit.error(0);
1557  Ea[1][0] = trk_fit.error(1);
1558  Ea[1][1] = trk_fit.error(2);
1559  Ea[2][0] = trk_fit.error(3);
1560  Ea[2][1] = trk_fit.error(4);
1561  Ea[2][2] = trk_fit.error(5);
1562  Ea[3][0] = trk_fit.error(6);
1563  Ea[3][1] = trk_fit.error(7);
1564  Ea[3][2] = trk_fit.error(8);
1565  Ea[3][3] = trk_fit.error(9);
1566  Ea[4][0] = trk_fit.error(10);
1567  Ea[4][1] = trk_fit.error(11);
1568  Ea[4][2] = trk_fit.error(12);
1569  Ea[4][3] = trk_fit.error(13);
1570  Ea[4][4] = trk_fit.error(14);
1571 
1572  Belle::Helix helix(pivot, a, Ea);
1573 
1574  if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
1575  helix.pivot(newPivot);
1576  } else {
1577  if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
1578  helix.pivot(HepPoint3D(0., 0., 0.));
1579  }
1580  }
1581 
1582  CLHEP::HepSymMatrix error5x5(5, 0);
1583  convertHelix(helix, helixParams, error5x5);
1584 
1585  unsigned int size = 5;
1586  unsigned int counter = 0;
1587  for (unsigned int i = 0; i < size; i++)
1588  for (unsigned int j = i; j < size; j++)
1589  helixError[counter++] = error5x5[i][j];
1590 }
1591 
1592 void B2BIIConvertMdstModule::convertHelix(Belle::Helix& helix, std::vector<float>& helixParams, CLHEP::HepSymMatrix& error5x5)
1593 {
1594  CLHEP::HepVector a(5);
1595  CLHEP::HepSymMatrix Ea(5, 0);
1596 
1597  a = helix.a();
1598  Ea = helix.Ea();
1599 
1600  // param 0: d_0 = d_rho
1601  helixParams[0] = a[0];
1602 
1603  // param 1: phi = phi_0 + pi/2
1604  helixParams[1] = adjustAngleRange(a[1] + TMath::Pi() / 2.0);
1605 
1606  // param 2: omega = Kappa * alpha = Kappa * B[Tesla] * speed_of_light[m/s] * 1e-11
1607  helixParams[2] = a[2] * KAPPA2OMEGA;
1608 
1609  // param 3: d_z = z0
1610  helixParams[3] = a[3];
1611 
1612  // param 4: tan(Lambda) = tanLambda
1613  helixParams[4] = a[4];
1614 
1615  unsigned int size = 5;
1616  for (unsigned int i = 0; i < size; i++) {
1617  for (unsigned int j = 0; j < size; j++) {
1618  error5x5[i][j] = Ea[i][j];
1619  if (i == 2)
1620  error5x5[i][j] *= KAPPA2OMEGA;
1621  if (j == 2)
1622  error5x5[i][j] *= KAPPA2OMEGA;
1623 
1624  if (std::isinf(error5x5[i][j])) {
1625  B2DEBUG(99, "Helix covariance matrix element found to be infinite. Setting value to DBL_MAX/2.0.");
1626  error5x5[i][j] = DBL_MAX / 2.0;
1627  }
1628  }
1629  }
1630 }
1631 
1632 void B2BIIConvertMdstModule::convertMdstChargedObject(const Belle::Mdst_charged& belleTrack, Track* track)
1633 {
1634  Belle::Mdst_trk& trk = belleTrack.trk();
1635 
1636  for (int mhyp = 0 ; mhyp < c_nHyp; ++mhyp) {
1638  double thisMass = pType.getMass();
1639 
1640  Belle::Mdst_trk_fit& trk_fit = trk.mhyp(mhyp);
1641 
1642  // Converted helix parameters
1643  std::vector<float> helixParam(5);
1644  // Converted 5x5 error matrix
1645  CLHEP::HepSymMatrix error5x5(5, 0);
1646  // 4-momentum
1647  CLHEP::HepLorentzVector momentum;
1648  // 7x7 (momentum, position) error matrix
1649  CLHEP::HepSymMatrix error7x7(7, 0);
1650  // position
1651  HepPoint3D position;
1652 
1653  getHelixParameters(trk_fit, thisMass, HepPoint3D(0., 0., 0.),
1654  helixParam, error5x5,
1655  momentum, position, error7x7, 0.0);
1656 
1657  std::vector<float> helixError(15);
1658  unsigned int size = 5;
1659  unsigned int counter = 0;
1660  for (unsigned int i = 0; i < size; i++)
1661  for (unsigned int j = i; j < size; j++)
1662  helixError[counter++] = error5x5[i][j];
1663 
1664  double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
1665 
1666  // Create an empty cdc hitpattern and set the number of total hits
1667  // use hits from 0: axial-wire, 1:stereo-wire, 2:cathode
1668  // the actual cdc hitpattern is not converted
1669 
1670  int cdcNHits = 0;
1671  for (unsigned int i = 0; i < 3; i++)
1672  cdcNHits += trk_fit.nhits(i);
1673 
1674  HitPatternCDC patternCdc;
1675  patternCdc.setNHits(cdcNHits);
1676 
1677  // conversion of track position in CDC layers
1678  if (m_convertTrkExtra) {
1679  double tof = 0;
1680  double path_length = 0;
1681  double tof_sigma = 0;
1682  short tof_qual = 0;
1683  int acc_ph = 0;
1684  short acc_qual = 0;
1685  double dedx = trk.dEdx();
1686  short dedx_qual = trk.quality_dedx();
1687 
1688  const Belle::Mdst_tof& tof_obj = belleTrack.tof();
1689  if (tof_obj) {
1690  tof = tof_obj.tof();
1691  path_length = tof_obj.path_length();
1692  tof_qual = tof_obj.quality();
1693  tof_sigma = tof_obj.sigma_tof();
1694  }
1695 
1696  const Belle::Mdst_acc& acc_obj = belleTrack.acc();
1697  if (acc_obj) {
1698  acc_ph = acc_obj.photo_electron();
1699  acc_qual = acc_obj.quality();
1700  }
1701 
1702 
1703  auto cdcExtraInfo = m_belleTrkExtra.appendNew(trk_fit.first_x(), trk_fit.first_y(), trk_fit.first_z(),
1704  trk_fit.last_x(), trk_fit.last_y(), trk_fit.last_z(),
1705  tof, path_length, tof_qual, tof_sigma,
1706  acc_ph, acc_qual, dedx, dedx_qual);
1707  track->addRelationTo(cdcExtraInfo);
1708  }
1709  // conversion of the SVD hit pattern
1710  int svdHitPattern = trk_fit.hit_svd();
1711  // use hits from 3: SVD-rphi, 4: SVD-z
1712  // int svdNHits = trk_fit.nhits(3) + trk_fit.nhits(4);
1713 
1714  std::bitset<32> svdBitSet(svdHitPattern);
1715 
1716  HitPatternVXD patternVxd;
1717 
1718  unsigned short svdLayers;
1719  // taken from: http://belle.kek.jp/group/indirectcp/cpfit/cpfit-festa/2004/talks/Apr.14/CPfesta-2005-Higuchi(3).pdf
1721  // mask for the rphi hits, first 6 (8) bits/ 2 bits per layer
1722  std::bitset<32> svdUMask(static_cast<std::string>("00000000000000000000000000000011"));
1723  // mask for the z hits, second 6 (8) bits/ 2 bits per layer
1724  std::bitset<32> svdVMask;
1725 
1726  // find out if the SVD has 3 (4) layers; if exp <= (>) exp 27
1727  if (event->getExperiment() <= 27) {
1728  svdVMask = svdUMask << 6;
1729  svdLayers = 3;
1730  } else {
1731  svdVMask = svdUMask << 8;
1732  svdLayers = 4;
1733  }
1734 
1735  // loop over all svd layers (layer index is shifted + 3 for basf2)
1736  for (unsigned short layerId = 0; layerId < svdLayers; layerId++) {
1737  unsigned short uHits = (svdBitSet & svdUMask).count();
1738  unsigned short vHits = (svdBitSet & svdVMask).count();
1739  patternVxd.setSVDLayer(layerId + 3, uHits, vHits);
1740  // shift masks to the left
1741  svdUMask <<= 2;
1742  svdVMask <<= 2;
1743  }
1744 
1745  TrackFitResult helixFromHelix(helixParam, helixError, pType, pValue, -1, patternVxd.getInteger(), 0);
1746 
1748  TMatrixDSym cartesianCovariance(6);
1749  for (unsigned i = 0; i < 7; i++) {
1750  if (i == 3)
1751  continue;
1752  for (unsigned j = 0; j < 7; j++) {
1753  if (j == 3)
1754  continue;
1755 
1756  cartesianCovariance(ERRMCONV[i], ERRMCONV[j]) = error7x7[i][j];
1757  }
1758  }
1759  UncertainHelix helixFromCartesian(helixFromHelix.getPosition(), helixFromHelix.getMomentum(), helixFromHelix.getChargeSign(),
1760  BFIELD, cartesianCovariance, pValue);
1761 
1762  TMatrixDSym helixCovariance = helixFromCartesian.getCovariance();
1763 
1764  counter = 0;
1765  for (unsigned int i = 0; i < 5; ++i)
1766  for (unsigned int j = i; j < 5; ++j)
1767  helixError[counter++] = helixCovariance(i, j);
1768  }
1769 
1770  auto trackFit = m_trackFitResults.appendNew(helixParam, helixError, pType, pValue, patternCdc.getInteger(),
1771  patternVxd.getInteger(), trk_fit.ndf());
1772  track->setTrackFitResultIndex(pType, trackFit->getArrayIndex());
1773  /*
1774  B2INFO("--- B1 Track: ");
1775  std::cout << "Momentum = " << momentum << std::endl;
1776  std::cout << "Position = " << position << std::endl;
1777  std::cout << "7x7 error matrix = " << error7x7 << std::endl;
1778  B2INFO("--- B2 Track: ");
1779  std::cout << "Momentum = " << std::endl;
1780  trackFit->get4Momentum().Print();
1781  std::cout << "Position = " << std::endl;
1782  trackFit->getPosition().Print();
1783  std::cout << "6x6 error matrix = " << std::endl;
1784  trackFit->getCovariance6().Print();
1785  TMatrixDSym b2Error7x7(7);
1786  fill7x7ErrorMatrix(trackFit, b2Error7x7, thisMass, 1.5);
1787  std::cout << "7x7 error matrix = " << std::endl;
1788  b2Error7x7.Print();
1789  */
1790  }
1791 }
1792 
1793 void B2BIIConvertMdstModule::convertGenHepevtObject(const Belle::Gen_hepevt& genHepevt, MCParticleGraph::GraphParticle* mcParticle)
1794 {
1795  //B2DEBUG(80, "Gen_ehepevt: idhep " << genHepevt.idhep() << " (" << genHepevt.isthep() << ") with ID = " << genHepevt.get_ID());
1796 
1797  // updating the GraphParticle information from the Gen_hepevt information
1798  const int idHep = recoverMoreThan24bitIDHEP(genHepevt.idhep());
1799  if (idHep == 0) {
1800  B2WARNING("Trying to convert Gen_hepevt with idhep = " << idHep <<
1801  ". This should never happen.");
1802  mcParticle->setPDG(Const::photon.getPDGCode());
1803  } else {
1804  mcParticle->setPDG(idHep);
1805  }
1806 
1807  if (genHepevt.isthep() > 0) {
1809  }
1810 
1811  mcParticle->setMass(genHepevt.M());
1812 
1813  ROOT::Math::PxPyPzEVector p4(genHepevt.PX(), genHepevt.PY(), genHepevt.PZ(), genHepevt.E());
1814  mcParticle->set4Vector(p4);
1815 
1816  mcParticle->setProductionVertex(genHepevt.VX()*Unit::mm, genHepevt.VY()*Unit::mm, genHepevt.VZ()*Unit::mm);
1817  mcParticle->setProductionTime(genHepevt.T()*Unit::mm / Const::speedOfLight);
1818 
1819  // decay time of this particle is production time of the daughter particle
1820  if (genHepevt.daFirst() > 0) {
1821  Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
1822  Belle::Gen_hepevt daughterParticle = genMgr(Belle::Panther_ID(genHepevt.daFirst()));
1823  mcParticle->setDecayTime(daughterParticle.T()*Unit::mm / Const::speedOfLight);
1824  mcParticle->setDecayVertex(daughterParticle.VX()*Unit::mm, daughterParticle.VY()*Unit::mm, daughterParticle.VZ()*Unit::mm);
1825  } else {
1826  //otherwise, assume it's stable
1827  mcParticle->setDecayTime(std::numeric_limits<float>::infinity());
1828  }
1829 
1830  mcParticle->setValidVertex(true);
1831 }
1832 
1833 void B2BIIConvertMdstModule::convertMdstECLObject(const Belle::Mdst_ecl& ecl, const Belle::Mdst_ecl_aux& eclAux,
1834  ECLCluster* eclCluster)
1835 {
1836  if (eclAux.e9oe25() < m_matchType2E9oE25Threshold)
1837  eclCluster->setIsTrack(ecl.match() > 0);
1838  else
1839  eclCluster->setIsTrack(ecl.match() == 1);
1840 
1841  eclCluster->setEnergy(ecl.energy()); //must happen before setCovarianceMatrix()!
1842  eclCluster->setPhi(ecl.phi());
1843  eclCluster->setTheta(ecl.theta());
1844  eclCluster->setR(ecl.r());
1845  eclCluster->setdeltaL(ecl.quality());
1846 
1847  double covarianceMatrix[6];
1848  covarianceMatrix[0] = ecl.error(0); // error on energy
1849  covarianceMatrix[1] = ecl.error(1);
1850  covarianceMatrix[2] = ecl.error(2); // error on phi
1851  covarianceMatrix[3] = ecl.error(3);
1852  covarianceMatrix[4] = ecl.error(4);
1853  covarianceMatrix[5] = ecl.error(5); // error on theta
1854  eclCluster->setCovarianceMatrix(covarianceMatrix);
1855 
1856  eclCluster->setLAT(eclAux.width());
1857  eclCluster->setEnergyRaw(eclAux.mass());
1858  eclCluster->setE9oE21(eclAux.e9oe25());
1859  eclCluster->setEnergyHighestCrystal(eclAux.seed());
1860  eclCluster->setTime(eclAux.property(0));
1861  eclCluster->setNumberOfCrystals(eclAux.nhits());
1862 }
1863 
1864 void B2BIIConvertMdstModule::convertMdstKLMObject(const Belle::Mdst_klm_cluster& klm_cluster, KLMCluster* klmCluster)
1865 {
1866  // note: Belle quality flag is not saved (no free int variable in Belle2 KLMCluster)
1867  klmCluster->setLayers(klm_cluster.layers());
1868  klmCluster->setInnermostLayer(klm_cluster.first_layer());
1869 }
1870 
1871 
1872 //-----------------------------------------------------------------------------
1873 // RELATIONS
1874 //-----------------------------------------------------------------------------
1876 {
1877  // Relations
1878  RelationArray tracksToECLClusters(m_tracks, m_eclClusters);
1879 
1880  Belle::Mdst_ecl_trk_Manager& m = Belle::Mdst_ecl_trk_Manager::get_manager();
1881  Belle::Mdst_charged_Manager& chgMg = Belle::Mdst_charged_Manager::get_manager();
1882 
1883  // We first insert relations to tracks which are directly matched (type == 1)
1884  // secondly we had CR matched tracks (connected region) (type == 2)
1885  // finally tracks which are geometrically matched (type == 0)
1886  std::vector<int> insert_order_types = {1, 2, 0};
1887  for (auto& insert_type : insert_order_types) {
1888  for (Belle::Mdst_ecl_trk_Manager::iterator ecltrkIterator = m.begin(); ecltrkIterator != m.end(); ++ecltrkIterator) {
1889  Belle::Mdst_ecl_trk mECLTRK = *ecltrkIterator;
1890 
1891  if (mECLTRK.type() != insert_type)
1892  continue;
1893 
1894  Belle::Mdst_ecl mdstEcl = mECLTRK.ecl();
1895  Belle::Mdst_trk mTRK = mECLTRK.trk();
1896 
1897  if (!mdstEcl)
1898  continue;
1899 
1900  // the numbering in mdst_charged
1901  // not necessarily the same as in mdst_trk
1902  // therfore have to find corresponding mdst_charged
1903  for (Belle::Mdst_charged_Manager::iterator chgIterator = chgMg.begin(); chgIterator != chgMg.end(); ++chgIterator) {
1904  Belle::Mdst_charged mChar = *chgIterator;
1905  Belle::Mdst_trk mTRK_in_charged = mChar.trk();
1906 
1907  if (mTRK_in_charged.get_ID() == mTRK.get_ID()) {
1908  // found the correct mdst_charged
1909  tracksToECLClusters.add(mChar.get_ID() - 1, mdstEcl.get_ID() - 1, 1.0);
1910  break;
1911  }
1912  }
1913  }
1914  }
1915 }
1916 
1917 
1919 {
1920  // Relations
1921  RelationArray klmClustersToTracks(m_klmClusters, m_tracks);
1922  RelationArray klmClustersToEclClusters(m_klmClusters, m_eclClusters);
1923 
1924  Belle::Mdst_klm_cluster_Manager& klm_cluster_manager = Belle::Mdst_klm_cluster_Manager::get_manager();
1925 
1926 
1927  for (Belle::Mdst_klm_cluster_Manager::iterator klmC_Ite = klm_cluster_manager.begin(); klmC_Ite != klm_cluster_manager.end();
1928  ++klmC_Ite) {
1929 
1930  Belle::Mdst_klm_cluster mdstKlm_cluster = *klmC_Ite;
1931  Belle::Mdst_trk mTRK = mdstKlm_cluster.trk();
1932  Belle::Mdst_ecl mECL = mdstKlm_cluster.ecl();
1933 
1934  if (mTRK) klmClustersToTracks.add(mdstKlm_cluster.get_ID() - 1, mTRK.get_ID() - 1);
1935  if (mECL) klmClustersToEclClusters.add(mdstKlm_cluster.get_ID() - 1, mECL.get_ID() - 1);
1936  }
1937 }
1938 
1939 
1940 //-----------------------------------------------------------------------------
1941 // MISC
1942 //-----------------------------------------------------------------------------
1943 
1945 {
1946  /*
1947  QUICK CHECK: most of the normal particles are smaller than
1948  0x100000, while all the corrupt id has some of the high bits on.
1949 
1950  This bit check has to be revised when the table below is updated.
1951  */
1952  const int mask = 0x00f00000;
1953  int high_bits = id & mask;
1954  if (high_bits == 0 || high_bits == mask) return id;
1955 
1956  switch (id) {
1957  case 7114363:
1958  return 91000443; // X(3940)
1959  case 6114363:
1960  return 90000443; // Y(3940)
1961  case 6114241:
1962  return 90000321; // K_0*(800)+
1963  case 6114231:
1964  return 90000311; // K_0*(800)0
1965  case -6865004:
1966  return 9912212; // p_diff+
1967  case -6865104:
1968  return 9912112; // n_diffr
1969  case -6866773:
1970  return 9910443; // psi_diff
1971  case -6866883:
1972  return 9910333; // phi_diff
1973  case -6866993:
1974  return 9910223; // omega_diff
1975  case -6867005:
1976  return 9910211; // pi_diff+
1977  case -6867103:
1978  return 9910113; // rho_diff0
1979  case -7746995:
1980  return 9030221; // f_0(1500)
1981  case -7756773:
1982  return 9020443; // psi(4415)
1983  case -7756995:
1984  return 9020221; // eta(1405)
1985  case -7766773:
1986  return 9010443; // psi(4160)
1987  case -7776663:
1988  return 9000553; // Upsilon(5S)
1989  case -7776773:
1990  return 9000443; // psi(4040)
1991  case -7776783:
1992  return 9000433; // D_sj(2700)+
1993  case -7776995:
1994  return 9000221; // f_0(600)
1995  case -6114241:
1996  return -90000321; // K_0*(800)-
1997  case -6114231:
1998  return -90000311; // anti-K_0*(800)0
1999  case 6865004:
2000  return -9912212; // anti-p_diff-
2001  case 6865104:
2002  return -9912112; // anti-n_diffr
2003  case 6867005:
2004  return -9910211; // pi_diff-
2005  case 7776783:
2006  return -9000433; // D_sj(2700)-
2007  default:
2008  return id;
2009  }
2010 }
2011 
2012 void B2BIIConvertMdstModule::testMCRelation(const Belle::Gen_hepevt& belleMC, const MCParticle* mcP, const std::string& objectName)
2013 {
2014  int bellePDGCode = belleMC.idhep();
2015  int belleIIPDGCode = mcP->getPDG();
2016 
2017  if (bellePDGCode == 0)
2018  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to Gen_hepevt with idhep = 0.");
2019 
2020  if (bellePDGCode != belleIIPDGCode)
2021  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to different MCParticle! " << bellePDGCode << " vs. " <<
2022  belleIIPDGCode);
2023 
2024  const double belleMomentum[] = { belleMC.PX(), belleMC.PY(), belleMC.PZ() };
2025  const double belle2Momentum[] = { mcP->get4Vector().Px(), mcP->get4Vector().Py(), mcP->get4Vector().Pz() };
2026 
2027  for (unsigned i = 0; i < 3; i++) {
2028  double relDev = (belle2Momentum[i] - belleMomentum[i]) / belleMomentum[i];
2029 
2030  if (relDev > 1e-3) {
2031  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to different MCParticle!");
2032  B2INFO(" - Gen_hepevt [" << bellePDGCode << "] px/py/pz = " << belleMC.PX() << "/" << belleMC.PY() << "/" << belleMC.PZ());
2033  B2INFO(" - TrackFitResult [" << belleIIPDGCode << "] px/py/pz = " << mcP->get4Vector().Px() << "/" << mcP->get4Vector().Py() << "/"
2034  << mcP->get4Vector().Pz());
2035  }
2036  }
2037 }
2038 
2039 void B2BIIConvertMdstModule::belleVeeDaughterToCartesian(const Belle::Mdst_vee2& vee, const int charge,
2040  const Const::ParticleType& pType,
2041  CLHEP::HepLorentzVector& momentum, HepPoint3D& position, CLHEP::HepSymMatrix& error)
2042 {
2043  const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2044  CLHEP::HepVector a(5);
2045  CLHEP::HepSymMatrix Ea(5, 0);
2046  if (charge > 0) {
2047  a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2048  a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2049  a[4] = vee.daut().helix_p(4);
2050  Ea[0][0] = vee.daut().error_p(0); Ea[1][0] = vee.daut().error_p(1);
2051  Ea[1][1] = vee.daut().error_p(2); Ea[2][0] = vee.daut().error_p(3);
2052  Ea[2][1] = vee.daut().error_p(4); Ea[2][2] = vee.daut().error_p(5);
2053  Ea[3][0] = vee.daut().error_p(6); Ea[3][1] = vee.daut().error_p(7);
2054  Ea[3][2] = vee.daut().error_p(8); Ea[3][3] = vee.daut().error_p(9);
2055  Ea[4][0] = vee.daut().error_p(10); Ea[4][1] = vee.daut().error_p(11);
2056  Ea[4][2] = vee.daut().error_p(12); Ea[4][3] = vee.daut().error_p(13);
2057  Ea[4][4] = vee.daut().error_p(14);
2058  } else {
2059  a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2060  a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2061  a[4] = vee.daut().helix_m(4);
2062  Ea[0][0] = vee.daut().error_m(0); Ea[1][0] = vee.daut().error_m(1);
2063  Ea[1][1] = vee.daut().error_m(2); Ea[2][0] = vee.daut().error_m(3);
2064  Ea[2][1] = vee.daut().error_m(4); Ea[2][2] = vee.daut().error_m(5);
2065  Ea[3][0] = vee.daut().error_m(6); Ea[3][1] = vee.daut().error_m(7);
2066  Ea[3][2] = vee.daut().error_m(8); Ea[3][3] = vee.daut().error_m(9);
2067  Ea[4][0] = vee.daut().error_m(10); Ea[4][1] = vee.daut().error_m(11);
2068  Ea[4][2] = vee.daut().error_m(12); Ea[4][3] = vee.daut().error_m(13);
2069  Ea[4][4] = vee.daut().error_m(14);
2070  }
2071 
2072  Belle::Helix helix(pivot, a, Ea);
2073 
2074  // this is Vee daughter momentum/position/error at pivot = V0 Decay Vertex
2075  momentum = helix.momentum(0., pType.getMass(), position, error);
2076 }
2077 
2078 void B2BIIConvertMdstModule::belleVeeDaughterHelix(const Belle::Mdst_vee2& vee, const int charge, std::vector<float>& helixParam,
2079  std::vector<float>& helixError)
2080 {
2081  const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2082  CLHEP::HepVector a(5);
2083  CLHEP::HepSymMatrix Ea(5, 0);
2084  if (charge > 0) {
2085  a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2086  a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2087  a[4] = vee.daut().helix_p(4);
2088  Ea[0][0] = vee.daut().error_p(0);
2089  Ea[1][0] = vee.daut().error_p(1);
2090  Ea[1][1] = vee.daut().error_p(2);
2091  Ea[2][0] = vee.daut().error_p(3);
2092  Ea[2][1] = vee.daut().error_p(4);
2093  Ea[2][2] = vee.daut().error_p(5);
2094  Ea[3][0] = vee.daut().error_p(6);
2095  Ea[3][1] = vee.daut().error_p(7);
2096  Ea[3][2] = vee.daut().error_p(8);
2097  Ea[3][3] = vee.daut().error_p(9);
2098  Ea[4][0] = vee.daut().error_p(10);
2099  Ea[4][1] = vee.daut().error_p(11);
2100  Ea[4][2] = vee.daut().error_p(12);
2101  Ea[4][3] = vee.daut().error_p(13);
2102  Ea[4][4] = vee.daut().error_p(14);
2103  } else {
2104  a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2105  a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2106  a[4] = vee.daut().helix_m(4);
2107  Ea[0][0] = vee.daut().error_m(0);
2108  Ea[1][0] = vee.daut().error_m(1);
2109  Ea[1][1] = vee.daut().error_m(2);
2110  Ea[2][0] = vee.daut().error_m(3);
2111  Ea[2][1] = vee.daut().error_m(4);
2112  Ea[2][2] = vee.daut().error_m(5);
2113  Ea[3][0] = vee.daut().error_m(6);
2114  Ea[3][1] = vee.daut().error_m(7);
2115  Ea[3][2] = vee.daut().error_m(8);
2116  Ea[3][3] = vee.daut().error_m(9);
2117  Ea[4][0] = vee.daut().error_m(10);
2118  Ea[4][1] = vee.daut().error_m(11);
2119  Ea[4][2] = vee.daut().error_m(12);
2120  Ea[4][3] = vee.daut().error_m(13);
2121  Ea[4][4] = vee.daut().error_m(14);
2122  }
2123 
2124  Belle::Helix helix(pivot, a, Ea);
2125 
2126  // go to the new pivot
2127  helix.pivot(HepPoint3D(0., 0., 0.));
2128 
2129  CLHEP::HepSymMatrix error5x5(5, 0);
2130  convertHelix(helix, helixParam, error5x5);
2131 
2132  unsigned int size = 5;
2133  unsigned int counter = 0;
2134  for (unsigned int i = 0; i < size; i++)
2135  for (unsigned int j = i; j < size; j++)
2136  helixError[counter++] = error5x5[i][j];
2137 }
2138 
2139 TrackFitResult B2BIIConvertMdstModule::createTrackFitResult(const CLHEP::HepLorentzVector& momentum,
2140  const HepPoint3D& position,
2141  const CLHEP::HepSymMatrix& error,
2142  const short int charge,
2143  const Const::ParticleType& pType,
2144  const float pValue,
2145  const uint64_t hitPatternCDCInitializer,
2146  const uint32_t hitPatternVXDInitializer,
2147  const uint16_t ndf)
2148 {
2149  ROOT::Math::XYZVector pos(position.x(), position.y(), position.z());
2150  ROOT::Math::XYZVector mom(momentum.px(), momentum.py(), momentum.pz());
2151 
2152  TMatrixDSym errMatrix(6);
2153  for (unsigned i = 0; i < 7; i++) {
2154  if (i == 3)
2155  continue;
2156  for (unsigned j = 0; j < 7; j++) {
2157  if (j == 3)
2158  continue;
2159 
2160  if (i == j)
2161  errMatrix(ERRMCONV[i], ERRMCONV[i]) = error[i][i];
2162  else
2163  errMatrix(ERRMCONV[i], ERRMCONV[j]) = errMatrix(ERRMCONV[j], ERRMCONV[i]) = error[i][j];
2164  }
2165  }
2166 
2167  return TrackFitResult(pos, mom, errMatrix, charge, pType, pValue, BFIELD, hitPatternCDCInitializer, hitPatternVXDInitializer, ndf);
2168 }
2169 
2170 double B2BIIConvertMdstModule::atcPID(const PIDLikelihood* pid, int sigHyp, int bkgHyp)
2171 {
2172  if (!pid) return 0.5;
2173 
2174  // ACC = ARICH
2175  Const::PIDDetectorSet set = Const::ARICH;
2176  double acc_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2177  double acc_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2178  double acc = 0.5;
2179  if (acc_sig + acc_bkg > 0.0)
2180  acc = acc_sig / (acc_sig + acc_bkg);
2181 
2182  // TOF = TOP
2183  set = Const::TOP;
2184  double tof_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2185  double tof_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2186  double tof = 0.5;
2187  double tof_all = tof_sig + tof_bkg;
2188  if (tof_all != 0) {
2189  tof = tof_sig / tof_all;
2190  if (tof < 0.001) tof = 0.001;
2191  if (tof > 0.999) tof = 0.999;
2192  }
2193 
2194  // dE/dx = CDC
2195  set = Const::CDC;
2196  double cdc_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2197  double cdc_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2198  double cdc = 0.5;
2199  double cdc_all = cdc_sig + cdc_bkg;
2200  if (cdc_all != 0) {
2201  cdc = cdc_sig / cdc_all;
2202  if (cdc < 0.001) cdc = 0.001;
2203  if (cdc > 0.999) cdc = 0.999;
2204  }
2205 
2206  // Combined
2207  double pid_sig = acc * tof * cdc;
2208  double pid_bkg = (1. - acc) * (1. - tof) * (1. - cdc);
2209 
2210  return pid_sig / (pid_sig + pid_bkg);
2211 }
2212 
2213 
2215 {
2216  B2DEBUG(99, "B2BIIConvertMdst: endRun done.");
2217 }
2218 
2219 
2221 {
2222  B2DEBUG(99, "B2BIIConvertMdst: terminate called");
2223 }
2224 
void belleVeeDaughterHelix(const Belle::Mdst_vee2 &vee, const int charge, std::vector< float > &helixParam, std::vector< float > &helixError)
obtains the helix parameters of the vee daughters
double atcPID(const PIDLikelihood *pid, int sigHyp, int bkgHyp)
calculates atc_pid(3,1,5,sigHyp,bkgHyp).prob() from converted PIDLikelihood
void convertMdstGammaTable()
Reads all entries of Mdst_Gamma Panther table, creates a particle list 'gamma:mdst' and adds them to ...
std::map< int, int > mdstKlmToKLMCluster
map of Mdst_klm Panther IDs and corresponding KLMCluster StoreArray indices
void convertMdstChargedTable()
Reads and converts all entries of Mdst_charged (Mdst_trk and Mdst_trk_fit) Panther table to Track (Tr...
void setTracksToECLClustersRelations()
Sets Track -> ECLCluster relations.
OptionalDBObjPtr< CollisionBoostVector > m_collisionBoostVectorDB
CollisionBoostVector for boost vector.
StoreObjPtr< EventExtraInfo > m_evtInfo
Event Extra Info.
const int ERRMCONV[7]
CONVERSION OF TRACK ERROR MATRIX ELEMENTS.
StoreArray< KLMCluster > m_klmClusters
KLM clusters.
void belleVeeDaughterToCartesian(const Belle::Mdst_vee2 &vee, const int charge, const Const::ParticleType &pType, CLHEP::HepLorentzVector &momentum, HepPoint3D &position, CLHEP::HepSymMatrix &error)
Fills 4-momentum, position and 7x7 error matrix from Belle Vee daughter.
@ c_GeneratorLevel
Match to generator-level particles.
virtual void initialize() override
Initialize the module.
void setLikelihoods(PIDLikelihood *pid, Const::EDetector det, double likelihoods[c_nHyp], bool discard_allzero=false)
Add given Belle likelihoods (not log-likelihoods, in Belle hypothesis order) for given detector to pi...
StoreArray< V0 > m_v0s
V0-particles.
virtual void event() override
Called for each event.
void convertGenHepEvtTable()
Reads and converts all entries of Gen_hepevt Panther table to MCParticle dataobjects and adds them to...
double cdc_pid(const Belle::Mdst_charged &chg, int idp)
Returns CDC likelihood for given hypothesis idp.
void convertRecTrgTable()
Reads and converts m_final from rectrg_summary3.
void setKLMClustersRelations()
Sets KLMCluster -> Track and ECLCluster relations.
bool m_convertEvtcls
Flag to switch on conversion of Evtcls table.
const double KAPPA2OMEGA
Conversion factor for Kappa -> Omega helix parameters.
virtual void endRun() override
Called when the current run is finished.
bool m_nisEnable
Flag to switch on conversion of nisKsFinder info.
StoreArray< TrackFitResult > m_trackFitResults
Track fir results.
BeamSpot m_beamSpot
Interaction Point of the beam.
StoreArray< Particle > m_particles
Particles.
virtual void terminate() override
Terminates the module.
OptionalDBObjPtr< CollisionInvariantMass > m_collisionInvMDB
CollisionInvariantMass for Invariant Mass of Beam.
void convertHelix(Belle::Helix &helix, std::vector< float > &helixParams, CLHEP::HepSymMatrix &error5x5)
Converts Belle's Helix parameters and it's covariance matrix to Belle II's version.
void convertMdstPi0Table()
Reads all entries of Mdst_Pi0 Panther table, creates a particle list 'pi0:mdst' and adds them to Stor...
std::map< int, int > genHepevtToMCParticle
map of Gen_hepevt Panther IDs and corresponding MCParticle StoreArray indices
void testMCRelation(const Belle::Gen_hepevt &belleMC, const MCParticle *mcP, const std::string &objectName)
Checks if the reconstructed object (Track, ECLCluster, ...) was matched to the same MCParticle.
void convertIPProfile(bool beginRun=false)
Stores the IPProfiles in BeamSpot (currently in DataStore)
void convertMdstECLObject(const Belle::Mdst_ecl &ecl, const Belle::Mdst_ecl_aux &eclAux, ECLCluster *eclCluster)
Converts Mdst_ecl(_aux) record to ECLCluster object.
CollisionInvariantMass m_collisionInvM
CollisionInvariantMass for the invariant mass of the beam.
bool m_convertBeamParameters
Convert beam parameters or use information stored in Belle II database.
static double acc_pid(const Belle::Mdst_charged &chg, int idp)
Returns ACC likelihood for given hypothesis idp.
virtual void beginRun() override
Module functions to be called from event process.
StoreArray< Track > m_tracks
Tracks.
void convertBeamEnergy()
Stores beam parameters (energy, angles) in CollisionInvariantMass and CollisionBoostVector (currently...
void convertPIDData(const Belle::Mdst_charged &belleTrack, const Track *track)
Get PID information for belleTrack and add it to PIDLikelihood (with relation from track).
void convertMdstKLMObject(const Belle::Mdst_klm_cluster &klm, KLMCluster *klmCluster)
Converts Mdst_klm_cluster record to KLMCluster object.
Belle2::MCParticleGraph m_particleGraph
MCParticle Graph to build Belle2 MC Particles.
void convertMdstKLongTable()
Reads all entries of Mdst_Klong Panther table, creates a particle list 'K_L0:mdst' and adds them to S...
OptionalDBObjPtr< BeamSpot > m_beamSpotDB
BeamSpot for IP.
bool m_realData
flag that tells whether given data sample is for real data or MC
StoreArray< ECLCluster > m_eclClusters
ECL clusters.
void convertGenHepevtObject(const Belle::Gen_hepevt &genHepevt, MCParticleGraph::GraphParticle *mcParticle)
Converts Gen_hepevt record to MCParticleGraph object.
void initializeDataStore()
Initializes Belle II DataStore.
void convertMdstECLTable()
Reads and converts all entries of Mdst_ecl(_aux) Panther table to ECLCluster dataobjects and adds the...
CollisionBoostVector m_collisionBoostVector
CollisionBoostVector for bosst vector of the beam.
std::string m_mcMatchingModeString
MC matching mode.
int m_lastIPProfileBin
variable to tell us which IPProfile bin was active last time we looked
void convertMdstKLMTable()
Reads and converts all entries of Mdst_klm_cluster Panther table to KLMCluster dataobjects and adds t...
bool m_convertRecTrg
Flag to switch on conversion of rectrg_summary3.
std::map< int, int > mdstGammaToParticle
map of gamma Panther IDs and corresponding Particle StoreArray indices
static const Const::ChargedStable c_belleHyp_to_chargedStable[c_nHyp]
maps Belle hypotheses to Const::ChargedStable (from http://belle.kek.jp/secured/wiki/doku....
std::map< int, int > mdstKlongToParticle
map of Klong Panther IDs and corresponding Particle StoreArray indices
double m_matchType2E9oE25Threshold
E9/E25 threshold value clusters with a value above this threshold are classified as neutral even if t...
void convertMdstChargedObject(const Belle::Mdst_charged &belleTrack, Track *track)
Converts Mdst_charged (Mdst_trk(_fit)) record to Track (TrackFitResult) object.
StoreArray< MCParticle > m_mcParticles
MC particles.
void convertMdstVee2Table()
Reads and converts all entries of Mdst_vee2 Panther table to V0 dataobjects and adds them to StoreArr...
StoreArray< PIDLikelihood > m_pidLikelihoods
output PIDLikelihood array.
std::map< int, int > mdstEclToECLCluster
map of Mdst_ecl Panther IDs and corresponding ECLCluster StoreArray indices
TrackFitResult createTrackFitResult(const CLHEP::HepLorentzVector &momentum, const HepPoint3D &position, const CLHEP::HepSymMatrix &error, const short int charge, const Const::ParticleType &pType, const float pValue, const uint64_t hitPatternCDCInitializer, const uint32_t hitPatternVXDInitializer, const uint16_t ndf)
Creates TrackFitResult and fills it.
virtual ~B2BIIConvertMdstModule() override
Destructor.
const double BFIELD
B filed in TESLA.
MCMatchingMode m_mcMatchingMode
C matching mode.
StoreArray< BelleTrkExtra > m_belleTrkExtra
Belle CDC extra information.
bool m_use6x6CovarianceMatrix4Tracks
flag that tells which form of covariance matrix should be used in the conversion of charged tracks
int recoverMoreThan24bitIDHEP(int id)
Helper function to recover falsely set idhep info in GenHepEvt list.
bool m_convertTrkExtra
Flag to switch on conversion of first(last)_{x,y,z} of mdst_trk_fit.
int getHelixParameters(const Belle::Mdst_trk_fit &trk_fit, const double mass, const HepPoint3D &newPivot, std::vector< float > &helixParams, CLHEP::HepSymMatrix &error5x5, CLHEP::HepLorentzVector &momentum, HepPoint3D &position, CLHEP::HepSymMatrix &error7x7, const double dPhi=0.0)
Fills Helix parameters (converted to Belle II version), 5x5 error matrix, 4-momentum,...
static const int c_nHyp
Number of Belle track hypotheses (see c_belleHyp_to_chargedStable).
void convertEvtclsTable()
Reads and converts all entries of evtcls Panther table.
This class contains the beam spot position and size modeled as a gaussian distribution in space.
Definition: BeamSpot.h:22
void setIP(const TVector3 &ipPosition, const TMatrixDSym &covariance)
Set the IP position and its error matrix.
Definition: BeamSpot.h:59
This class contains the measured average boost vector vec(beta) = (beta_x, beta_y,...
void setBoost(const TVector3 &boost, const TMatrixDSym &covariance)
Set the boost vector and its error matrix.
This class contains the measured average center-of-mass energy, which is equal to the invariant mass ...
void setMass(double mass, double error, double spread)
Set the CMS energy and its uncertainty.
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:580
The ParticleType class for identifying different particle types.
Definition: Const.h:399
int getPDGCode() const
PDG code.
Definition: Const.h:464
double getMass() const
Particle mass.
Definition: UnitConst.cc:356
A class for sets of detector IDs whose content is limited to restricted set of valid detector IDs.
Definition: Const.h:287
static const ChargedStable muon
muon particle
Definition: Const.h:651
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:42
static const ChargedStable pion
charged pion particle
Definition: Const.h:652
static const ParticleType Klong
K^0_L particle.
Definition: Const.h:668
static const double speedOfLight
[cm/ns]
Definition: Const.h:685
static const double electronMass
electron mass
Definition: Const.h:675
static const ChargedStable proton
proton particle
Definition: Const.h:654
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:653
static const ParticleType photon
photon particle
Definition: Const.h:664
static const ChargedStable electron
electron particle
Definition: Const.h:650
static const ChargedStable deuteron
deuteron particle
Definition: Const.h:655
ECL cluster data.
Definition: ECLCluster.h:27
void setE9oE21(double E9oE21)
Set E9/E21 energy ratio.
Definition: ECLCluster.h:193
void setTheta(double theta)
Set Theta of Shower (radian).
Definition: ECLCluster.h:217
void setPhi(double phi)
Set Phi of Shower (radian).
Definition: ECLCluster.h:220
void setdeltaL(double deltaL)
Set deltaL for shower shape.
Definition: ECLCluster.h:172
void setEnergyHighestCrystal(double energyhighestcrystal)
Set energy of highest energetic crystal (GeV).
Definition: ECLCluster.h:232
void setEnergyRaw(double energyraw)
Set Uncorrect Energy deposited (GeV).
Definition: ECLCluster.h:229
void setTime(double time)
Set time information.
Definition: ECLCluster.h:211
void setCovarianceMatrix(double covArray[6])
Set covariance matrix (3x3), i.e.
Definition: ECLCluster.h:152
void setLAT(double LAT)
Set Lateral distribution parameter.
Definition: ECLCluster.h:205
void setNumberOfCrystals(double noc)
Set number of crystals (sum of weights).
Definition: ECLCluster.h:208
void setEnergy(double energy)
Set Corrected Energy (GeV).
Definition: ECLCluster.h:226
void setIsTrack(bool istrack)
Set m_isTrack true if the cluster matches with a track.
Definition: ECLCluster.h:104
void setR(double r)
Set R (in cm).
Definition: ECLCluster.h:223
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Definition: Helix.cc:109
Hit pattern of CDC hits within a track.
Definition: HitPatternCDC.h:35
void setNHits(unsigned short nHits)
Sets the 8 MSBs to the total number of hits in the CDC.
ULong64_t getInteger() const
Getter for underlying integer type.
Hit pattern of the VXD within a track.
Definition: HitPatternVXD.h:37
unsigned int getInteger() const
Getter for the underlying integer.
Definition: HitPatternVXD.h:58
void setSVDLayer(const unsigned short layerId, unsigned short uHits, unsigned short vHits)
Set the number of hits in a specific layer of the SVD.
A class that describes the interval of experiments/runs for which an object in the database is valid.
KLM cluster data.
Definition: KLMCluster.h:28
void setLayers(int layers)
Set number of layers with hits.
Definition: KLMCluster.h:145
void setClusterPosition(float globalX, float globalY, float globalZ)
Set global position.
Definition: KLMCluster.h:161
void setInnermostLayer(int innermostLayer)
Set number of the innermost layer with hits.
Definition: KLMCluster.h:152
Class to represent Particle data in graph.
void decaysInto(GraphParticle &daughter)
Tells the graph that this particle decays into daughter.
size_t size() const
Return the number of particles in the graph.
void generateList(const std::string &name="", int options=c_setNothing)
Generates the MCParticle list and stores it in the StoreArray with the given name.
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:32
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:47
void setDecayTime(float time)
Set decay time.
Definition: MCParticle.h:390
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:366
void setDecayVertex(const ROOT::Math::XYZVector &vertex)
Set decay vertex.
Definition: MCParticle.h:447
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:378
void setProductionVertex(const ROOT::Math::XYZVector &vertex)
Set production vertex position.
Definition: MCParticle.h:396
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:207
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:335
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:112
void set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets the 4Vector of particle.
Definition: MCParticle.h:438
void setStatus(unsigned short int status)
Set Status code for the particle.
Definition: MCParticle.h:346
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:384
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 collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:26
Class to store reconstructed particles.
Definition: Particle.h:74
void appendDaughter(const Particle *daughter, const bool updateType=true, const int daughterProperty=c_Ordinary)
Appends index of daughter to daughters index array.
Definition: Particle.cc:708
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
Definition: Particle.h:310
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1363
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:353
void updateMomentum(const ROOT::Math::PxPyPzEVector &p4, const ROOT::Math::XYZVector &vertex, const TMatrixFSym &errMatrix, double pValue)
Sets Lorentz vector, position, 7x7 error matrix and p-value.
Definition: Particle.h:373
static bool parallelProcessingUsed()
Returns true if multiple processes have been spawned, false in single-core mode.
Definition: ProcHandler.cc:226
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:62
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
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.
const std::string & getName() const
Return name under which the object is saved in the DataStore.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
bool create(bool replace=false)
Create a default object in the data store.
Type-safe access to single objects in the data store.
Definition: StoreObjPtr.h:95
Values of the result of a track fit with a given particle hypothesis.
Helix getHelix() const
Conversion to framework Helix (without covariance).
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
short getChargeSign() const
Return track charge (1 or -1).
double getOmega() const
Getter for omega.
double getD0() const
Getter for d0.
double getTanLambda() const
Getter for tanLambda.
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.
double getPhi0() const
Getter for phi0.
Class that bundles various TrackFitResults.
Definition: Track.h:25
This class represents an ideal helix in perigee parameterization including the covariance matrix of t...
const TMatrixDSym & getCovariance() const
Getter for covariance matrix of perigee parameters in matrix form.
static const double mm
[millimeters]
Definition: Unit.h:70
REG_MODULE(B2BIIConvertBeamParams)
Register the module.
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
static Database & Instance()
Instance of a singleton Database.
Definition: Database.cc:42
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:28
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
Definition: Database.cc:141
void addConstantOverride(const std::string &name, TObject *obj, bool oneRun=false)
Add constant override payload.
Definition: DBStore.cc:196
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
void clear()
Reset particles and decay information to make the class reusable.
GraphParticle & addParticle()
Add new particle to the graph.
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:23