Belle II Software  release-05-02-19
B2BIIConvertMdstModule.cc
1 //+
2 // File : PantherInutModule.cc
3 // Description : A module to read panther records in basf2
4 //
5 // Author : Ryosuke Itoh, IPNS, KEK
6 // Date : 16 - Feb - 2015
7 //
8 // Contributors: Anze Zupanc, Matic Lubej
9 //-
10 
11 #include <b2bii/modules/B2BIIMdstInput/B2BIIConvertMdstModule.h>
12 
13 #include <framework/datastore/StoreObjPtr.h>
14 #include <framework/datastore/RelationArray.h>
15 #include <framework/database/Database.h>
16 #include <framework/pcore/ProcHandler.h>
17 
18 #include <mdst/dataobjects/HitPatternVXD.h>
19 #include <mdst/dataobjects/HitPatternCDC.h>
20 
21 // Belle II utilities
22 #include <framework/gearbox/Unit.h>
23 #include <framework/gearbox/Const.h>
24 #include <analysis/dataobjects/ParticleExtraInfoMap.h>
25 
26 // Belle II dataobjects
27 #include <framework/dataobjects/EventMetaData.h>
28 #include <framework/dataobjects/Helix.h>
29 #include <framework/dataobjects/UncertainHelix.h>
30 
31 // Belle utilities
32 #include <b2bii/utility/BelleMdstToGenHepevt.h>
33 
34 // ROOT
35 #include <TVector3.h>
36 #include <TLorentzVector.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 
75 const Const::ChargedStable B2BIIConvertMdstModule::c_belleHyp_to_chargedStable[c_nHyp] = { Const::electron, Const::muon, Const::pion, Const::kaon, Const::proton };
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 
160 B2BIIConvertMdstModule::B2BIIConvertMdstModule() : Module(),
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 
189 B2BIIConvertMdstModule::~B2BIIConvertMdstModule()
190 {
191 }
192 
193 void B2BIIConvertMdstModule::initialize()
194 {
195  // Initialize Belle II DataStore
196  initializeDataStore();
197  if (m_mcMatchingModeString == "Direct")
198  m_mcMatchingMode = c_Direct;
199  else if (m_mcMatchingModeString == "GeneratorLevel")
200  m_mcMatchingMode = c_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 
208 void B2BIIConvertMdstModule::initializeDataStore()
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 
220  StoreObjPtr<ParticleExtraInfoMap> extraInfoMap;
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 
262 void B2BIIConvertMdstModule::beginRun()
263 {
264  B2DEBUG(99, "B2BIIConvertMdst: beginRun called.");
265 
266  if (m_convertBeamParameters) {
267  //BeamEnergy class updated by fixmdst module in beginRun()
268  Belle::BeamEnergy::begin_run();
269  convertBeamEnergy();
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 
288 void B2BIIConvertMdstModule::event()
289 {
290  if (m_convertBeamParameters) {
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
301  convertIPProfile();
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) ||
306  !m_collisionBoostVectorDB || !m_collisionInvMDB) {
307  if ((!m_beamSpotDB || !m_collisionBoostVectorDB || !m_collisionInvMDB) && !m_realData) {
308  B2INFO("No database entry for this run yet, create one");
309  StoreObjPtr<EventMetaData> event;
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  }
322  if (ProcHandler::parallelProcessingUsed()) {
323  B2FATAL("Cannot reliably override the Database content in parallel processing "
324  "mode, please run the conversion in single processing mode");
325  }
326  DBStore::Instance().addConstantOverride("CollisionBoostVector", new CollisionBoostVector(m_collisionBoostVector), true);
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
333  convertGenHepEvtTable();
334 
335  // 2. Convert ECL information
336  convertMdstECLTable();
337 
338  // 3. Convert KLM information
339  convertMdstKLMTable();
340 
341  // 4. Convert Tracking information
342  convertMdstChargedTable();
343 
344  // 5. Set Track -> ECLCluster relations
345  setTracksToECLClustersRelations();
346 
347  // 6. Set KLMCluster -> Track, ECLCluster relations
348  setKLMClustersRelations();
349 
350  // 7. Convert Gamma information
351  convertMdstGammaTable();
352 
353  // 8. Convert Pi0 information
354  convertMdstPi0Table();
355 
356  // 9. Convert V0s
357  convertMdstVee2Table();
358 
359  // 10. Convert KLong information
360  convertMdstKLongTable();
361 
362  // 11. Convert Evtcls panther table information
363  if (m_convertEvtcls) convertEvtclsTable();
364 
365  // 12. Convert trigger information from rectrg_summary3
366  if (m_convertRecTrg) convertRecTrgTable();
367 
368 }
369 
370 
371 //-----------------------------------------------------------------------------
372 // CONVERT TABLES
373 //-----------------------------------------------------------------------------
374 void B2BIIConvertMdstModule::convertBeamEnergy()
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  TLorentzVector P_her(0.0, 0.0, TMath::Sqrt(Eher * Eher - mass_e * mass_e), Eher);
387  P_her.RotateY(angleHer);
388  TLorentzVector P_ler(0.0, 0.0, TMath::Sqrt(Eler * Eler - mass_e * mass_e), Eler);
389  P_ler.RotateY(angleLer);
390 
391  // Get four momentum of beam
392  TLorentzVector P_beam = P_her + P_ler;
393 
394  m_collisionBoostVector.setBoost(P_beam.BoostVector(), covariance);
395  m_collisionInvM.setMass(P_beam.M(), 0.0 , 0.0);
396 
397  B2DEBUG(99, "Beam Energy: E_HER = " << Eher << "; E_LER = " << Eler << "; angle = " << crossingAngle);
398  B2DEBUG(99, "Beam Momentum (pre-convert) : P_X = " << p_beam.px() << "; P_Y = " << p_beam.py() << "; P_Z = " << p_beam.pz());
399  B2DEBUG(99, "Beam Momentum (post-convert) : P_X = " << P_beam.Px() << "; P_Y = " << P_beam.Py() << "; P_Z = " << P_beam.Pz());
400 }
401 
402 void B2BIIConvertMdstModule::convertIPProfile(bool beginRun)
403 {
404  if (!Belle::IpProfile::usable()) {
405  // No IPProfile for this run ...
406  if (beginRun) {
407  // no IPProfile, set vertex to NaN without errors for the full run
408  m_beamSpot.setIP(
409  TVector3(std::numeric_limits<double>::quiet_NaN(),
410  std::numeric_limits<double>::quiet_NaN(),
411  std::numeric_limits<double>::quiet_NaN()
412  ), TMatrixTSym<double>()
413  );
414  }
415  return;
416  }
417  HepPoint3D ip;
418  CLHEP::HepSymMatrix ipErr;
419  if (beginRun) {
420  // use event independent average in begin run
421  ip = Belle::IpProfile::position();
422  ipErr = Belle::IpProfile::position_err();
423  } else {
424  // update evtbin
425  Belle::IpProfile::set_evtbin_number();
426  // check if it changed, if not there's nothing to do
427  if (Belle::IpProfile::EvtBinNo() == m_lastIPProfileBin) return;
428  // get event dependent position and error
429  ip = Belle::IpProfile::e_position();
430  ipErr = Belle::IpProfile::e_position_err();
431  }
432  // reset last ipprofile bin
433  m_lastIPProfileBin = Belle::IpProfile::EvtBinNo();
434 
435  TMatrixDSym cov(ipErr.num_col());
436  for (int i = 0; i < ipErr.num_row(); ++i) {
437  for (int j = 0; j < ipErr.num_col(); ++j) {
438  cov(i, j) = ipErr(i + 1, j + 1);
439  }
440  }
441  m_beamSpot.setIP(TVector3(ip.x(), ip.y(), ip.z()), cov);
442 }
443 
444 void B2BIIConvertMdstModule::convertMdstChargedTable()
445 {
446  // Relations
447  RelationArray tracksToMCParticles(m_tracks, m_mcParticles);
448 
449  // Loop over all Belle charged tracks
450  Belle::Mdst_charged_Manager& m = Belle::Mdst_charged_Manager::get_manager();
451  for (Belle::Mdst_charged_Manager::iterator chargedIterator = m.begin(); chargedIterator != m.end(); ++chargedIterator) {
452  Belle::Mdst_charged belleTrack = *chargedIterator;
453 
454  auto track = m_tracks.appendNew();
455 
456  // convert MDST_Charged -> Track
457  convertMdstChargedObject(belleTrack, track);
458 
459  convertPIDData(belleTrack, track);
460 
461  if (m_realData)
462  continue;
463 
464  // create Track -> MCParticle relation
465  // step 1: MDSTCharged -> Gen_hepevt
466  const Belle::Gen_hepevt& hep0 = get_hepevt(belleTrack);
467  if (hep0 == 0)
468  continue;
469  const Belle::Gen_hepevt* hep = nullptr;
470  switch (m_mcMatchingMode) {
471  case c_Direct:
472  hep = &hep0;
473  break;
474  case c_GeneratorLevel:
475  hep = &gen_level(hep0);
476  break;
477  }
478  // step 2: Gen_hepevt -> MCParticle
479  if (genHepevtToMCParticle.count(hep->get_ID()) > 0) {
480  int matchedMCParticle = genHepevtToMCParticle[hep->get_ID()];
481 
482  // step 3: set the relation
483  tracksToMCParticles.add(track->getArrayIndex(), matchedMCParticle);
484 
485  testMCRelation(*hep, m_mcParticles[matchedMCParticle], "Track");
486  } else {
487  B2DEBUG(99, "Can not find MCParticle corresponding to this gen_hepevt (Panther ID = " << hep->get_ID() << ")");
488  B2DEBUG(99, "Gen_hepevt: Panther ID = " << hep->get_ID() << "; idhep = " << hep->idhep() << "; isthep = " << hep->isthep());
489  }
490  }
491 }
492 
493 void B2BIIConvertMdstModule::convertMdstVee2Table()
494 {
495  // Create and initialize K_S0 particle list
496  StoreObjPtr<ParticleList> ksPList("K_S0:mdst");
497  ksPList.create();
498  ksPList->initialize(310, ksPList.getName());
499 
500  // Create and initialize Lambda0 and anti-Lamda0 particle list
501  StoreObjPtr<ParticleList> lambda0PList("Lambda0:mdst");
502  lambda0PList.create();
503  lambda0PList->initialize(3122, lambda0PList.getName());
504 
505  StoreObjPtr<ParticleList> antiLambda0PList("anti-Lambda0:mdst");
506  antiLambda0PList.create();
507  antiLambda0PList->initialize(-3122, antiLambda0PList.getName());
508 
509  antiLambda0PList->bindAntiParticleList(*lambda0PList);
510 
511  // Create and initialize converted gamma particle list
512  StoreObjPtr<ParticleList> convGammaPList("gamma:v0mdst");
513  convGammaPList.create();
514  convGammaPList->initialize(22, convGammaPList.getName());
515 
516  // Loop over all Belle Vee2 candidates
517  Belle::Mdst_vee2_Manager& m = Belle::Mdst_vee2_Manager::get_manager();
518  for (Belle::Mdst_vee2_Manager::iterator vee2Iterator = m.begin(); vee2Iterator != m.end(); ++vee2Iterator) {
519  Belle::Mdst_vee2 belleV0 = *vee2Iterator;
520 
521  // +ve track
522  Belle::Mdst_charged belleTrackP = belleV0.chgd(0);
523  // -ve track
524  Belle::Mdst_charged belleTrackM = belleV0.chgd(1);
525 
526  // type of V0
527  Const::ChargedStable pTypeP(Const::pion);
528  Const::ChargedStable pTypeM(Const::pion);
529  int belleHypP = -1;
530  int belleHypM = -1;
531 
532  switch (belleV0.kind()) {
533  case 1 : // K0s -> pi+ pi-
534  pTypeP = Const::pion;
535  pTypeM = Const::pion;
536  belleHypP = 2;
537  belleHypM = 2;
538  break;
539  case 2 : // Lambda -> p+ pi-
540  pTypeP = Const::proton;
541  pTypeM = Const::pion;
542  belleHypP = 4;
543  belleHypM = 2;
544  break;
545  case 3 : // anti-Lambda -> pi+ anti-p-
546  pTypeP = Const::pion;
547  pTypeM = Const::proton;
548  belleHypP = 2;
549  belleHypM = 4;
550  break;
551  case 4 : // gamma -> e+ e-
552  pTypeP = Const::electron;
553  pTypeM = Const::electron;
554  belleHypP = 0;
555  belleHypM = 0;
556  break;
557  default :
558  B2WARNING("Conversion of vee2 candidate of unknown kind! kind = " << belleV0.kind());
559  }
560 
561  // This part is copied from Relation.cc in BASF
562  int trackID[2] = {0, 0};
563  unsigned nTrack = 0;
564  Belle::Mdst_charged_Manager& charged_mag = Belle::Mdst_charged_Manager::get_manager();
565  for (std::vector<Belle::Mdst_charged>::iterator chgIterator = charged_mag.begin(); chgIterator != charged_mag.end();
566  ++chgIterator) {
567  if (belleV0.chgd(0).get_ID() >= 1 && trackID[0] == 0 && belleV0.chgd(0).get_ID() == chgIterator->get_ID()) {
568  trackID[0] = (int)(chgIterator->get_ID()); //+ve trac
569  ++nTrack;
570  }
571  if (belleV0.chgd(1).get_ID() >= 1 && trackID[1] == 0 && belleV0.chgd(1).get_ID() == chgIterator->get_ID()) {
572  trackID[1] = (int)(chgIterator->get_ID()); //-ve trac
573  ++nTrack;
574  }
575  if (nTrack == 2)
576  break;
577  }
578 
579  HepPoint3D dauPivot(belleV0.vx(), belleV0.vy(), belleV0.vz());
580  int trackFitPIndex = -1;
581  int trackFitMIndex = -1;
582  Particle daughterP, daughterM;
583  CLHEP::HepLorentzVector momentumP;
584  CLHEP::HepSymMatrix error7x7P(7, 0);
585  HepPoint3D positionP;
586  TMatrixFSym errMatrixP(7);
587  CLHEP::HepLorentzVector momentumM;
588  CLHEP::HepSymMatrix error7x7M(7, 0);
589  HepPoint3D positionM;
590  TMatrixFSym errMatrixM(7);
591  CLHEP::HepSymMatrix error5x5(5, 0);
592  if (trackID[0] >= 1) {
593  if (belleV0.daut()) {
594  std::vector<float> helixParam(5);
595  std::vector<float> helixError(15);
596  belleVeeDaughterHelix(belleV0, 1, helixParam, helixError);
597 
598  auto trackFitP = m_trackFitResults.appendNew(helixParam, helixError, pTypeP, 0.5, -1, -1, 0);
599  trackFitPIndex = trackFitP->getArrayIndex();
600 
601  belleVeeDaughterToCartesian(belleV0, 1, pTypeP, momentumP, positionP, error7x7P);
602  TrackFitResult* tmpTFR = new TrackFitResult(createTrackFitResult(momentumP, positionP, error7x7P, 1, pTypeP, 0.5, -1, -1, 0));
603  // TrackFitResult internaly stores helix parameters at pivot = (0,0,0) so the momentum of the Particle will be wrong again.
604  // Overwrite it.
605 
606  for (unsigned i = 0; i < 7; i++)
607  for (unsigned j = 0; j < 7; j++)
608  errMatrixP(i, j) = error7x7P[i][j];
609 
610  daughterP = Particle(trackID[0] - 1, tmpTFR, pTypeP);
611  daughterP.updateMomentum(TLorentzVector(momentumP.px(), momentumP.py(), momentumP.pz(), momentumP.e()),
612  TVector3(positionP.x(), positionP.y(), positionP.z()),
613  errMatrixP, 0.5);
614  delete tmpTFR;
615  } else {
616  Belle::Mdst_trk_fit& trk_fit = charged_mag[trackID[0] - 1].trk().mhyp(belleHypP);
617  double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
618 
619  std::vector<float> helixParam(5);
620  std::vector<float> helixError(15);
621  convertHelix(trk_fit, HepPoint3D(0., 0., 0.), helixParam, helixError);
622 
623  // Checking for invalid helix curvature with parameter 2 equal to 0:
624  if (helixParam[2] == 0) {
625  B2WARNING("Helix parameter for curvature == 0. Skipping Track! The parameter is: " << helixParam[2] << "...");
626  continue;
627  }
628 
629  auto trackFitP = m_trackFitResults.appendNew(helixParam, helixError, pTypeP, pValue, -1, -1, 0);
630 
631  trackFitPIndex = trackFitP->getArrayIndex();
632 
633  daughterP = Particle(trackID[0] - 1, trackFitP, pTypeP);
634  // set momentum/positions at pivot = V0 decay vertex
635  getHelixParameters(trk_fit, pTypeP.getMass(), dauPivot,
636  helixParam, error5x5,
637  momentumP, positionP, error7x7P);
638 
639  for (unsigned i = 0; i < 7; i++)
640  for (unsigned j = 0; j < 7; j++)
641  errMatrixP(i, j) = error7x7P[i][j];
642 
643  daughterP.updateMomentum(TLorentzVector(momentumP.px(), momentumP.py(), momentumP.pz(), momentumP.e()),
644  TVector3(positionP.x(), positionP.y(), positionP.z()),
645  errMatrixP, pValue);
646  }
647  }
648  if (trackID[1] >= 1) {
649  if (belleV0.daut()) {
650  std::vector<float> helixParam(5);
651  std::vector<float> helixError(15);
652  belleVeeDaughterHelix(belleV0, -1, helixParam, helixError);
653 
654  auto trackFitM = m_trackFitResults.appendNew(helixParam, helixError, pTypeM, 0.5, -1, -1, 0);
655  trackFitMIndex = trackFitM->getArrayIndex();
656 
657  belleVeeDaughterToCartesian(belleV0, -1, pTypeM, momentumM, positionM, error7x7M);
658  TrackFitResult* tmpTFR = new TrackFitResult(createTrackFitResult(momentumM, positionM, error7x7M, -1, pTypeM, 0.5, -1, -1, 0));
659  // TrackFitResult internaly stores helix parameters at pivot = (0,0,0) so the momentum of the Particle will be wrong again.
660  // Overwrite it.
661  for (unsigned i = 0; i < 7; i++)
662  for (unsigned j = 0; j < 7; j++)
663  errMatrixM(i, j) = error7x7M[i][j];
664 
665  daughterM = Particle(trackID[1] - 1, tmpTFR, pTypeM);
666  daughterM.updateMomentum(TLorentzVector(momentumM.px(), momentumM.py(), momentumM.pz(), momentumM.e()),
667  TVector3(positionM.x(), positionM.y(), positionM.z()),
668  errMatrixM, 0.5);
669  delete tmpTFR;
670  } else {
671  Belle::Mdst_trk_fit& trk_fit = charged_mag[trackID[1] - 1].trk().mhyp(belleHypM);
672  double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
673 
674  std::vector<float> helixParam(5);
675  std::vector<float> helixError(15);
676  convertHelix(trk_fit, HepPoint3D(0., 0., 0.), helixParam, helixError);
677 
678  // Checking for invalid helix curvature with parameter 2 equal to 0:
679  if (helixParam[2] == 0) {
680  B2WARNING("Helix parameter for curvature == 0. Skipping Track! The parameter is: " << helixParam[2] << "...");
681  continue;
682  }
683 
684  auto trackFitM = m_trackFitResults.appendNew(helixParam, helixError, pTypeM, pValue, -1, -1, 0);
685 
686  trackFitMIndex = trackFitM->getArrayIndex();
687 
688  daughterM = Particle(trackID[1] - 1, trackFitM, pTypeM);
689  // set momentum/positions at pivot = V0 decay vertex
690  getHelixParameters(trk_fit, pTypeM.getMass(), dauPivot,
691  helixParam, error5x5,
692  momentumM, positionM, error7x7M);
693 
694  for (unsigned i = 0; i < 7; i++)
695  for (unsigned j = 0; j < 7; j++)
696  errMatrixM(i, j) = error7x7M[i][j];
697 
698  daughterM.updateMomentum(TLorentzVector(momentumM.px(), momentumM.py(), momentumM.pz(), momentumM.e()),
699  TVector3(positionM.x(), positionM.y(), positionM.z()),
700  errMatrixM, pValue);
701  }
702  }
703 
704  Track* trackP = m_tracks[trackID[0] - 1];
705  Track* trackM = m_tracks[trackID[1] - 1];
706 
707  TrackFitResult* trackFitP = m_trackFitResults[trackFitPIndex];
708  TrackFitResult* trackFitM = m_trackFitResults[trackFitMIndex];
709 
710  m_v0s.appendNew(std::make_pair(trackP, trackFitP), std::make_pair(trackM, trackFitM));
711 
712  // create Ks Particle and add it to the 'K_S0:mdst' ParticleList
713  const PIDLikelihood* pidP = trackP->getRelated<PIDLikelihood>();
714  const PIDLikelihood* pidM = trackM->getRelated<PIDLikelihood>();
715  const MCParticle* mcParticleP = trackP->getRelated<MCParticle>();
716  const MCParticle* mcParticleM = trackM->getRelated<MCParticle>();
717 
718  Particle* newDaugP = m_particles.appendNew(daughterP);
719  if (pidP)
720  newDaugP->addRelationTo(pidP);
721  if (mcParticleP)
722  newDaugP->addRelationTo(mcParticleP);
723  Particle* newDaugM = m_particles.appendNew(daughterM);
724  if (pidM)
725  newDaugM->addRelationTo(pidM);
726  if (mcParticleM)
727  newDaugM->addRelationTo(mcParticleM);
728 
729  TLorentzVector v0Momentum(belleV0.px(), belleV0.py(), belleV0.pz(), belleV0.energy());
730  TVector3 v0Vertex(belleV0.vx(), belleV0.vy(), belleV0.vz());
731 
732  /*
733  * Documentation of Mdst_vee2 vertex fit:
734  * /sw/belle/belle/b20090127_0910/share/tables/mdst.tdf (L96-L125)
735  */
736  auto appendVertexFitInfo = [](Belle::Mdst_vee2 & _belle_V0, Particle & _belle2_V0) {
737  // Add chisq of vertex fit. chiSquared=10^10 means the fit fails.
738  _belle2_V0.addExtraInfo("chiSquared", _belle_V0.chisq());
739  // Ndf of the vertex kinematic fit is 1
740  _belle2_V0.addExtraInfo("ndf", 1);
741  // Add p-value to extra Info
742  double prob = TMath::Prob(_belle_V0.chisq(), 1);
743  _belle2_V0.setPValue(prob);
744  };
745 
746  Particle* newV0 = nullptr;
747  if (belleV0.kind() == 1) { // K0s -> pi+ pi-
748  Particle KS(v0Momentum, 310);
749  KS.appendDaughter(newDaugP);
750  KS.appendDaughter(newDaugM);
751  KS.setVertex(v0Vertex);
752  appendVertexFitInfo(belleV0, KS);
753  newV0 = m_particles.appendNew(KS);
754  ksPList->addParticle(newV0);
755 
756  // append extra info: goodKs flag
757  Belle::FindKs belleKSFinder;
758  belleKSFinder.candidates(belleV0, Belle::IpProfile::position(1));
759  newV0->addExtraInfo("goodKs", belleKSFinder.goodKs());
760 
761  /*
762  std::cout << " ---- B1 Ks ---- " << std::endl;
763  std::cout << " momentum = " << std::endl;
764  v0Momentum.Print();
765  std::cout << " position = " << std::endl;
766  v0Vertex.Print();
767  std::cout << " ---- B2 Ks ---- " << std::endl;
768  std::cout << " momentum = " << std::endl;
769  newKS->get4Vector().Print();
770  std::cout << " position = " << std::endl;
771  newKS->getVertex().Print();
772  std::cout << " ---- B1 Ks.child(0) ---- " << std::endl;
773  std::cout << " momentum = " << momentumP << std::endl;
774  std::cout << " position = " << positionP << std::endl;
775  std::cout << " error7x7 = " << error7x7P << std::endl;
776  std::cout << " ---- B2 Ks.child(0) ---- " << std::endl;
777  std::cout << " momentum = " << std::endl;
778  newKS->getDaughter(0)->get4Vector().Print();
779  std::cout << " position = " << std::endl;
780  newKS->getDaughter(0)->getVertex().Print();
781  std::cout << " error7x7 = " << std::endl;
782  newKS->getDaughter(0)->getMomentumVertexErrorMatrix().Print();
783  std::cout << " ---- B1 Ks.child(1) ---- " << std::endl;
784  std::cout << " momentum = " << momentumM << std::endl;
785  std::cout << " position = " << positionM << std::endl;
786  std::cout << " error7x7 = " << error7x7M << std::endl;
787  std::cout << " ---- B2 Ks.child(1) ---- " << std::endl;
788  std::cout << " momentum = " << std::endl;
789  newKS->getDaughter(1)->get4Vector().Print();
790  std::cout << " position = " << std::endl;
791  newKS->getDaughter(1)->getVertex().Print();
792  std::cout << " error7x7 = " << std::endl;
793  newKS->getDaughter(1)->getMomentumVertexErrorMatrix().Print();
794  */
795  } else if (belleV0.kind() == 2) { // Lambda -> p+ pi-
796  Particle Lambda0(v0Momentum, 3122);
797  Lambda0.appendDaughter(newDaugP);
798  Lambda0.appendDaughter(newDaugM);
799  Lambda0.setVertex(v0Vertex);
800  appendVertexFitInfo(belleV0, Lambda0);
801  newV0 = m_particles.appendNew(Lambda0);
802  lambda0PList->addParticle(newV0);
803 
804  // GoodLambda flag as extra info
805  Belle::FindLambda lambdaFinder;
806  lambdaFinder.candidates(belleV0, Belle::IpProfile::position(1));
807  newV0->addExtraInfo("goodLambda", lambdaFinder.goodLambda());
808  } else if (belleV0.kind() == 3) { // anti-Lambda -> pi+ anti-p
809  Particle antiLambda0(v0Momentum, -3122);
810  antiLambda0.appendDaughter(newDaugM);
811  antiLambda0.appendDaughter(newDaugP);
812  antiLambda0.setVertex(v0Vertex);
813  appendVertexFitInfo(belleV0, antiLambda0);
814  newV0 = m_particles.appendNew(antiLambda0);
815  antiLambda0PList->addParticle(newV0);
816 
817  // GoodLambda flag as extra info
818  Belle::FindLambda lambdaFinder;
819  lambdaFinder.candidates(belleV0, Belle::IpProfile::position(1));
820  newV0->addExtraInfo("goodLambda", lambdaFinder.goodLambda());
821  } else if (belleV0.kind() == 4) { // gamma -> e+ e-
822  Particle gamma(v0Momentum, 22);
823  gamma.appendDaughter(newDaugP);
824  gamma.appendDaughter(newDaugM);
825  gamma.setVertex(v0Vertex);
826  appendVertexFitInfo(belleV0, gamma);
827  newV0 = m_particles.appendNew(gamma);
828  convGammaPList->addParticle(newV0);
829  }
830  // append extra info: nisKsFinder quality indicators
831  if (m_nisEnable) {
832  if (belleV0.kind() <= 3) { // K_S0, Lambda, anti-Lambda
833  Belle::nisKsFinder ksnb;
834  double protIDP = atcPID(pidP, 2, 4);
835  double protIDM = atcPID(pidM, 2, 4);
836  ksnb.candidates(belleV0, Belle::IpProfile::position(1), momentumP, protIDP, protIDM);
837  // K_S0 and Lambda (inverse cut on ksnbNoLam for Lambda selection).
838  newV0->addExtraInfo("ksnbVLike", ksnb.nb_vlike());
839  newV0->addExtraInfo("ksnbNoLam", ksnb.nb_nolam());
840  // K_S0 only
841  if (belleV0.kind() == 1)
842  newV0->addExtraInfo("ksnbStandard", ksnb.standard());
843  }
844  }
845  }
846 }
847 
848 void B2BIIConvertMdstModule::convertGenHepEvtTable()
849 {
850  if (m_realData)
851  return;
852 
853  // clear the Gen_hepevt_ID <-> MCParticleGraphPosition map
854  genHepevtToMCParticle.clear();
855 
856  // check if the Gen_hepevt table has any entries
857  Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
858  if (genMgr.count() == 0)
859  return;
860 
861  typedef std::pair<MCParticleGraph::GraphParticle*, Belle::Gen_hepevt> halfFamily;
862  halfFamily currFamily;
863  halfFamily family;
864  std::queue < halfFamily > heritancesQueue;
865 
866  // Add motherless particles. The root particle is often the first one,
867  // but this is not correct in general case; thus, all particles are added,
868  // including the beam-background ones.
869  m_particleGraph.clear();
870  for (Belle::Gen_hepevt_Manager::iterator genIterator = genMgr.begin();
871  genIterator != genMgr.end(); ++genIterator) {
872  Belle::Gen_hepevt hep = *genIterator;
873  // Select particles without mother.
874  if (!(hep.moFirst() == 0 && hep.moLast() == 0))
875  continue;
876  // Ignore particles with code 911, they are used for CDC data.
877  if (hep.idhep() == 911)
878  continue;
879  int position = m_particleGraph.size();
880  m_particleGraph.addParticle();
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 or 911 (CDC data).
905  if (currDaughter.idhep() == 0 || currDaughter.idhep() == 911)
906  continue;
907 
908  //putting the daughter in the graph:
909  int position = m_particleGraph.size();
910  m_particleGraph.addParticle();
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 
935  m_particleGraph.generateList();
936 }
937 
938 void B2BIIConvertMdstModule::convertMdstECLTable()
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 
998 void B2BIIConvertMdstModule::convertMdstKLMTable()
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 
1024 void B2BIIConvertMdstModule::convertMdstGammaTable()
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 
1071 void B2BIIConvertMdstModule::convertMdstPi0Table()
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  TLorentzVector p4(mdstPi0.px(), mdstPi0.py(), mdstPi0.pz(), mdstPi0.energy());
1090 
1091  // Create Particle from TLorentzVector 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 
1119 void B2BIIConvertMdstModule::convertMdstKLongTable()
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(130, "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()) == 130 && 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 
1204 void B2BIIConvertMdstModule::convertEvtclsTable()
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 
1249 void B2BIIConvertMdstModule::convertRecTrgTable()
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
1257  StoreObjPtr<EventMetaData> event;
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;
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);
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);
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);
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) {
1637  const Const::ChargedStable& pType = c_belleHyp_to_chargedStable[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  auto cdcExtraInfo = m_belleTrkExtra.appendNew(trk_fit.first_x(), trk_fit.first_y(), trk_fit.first_z(),
1680  trk_fit.last_x(), trk_fit.last_y(), trk_fit.last_z());
1681  track->addRelationTo(cdcExtraInfo);
1682  }
1683  // conversion of the SVD hit pattern
1684  int svdHitPattern = trk_fit.hit_svd();
1685  // use hits from 3: SVD-rphi, 4: SVD-z
1686  // int svdNHits = trk_fit.nhits(3) + trk_fit.nhits(4);
1687 
1688  std::bitset<32> svdBitSet(svdHitPattern);
1689 
1690  HitPatternVXD patternVxd;
1691 
1692  unsigned short svdLayers;
1693  // taken from: http://belle.kek.jp/group/indirectcp/cpfit/cpfit-festa/2004/talks/Apr.14/CPfesta-2005-Higuchi(3).pdf
1694  StoreObjPtr<EventMetaData> event;
1695  // mask for the rphi hits, first 6 (8) bits/ 2 bits per layer
1696  std::bitset<32> svdUMask(static_cast<std::string>("00000000000000000000000000000011"));
1697  // mask for the z hits, second 6 (8) bits/ 2 bits per layer
1698  std::bitset<32> svdVMask;
1699 
1700  // find out if the SVD has 3 (4) layers; if exp <= (>) exp 27
1701  if (event->getExperiment() <= 27) {
1702  svdVMask = svdUMask << 6;
1703  svdLayers = 3;
1704  } else {
1705  svdVMask = svdUMask << 8;
1706  svdLayers = 4;
1707  }
1708 
1709  // loop over all svd layers (layer index is shifted + 3 for basf2)
1710  for (unsigned short layerId = 0; layerId < svdLayers; layerId++) {
1711  unsigned short uHits = (svdBitSet & svdUMask).count();
1712  unsigned short vHits = (svdBitSet & svdVMask).count();
1713  patternVxd.setSVDLayer(layerId + 3, uHits, vHits);
1714  // shift masks to the left
1715  svdUMask <<= 2;
1716  svdVMask <<= 2;
1717  }
1718 
1719  TrackFitResult helixFromHelix(helixParam, helixError, pType, pValue, -1, patternVxd.getInteger(), 0);
1720 
1721  if (m_use6x6CovarianceMatrix4Tracks) {
1722  TMatrixDSym cartesianCovariance(6);
1723  for (unsigned i = 0; i < 7; i++) {
1724  if (i == 3)
1725  continue;
1726  for (unsigned j = 0; j < 7; j++) {
1727  if (j == 3)
1728  continue;
1729 
1730  cartesianCovariance(ERRMCONV[i], ERRMCONV[j]) = error7x7[i][j];
1731  }
1732  }
1733  UncertainHelix helixFromCartesian(helixFromHelix.getPosition(), helixFromHelix.getMomentum(), helixFromHelix.getChargeSign(),
1734  BFIELD, cartesianCovariance, pValue);
1735 
1736  TMatrixDSym helixCovariance = helixFromCartesian.getCovariance();
1737 
1738  counter = 0;
1739  for (unsigned int i = 0; i < 5; ++i)
1740  for (unsigned int j = i; j < 5; ++j)
1741  helixError[counter++] = helixCovariance(i, j);
1742  }
1743 
1744  auto trackFit = m_trackFitResults.appendNew(helixParam, helixError, pType, pValue, patternCdc.getInteger(),
1745  patternVxd.getInteger(), trk_fit.ndf());
1746  track->setTrackFitResultIndex(pType, trackFit->getArrayIndex());
1747  /*
1748  B2INFO("--- B1 Track: ");
1749  std::cout << "Momentum = " << momentum << std::endl;
1750  std::cout << "Position = " << position << std::endl;
1751  std::cout << "7x7 error matrix = " << error7x7 << std::endl;
1752  B2INFO("--- B2 Track: ");
1753  std::cout << "Momentum = " << std::endl;
1754  trackFit->get4Momentum().Print();
1755  std::cout << "Position = " << std::endl;
1756  trackFit->getPosition().Print();
1757  std::cout << "6x6 error matrix = " << std::endl;
1758  trackFit->getCovariance6().Print();
1759  TMatrixDSym b2Error7x7(7);
1760  fill7x7ErrorMatrix(trackFit, b2Error7x7, thisMass, 1.5);
1761  std::cout << "7x7 error matrix = " << std::endl;
1762  b2Error7x7.Print();
1763  */
1764  }
1765 }
1766 
1767 void B2BIIConvertMdstModule::convertGenHepevtObject(const Belle::Gen_hepevt& genHepevt, MCParticleGraph::GraphParticle* mcParticle)
1768 {
1769  //B2DEBUG(80, "Gen_ehepevt: idhep " << genHepevt.idhep() << " (" << genHepevt.isthep() << ") with ID = " << genHepevt.get_ID());
1770 
1771  // updating the GraphParticle information from the Gen_hepevt information
1772  const int idHep = recoverMoreThan24bitIDHEP(genHepevt.idhep());
1773  if (idHep == 0) {
1774  B2WARNING("Trying to convert Gen_hepevt with idhep = " << idHep <<
1775  ". This should never happen.");
1776  mcParticle->setPDG(22);
1777  } else {
1778  mcParticle->setPDG(idHep);
1779  }
1780 
1781  if (genHepevt.isthep() > 0) {
1782  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1783  }
1784 
1785  mcParticle->setMass(genHepevt.M());
1786 
1787  TLorentzVector p4(genHepevt.PX(), genHepevt.PY(), genHepevt.PZ(), genHepevt.E());
1788  mcParticle->set4Vector(p4);
1789 
1790  mcParticle->setProductionVertex(genHepevt.VX()*Unit::mm, genHepevt.VY()*Unit::mm, genHepevt.VZ()*Unit::mm);
1791  mcParticle->setProductionTime(genHepevt.T()*Unit::mm / Const::speedOfLight);
1792 
1793  // decay time of this particle is production time of the daughter particle
1794  if (genHepevt.daFirst() > 0) {
1795  Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
1796  Belle::Gen_hepevt daughterParticle = genMgr(Belle::Panther_ID(genHepevt.daFirst()));
1797  mcParticle->setDecayTime(daughterParticle.T()*Unit::mm / Const::speedOfLight);
1798  mcParticle->setDecayVertex(daughterParticle.VX()*Unit::mm, daughterParticle.VY()*Unit::mm, daughterParticle.VZ()*Unit::mm);
1799  } else {
1800  //otherwise, assume it's stable
1801  mcParticle->setDecayTime(std::numeric_limits<float>::infinity());
1802  }
1803 
1804  mcParticle->setValidVertex(true);
1805 }
1806 
1807 void B2BIIConvertMdstModule::convertMdstECLObject(const Belle::Mdst_ecl& ecl, const Belle::Mdst_ecl_aux& eclAux,
1808  ECLCluster* eclCluster)
1809 {
1810  if (eclAux.e9oe25() < m_matchType2E9oE25Threshold)
1811  eclCluster->setIsTrack(ecl.match() > 0);
1812  else
1813  eclCluster->setIsTrack(ecl.match() == 1);
1814 
1815  eclCluster->setEnergy(ecl.energy()); //must happen before setCovarianceMatrix()!
1816  eclCluster->setPhi(ecl.phi());
1817  eclCluster->setTheta(ecl.theta());
1818  eclCluster->setR(ecl.r());
1819  eclCluster->setdeltaL(ecl.quality());
1820 
1821  double covarianceMatrix[6];
1822  covarianceMatrix[0] = ecl.error(0); // error on energy
1823  covarianceMatrix[1] = ecl.error(1);
1824  covarianceMatrix[2] = ecl.error(2); // error on phi
1825  covarianceMatrix[3] = ecl.error(3);
1826  covarianceMatrix[4] = ecl.error(4);
1827  covarianceMatrix[5] = ecl.error(5); // error on theta
1828  eclCluster->setCovarianceMatrix(covarianceMatrix);
1829 
1830  eclCluster->setLAT(eclAux.width());
1831  eclCluster->setEnergyRaw(eclAux.mass());
1832  eclCluster->setE9oE21(eclAux.e9oe25());
1833  eclCluster->setEnergyHighestCrystal(eclAux.seed());
1834  eclCluster->setTime(eclAux.property(0));
1835  eclCluster->setNumberOfCrystals(eclAux.nhits());
1836 }
1837 
1838 void B2BIIConvertMdstModule::convertMdstKLMObject(const Belle::Mdst_klm_cluster& klm_cluster, KLMCluster* klmCluster)
1839 {
1840  // note: Belle quality flag is not saved (no free int variable in Belle2 KLMCluster)
1841  klmCluster->setLayers(klm_cluster.layers());
1842  klmCluster->setInnermostLayer(klm_cluster.first_layer());
1843 }
1844 
1845 
1846 //-----------------------------------------------------------------------------
1847 // RELATIONS
1848 //-----------------------------------------------------------------------------
1849 void B2BIIConvertMdstModule::setTracksToECLClustersRelations()
1850 {
1851  // Relations
1852  RelationArray tracksToECLClusters(m_tracks, m_eclClusters);
1853 
1854  Belle::Mdst_ecl_trk_Manager& m = Belle::Mdst_ecl_trk_Manager::get_manager();
1855  Belle::Mdst_charged_Manager& chgMg = Belle::Mdst_charged_Manager::get_manager();
1856 
1857  // We first insert relations to tracks which are directly matched (type == 1)
1858  // secondly we had CR matched tracks (connected region) (type == 2)
1859  // finally tracks which are geometrically matched (type == 0)
1860  std::vector<int> insert_order_types = {1, 2, 0};
1861  for (auto& insert_type : insert_order_types) {
1862  for (Belle::Mdst_ecl_trk_Manager::iterator ecltrkIterator = m.begin(); ecltrkIterator != m.end(); ++ecltrkIterator) {
1863  Belle::Mdst_ecl_trk mECLTRK = *ecltrkIterator;
1864 
1865  if (mECLTRK.type() != insert_type)
1866  continue;
1867 
1868  Belle::Mdst_ecl mdstEcl = mECLTRK.ecl();
1869  Belle::Mdst_trk mTRK = mECLTRK.trk();
1870 
1871  if (!mdstEcl)
1872  continue;
1873 
1874  // the numbering in mdst_charged
1875  // not necessarily the same as in mdst_trk
1876  // therfore have to find corresponding mdst_charged
1877  for (Belle::Mdst_charged_Manager::iterator chgIterator = chgMg.begin(); chgIterator != chgMg.end(); ++chgIterator) {
1878  Belle::Mdst_charged mChar = *chgIterator;
1879  Belle::Mdst_trk mTRK_in_charged = mChar.trk();
1880 
1881  if (mTRK_in_charged.get_ID() == mTRK.get_ID()) {
1882  // found the correct mdst_charged
1883  tracksToECLClusters.add(mChar.get_ID() - 1, mdstEcl.get_ID() - 1, 1.0);
1884  break;
1885  }
1886  }
1887  }
1888  }
1889 }
1890 
1891 
1892 void B2BIIConvertMdstModule::setKLMClustersRelations()
1893 {
1894  // Relations
1895  RelationArray klmClustersToTracks(m_klmClusters, m_tracks);
1896  RelationArray klmClustersToEclClusters(m_klmClusters, m_eclClusters);
1897 
1898  Belle::Mdst_klm_cluster_Manager& klm_cluster_manager = Belle::Mdst_klm_cluster_Manager::get_manager();
1899 
1900 
1901  for (Belle::Mdst_klm_cluster_Manager::iterator klmC_Ite = klm_cluster_manager.begin(); klmC_Ite != klm_cluster_manager.end();
1902  ++klmC_Ite) {
1903 
1904  Belle::Mdst_klm_cluster mdstKlm_cluster = *klmC_Ite;
1905  Belle::Mdst_trk mTRK = mdstKlm_cluster.trk();
1906  Belle::Mdst_ecl mECL = mdstKlm_cluster.ecl();
1907 
1908  if (mTRK) klmClustersToTracks.add(mdstKlm_cluster.get_ID() - 1, mTRK.get_ID() - 1);
1909  if (mECL) klmClustersToEclClusters.add(mdstKlm_cluster.get_ID() - 1, mECL.get_ID() - 1);
1910  }
1911 }
1912 
1913 
1914 //-----------------------------------------------------------------------------
1915 // MISC
1916 //-----------------------------------------------------------------------------
1917 
1918 int B2BIIConvertMdstModule::recoverMoreThan24bitIDHEP(int id)
1919 {
1920  /*
1921  QUICK CHECK: most of the normal particles are smaller than
1922  0x100000, while all the corrupt id has some of the high bits on.
1923 
1924  This bit check has to be revised when the table below is updated.
1925  */
1926  const int mask = 0x00f00000;
1927  int high_bits = id & mask;
1928  if (high_bits == 0 || high_bits == mask) return id;
1929 
1930  switch (id) {
1931  case 7114363:
1932  return 91000443; // X(3940)
1933  case 6114363:
1934  return 90000443; // Y(3940)
1935  case 6114241:
1936  return 90000321; // K_0*(800)+
1937  case 6114231:
1938  return 90000311; // K_0*(800)0
1939  case -6865004:
1940  return 9912212; // p_diff+
1941  case -6865104:
1942  return 9912112; // n_diffr
1943  case -6866773:
1944  return 9910443; // psi_diff
1945  case -6866883:
1946  return 9910333; // phi_diff
1947  case -6866993:
1948  return 9910223; // omega_diff
1949  case -6867005:
1950  return 9910211; // pi_diff+
1951  case -6867103:
1952  return 9910113; // rho_diff0
1953  case -7746995:
1954  return 9030221; // f_0(1500)
1955  case -7756773:
1956  return 9020443; // psi(4415)
1957  case -7756995:
1958  return 9020221; // eta(1405)
1959  case -7766773:
1960  return 9010443; // psi(4160)
1961  case -7776663:
1962  return 9000553; // Upsilon(5S)
1963  case -7776773:
1964  return 9000443; // psi(4040)
1965  case -7776783:
1966  return 9000433; // D_sj(2700)+
1967  case -7776995:
1968  return 9000221; // f_0(600)
1969  case -6114241:
1970  return -90000321; // K_0*(800)-
1971  case -6114231:
1972  return -90000311; // anti-K_0*(800)0
1973  case 6865004:
1974  return -9912212; // anti-p_diff-
1975  case 6865104:
1976  return -9912112; // anti-n_diffr
1977  case 6867005:
1978  return -9910211; // pi_diff-
1979  case 7776783:
1980  return -9000433; // D_sj(2700)-
1981  default:
1982  return id;
1983  }
1984 }
1985 
1986 void B2BIIConvertMdstModule::testMCRelation(const Belle::Gen_hepevt& belleMC, const MCParticle* mcP, const std::string& objectName)
1987 {
1988  int bellePDGCode = belleMC.idhep();
1989  int belleIIPDGCode = mcP->getPDG();
1990 
1991  if (bellePDGCode == 0)
1992  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to Gen_hepevt with idhep = 0.");
1993 
1994  if (bellePDGCode != belleIIPDGCode)
1995  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to different MCParticle! " << bellePDGCode << " vs. " <<
1996  belleIIPDGCode);
1997 
1998  double belleMomentum[] = { belleMC.PX(), belleMC.PY(), belleMC.PZ() };
1999  double belle2Momentum[] = { mcP->get4Vector().Px(), mcP->get4Vector().Py(), mcP->get4Vector().Pz() };
2000 
2001  for (unsigned i = 0; i < 3; i++) {
2002  double relDev = (belle2Momentum[i] - belleMomentum[i]) / belleMomentum[i];
2003 
2004  if (relDev > 1e-3) {
2005  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to different MCParticle!");
2006  B2INFO(" - Gen_hepevt [" << bellePDGCode << "] px/py/pz = " << belleMC.PX() << "/" << belleMC.PY() << "/" << belleMC.PZ());
2007  B2INFO(" - TrackFitResult [" << belleIIPDGCode << "] px/py/pz = " << mcP->get4Vector().Px() << "/" << mcP->get4Vector().Py() << "/"
2008  << mcP->get4Vector().Pz());
2009  }
2010  }
2011 }
2012 
2013 void B2BIIConvertMdstModule::belleVeeDaughterToCartesian(const Belle::Mdst_vee2& vee, const int charge,
2014  const Const::ParticleType& pType,
2015  CLHEP::HepLorentzVector& momentum, HepPoint3D& position, CLHEP::HepSymMatrix& error)
2016 {
2017  const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2018  CLHEP::HepVector a(5);
2019  CLHEP::HepSymMatrix Ea(5, 0);
2020  if (charge > 0) {
2021  a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2022  a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2023  a[4] = vee.daut().helix_p(4);
2024  Ea[0][0] = vee.daut().error_p(0); Ea[1][0] = vee.daut().error_p(1);
2025  Ea[1][1] = vee.daut().error_p(2); Ea[2][0] = vee.daut().error_p(3);
2026  Ea[2][1] = vee.daut().error_p(4); Ea[2][2] = vee.daut().error_p(5);
2027  Ea[3][0] = vee.daut().error_p(6); Ea[3][1] = vee.daut().error_p(7);
2028  Ea[3][2] = vee.daut().error_p(8); Ea[3][3] = vee.daut().error_p(9);
2029  Ea[4][0] = vee.daut().error_p(10); Ea[4][1] = vee.daut().error_p(11);
2030  Ea[4][2] = vee.daut().error_p(12); Ea[4][3] = vee.daut().error_p(13);
2031  Ea[4][4] = vee.daut().error_p(14);
2032  } else {
2033  a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2034  a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2035  a[4] = vee.daut().helix_m(4);
2036  Ea[0][0] = vee.daut().error_m(0); Ea[1][0] = vee.daut().error_m(1);
2037  Ea[1][1] = vee.daut().error_m(2); Ea[2][0] = vee.daut().error_m(3);
2038  Ea[2][1] = vee.daut().error_m(4); Ea[2][2] = vee.daut().error_m(5);
2039  Ea[3][0] = vee.daut().error_m(6); Ea[3][1] = vee.daut().error_m(7);
2040  Ea[3][2] = vee.daut().error_m(8); Ea[3][3] = vee.daut().error_m(9);
2041  Ea[4][0] = vee.daut().error_m(10); Ea[4][1] = vee.daut().error_m(11);
2042  Ea[4][2] = vee.daut().error_m(12); Ea[4][3] = vee.daut().error_m(13);
2043  Ea[4][4] = vee.daut().error_m(14);
2044  }
2045 
2046  Belle::Helix helix(pivot, a, Ea);
2047 
2048  // this is Vee daughter momentum/position/error at pivot = V0 Decay Vertex
2049  momentum = helix.momentum(0., pType.getMass(), position, error);
2050 }
2051 
2052 void B2BIIConvertMdstModule::belleVeeDaughterHelix(const Belle::Mdst_vee2& vee, const int charge, std::vector<float>& helixParam,
2053  std::vector<float>& helixError)
2054 {
2055  const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2056  CLHEP::HepVector a(5);
2057  CLHEP::HepSymMatrix Ea(5, 0);
2058  if (charge > 0) {
2059  a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2060  a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2061  a[4] = vee.daut().helix_p(4);
2062  Ea[0][0] = vee.daut().error_p(0);
2063  Ea[1][0] = vee.daut().error_p(1);
2064  Ea[1][1] = vee.daut().error_p(2);
2065  Ea[2][0] = vee.daut().error_p(3);
2066  Ea[2][1] = vee.daut().error_p(4);
2067  Ea[2][2] = vee.daut().error_p(5);
2068  Ea[3][0] = vee.daut().error_p(6);
2069  Ea[3][1] = vee.daut().error_p(7);
2070  Ea[3][2] = vee.daut().error_p(8);
2071  Ea[3][3] = vee.daut().error_p(9);
2072  Ea[4][0] = vee.daut().error_p(10);
2073  Ea[4][1] = vee.daut().error_p(11);
2074  Ea[4][2] = vee.daut().error_p(12);
2075  Ea[4][3] = vee.daut().error_p(13);
2076  Ea[4][4] = vee.daut().error_p(14);
2077  } else {
2078  a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2079  a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2080  a[4] = vee.daut().helix_m(4);
2081  Ea[0][0] = vee.daut().error_m(0);
2082  Ea[1][0] = vee.daut().error_m(1);
2083  Ea[1][1] = vee.daut().error_m(2);
2084  Ea[2][0] = vee.daut().error_m(3);
2085  Ea[2][1] = vee.daut().error_m(4);
2086  Ea[2][2] = vee.daut().error_m(5);
2087  Ea[3][0] = vee.daut().error_m(6);
2088  Ea[3][1] = vee.daut().error_m(7);
2089  Ea[3][2] = vee.daut().error_m(8);
2090  Ea[3][3] = vee.daut().error_m(9);
2091  Ea[4][0] = vee.daut().error_m(10);
2092  Ea[4][1] = vee.daut().error_m(11);
2093  Ea[4][2] = vee.daut().error_m(12);
2094  Ea[4][3] = vee.daut().error_m(13);
2095  Ea[4][4] = vee.daut().error_m(14);
2096  }
2097 
2098  Belle::Helix helix(pivot, a, Ea);
2099 
2100  // go to the new pivot
2101  helix.pivot(HepPoint3D(0., 0., 0.));
2102 
2103  CLHEP::HepSymMatrix error5x5(5, 0);
2104  convertHelix(helix, helixParam, error5x5);
2105 
2106  unsigned int size = 5;
2107  unsigned int counter = 0;
2108  for (unsigned int i = 0; i < size; i++)
2109  for (unsigned int j = i; j < size; j++)
2110  helixError[counter++] = error5x5[i][j];
2111 }
2112 
2113 void B2BIIConvertMdstModule::belleHelixToHelix(const Belle::Mdst_trk_fit& trk_fit,
2114  std::vector<float>& helixParam, std::vector<float>& helixError)
2115 {
2116  const HepPoint3D pivot(trk_fit.pivot_x(),
2117  trk_fit.pivot_y(),
2118  trk_fit.pivot_z());
2119 
2120  CLHEP::HepVector a(5);
2121  a[0] = trk_fit.helix(0);
2122  a[1] = trk_fit.helix(1);
2123  a[2] = trk_fit.helix(2);
2124  a[3] = trk_fit.helix(3);
2125  a[4] = trk_fit.helix(4);
2126  CLHEP::HepSymMatrix Ea(5, 0);
2127  Ea[0][0] = trk_fit.error(0);
2128  Ea[1][0] = trk_fit.error(1);
2129  Ea[1][1] = trk_fit.error(2);
2130  Ea[2][0] = trk_fit.error(3);
2131  Ea[2][1] = trk_fit.error(4);
2132  Ea[2][2] = trk_fit.error(5);
2133  Ea[3][0] = trk_fit.error(6);
2134  Ea[3][1] = trk_fit.error(7);
2135  Ea[3][2] = trk_fit.error(8);
2136  Ea[3][3] = trk_fit.error(9);
2137  Ea[4][0] = trk_fit.error(10);
2138  Ea[4][1] = trk_fit.error(11);
2139  Ea[4][2] = trk_fit.error(12);
2140  Ea[4][3] = trk_fit.error(13);
2141  Ea[4][4] = trk_fit.error(14);
2142  Belle::Helix helix(pivot, a, Ea);
2143 
2144  helix.pivot(HepPoint3D(0., 0., 0.));
2145 
2146  // convert to Belle II helix parameters
2147  // param 0: d_0 = d_rho
2148  helixParam[0] = helix.a()[0];
2149 
2150  // param 1: phi = phi_0 + pi/2
2151  helixParam[1] = helix.a()[1] + TMath::Pi() / 2.0;
2152 
2153  // param 2: omega = Kappa * alpha = Kappa * B[Tesla] * speed_of_light[m/s] * 1e-11
2154  helixParam[2] = helix.a()[2] * KAPPA2OMEGA;
2155 
2156  // param 3: d_z = z0
2157  helixParam[3] = helix.a()[3];
2158 
2159  // param 4: tan(Lambda) = cotTheta
2160  helixParam[4] = helix.a()[4];
2161 
2162  Ea = helix.Ea();
2163 
2164  helixError[0] = Ea[0][0];
2165  helixError[1] = Ea[1][0];
2166  helixError[2] = Ea[2][0] * KAPPA2OMEGA;
2167  helixError[3] = Ea[3][0];
2168  helixError[4] = Ea[4][0];
2169  helixError[5] = Ea[1][1];
2170  helixError[6] = Ea[2][1] * KAPPA2OMEGA;
2171  helixError[7] = Ea[3][1];
2172  helixError[8] = Ea[4][1];
2173  helixError[9] = Ea[2][2] * KAPPA2OMEGA * KAPPA2OMEGA;
2174  helixError[10] = Ea[3][2] * KAPPA2OMEGA;
2175  helixError[11] = Ea[4][2] * KAPPA2OMEGA;
2176  helixError[12] = Ea[3][3];
2177  helixError[13] = Ea[4][3];
2178  helixError[14] = Ea[4][4];
2179 }
2180 
2181 int B2BIIConvertMdstModule::belleHelixToCartesian(const Belle::Mdst_trk_fit& trk_fit, const double mass,
2182  const HepPoint3D& newPivot,
2183  CLHEP::HepLorentzVector& momentum, HepPoint3D& position, CLHEP::HepSymMatrix& error, double dPhi)
2184 {
2185  const HepPoint3D pivot(trk_fit.pivot_x(),
2186  trk_fit.pivot_y(),
2187  trk_fit.pivot_z());
2188 
2189  CLHEP::HepVector a(5);
2190  a[0] = trk_fit.helix(0);
2191  a[1] = trk_fit.helix(1);
2192  a[2] = trk_fit.helix(2);
2193  a[3] = trk_fit.helix(3);
2194  a[4] = trk_fit.helix(4);
2195  CLHEP::HepSymMatrix Ea(5, 0);
2196  Ea[0][0] = trk_fit.error(0);
2197  Ea[1][0] = trk_fit.error(1);
2198  Ea[1][1] = trk_fit.error(2);
2199  Ea[2][0] = trk_fit.error(3);
2200  Ea[2][1] = trk_fit.error(4);
2201  Ea[2][2] = trk_fit.error(5);
2202  Ea[3][0] = trk_fit.error(6);
2203  Ea[3][1] = trk_fit.error(7);
2204  Ea[3][2] = trk_fit.error(8);
2205  Ea[3][3] = trk_fit.error(9);
2206  Ea[4][0] = trk_fit.error(10);
2207  Ea[4][1] = trk_fit.error(11);
2208  Ea[4][2] = trk_fit.error(12);
2209  Ea[4][3] = trk_fit.error(13);
2210  Ea[4][4] = trk_fit.error(14);
2211 
2212  Belle::Helix helix(pivot, a, Ea);
2213 
2214  int charge = 0;
2215  if (helix.kappa() > 0)
2216  charge = 1;
2217  else
2218  charge = -1;
2219 
2220  if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
2221  helix.pivot(newPivot);
2222  momentum = helix.momentum(dPhi, mass, position, error);
2223  } else {
2224  if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
2225  helix.pivot(HepPoint3D(0., 0., 0.));
2226  momentum = helix.momentum(dPhi, mass, position, error);
2227  } else {
2228  momentum = helix.momentum(dPhi, mass, position, error);
2229  }
2230  }
2231  return charge;
2232 }
2233 
2234 TrackFitResult B2BIIConvertMdstModule::createTrackFitResult(const CLHEP::HepLorentzVector& momentum,
2235  const HepPoint3D& position,
2236  const CLHEP::HepSymMatrix& error,
2237  const short int charge,
2238  const Const::ParticleType& pType,
2239  const float pValue,
2240  const uint64_t hitPatternCDCInitializer,
2241  const uint32_t hitPatternVXDInitializer,
2242  const uint16_t ndf)
2243 {
2244  TVector3 pos(position.x(), position.y(), position.z());
2245  TVector3 mom(momentum.px(), momentum.py(), momentum.pz());
2246 
2247  TMatrixDSym errMatrix(6);
2248  for (unsigned i = 0; i < 7; i++) {
2249  if (i == 3)
2250  continue;
2251  for (unsigned j = 0; j < 7; j++) {
2252  if (j == 3)
2253  continue;
2254 
2255  if (i == j)
2256  errMatrix(ERRMCONV[i], ERRMCONV[i]) = error[i][i];
2257  else
2258  errMatrix(ERRMCONV[i], ERRMCONV[j]) = errMatrix(ERRMCONV[j], ERRMCONV[i]) = error[i][j];
2259  }
2260  }
2261 
2262  return TrackFitResult(pos, mom, errMatrix, charge, pType, pValue, BFIELD, hitPatternCDCInitializer, hitPatternVXDInitializer, ndf);
2263 }
2264 
2265 double B2BIIConvertMdstModule::atcPID(const PIDLikelihood* pid, int sigHyp, int bkgHyp)
2266 {
2267  if (!pid) return 0.5;
2268 
2269  // ACC = ARICH
2270  Const::PIDDetectorSet set = Const::ARICH;
2271  double acc_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2272  double acc_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2273  double acc = 0.5;
2274  if (acc_sig + acc_bkg > 0.0)
2275  acc = acc_sig / (acc_sig + acc_bkg);
2276 
2277  // TOF = TOP
2278  set = Const::TOP;
2279  double tof_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2280  double tof_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2281  double tof = 0.5;
2282  double tof_all = tof_sig + tof_bkg;
2283  if (tof_all != 0) {
2284  tof = tof_sig / tof_all;
2285  if (tof < 0.001) tof = 0.001;
2286  if (tof > 0.999) tof = 0.999;
2287  }
2288 
2289  // dE/dx = CDC
2290  set = Const::CDC;
2291  double cdc_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2292  double cdc_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2293  double cdc = 0.5;
2294  double cdc_all = cdc_sig + cdc_bkg;
2295  if (cdc_all != 0) {
2296  cdc = cdc_sig / cdc_all;
2297  if (cdc < 0.001) cdc = 0.001;
2298  if (cdc > 0.999) cdc = 0.999;
2299  }
2300 
2301  // Combined
2302  double pid_sig = acc * tof * cdc;
2303  double pid_bkg = (1. - acc) * (1. - tof) * (1. - cdc);
2304 
2305  return pid_sig / (pid_sig + pid_bkg);
2306 }
2307 
2308 
2309 void B2BIIConvertMdstModule::endRun()
2310 {
2311  B2DEBUG(99, "B2BIIConvertMdst: endRun done.");
2312 }
2313 
2314 
2315 void B2BIIConvertMdstModule::terminate()
2316 {
2317  B2DEBUG(99, "B2BIIConvertMdst: terminate called");
2318 }
2319 
Belle2::EvtPDLUtil::charge
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:46
Belle2::StoreArray::appendNew
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:256
Belle2::B2BIIConvertMdstModule::convertMdstVee2Table
void convertMdstVee2Table()
Reads and converts all entries of Mdst_vee2 Panther table to V0 dataobjects and adds them to StoreArr...
Definition: B2BIIConvertMdstModule.cc:493
Belle2::IntervalOfValidity
A class that describes the interval of experiments/runs for which an object in the database is valid.
Definition: IntervalOfValidity.h:35
Belle2::B2BIIConvertMdstModule::m_eclClusters
StoreArray< ECLCluster > m_eclClusters
ECL clusters.
Definition: B2BIIConvertMdstModule.h:390
Belle2::MCParticleGraph::generateList
void generateList(const std::string &name="", int options=c_setNothing)
Generates the MCParticle list and stores it in the StoreArray with the given name.
Definition: MCParticleGraph.cc:211
Belle2::B2BIIConvertMdstModule::m_convertEvtcls
bool m_convertEvtcls
Flag to switch on conversion of Evtcls table.
Definition: B2BIIConvertMdstModule.h:156
Belle2::B2BIIConvertMdstModule::m_nisEnable
bool m_nisEnable
Flag to switch on conversion of nisKsFinder info.
Definition: B2BIIConvertMdstModule.h:158
Belle2::B2BIIConvertMdstModule::m_convertTrkExtra
bool m_convertTrkExtra
Flag to switch on conversion of first(last)_{x,y,z} of mdst_trk_fit.
Definition: B2BIIConvertMdstModule.h:162
Belle2::ECLCluster::setTime
void setTime(double time)
Set time information.
Definition: ECLCluster.h:223
Belle2::B2BIIConvertMdstModule::testMCRelation
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.
Definition: B2BIIConvertMdstModule.cc:1986
Belle2::RelationArray
Low-level class to create/modify relations between StoreArrays.
Definition: RelationArray.h:72
Belle2::ECLCluster::setIsTrack
void setIsTrack(bool istrack)
Set m_isTrack true if the cluster matches with a track.
Definition: ECLCluster.h:115
Belle2::B2BIIConvertMdstModule::mdstGammaToParticle
std::map< int, int > mdstGammaToParticle
map of gamma Panther IDs and corresponding Particle StoreArray indices
Definition: B2BIIConvertMdstModule.h:380
Belle2::MCParticleGraph::size
size_t size() const
Return the number of particles in the graph.
Definition: MCParticleGraph.h:253
Belle2::B2BIIConvertMdstModule::beginRun
virtual void beginRun() override
Module functions to be called from event process.
Definition: B2BIIConvertMdstModule.cc:262
Belle2::UncertainHelix
This class represents an ideal helix in perigee parameterization including the covariance matrix of t...
Definition: UncertainHelix.h:40
Belle2::B2BIIConvertMdstModule::convertMdstECLTable
void convertMdstECLTable()
Reads and converts all entries of Mdst_ecl(_aux) Panther table to ECLCluster dataobjects and adds the...
Definition: B2BIIConvertMdstModule.cc:938
Belle2::CollisionBoostVector::setBoost
void setBoost(const TVector3 &boost, const TMatrixDSym &covariance)
Set the boost vector and its error matrix.
Definition: CollisionBoostVector.h:46
Belle2::StoreArray::registerRelationTo
bool registerRelationTo(const StoreArray< TO > &toArray, DataStore::EDurability durability=DataStore::c_Event, DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut, const std::string &namedRelation="") const
Register a relation to the given StoreArray.
Definition: StoreArray.h:150
Belle2::B2BIIConvertMdstModule::convertIPProfile
void convertIPProfile(bool beginRun=false)
Stores the IPProfiles in BeamSpot (currently in DataStore)
Definition: B2BIIConvertMdstModule.cc:402
Belle2::B2BIIConvertMdstModule::createTrackFitResult
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.
Definition: B2BIIConvertMdstModule.cc:2234
Belle2::B2BIIConvertMdstModule::m_collisionBoostVector
CollisionBoostVector m_collisionBoostVector
CollisionBoostVector for bosst vector of the beam.
Definition: B2BIIConvertMdstModule.h:422
Belle2::DBStore::addConstantOverride
void addConstantOverride(const std::string &name, TObject *obj, bool oneRun=false)
Add constant override payload.
Definition: DBStore.cc:204
prepareAsicCrosstalkSimDB.e
e
aux.
Definition: prepareAsicCrosstalkSimDB.py:53
Belle2::B2BIIConvertMdstModule::convertHelix
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.
Definition: B2BIIConvertMdstModule.cc:1592
Belle2::B2BIIConvertMdstModule::m_belleTrkExtra
StoreArray< BelleTrkExtra > m_belleTrkExtra
Belle CDC extra information.
Definition: B2BIIConvertMdstModule.h:408
Belle2::B2BIIConvertMdstModule::mdstKlongToParticle
std::map< int, int > mdstKlongToParticle
map of Klong Panther IDs and corresponding Particle StoreArray indices
Definition: B2BIIConvertMdstModule.h:384
Belle2::B2BIIConvertMdstModule::convertEvtclsTable
void convertEvtclsTable()
Reads and converts all entries of evtcls Panther table.
Definition: B2BIIConvertMdstModule.cc:1204
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:116
Belle2::B2BIIConvertMdstModule::m_particles
StoreArray< Particle > m_particles
Particles.
Definition: B2BIIConvertMdstModule.h:405
Belle2::ECLCluster
ECL cluster data.
Definition: ECLCluster.h:39
Belle2::Const::electron
static const ChargedStable electron
electron particle
Definition: Const.h:533
REG_MODULE
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:652
Belle2::B2BIIConvertMdstModule::convertMdstPi0Table
void convertMdstPi0Table()
Reads all entries of Mdst_Pi0 Panther table, creates a particle list 'pi0:mdst' and adds them to Stor...
Definition: B2BIIConvertMdstModule.cc:1071
Belle2::B2BIIConvertMdstModule::setTracksToECLClustersRelations
void setTracksToECLClustersRelations()
Sets Track -> ECLCluster relations.
Definition: B2BIIConvertMdstModule.cc:1849
Belle2::B2BIIConvertMdstModule::terminate
virtual void terminate() override
Terminates the module.
Definition: B2BIIConvertMdstModule.cc:2315
Belle2::TrackFitResult::getOmega
double getOmega() const
Getter for omega.
Definition: TrackFitResult.h:194
Belle2::B2BIIConvertMdstModule::event
virtual void event() override
Called for each event.
Definition: B2BIIConvertMdstModule.cc:288
Belle2::B2BIIConvertMdstModule::setKLMClustersRelations
void setKLMClustersRelations()
Sets KLMCluster -> Track and ECLCluster relations.
Definition: B2BIIConvertMdstModule.cc:1892
Belle2::B2BIIConvertMdstModule::m_collisionInvM
CollisionInvariantMass m_collisionInvM
CollisionInvariantMass for the invariant mass of the beam.
Definition: B2BIIConvertMdstModule.h:426
Belle2::ECLCluster::setNumberOfCrystals
void setNumberOfCrystals(double noc)
Set number of crystals (sum of weights).
Definition: ECLCluster.h:220
Belle2::Particle::addExtraInfo
void addExtraInfo(const std::string &name, float value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1224
Belle2::Database::storeData
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
Definition: Database.cc:152
Belle2::B2BIIConvertMdstModule::m_convertRecTrg
bool m_convertRecTrg
Flag to switch on conversion of rectrg_summary3.
Definition: B2BIIConvertMdstModule.h:160
Belle2::B2BIIConvertMdstModule::acc_pid
static double acc_pid(const Belle::Mdst_charged &chg, int idp)
Returns ACC likelihood for given hypothesis idp.
Belle2::B2BIIConvertMdstModule::m_lastIPProfileBin
int m_lastIPProfileBin
variable to tell us which IPProfile bin was active last time we looked
Definition: B2BIIConvertMdstModule.h:172
Belle2::MCParticle::setDecayVertex
void setDecayVertex(const TVector3 &vertex)
Set decay vertex.
Definition: MCParticle.h:452
Belle2::B2BIIConvertMdstModule::belleHelixToCartesian
int belleHelixToCartesian(const Belle::Mdst_trk_fit &trk_fit, const double mass, const HepPoint3D &newPivot, CLHEP::HepLorentzVector &momentum, HepPoint3D &position, CLHEP::HepSymMatrix &error, double dPhi=0.0)
Fills 4-momentum, position and 7x7 error matrix from Belle Helix stored in Mdst_trk_fit.
Definition: B2BIIConvertMdstModule.cc:2181
Belle2::RelationsInterface::addRelationTo
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).
Definition: RelationsObject.h:144
Belle2::B2BIIConvertMdstModule::m_mcParticles
StoreArray< MCParticle > m_mcParticles
MC particles.
Definition: B2BIIConvertMdstModule.h:387
Belle2::MCParticle::setStatus
void setStatus(unsigned short int status)
Set Status code for the particle.
Definition: MCParticle.h:354
Belle2::RelationsInterface::getRelated
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
Definition: RelationsObject.h:280
Belle2::PIDLikelihood
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Definition: PIDLikelihood.h:37
Belle2::B2BIIConvertMdstModule::belleVeeDaughterToCartesian
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.
Definition: B2BIIConvertMdstModule.cc:2013
Belle2::TrackFitResult::getTanLambda
double getTanLambda() const
Getter for tanLambda.
Definition: TrackFitResult.h:206
Belle2::B2BIIConvertMdstModule::initializeDataStore
void initializeDataStore()
Initializes Belle II DataStore.
Definition: B2BIIConvertMdstModule.cc:208
Belle2::Particle::updateMomentum
void updateMomentum(const TLorentzVector &p4, const TVector3 &vertex, const TMatrixFSym &errMatrix, float pValue)
Sets Lorentz vector, position, 7x7 error matrix and p-value.
Definition: Particle.h:328
Belle2::B2BIIConvertMdstModule::convertPIDData
void convertPIDData(const Belle::Mdst_charged &belleTrack, const Track *track)
Get PID information for belleTrack and add it to PIDLikelihood (with relation from track).
Definition: B2BIIConvertMdstModule.cc:1352
Belle2::ECLCluster::setEnergy
void setEnergy(double energy)
Set Corrected Energy (GeV).
Definition: ECLCluster.h:238
Belle2::B2BIIConvertMdstModule::cdc_pid
double cdc_pid(const Belle::Mdst_charged &chg, int idp)
Returns CDC likelihood for given hypothesis idp.
Belle2::Const::EDetector
EDetector
Enum for identifying the detector components (detector and subdetector).
Definition: Const.h:44
Belle2::B2BIIConvertMdstModule::setLikelihoods
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...
Definition: B2BIIConvertMdstModule.cc:1334
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::B2BIIConvertMdstModule::KAPPA2OMEGA
const double KAPPA2OMEGA
Conversion factor for Kappa -> Omega helix parameters.
Definition: B2BIIConvertMdstModule.h:443
Belle2::BeamSpot
This class contains the beam spot position and size modeled as a gaussian distribution in space.
Definition: BeamSpot.h:32
Belle2::B2BIIConvertMdstModule::m_collisionBoostVectorDB
OptionalDBObjPtr< CollisionBoostVector > m_collisionBoostVectorDB
CollisionBoostVector for boost vector.
Definition: B2BIIConvertMdstModule.h:421
Belle2::Const::electronMass
static const double electronMass
electron mass
Definition: Const.h:558
Belle2::TrackFitResult::getPosition
TVector3 getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:109
Belle2::B2BIIConvertMdstModule::convertMdstChargedTable
void convertMdstChargedTable()
Reads and converts all entries of Mdst_charged (Mdst_trk and Mdst_trk_fit) Panther table to Track (Tr...
Definition: B2BIIConvertMdstModule.cc:444
Belle2::B2BIIConvertMdstModule::c_GeneratorLevel
@ c_GeneratorLevel
Match to generator-level particles.
Definition: B2BIIConvertMdstModule.h:108
Belle2::B2BIIConvertMdstModule::convertGenHepevtObject
void convertGenHepevtObject(const Belle::Gen_hepevt &genHepevt, MCParticleGraph::GraphParticle *mcParticle)
Converts Gen_hepevt record to MCParticleGraph object.
Definition: B2BIIConvertMdstModule.cc:1767
Belle2::B2BIIConvertMdstModule::m_beamSpot
BeamSpot m_beamSpot
Interaction Point of the beam.
Definition: B2BIIConvertMdstModule.h:418
Belle2::Const::kaon
static const ChargedStable kaon
charged kaon particle
Definition: Const.h:536
Belle2::B2BIIConvertMdstModule::convertBeamEnergy
void convertBeamEnergy()
Stores beam parameters (energy, angles) in CollisionInvariantMass and CollisionBoostVector (currently...
Definition: B2BIIConvertMdstModule.cc:374
Belle2::Module
Base class for Modules.
Definition: Module.h:74
Belle2::B2BIIConvertMdstModule::convertMdstKLongTable
void convertMdstKLongTable()
Reads all entries of Mdst_Klong Panther table, creates a particle list 'K_L0:mdst' and adds them to S...
Definition: B2BIIConvertMdstModule.cc:1119
Belle2::B2BIIConvertMdstModule::BFIELD
const double BFIELD
B filed in TESLA.
Definition: B2BIIConvertMdstModule.h:440
Belle2::B2BIIConvertMdstModule::mdstKlmToKLMCluster
std::map< int, int > mdstKlmToKLMCluster
map of Mdst_klm Panther IDs and corresponding KLMCluster StoreArray indices
Definition: B2BIIConvertMdstModule.h:382
Belle2::B2BIIConvertMdstModule::m_klmClusters
StoreArray< KLMCluster > m_klmClusters
KLM clusters.
Definition: B2BIIConvertMdstModule.h:393
Belle2::B2BIIConvertMdstModule::m_v0s
StoreArray< V0 > m_v0s
V0-particles.
Definition: B2BIIConvertMdstModule.h:402
Belle2::Const::RestrictedDetectorSet
A class for sets of detector IDs whose content is limited to restricted set of valid detector IDs.
Definition: Const.h:172
Belle2::Const::speedOfLight
static const double speedOfLight
[cm/ns]
Definition: Const.h:568
Belle2::B2BIIConvertMdstModule::m_matchType2E9oE25Threshold
double m_matchType2E9oE25Threshold
E9/E25 threshold value clusters with a value above this threshold are classified as neutral even if t...
Definition: B2BIIConvertMdstModule.h:169
Belle2::Const::pion
static const ChargedStable pion
charged pion particle
Definition: Const.h:535
Belle2::B2BIIConvertMdstModule::m_trackFitResults
StoreArray< TrackFitResult > m_trackFitResults
Track fir results.
Definition: B2BIIConvertMdstModule.h:399
Belle2::MCParticle::setProductionTime
void setProductionTime(float time)
Set production time.
Definition: MCParticle.h:392
Belle2::B2BIIConvertMdstModule::c_nHyp
const static int c_nHyp
Number of Belle track hypotheses (see c_belleHyp_to_chargedStable).
Definition: B2BIIConvertMdstModule.h:327
Belle2::B2BIIConvertMdstModule::convertRecTrgTable
void convertRecTrgTable()
Reads and converts m_final from rectrg_summary3.
Definition: B2BIIConvertMdstModule.cc:1249
Belle2::B2BIIConvertMdstModule::~B2BIIConvertMdstModule
virtual ~B2BIIConvertMdstModule() override
Destructor.
Definition: B2BIIConvertMdstModule.cc:189
Belle2::ECLCluster::setPhi
void setPhi(double phi)
Set Phi of Shower (radian).
Definition: ECLCluster.h:232
Belle2::ECLCluster::setLAT
void setLAT(double LAT)
Set Lateral distribution parameter.
Definition: ECLCluster.h:217
Belle2::B2BIIConvertMdstModule::convertGenHepEvtTable
void convertGenHepEvtTable()
Reads and converts all entries of Gen_hepevt Panther table to MCParticle dataobjects and adds them to...
Definition: B2BIIConvertMdstModule.cc:848
Belle2::KLMCluster::setClusterPosition
void setClusterPosition(float globalX, float globalY, float globalZ)
Set global position.
Definition: KLMCluster.h:171
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::StoreObjPtr
Type-safe access to single objects in the data store.
Definition: ParticleList.h:33
Belle2::BeamSpot::setIP
void setIP(const TVector3 &ipPosition, const TMatrixDSym &covariance)
Set the IP position and its error matrix.
Definition: BeamSpot.h:69
Belle2::MCParticle::setPDG
void setPDG(int pdg)
Set PDG code of the particle.
Definition: MCParticle.h:343
Belle2::Const::deuteron
static const ChargedStable deuteron
deuteron particle
Definition: Const.h:538
Belle2::ECLCluster::setdeltaL
void setdeltaL(double deltaL)
Set deltaL for shower shape.
Definition: ECLCluster.h:183
Belle2::Particle::appendDaughter
void appendDaughter(const Particle *daughter, const bool updateType=true)
Appends index of daughter to daughters index array.
Definition: Particle.cc:633
Belle2::HitPatternCDC::setNHits
void setNHits(unsigned short nHits)
Sets the 8 MSBs to the total number of hits in the CDC.
Definition: HitPatternCDC.cc:55
Belle2::B2BIIConvertMdstModule::atcPID
double atcPID(const PIDLikelihood *pid, int sigHyp, int bkgHyp)
calculates atc_pid(3,1,5,sigHyp,bkgHyp).prob() from converted PIDLikelihood
Definition: B2BIIConvertMdstModule.cc:2265
Belle2::B2BIIConvertMdstModule::m_mcMatchingModeString
std::string m_mcMatchingModeString
MC matching mode.
Definition: B2BIIConvertMdstModule.h:151
Belle2::ECLCluster::setEnergyRaw
void setEnergyRaw(double energyraw)
Set Uncorrect Energy deposited (GeV).
Definition: ECLCluster.h:241
Belle2::HitPatternVXD::setSVDLayer
void setSVDLayer(const unsigned short layerId, unsigned short uHits, unsigned short vHits)
Set the number of hits in a specific layer of the SVD.
Definition: HitPatternVXD.cc:51
Belle2::B2BIIConvertMdstModule
Module converts Belle MDST objects (Panther records) to Belle II MDST objects.
Definition: B2BIIConvertMdstModule.h:97
Belle2::B2BIIConvertMdstModule::m_tracks
StoreArray< Track > m_tracks
Tracks.
Definition: B2BIIConvertMdstModule.h:396
Belle2::MCParticle::getPDG
int getPDG() const
Return PDG code of particle.
Definition: MCParticle.h:123
Belle2::MCParticleGraph::addParticle
GraphParticle & addParticle()
Add new particle to the graph.
Definition: MCParticleGraph.h:305
Belle2::B2BIIConvertMdstModule::convertMdstGammaTable
void convertMdstGammaTable()
Reads all entries of Mdst_Gamma Panther table, creates a particle list 'gamma:mdst' and adds them to ...
Definition: B2BIIConvertMdstModule.cc:1024
Belle2::UncertainHelix::getCovariance
const TMatrixDSym & getCovariance() const
Getter for covariance matrix of perigee parameters in matrix form.
Definition: UncertainHelix.h:104
Belle2::DBStore::Instance
static DBStore & Instance()
Instance of a singleton DBStore.
Definition: DBStore.cc:36
Belle2::CollisionInvariantMass
This class contains the measured average center-of-mass energy, which is equal to the invariant mass ...
Definition: CollisionInvariantMass.h:31
Belle2::RelationArray::add
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
Definition: RelationArray.h:291
Belle2::MCParticle::setProductionVertex
void setProductionVertex(const TVector3 &vertex)
Set production vertex position.
Definition: MCParticle.h:404
Belle2::TrackFitResult::getPhi0
double getPhi0() const
Getter for phi0.
Definition: TrackFitResult.h:184
Belle2::Particle::setVertex
void setVertex(const TVector3 &vertex)
Sets position (decay vertex)
Definition: Particle.h:291
Belle2::B2BIIConvertMdstModule::m_convertBeamParameters
bool m_convertBeamParameters
Convert beam parameters or use information stored in Belle II database.
Definition: B2BIIConvertMdstModule.h:143
Belle2::B2BIIConvertMdstModule::convertMdstChargedObject
void convertMdstChargedObject(const Belle::Mdst_charged &belleTrack, Track *track)
Converts Mdst_charged (Mdst_trk(_fit)) record to Track (TrackFitResult) object.
Definition: B2BIIConvertMdstModule.cc:1632
Belle2::B2BIIConvertMdstModule::m_beamSpotDB
OptionalDBObjPtr< BeamSpot > m_beamSpotDB
BeamSpot for IP.
Definition: B2BIIConvertMdstModule.h:417
Belle2::B2BIIConvertMdstModule::convertMdstECLObject
void convertMdstECLObject(const Belle::Mdst_ecl &ecl, const Belle::Mdst_ecl_aux &eclAux, ECLCluster *eclCluster)
Converts Mdst_ecl(_aux) record to ECLCluster object.
Definition: B2BIIConvertMdstModule.cc:1807
Belle2::RelationsInterface::getArrayIndex
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Definition: RelationsObject.h:387
Belle2::KLMCluster
KLM cluster data.
Definition: KLMCluster.h:38
Belle2::ECLCluster::setR
void setR(double r)
Set R (in cm).
Definition: ECLCluster.h:235
Belle2::B2BIIConvertMdstModule::m_realData
bool m_realData
flag that tells whether given data sample is for real data or MC
Definition: B2BIIConvertMdstModule.h:140
Belle2::B2BIIConvertMdstModule::convertMdstKLMTable
void convertMdstKLMTable()
Reads and converts all entries of Mdst_klm_cluster Panther table to KLMCluster dataobjects and adds t...
Definition: B2BIIConvertMdstModule.cc:998
Belle2::ECLCluster::setE9oE21
void setE9oE21(double E9oE21)
Set E9/E21 energy ratio.
Definition: ECLCluster.h:201
Belle2::B2BIIConvertMdstModule::getHelixParameters
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,...
Definition: B2BIIConvertMdstModule.cc:1480
Belle2::B2BIIConvertMdstModule::belleHelixToHelix
void belleHelixToHelix(const Belle::Mdst_trk_fit &trk_fit, std::vector< float > &helixParam, std::vector< float > &helixError)
obtains the helix parameters from trk_fit and moves the pivot to 0,0,0
Definition: B2BIIConvertMdstModule.cc:2113
Belle2::B2BIIConvertMdstModule::m_use6x6CovarianceMatrix4Tracks
bool m_use6x6CovarianceMatrix4Tracks
flag that tells which form of covariance matrix should be used in the conversion of charged tracks
Definition: B2BIIConvertMdstModule.h:148
Belle2::B2BIIConvertMdstModule::m_pidLikelihoods
StoreArray< PIDLikelihood > m_pidLikelihoods
output PIDLikelihood array.
Definition: B2BIIConvertMdstModule.h:411
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::KLMCluster::setLayers
void setLayers(int layers)
Set number of layers with hits.
Definition: KLMCluster.h:155
Belle2::Unit::mm
static const double mm
[millimeters]
Definition: Unit.h:80
Belle2::TrackFitResult::getCovariance5
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
Definition: TrackFitResult.cc:106
Belle2::MCParticle::get4Vector
TLorentzVector get4Vector() const
Return 4Vector of particle.
Definition: MCParticle.h:218
Belle2::Const::ParticleType
The ParticleType class for identifying different particle types.
Definition: Const.h:284
Belle2::Const::proton
static const ChargedStable proton
proton particle
Definition: Const.h:537
Belle2::Database::Instance
static Database & Instance()
Instance of a singleton Database.
Definition: Database.cc:54
Belle2::ECLCluster::setEnergyHighestCrystal
void setEnergyHighestCrystal(double energyhighestcrystal)
Set energy of highest energetic crystal (GeV).
Definition: ECLCluster.h:244
Belle2::B2BIIConvertMdstModule::genHepevtToMCParticle
std::map< int, int > genHepevtToMCParticle
map of Gen_hepevt Panther IDs and corresponding MCParticle StoreArray indices
Definition: B2BIIConvertMdstModule.h:376
Belle2::HitPatternCDC::getInteger
ULong64_t getInteger() const
Getter for underlying integer type.
Definition: HitPatternCDC.cc:71
Belle2::Const::ChargedStable
Provides a type-safe way to pass members of the chargedStableSet set.
Definition: Const.h:465
Belle2::MCParticle::setMass
void setMass(float mass)
Set particle mass.
Definition: MCParticle.h:374
Belle2::MCParticle::set4Vector
void set4Vector(const TLorentzVector &p4)
Sets the 4Vector of particle.
Definition: MCParticle.h:446
Belle2::B2BIIConvertMdstModule::m_particleGraph
Belle2::MCParticleGraph m_particleGraph
MCParticle Graph to build Belle2 MC Particles.
Definition: B2BIIConvertMdstModule.h:373
Belle2::B2BIIConvertMdstModule::mdstEclToECLCluster
std::map< int, int > mdstEclToECLCluster
map of Mdst_ecl Panther IDs and corresponding ECLCluster StoreArray indices
Definition: B2BIIConvertMdstModule.h:378
Belle2::Const::muon
static const ChargedStable muon
muon particle
Definition: Const.h:534
Belle2::ProcHandler::parallelProcessingUsed
static bool parallelProcessingUsed()
Returns true if multiple processes have been spawned, false in single-core mode.
Definition: ProcHandler.cc:221
Belle2::Track
Class that bundles various TrackFitResults.
Definition: Track.h:35
HepGeom::Point3D< double >
Belle2::MCParticle
A Class to store the Monte Carlo particle information.
Definition: MCParticle.h:43
Belle2::B2BIIConvertMdstModule::convertMdstKLMObject
void convertMdstKLMObject(const Belle::Mdst_klm_cluster &klm, KLMCluster *klmCluster)
Converts Mdst_klm_cluster record to KLMCluster object.
Definition: B2BIIConvertMdstModule.cc:1838
Belle2::MCParticleGraph::GraphParticle::decaysInto
void decaysInto(GraphParticle &daughter)
Tells the graph that this particle decays into daughter.
Definition: MCParticleGraph.h:110
Belle2::B2BIIConvertMdstModule::c_Direct
@ c_Direct
Direct matching.
Definition: B2BIIConvertMdstModule.h:105
Belle2::TrackFitResult::getHelix
Helix getHelix() const
Conversion to framework Helix (without covariance).
Definition: TrackFitResult.h:230
Belle2::KLMCluster::setInnermostLayer
void setInnermostLayer(int innermostLayer)
Set number of the innermost layer with hits.
Definition: KLMCluster.h:162
Belle2::MCParticle::setDecayTime
void setDecayTime(float time)
Set decay time.
Definition: MCParticle.h:398
Belle2::MCParticleGraph::clear
void clear()
Reset particles and decay information to make the class reusable.
Definition: MCParticleGraph.h:298
Belle2::CollisionBoostVector
This class contains the measured average boost vector vec(beta) = (beta_x, beta_y,...
Definition: CollisionBoostVector.h:33
Belle2::B2BIIConvertMdstModule::recoverMoreThan24bitIDHEP
int recoverMoreThan24bitIDHEP(int id)
Helper function to recover falsely set idhep info in GenHepEvt list.
Definition: B2BIIConvertMdstModule.cc:1918
Belle2::B2BIIConvertMdstModule::c_belleHyp_to_chargedStable
const static Const::ChargedStable c_belleHyp_to_chargedStable[c_nHyp]
maps Belle hypotheses to Const::ChargedStable (from http://belle.kek.jp/secured/wiki/doku....
Definition: B2BIIConvertMdstModule.h:329
Belle2::HitPatternVXD::getInteger
unsigned int getInteger() const
Getter for the underlying integer.
Definition: HitPatternVXD.h:65
Belle2::CollisionInvariantMass::setMass
void setMass(double mass, double error, double spread)
Set the CMS energy and its uncertainty.
Definition: CollisionInvariantMass.h:45
Belle2::MCParticle::setValidVertex
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
Definition: MCParticle.h:386
Belle2::B2BIIConvertMdstModule::m_mcMatchingMode
MCMatchingMode m_mcMatchingMode
C matching mode.
Definition: B2BIIConvertMdstModule.h:154
Belle2::B2BIIConvertMdstModule::ERRMCONV
const int ERRMCONV[7]
CONVERSION OF TRACK ERROR MATRIX ELEMENTS.
Definition: B2BIIConvertMdstModule.h:437
Belle2::MCParticleGraph::GraphParticle
Class to represent Particle data in graph.
Definition: MCParticleGraph.h:86
Belle2::HitPatternVXD
Hit pattern of the VXD within a track.
Definition: HitPatternVXD.h:44
Belle2::TrackFitResult::getChargeSign
short getChargeSign() const
Return track charge (1 or -1).
Definition: TrackFitResult.h:160
Belle2::MCParticle::c_PrimaryParticle
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition: MCParticle.h:58
Belle2::Particle::setPValue
void setPValue(float pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:308
Belle2::ECLCluster::setCovarianceMatrix
void setCovarianceMatrix(double covArray[6])
Set covariance matrix (3x3), i.e.
Definition: ECLCluster.h:163
Belle2::HitPatternCDC
Hit pattern of CDC hits within a track.
Definition: HitPatternCDC.h:45
Belle2::ECLCluster::setTheta
void setTheta(double theta)
Set Theta of Shower (radian).
Definition: ECLCluster.h:229
Belle2::B2BIIConvertMdstModule::belleVeeDaughterHelix
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
Definition: B2BIIConvertMdstModule.cc:2052
Belle2::TrackFitResult::getD0
double getD0() const
Getter for d0.
Definition: TrackFitResult.h:178
Belle2::B2BIIConvertMdstModule::endRun
virtual void endRun() override
Called when the current run is finished.
Definition: B2BIIConvertMdstModule.cc:2309
Belle2::B2BIIConvertMdstModule::initialize
virtual void initialize() override
Initialize the module.
Definition: B2BIIConvertMdstModule.cc:193
Belle2::B2BIIConvertMdstModule::m_collisionInvMDB
OptionalDBObjPtr< CollisionInvariantMass > m_collisionInvMDB
CollisionInvariantMass for Invariant Mass of Beam.
Definition: B2BIIConvertMdstModule.h:425
Belle2::B2BIIConvertMdstModule::m_evtInfo
StoreObjPtr< EventExtraInfo > m_evtInfo
Event Extra Info.
Definition: B2BIIConvertMdstModule.h:414
Belle2::Const::ParticleType::getMass
double getMass() const
Particle mass.
Definition: UnitConst.cc:310