Belle II Software  release-05-01-25
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 
1257  // Pull rectrg_summary3 from manager
1258  Belle::Rectrg_summary3_Manager& RecTrgSummary3Mgr = Belle::Rectrg_summary3_Manager::get_manager();
1259 
1260  std::string name = "rectrg_summary3_m_final";
1261  // Only one entry in each event
1262  std::vector<Belle::Rectrg_summary3>::iterator eflagIterator = RecTrgSummary3Mgr.begin();
1263 
1264  // Converting m_final(3)
1265  for (int index = 0; index < 3; ++index) {
1266  std::string iVar = name + std::to_string(index);
1267  m_evtInfo->addExtraInfo(iVar, (*eflagIterator).final(index));
1268  B2DEBUG(99, "m_final(" << index << ") = " << m_evtInfo->getExtraInfo(iVar));
1269  }
1270 
1271 }
1272 
1273 
1274 //-----------------------------------------------------------------------------
1275 // CONVERT OBJECTS
1276 //-----------------------------------------------------------------------------
1277 
1278 #ifdef HAVE_KID_ACC
1279 double B2BIIConvertMdstModule::acc_pid(const Belle::Mdst_charged& chg, int idp)
1280 {
1281  static Belle::kid_acc acc_pdf(0);
1282  //static kid_acc acc_pdf(1);
1283 
1284  const double pmass[5] = { 0.00051099907, 0.105658389, 0.13956995, 0.493677, 0.93827231 };
1285 
1286  CLHEP::Hep3Vector mom(chg.px(), chg.py(), chg.pz());
1287  double cos_theta = mom.cosTheta();
1288  double pval = mom.mag();
1289 
1290  double npe = chg.acc().photo_electron();
1291  double beta = pval / sqrt(pval * pval + pmass[idp] * pmass[idp]);
1292  double pdfval = acc_pdf.npe2pdf(cos_theta, beta, npe);
1293 
1294  return pdfval;
1295 }
1296 
1297 // this is CDC_prob5
1298 double B2BIIConvertMdstModule::cdc_pid(const Belle::Mdst_charged& chg, int idp)
1299 {
1300  CLHEP::Hep3Vector mom(chg.px(), chg.py(), chg.pz());
1301  double pval = mom.mag();
1302 
1303  Belle::kid_cdc kidCdc(5);
1304  float factor0 = kidCdc.factor0();
1305  float factor1 = kidCdc.factor1(idp, pval);
1306 
1307  if (factor0 == 1.0 && factor1 == 1.0) return chg.trk().pid(idp);
1308  //
1309  double m = chg.trk().dEdx() / factor0;
1310  double e = chg.trk().dEdx_exp(idp) * factor1;
1311  double s = chg.trk().sigma_dEdx(idp);
1312  double val = 1. / sqrt(2.*M_PI) / s * exp(-0.5 * (m - e) * (m - e) / s / s);
1313 
1314  return val;
1315 }
1316 #endif
1317 
1318 void B2BIIConvertMdstModule::setLikelihoods(PIDLikelihood* pid, Const::EDetector det, double likelihoods[c_nHyp],
1319  bool discard_allzero)
1320 {
1321  if (discard_allzero) {
1322  const double max_l = *std::max_element(likelihoods, likelihoods + c_nHyp);
1323  if (max_l <= 0.0) {
1324  return; //likelihoods broken, ignore
1325  }
1326  }
1327 
1328  for (int i = 0; i < c_nHyp; i++) {
1329  float logl = log(likelihoods[i]);
1330  pid->setLogLikelihood(det, c_belleHyp_to_chargedStable[i], logl);
1331  }
1332  //copy proton likelihood to deuterons
1333  pid->setLogLikelihood(det, Const::deuteron, pid->getLogL(Const::proton, det));
1334 }
1335 
1336 void B2BIIConvertMdstModule::convertPIDData(const Belle::Mdst_charged& belleTrack, const Track* track)
1337 {
1338  PIDLikelihood* pid = m_pidLikelihoods.appendNew();
1339  track->addRelationTo(pid);
1340 
1341  //convert data handled by atc_pid: dE/dx (-> CDC), TOF (-> TOP), ACC ( -> ARICH)
1342  //this should result in the same likelihoods used when creating atc_pid(3, 1, 5, ..., ...)
1343  //and calling prob(const Mdst_charged & chg).
1344 
1345  double likelihoods[c_nHyp];
1346  double accL[c_nHyp];
1347  double tofL[c_nHyp];
1348  double cdcL[c_nHyp];
1349  for (int i = 0; i < c_nHyp; i++) {
1350  accL[i] = tofL[i] = cdcL[i] = 1.0;
1351  }
1352 #ifdef HAVE_KID_ACC
1353  //accq0 = 3, as implemented in acc_prob3()
1354  const auto& acc = belleTrack.acc();
1355  if (acc and acc.quality() == 0) {
1356  for (int i = 0; i < c_nHyp; i++)
1357  accL[i] = likelihoods[i] = acc_pid(belleTrack, i);
1358  setLikelihoods(pid, Const::ARICH, likelihoods, true);
1359  }
1360 #endif
1361 
1362  //tofq0 = 1, as implemented in tof_prob1()
1363  //uses p1 / (p1 + p2) to create probability, so this should map directly to likelihoods
1364  const Belle::Mdst_tof& tof = belleTrack.tof();
1365  if (tof and tof.quality() == 0) {
1366  for (int i = 0; i < c_nHyp; i++)
1367  tofL[i] = likelihoods[i] = tof.pid(i);
1368  setLikelihoods(pid, Const::TOP, likelihoods, true);
1369  }
1370 
1371  // cdcq0 = 5, as implemented in cdc_prob0() (which is used for all values of cdcq0!)
1372  //uses p1 / (p1 + p2) to create probability, so this should map directly to likelihoods
1373  // eID actually uses cdc_pid (cdc_prob5)
1374  const Belle::Mdst_trk& trk = belleTrack.trk();
1375  if (trk.dEdx() > 0) {
1376  for (int i = 0; i < c_nHyp; i++) {
1377  likelihoods[i] = trk.pid(i);
1378  cdcL[i] = cdc_pid(belleTrack, i);
1379  }
1380  setLikelihoods(pid, Const::CDC, likelihoods, true);
1381  }
1382 
1383 
1384  // eid
1385  // eid is combination of atc_pid and ecl related information
1386  // since atc_pid part is already converted above only the ECL part
1387  // is converted
1388  // ECL pdfs are available only for electrons and hadrons (assumed to be pions)
1389  // likelihoods for others are set to 0
1390 
1391 #ifdef HAVE_EID
1392  Belle::eid electronID(belleTrack);
1393  float eclID_e_pdf = electronID.pdf_e_ecl();
1394  float eclID_h_pdf = electronID.pdf_h_ecl();
1395  float atcID_e_pdf = electronID.atc_pid_pdf(true, accL, tofL, cdcL);
1396  float atcID_h_pdf = electronID.atc_pid_pdf(false, accL, tofL, cdcL);
1397 
1398  // eID
1399  float eclProb = eclID_e_pdf / (eclID_e_pdf + eclID_h_pdf);
1400  float atcProb = atcID_e_pdf / (atcID_e_pdf + atcID_h_pdf);
1401 
1402  if (atcProb > 0.999999) atcProb = 0.999999;
1403  // combine the two probabilities.
1404  double eidCombinedSig = eclProb * atcProb;
1405  double eidCombinedBkg = (1. - eclProb) * (1. - atcProb);
1406 
1407  likelihoods[0] = eidCombinedSig;
1408  likelihoods[1] = 0; // no muons
1409  likelihoods[2] = eidCombinedBkg;
1410  likelihoods[3] = 0; // no kaons
1411  likelihoods[4] = 0; // no protons
1412 
1413  setLikelihoods(pid, Const::ECL, likelihoods, true);
1414 
1415  //Hep3Vector mom(belleTrack.px(), belleTrack.py(), belleTrack.pz());
1416  //B2INFO(" p = " << mom.mag() << " le_ecl = " << electronID.le_ecl());
1417 #endif
1418 
1419  //muid
1420  //Note that though it says "_likelihood()" on the label, those are
1421  //actually likelihood ratios of the type L(hyp) / (L(mu) + L(pi) + L(K)),
1422  //which are set in the FixMdst module.
1423  int muid_trackid = belleTrack.muid_ID();
1424  if (muid_trackid) {
1425  //Using approach 2. from http://belle.kek.jp/secured/muid/usage_muid.html since
1426  //it's much simpler than what Muid_mdst does.
1427  Belle::Mdst_klm_mu_ex_Manager& ex_mgr = Belle::Mdst_klm_mu_ex_Manager::get_manager();
1428  Belle::Mdst_klm_mu_ex& ex = ex_mgr(Belle::Panther_ID(muid_trackid));
1429 
1430  //filter out tracks with insufficient #hits (equal to cut on Muid_mdst::Chi_2())
1431  if (ex.Chi_2() > 0) {
1432  likelihoods[0] = 0; //no electrons
1433  likelihoods[1] = ex.Muon_likelihood();
1434  likelihoods[2] = ex.Pion_likelihood();
1435  likelihoods[3] = ex.Kaon_likelihood();
1436  likelihoods[4] = 0; //no protons
1437  //Miss_likelihood should only be != 0 for tracks that do not pass the Chi_2 cut.
1438 
1439  // in some cases the ex.XYZ_likelihood() < 0; Set it to 0 in these cases.
1440  for (int i = 0; i < 5; i++)
1441  if (likelihoods[i] < 0)
1442  likelihoods[i] = 0;
1443 
1444  //note: discard_allzero = false since all likelihoods = 0 usually means that Junk_likelihood is 1
1445  // PIDLikelihood::getProbability(hyp) will correctly return 0 then.
1446  setLikelihoods(pid, Const::KLM, likelihoods);
1447 
1448  /*
1449  const double tolerance = 1e-7;
1450  if (fabs(pid->getProbability(Const::muon, nullptr, Const::KLM) - ex.Muon_likelihood()) > tolerance ||
1451  fabs(pid->getProbability(Const::pion, nullptr, Const::KLM) - ex.Pion_likelihood()) > tolerance ||
1452  fabs(pid->getProbability(Const::kaon, nullptr, Const::KLM) - ex.Kaon_likelihood()) > tolerance) {
1453 
1454  B2INFO("muons: " << pid->getProbability(Const::muon, nullptr, Const::KLM) << " " << ex.Muon_likelihood());
1455  B2INFO("pion: " << pid->getProbability(Const::pion, nullptr, Const::KLM) << " " << ex.Pion_likelihood());
1456  B2INFO("kaon: " << pid->getProbability(Const::kaon, nullptr, Const::KLM) << " " << ex.Kaon_likelihood());
1457  B2INFO("miss/junk: " << ex.Miss_likelihood() << " " << ex.Junk_likelihood());
1458  }
1459  */
1460  }
1461  }
1462 }
1463 
1464 int B2BIIConvertMdstModule::getHelixParameters(const Belle::Mdst_trk_fit& trk_fit,
1465  const double mass,
1466  const HepPoint3D& newPivot,
1467  std::vector<float>& helixParams,
1468  CLHEP::HepSymMatrix& error5x5,
1469  CLHEP::HepLorentzVector& momentum,
1470  HepPoint3D& position,
1471  CLHEP::HepSymMatrix& error7x7, const double dPhi)
1472 {
1473  const HepPoint3D pivot(trk_fit.pivot_x(),
1474  trk_fit.pivot_y(),
1475  trk_fit.pivot_z());
1476 
1477  CLHEP::HepVector a(5);
1478  a[0] = trk_fit.helix(0);
1479  a[1] = trk_fit.helix(1);
1480  a[2] = trk_fit.helix(2);
1481  a[3] = trk_fit.helix(3);
1482  a[4] = trk_fit.helix(4);
1483  CLHEP::HepSymMatrix Ea(5, 0);
1484  Ea[0][0] = trk_fit.error(0);
1485  Ea[1][0] = trk_fit.error(1);
1486  Ea[1][1] = trk_fit.error(2);
1487  Ea[2][0] = trk_fit.error(3);
1488  Ea[2][1] = trk_fit.error(4);
1489  Ea[2][2] = trk_fit.error(5);
1490  Ea[3][0] = trk_fit.error(6);
1491  Ea[3][1] = trk_fit.error(7);
1492  Ea[3][2] = trk_fit.error(8);
1493  Ea[3][3] = trk_fit.error(9);
1494  Ea[4][0] = trk_fit.error(10);
1495  Ea[4][1] = trk_fit.error(11);
1496  Ea[4][2] = trk_fit.error(12);
1497  Ea[4][3] = trk_fit.error(13);
1498  Ea[4][4] = trk_fit.error(14);
1499 
1500  Belle::Helix helix(pivot, a, Ea);
1501 
1502  int charge = 0;
1503  if (helix.kappa() > 0)
1504  charge = 1;
1505  else
1506  charge = -1;
1507 
1508  if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
1509  helix.pivot(newPivot);
1510  momentum = helix.momentum(dPhi, mass, position, error7x7);
1511  } else {
1512  if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
1513  helix.pivot(HepPoint3D(0., 0., 0.));
1514  momentum = helix.momentum(dPhi, mass, position, error7x7);
1515  } else {
1516  momentum = helix.momentum(dPhi, mass, position, error7x7);
1517  }
1518  }
1519 
1520  convertHelix(helix, helixParams, error5x5);
1521 
1522  return charge;
1523 }
1524 
1525 void B2BIIConvertMdstModule::convertHelix(const Belle::Mdst_trk_fit& trk_fit,
1526  const HepPoint3D& newPivot,
1527  std::vector<float>& helixParams, std::vector<float>& helixError)
1528 {
1529  const HepPoint3D pivot(trk_fit.pivot_x(),
1530  trk_fit.pivot_y(),
1531  trk_fit.pivot_z());
1532 
1533  CLHEP::HepVector a(5);
1534  a[0] = trk_fit.helix(0);
1535  a[1] = trk_fit.helix(1);
1536  a[2] = trk_fit.helix(2);
1537  a[3] = trk_fit.helix(3);
1538  a[4] = trk_fit.helix(4);
1539  CLHEP::HepSymMatrix Ea(5, 0);
1540  Ea[0][0] = trk_fit.error(0);
1541  Ea[1][0] = trk_fit.error(1);
1542  Ea[1][1] = trk_fit.error(2);
1543  Ea[2][0] = trk_fit.error(3);
1544  Ea[2][1] = trk_fit.error(4);
1545  Ea[2][2] = trk_fit.error(5);
1546  Ea[3][0] = trk_fit.error(6);
1547  Ea[3][1] = trk_fit.error(7);
1548  Ea[3][2] = trk_fit.error(8);
1549  Ea[3][3] = trk_fit.error(9);
1550  Ea[4][0] = trk_fit.error(10);
1551  Ea[4][1] = trk_fit.error(11);
1552  Ea[4][2] = trk_fit.error(12);
1553  Ea[4][3] = trk_fit.error(13);
1554  Ea[4][4] = trk_fit.error(14);
1555 
1556  Belle::Helix helix(pivot, a, Ea);
1557 
1558  if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
1559  helix.pivot(newPivot);
1560  } else {
1561  if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
1562  helix.pivot(HepPoint3D(0., 0., 0.));
1563  }
1564  }
1565 
1566  CLHEP::HepSymMatrix error5x5(5, 0);
1567  convertHelix(helix, helixParams, error5x5);
1568 
1569  unsigned int size = 5;
1570  unsigned int counter = 0;
1571  for (unsigned int i = 0; i < size; i++)
1572  for (unsigned int j = i; j < size; j++)
1573  helixError[counter++] = error5x5[i][j];
1574 }
1575 
1576 void B2BIIConvertMdstModule::convertHelix(Belle::Helix& helix, std::vector<float>& helixParams, CLHEP::HepSymMatrix& error5x5)
1577 {
1578  CLHEP::HepVector a(5);
1579  CLHEP::HepSymMatrix Ea(5, 0);
1580 
1581  a = helix.a();
1582  Ea = helix.Ea();
1583 
1584  // param 0: d_0 = d_rho
1585  helixParams[0] = a[0];
1586 
1587  // param 1: phi = phi_0 + pi/2
1588  helixParams[1] = adjustAngleRange(a[1] + TMath::Pi() / 2.0);
1589 
1590  // param 2: omega = Kappa * alpha = Kappa * B[Tesla] * speed_of_light[m/s] * 1e-11
1591  helixParams[2] = a[2] * KAPPA2OMEGA;
1592 
1593  // param 3: d_z = z0
1594  helixParams[3] = a[3];
1595 
1596  // param 4: tan(Lambda) = tanLambda
1597  helixParams[4] = a[4];
1598 
1599  unsigned int size = 5;
1600  for (unsigned int i = 0; i < size; i++) {
1601  for (unsigned int j = 0; j < size; j++) {
1602  error5x5[i][j] = Ea[i][j];
1603  if (i == 2)
1604  error5x5[i][j] *= KAPPA2OMEGA;
1605  if (j == 2)
1606  error5x5[i][j] *= KAPPA2OMEGA;
1607 
1608  if (std::isinf(error5x5[i][j])) {
1609  B2DEBUG(99, "Helix covariance matrix element found to be infinite. Setting value to DBL_MAX/2.0.");
1610  error5x5[i][j] = DBL_MAX / 2.0;
1611  }
1612  }
1613  }
1614 }
1615 
1616 void B2BIIConvertMdstModule::convertMdstChargedObject(const Belle::Mdst_charged& belleTrack, Track* track)
1617 {
1618  Belle::Mdst_trk& trk = belleTrack.trk();
1619 
1620  for (int mhyp = 0 ; mhyp < c_nHyp; ++mhyp) {
1621  const Const::ChargedStable& pType = c_belleHyp_to_chargedStable[mhyp];
1622  double thisMass = pType.getMass();
1623 
1624  Belle::Mdst_trk_fit& trk_fit = trk.mhyp(mhyp);
1625 
1626  // Converted helix parameters
1627  std::vector<float> helixParam(5);
1628  // Converted 5x5 error matrix
1629  CLHEP::HepSymMatrix error5x5(5, 0);
1630  // 4-momentum
1631  CLHEP::HepLorentzVector momentum;
1632  // 7x7 (momentum, position) error matrix
1633  CLHEP::HepSymMatrix error7x7(7, 0);
1634  // position
1635  HepPoint3D position;
1636 
1637  getHelixParameters(trk_fit, thisMass, HepPoint3D(0., 0., 0.),
1638  helixParam, error5x5,
1639  momentum, position, error7x7, 0.0);
1640 
1641  std::vector<float> helixError(15);
1642  unsigned int size = 5;
1643  unsigned int counter = 0;
1644  for (unsigned int i = 0; i < size; i++)
1645  for (unsigned int j = i; j < size; j++)
1646  helixError[counter++] = error5x5[i][j];
1647 
1648  double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
1649 
1650  // Create an empty cdc hitpattern and set the number of total hits
1651  // use hits from 0: axial-wire, 1:stereo-wire, 2:cathode
1652  // the actual cdc hitpattern is not converted
1653 
1654  int cdcNHits = 0;
1655  for (unsigned int i = 0; i < 3; i++)
1656  cdcNHits += trk_fit.nhits(i);
1657 
1658  HitPatternCDC patternCdc;
1659  patternCdc.setNHits(cdcNHits);
1660 
1661  // conversion of track position in CDC layers
1662  if (m_convertTrkExtra) {
1663  auto cdcExtraInfo = m_belleTrkExtra.appendNew(trk_fit.first_x(), trk_fit.first_y(), trk_fit.first_z(),
1664  trk_fit.last_x(), trk_fit.last_y(), trk_fit.last_z());
1665  track->addRelationTo(cdcExtraInfo);
1666  }
1667  // conversion of the SVD hit pattern
1668  int svdHitPattern = trk_fit.hit_svd();
1669  // use hits from 3: SVD-rphi, 4: SVD-z
1670  // int svdNHits = trk_fit.nhits(3) + trk_fit.nhits(4);
1671 
1672  std::bitset<32> svdBitSet(svdHitPattern);
1673 
1674  HitPatternVXD patternVxd;
1675 
1676  unsigned short svdLayers;
1677  // taken from: http://belle.kek.jp/group/indirectcp/cpfit/cpfit-festa/2004/talks/Apr.14/CPfesta-2005-Higuchi(3).pdf
1678  StoreObjPtr<EventMetaData> event;
1679  // mask for the rphi hits, first 6 (8) bits/ 2 bits per layer
1680  std::bitset<32> svdUMask(static_cast<std::string>("00000000000000000000000000000011"));
1681  // mask for the z hits, second 6 (8) bits/ 2 bits per layer
1682  std::bitset<32> svdVMask;
1683 
1684  // find out if the SVD has 3 (4) layers; if exp <= (>) exp 27
1685  if (event->getExperiment() <= 27) {
1686  svdVMask = svdUMask << 6;
1687  svdLayers = 3;
1688  } else {
1689  svdVMask = svdUMask << 8;
1690  svdLayers = 4;
1691  }
1692 
1693  // loop over all svd layers (layer index is shifted + 3 for basf2)
1694  for (unsigned short layerId = 0; layerId < svdLayers; layerId++) {
1695  unsigned short uHits = (svdBitSet & svdUMask).count();
1696  unsigned short vHits = (svdBitSet & svdVMask).count();
1697  patternVxd.setSVDLayer(layerId + 3, uHits, vHits);
1698  // shift masks to the left
1699  svdUMask <<= 2;
1700  svdVMask <<= 2;
1701  }
1702 
1703  TrackFitResult helixFromHelix(helixParam, helixError, pType, pValue, -1, patternVxd.getInteger(), 0);
1704 
1705  if (m_use6x6CovarianceMatrix4Tracks) {
1706  TMatrixDSym cartesianCovariance(6);
1707  for (unsigned i = 0; i < 7; i++) {
1708  if (i == 3)
1709  continue;
1710  for (unsigned j = 0; j < 7; j++) {
1711  if (j == 3)
1712  continue;
1713 
1714  cartesianCovariance(ERRMCONV[i], ERRMCONV[j]) = error7x7[i][j];
1715  }
1716  }
1717  UncertainHelix helixFromCartesian(helixFromHelix.getPosition(), helixFromHelix.getMomentum(), helixFromHelix.getChargeSign(),
1718  BFIELD, cartesianCovariance, pValue);
1719 
1720  TMatrixDSym helixCovariance = helixFromCartesian.getCovariance();
1721 
1722  counter = 0;
1723  for (unsigned int i = 0; i < 5; ++i)
1724  for (unsigned int j = i; j < 5; ++j)
1725  helixError[counter++] = helixCovariance(i, j);
1726  }
1727 
1728  auto trackFit = m_trackFitResults.appendNew(helixParam, helixError, pType, pValue, patternCdc.getInteger(),
1729  patternVxd.getInteger(), trk_fit.ndf());
1730  track->setTrackFitResultIndex(pType, trackFit->getArrayIndex());
1731  /*
1732  B2INFO("--- B1 Track: ");
1733  std::cout << "Momentum = " << momentum << std::endl;
1734  std::cout << "Position = " << position << std::endl;
1735  std::cout << "7x7 error matrix = " << error7x7 << std::endl;
1736  B2INFO("--- B2 Track: ");
1737  std::cout << "Momentum = " << std::endl;
1738  trackFit->get4Momentum().Print();
1739  std::cout << "Position = " << std::endl;
1740  trackFit->getPosition().Print();
1741  std::cout << "6x6 error matrix = " << std::endl;
1742  trackFit->getCovariance6().Print();
1743  TMatrixDSym b2Error7x7(7);
1744  fill7x7ErrorMatrix(trackFit, b2Error7x7, thisMass, 1.5);
1745  std::cout << "7x7 error matrix = " << std::endl;
1746  b2Error7x7.Print();
1747  */
1748  }
1749 }
1750 
1751 void B2BIIConvertMdstModule::convertGenHepevtObject(const Belle::Gen_hepevt& genHepevt, MCParticleGraph::GraphParticle* mcParticle)
1752 {
1753  //B2DEBUG(80, "Gen_ehepevt: idhep " << genHepevt.idhep() << " (" << genHepevt.isthep() << ") with ID = " << genHepevt.get_ID());
1754 
1755  // updating the GraphParticle information from the Gen_hepevt information
1756  const int idHep = recoverMoreThan24bitIDHEP(genHepevt.idhep());
1757  if (idHep == 0) {
1758  B2WARNING("Trying to convert Gen_hepevt with idhep = " << idHep <<
1759  ". This should never happen.");
1760  mcParticle->setPDG(22);
1761  } else {
1762  mcParticle->setPDG(idHep);
1763  }
1764 
1765  if (genHepevt.isthep() > 0) {
1766  mcParticle->setStatus(MCParticle::c_PrimaryParticle);
1767  }
1768 
1769  mcParticle->setMass(genHepevt.M());
1770 
1771  TLorentzVector p4(genHepevt.PX(), genHepevt.PY(), genHepevt.PZ(), genHepevt.E());
1772  mcParticle->set4Vector(p4);
1773 
1774  mcParticle->setProductionVertex(genHepevt.VX()*Unit::mm, genHepevt.VY()*Unit::mm, genHepevt.VZ()*Unit::mm);
1775  mcParticle->setProductionTime(genHepevt.T()*Unit::mm / Const::speedOfLight);
1776 
1777  // decay time of this particle is production time of the daughter particle
1778  if (genHepevt.daFirst() > 0) {
1779  Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
1780  Belle::Gen_hepevt daughterParticle = genMgr(Belle::Panther_ID(genHepevt.daFirst()));
1781  mcParticle->setDecayTime(daughterParticle.T()*Unit::mm / Const::speedOfLight);
1782  mcParticle->setDecayVertex(daughterParticle.VX()*Unit::mm, daughterParticle.VY()*Unit::mm, daughterParticle.VZ()*Unit::mm);
1783  } else {
1784  //otherwise, assume it's stable
1785  mcParticle->setDecayTime(std::numeric_limits<float>::infinity());
1786  }
1787 
1788  mcParticle->setValidVertex(true);
1789 }
1790 
1791 void B2BIIConvertMdstModule::convertMdstECLObject(const Belle::Mdst_ecl& ecl, const Belle::Mdst_ecl_aux& eclAux,
1792  ECLCluster* eclCluster)
1793 {
1794  if (eclAux.e9oe25() < m_matchType2E9oE25Threshold)
1795  eclCluster->setIsTrack(ecl.match() > 0);
1796  else
1797  eclCluster->setIsTrack(ecl.match() == 1);
1798 
1799  eclCluster->setEnergy(ecl.energy()); //must happen before setCovarianceMatrix()!
1800  eclCluster->setPhi(ecl.phi());
1801  eclCluster->setTheta(ecl.theta());
1802  eclCluster->setR(ecl.r());
1803  eclCluster->setdeltaL(ecl.quality());
1804 
1805  double covarianceMatrix[6];
1806  covarianceMatrix[0] = ecl.error(0); // error on energy
1807  covarianceMatrix[1] = ecl.error(1);
1808  covarianceMatrix[2] = ecl.error(2); // error on phi
1809  covarianceMatrix[3] = ecl.error(3);
1810  covarianceMatrix[4] = ecl.error(4);
1811  covarianceMatrix[5] = ecl.error(5); // error on theta
1812  eclCluster->setCovarianceMatrix(covarianceMatrix);
1813 
1814  eclCluster->setLAT(eclAux.width());
1815  eclCluster->setEnergyRaw(eclAux.mass());
1816  eclCluster->setE9oE21(eclAux.e9oe25());
1817  eclCluster->setEnergyHighestCrystal(eclAux.seed());
1818  eclCluster->setTime(eclAux.property(0));
1819  eclCluster->setNumberOfCrystals(eclAux.nhits());
1820 }
1821 
1822 void B2BIIConvertMdstModule::convertMdstKLMObject(const Belle::Mdst_klm_cluster& klm_cluster, KLMCluster* klmCluster)
1823 {
1824  // note: Belle quality flag is not saved (no free int variable in Belle2 KLMCluster)
1825  klmCluster->setLayers(klm_cluster.layers());
1826  klmCluster->setInnermostLayer(klm_cluster.first_layer());
1827 }
1828 
1829 
1830 //-----------------------------------------------------------------------------
1831 // RELATIONS
1832 //-----------------------------------------------------------------------------
1833 void B2BIIConvertMdstModule::setTracksToECLClustersRelations()
1834 {
1835  // Relations
1836  RelationArray tracksToECLClusters(m_tracks, m_eclClusters);
1837 
1838  Belle::Mdst_ecl_trk_Manager& m = Belle::Mdst_ecl_trk_Manager::get_manager();
1839  Belle::Mdst_charged_Manager& chgMg = Belle::Mdst_charged_Manager::get_manager();
1840 
1841  // We first insert relations to tracks which are directly matched (type == 1)
1842  // secondly we had CR matched tracks (connected region) (type == 2)
1843  // finally tracks which are geometrically matched (type == 0)
1844  std::vector<int> insert_order_types = {1, 2, 0};
1845  for (auto& insert_type : insert_order_types) {
1846  for (Belle::Mdst_ecl_trk_Manager::iterator ecltrkIterator = m.begin(); ecltrkIterator != m.end(); ++ecltrkIterator) {
1847  Belle::Mdst_ecl_trk mECLTRK = *ecltrkIterator;
1848 
1849  if (mECLTRK.type() != insert_type)
1850  continue;
1851 
1852  Belle::Mdst_ecl mdstEcl = mECLTRK.ecl();
1853  Belle::Mdst_trk mTRK = mECLTRK.trk();
1854 
1855  if (!mdstEcl)
1856  continue;
1857 
1858  // the numbering in mdst_charged
1859  // not necessarily the same as in mdst_trk
1860  // therfore have to find corresponding mdst_charged
1861  for (Belle::Mdst_charged_Manager::iterator chgIterator = chgMg.begin(); chgIterator != chgMg.end(); ++chgIterator) {
1862  Belle::Mdst_charged mChar = *chgIterator;
1863  Belle::Mdst_trk mTRK_in_charged = mChar.trk();
1864 
1865  if (mTRK_in_charged.get_ID() == mTRK.get_ID()) {
1866  // found the correct mdst_charged
1867  tracksToECLClusters.add(mChar.get_ID() - 1, mdstEcl.get_ID() - 1, 1.0);
1868  break;
1869  }
1870  }
1871  }
1872  }
1873 }
1874 
1875 
1876 void B2BIIConvertMdstModule::setKLMClustersRelations()
1877 {
1878  // Relations
1879  RelationArray klmClustersToTracks(m_klmClusters, m_tracks);
1880  RelationArray klmClustersToEclClusters(m_klmClusters, m_eclClusters);
1881 
1882  Belle::Mdst_klm_cluster_Manager& klm_cluster_manager = Belle::Mdst_klm_cluster_Manager::get_manager();
1883 
1884 
1885  for (Belle::Mdst_klm_cluster_Manager::iterator klmC_Ite = klm_cluster_manager.begin(); klmC_Ite != klm_cluster_manager.end();
1886  ++klmC_Ite) {
1887 
1888  Belle::Mdst_klm_cluster mdstKlm_cluster = *klmC_Ite;
1889  Belle::Mdst_trk mTRK = mdstKlm_cluster.trk();
1890  Belle::Mdst_ecl mECL = mdstKlm_cluster.ecl();
1891 
1892  if (mTRK) klmClustersToTracks.add(mdstKlm_cluster.get_ID() - 1, mTRK.get_ID() - 1);
1893  if (mECL) klmClustersToEclClusters.add(mdstKlm_cluster.get_ID() - 1, mECL.get_ID() - 1);
1894  }
1895 }
1896 
1897 
1898 //-----------------------------------------------------------------------------
1899 // MISC
1900 //-----------------------------------------------------------------------------
1901 
1902 int B2BIIConvertMdstModule::recoverMoreThan24bitIDHEP(int id)
1903 {
1904  /*
1905  QUICK CHECK: most of the normal particles are smaller than
1906  0x100000, while all the corrupt id has some of the high bits on.
1907 
1908  This bit check has to be revised when the table below is updated.
1909  */
1910  const int mask = 0x00f00000;
1911  int high_bits = id & mask;
1912  if (high_bits == 0 || high_bits == mask) return id;
1913 
1914  switch (id) {
1915  case 7114363:
1916  return 91000443; // X(3940)
1917  case 6114363:
1918  return 90000443; // Y(3940)
1919  case 6114241:
1920  return 90000321; // K_0*(800)+
1921  case 6114231:
1922  return 90000311; // K_0*(800)0
1923  case -6865004:
1924  return 9912212; // p_diff+
1925  case -6865104:
1926  return 9912112; // n_diffr
1927  case -6866773:
1928  return 9910443; // psi_diff
1929  case -6866883:
1930  return 9910333; // phi_diff
1931  case -6866993:
1932  return 9910223; // omega_diff
1933  case -6867005:
1934  return 9910211; // pi_diff+
1935  case -6867103:
1936  return 9910113; // rho_diff0
1937  case -7746995:
1938  return 9030221; // f_0(1500)
1939  case -7756773:
1940  return 9020443; // psi(4415)
1941  case -7756995:
1942  return 9020221; // eta(1405)
1943  case -7766773:
1944  return 9010443; // psi(4160)
1945  case -7776663:
1946  return 9000553; // Upsilon(5S)
1947  case -7776773:
1948  return 9000443; // psi(4040)
1949  case -7776783:
1950  return 9000433; // D_sj(2700)+
1951  case -7776995:
1952  return 9000221; // f_0(600)
1953  case -6114241:
1954  return -90000321; // K_0*(800)-
1955  case -6114231:
1956  return -90000311; // anti-K_0*(800)0
1957  case 6865004:
1958  return -9912212; // anti-p_diff-
1959  case 6865104:
1960  return -9912112; // anti-n_diffr
1961  case 6867005:
1962  return -9910211; // pi_diff-
1963  case 7776783:
1964  return -9000433; // D_sj(2700)-
1965  default:
1966  return id;
1967  }
1968 }
1969 
1970 void B2BIIConvertMdstModule::testMCRelation(const Belle::Gen_hepevt& belleMC, const MCParticle* mcP, const std::string& objectName)
1971 {
1972  int bellePDGCode = belleMC.idhep();
1973  int belleIIPDGCode = mcP->getPDG();
1974 
1975  if (bellePDGCode == 0)
1976  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to Gen_hepevt with idhep = 0.");
1977 
1978  if (bellePDGCode != belleIIPDGCode)
1979  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to different MCParticle! " << bellePDGCode << " vs. " <<
1980  belleIIPDGCode);
1981 
1982  double belleMomentum[] = { belleMC.PX(), belleMC.PY(), belleMC.PZ() };
1983  double belle2Momentum[] = { mcP->get4Vector().Px(), mcP->get4Vector().Py(), mcP->get4Vector().Pz() };
1984 
1985  for (unsigned i = 0; i < 3; i++) {
1986  double relDev = (belle2Momentum[i] - belleMomentum[i]) / belleMomentum[i];
1987 
1988  if (relDev > 1e-3) {
1989  B2WARNING("[B2BIIConvertMdstModule] " << objectName << " matched to different MCParticle!");
1990  B2INFO(" - Gen_hepevt [" << bellePDGCode << "] px/py/pz = " << belleMC.PX() << "/" << belleMC.PY() << "/" << belleMC.PZ());
1991  B2INFO(" - TrackFitResult [" << belleIIPDGCode << "] px/py/pz = " << mcP->get4Vector().Px() << "/" << mcP->get4Vector().Py() << "/"
1992  << mcP->get4Vector().Pz());
1993  }
1994  }
1995 }
1996 
1997 void B2BIIConvertMdstModule::belleVeeDaughterToCartesian(const Belle::Mdst_vee2& vee, const int charge,
1998  const Const::ParticleType& pType,
1999  CLHEP::HepLorentzVector& momentum, HepPoint3D& position, CLHEP::HepSymMatrix& error)
2000 {
2001  const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2002  CLHEP::HepVector a(5);
2003  CLHEP::HepSymMatrix Ea(5, 0);
2004  if (charge > 0) {
2005  a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2006  a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2007  a[4] = vee.daut().helix_p(4);
2008  Ea[0][0] = vee.daut().error_p(0); Ea[1][0] = vee.daut().error_p(1);
2009  Ea[1][1] = vee.daut().error_p(2); Ea[2][0] = vee.daut().error_p(3);
2010  Ea[2][1] = vee.daut().error_p(4); Ea[2][2] = vee.daut().error_p(5);
2011  Ea[3][0] = vee.daut().error_p(6); Ea[3][1] = vee.daut().error_p(7);
2012  Ea[3][2] = vee.daut().error_p(8); Ea[3][3] = vee.daut().error_p(9);
2013  Ea[4][0] = vee.daut().error_p(10); Ea[4][1] = vee.daut().error_p(11);
2014  Ea[4][2] = vee.daut().error_p(12); Ea[4][3] = vee.daut().error_p(13);
2015  Ea[4][4] = vee.daut().error_p(14);
2016  } else {
2017  a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2018  a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2019  a[4] = vee.daut().helix_m(4);
2020  Ea[0][0] = vee.daut().error_m(0); Ea[1][0] = vee.daut().error_m(1);
2021  Ea[1][1] = vee.daut().error_m(2); Ea[2][0] = vee.daut().error_m(3);
2022  Ea[2][1] = vee.daut().error_m(4); Ea[2][2] = vee.daut().error_m(5);
2023  Ea[3][0] = vee.daut().error_m(6); Ea[3][1] = vee.daut().error_m(7);
2024  Ea[3][2] = vee.daut().error_m(8); Ea[3][3] = vee.daut().error_m(9);
2025  Ea[4][0] = vee.daut().error_m(10); Ea[4][1] = vee.daut().error_m(11);
2026  Ea[4][2] = vee.daut().error_m(12); Ea[4][3] = vee.daut().error_m(13);
2027  Ea[4][4] = vee.daut().error_m(14);
2028  }
2029 
2030  Belle::Helix helix(pivot, a, Ea);
2031 
2032  // this is Vee daughter momentum/position/error at pivot = V0 Decay Vertex
2033  momentum = helix.momentum(0., pType.getMass(), position, error);
2034 }
2035 
2036 void B2BIIConvertMdstModule::belleVeeDaughterHelix(const Belle::Mdst_vee2& vee, const int charge, std::vector<float>& helixParam,
2037  std::vector<float>& helixError)
2038 {
2039  const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2040  CLHEP::HepVector a(5);
2041  CLHEP::HepSymMatrix Ea(5, 0);
2042  if (charge > 0) {
2043  a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2044  a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2045  a[4] = vee.daut().helix_p(4);
2046  Ea[0][0] = vee.daut().error_p(0);
2047  Ea[1][0] = vee.daut().error_p(1);
2048  Ea[1][1] = vee.daut().error_p(2);
2049  Ea[2][0] = vee.daut().error_p(3);
2050  Ea[2][1] = vee.daut().error_p(4);
2051  Ea[2][2] = vee.daut().error_p(5);
2052  Ea[3][0] = vee.daut().error_p(6);
2053  Ea[3][1] = vee.daut().error_p(7);
2054  Ea[3][2] = vee.daut().error_p(8);
2055  Ea[3][3] = vee.daut().error_p(9);
2056  Ea[4][0] = vee.daut().error_p(10);
2057  Ea[4][1] = vee.daut().error_p(11);
2058  Ea[4][2] = vee.daut().error_p(12);
2059  Ea[4][3] = vee.daut().error_p(13);
2060  Ea[4][4] = vee.daut().error_p(14);
2061  } else {
2062  a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2063  a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2064  a[4] = vee.daut().helix_m(4);
2065  Ea[0][0] = vee.daut().error_m(0);
2066  Ea[1][0] = vee.daut().error_m(1);
2067  Ea[1][1] = vee.daut().error_m(2);
2068  Ea[2][0] = vee.daut().error_m(3);
2069  Ea[2][1] = vee.daut().error_m(4);
2070  Ea[2][2] = vee.daut().error_m(5);
2071  Ea[3][0] = vee.daut().error_m(6);
2072  Ea[3][1] = vee.daut().error_m(7);
2073  Ea[3][2] = vee.daut().error_m(8);
2074  Ea[3][3] = vee.daut().error_m(9);
2075  Ea[4][0] = vee.daut().error_m(10);
2076  Ea[4][1] = vee.daut().error_m(11);
2077  Ea[4][2] = vee.daut().error_m(12);
2078  Ea[4][3] = vee.daut().error_m(13);
2079  Ea[4][4] = vee.daut().error_m(14);
2080  }
2081 
2082  Belle::Helix helix(pivot, a, Ea);
2083 
2084  // go to the new pivot
2085  helix.pivot(HepPoint3D(0., 0., 0.));
2086 
2087  CLHEP::HepSymMatrix error5x5(5, 0);
2088  convertHelix(helix, helixParam, error5x5);
2089 
2090  unsigned int size = 5;
2091  unsigned int counter = 0;
2092  for (unsigned int i = 0; i < size; i++)
2093  for (unsigned int j = i; j < size; j++)
2094  helixError[counter++] = error5x5[i][j];
2095 }
2096 
2097 void B2BIIConvertMdstModule::belleHelixToHelix(const Belle::Mdst_trk_fit& trk_fit,
2098  std::vector<float>& helixParam, std::vector<float>& helixError)
2099 {
2100  const HepPoint3D pivot(trk_fit.pivot_x(),
2101  trk_fit.pivot_y(),
2102  trk_fit.pivot_z());
2103 
2104  CLHEP::HepVector a(5);
2105  a[0] = trk_fit.helix(0);
2106  a[1] = trk_fit.helix(1);
2107  a[2] = trk_fit.helix(2);
2108  a[3] = trk_fit.helix(3);
2109  a[4] = trk_fit.helix(4);
2110  CLHEP::HepSymMatrix Ea(5, 0);
2111  Ea[0][0] = trk_fit.error(0);
2112  Ea[1][0] = trk_fit.error(1);
2113  Ea[1][1] = trk_fit.error(2);
2114  Ea[2][0] = trk_fit.error(3);
2115  Ea[2][1] = trk_fit.error(4);
2116  Ea[2][2] = trk_fit.error(5);
2117  Ea[3][0] = trk_fit.error(6);
2118  Ea[3][1] = trk_fit.error(7);
2119  Ea[3][2] = trk_fit.error(8);
2120  Ea[3][3] = trk_fit.error(9);
2121  Ea[4][0] = trk_fit.error(10);
2122  Ea[4][1] = trk_fit.error(11);
2123  Ea[4][2] = trk_fit.error(12);
2124  Ea[4][3] = trk_fit.error(13);
2125  Ea[4][4] = trk_fit.error(14);
2126  Belle::Helix helix(pivot, a, Ea);
2127 
2128  helix.pivot(HepPoint3D(0., 0., 0.));
2129 
2130  // convert to Belle II helix parameters
2131  // param 0: d_0 = d_rho
2132  helixParam[0] = helix.a()[0];
2133 
2134  // param 1: phi = phi_0 + pi/2
2135  helixParam[1] = helix.a()[1] + TMath::Pi() / 2.0;
2136 
2137  // param 2: omega = Kappa * alpha = Kappa * B[Tesla] * speed_of_light[m/s] * 1e-11
2138  helixParam[2] = helix.a()[2] * KAPPA2OMEGA;
2139 
2140  // param 3: d_z = z0
2141  helixParam[3] = helix.a()[3];
2142 
2143  // param 4: tan(Lambda) = cotTheta
2144  helixParam[4] = helix.a()[4];
2145 
2146  Ea = helix.Ea();
2147 
2148  helixError[0] = Ea[0][0];
2149  helixError[1] = Ea[1][0];
2150  helixError[2] = Ea[2][0] * KAPPA2OMEGA;
2151  helixError[3] = Ea[3][0];
2152  helixError[4] = Ea[4][0];
2153  helixError[5] = Ea[1][1];
2154  helixError[6] = Ea[2][1] * KAPPA2OMEGA;
2155  helixError[7] = Ea[3][1];
2156  helixError[8] = Ea[4][1];
2157  helixError[9] = Ea[2][2] * KAPPA2OMEGA * KAPPA2OMEGA;
2158  helixError[10] = Ea[3][2] * KAPPA2OMEGA;
2159  helixError[11] = Ea[4][2] * KAPPA2OMEGA;
2160  helixError[12] = Ea[3][3];
2161  helixError[13] = Ea[4][3];
2162  helixError[14] = Ea[4][4];
2163 }
2164 
2165 int B2BIIConvertMdstModule::belleHelixToCartesian(const Belle::Mdst_trk_fit& trk_fit, const double mass,
2166  const HepPoint3D& newPivot,
2167  CLHEP::HepLorentzVector& momentum, HepPoint3D& position, CLHEP::HepSymMatrix& error, double dPhi)
2168 {
2169  const HepPoint3D pivot(trk_fit.pivot_x(),
2170  trk_fit.pivot_y(),
2171  trk_fit.pivot_z());
2172 
2173  CLHEP::HepVector a(5);
2174  a[0] = trk_fit.helix(0);
2175  a[1] = trk_fit.helix(1);
2176  a[2] = trk_fit.helix(2);
2177  a[3] = trk_fit.helix(3);
2178  a[4] = trk_fit.helix(4);
2179  CLHEP::HepSymMatrix Ea(5, 0);
2180  Ea[0][0] = trk_fit.error(0);
2181  Ea[1][0] = trk_fit.error(1);
2182  Ea[1][1] = trk_fit.error(2);
2183  Ea[2][0] = trk_fit.error(3);
2184  Ea[2][1] = trk_fit.error(4);
2185  Ea[2][2] = trk_fit.error(5);
2186  Ea[3][0] = trk_fit.error(6);
2187  Ea[3][1] = trk_fit.error(7);
2188  Ea[3][2] = trk_fit.error(8);
2189  Ea[3][3] = trk_fit.error(9);
2190  Ea[4][0] = trk_fit.error(10);
2191  Ea[4][1] = trk_fit.error(11);
2192  Ea[4][2] = trk_fit.error(12);
2193  Ea[4][3] = trk_fit.error(13);
2194  Ea[4][4] = trk_fit.error(14);
2195 
2196  Belle::Helix helix(pivot, a, Ea);
2197 
2198  int charge = 0;
2199  if (helix.kappa() > 0)
2200  charge = 1;
2201  else
2202  charge = -1;
2203 
2204  if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
2205  helix.pivot(newPivot);
2206  momentum = helix.momentum(dPhi, mass, position, error);
2207  } else {
2208  if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
2209  helix.pivot(HepPoint3D(0., 0., 0.));
2210  momentum = helix.momentum(dPhi, mass, position, error);
2211  } else {
2212  momentum = helix.momentum(dPhi, mass, position, error);
2213  }
2214  }
2215  return charge;
2216 }
2217 
2218 TrackFitResult B2BIIConvertMdstModule::createTrackFitResult(const CLHEP::HepLorentzVector& momentum,
2219  const HepPoint3D& position,
2220  const CLHEP::HepSymMatrix& error,
2221  const short int charge,
2222  const Const::ParticleType& pType,
2223  const float pValue,
2224  const uint64_t hitPatternCDCInitializer,
2225  const uint32_t hitPatternVXDInitializer,
2226  const uint16_t ndf)
2227 {
2228  TVector3 pos(position.x(), position.y(), position.z());
2229  TVector3 mom(momentum.px(), momentum.py(), momentum.pz());
2230 
2231  TMatrixDSym errMatrix(6);
2232  for (unsigned i = 0; i < 7; i++) {
2233  if (i == 3)
2234  continue;
2235  for (unsigned j = 0; j < 7; j++) {
2236  if (j == 3)
2237  continue;
2238 
2239  if (i == j)
2240  errMatrix(ERRMCONV[i], ERRMCONV[i]) = error[i][i];
2241  else
2242  errMatrix(ERRMCONV[i], ERRMCONV[j]) = errMatrix(ERRMCONV[j], ERRMCONV[i]) = error[i][j];
2243  }
2244  }
2245 
2246  return TrackFitResult(pos, mom, errMatrix, charge, pType, pValue, BFIELD, hitPatternCDCInitializer, hitPatternVXDInitializer, ndf);
2247 }
2248 
2249 double B2BIIConvertMdstModule::atcPID(const PIDLikelihood* pid, int sigHyp, int bkgHyp)
2250 {
2251  if (!pid) return 0.5;
2252 
2253  // ACC = ARICH
2254  Const::PIDDetectorSet set = Const::ARICH;
2255  double acc_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2256  double acc_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2257  double acc = 0.5;
2258  if (acc_sig + acc_bkg > 0.0)
2259  acc = acc_sig / (acc_sig + acc_bkg);
2260 
2261  // TOF = TOP
2262  set = Const::TOP;
2263  double tof_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2264  double tof_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2265  double tof = 0.5;
2266  double tof_all = tof_sig + tof_bkg;
2267  if (tof_all != 0) {
2268  tof = tof_sig / tof_all;
2269  if (tof < 0.001) tof = 0.001;
2270  if (tof > 0.999) tof = 0.999;
2271  }
2272 
2273  // dE/dx = CDC
2274  set = Const::CDC;
2275  double cdc_sig = exp(pid->getLogL(c_belleHyp_to_chargedStable[sigHyp], set));
2276  double cdc_bkg = exp(pid->getLogL(c_belleHyp_to_chargedStable[bkgHyp], set));
2277  double cdc = 0.5;
2278  double cdc_all = cdc_sig + cdc_bkg;
2279  if (cdc_all != 0) {
2280  cdc = cdc_sig / cdc_all;
2281  if (cdc < 0.001) cdc = 0.001;
2282  if (cdc > 0.999) cdc = 0.999;
2283  }
2284 
2285  // Combined
2286  double pid_sig = acc * tof * cdc;
2287  double pid_bkg = (1. - acc) * (1. - tof) * (1. - cdc);
2288 
2289  return pid_sig / (pid_sig + pid_bkg);
2290 }
2291 
2292 
2293 void B2BIIConvertMdstModule::endRun()
2294 {
2295  B2DEBUG(99, "B2BIIConvertMdst: endRun done.");
2296 }
2297 
2298 
2299 void B2BIIConvertMdstModule::terminate()
2300 {
2301  B2DEBUG(99, "B2BIIConvertMdst: terminate called");
2302 }
2303 
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:1970
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:2218
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:1576
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:1833
Belle2::B2BIIConvertMdstModule::terminate
virtual void terminate() override
Terminates the module.
Definition: B2BIIConvertMdstModule.cc:2299
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:1876
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:2165
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:1997
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:1336
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:1318
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:1751
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:2249
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:1616
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:1791
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:1464
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:2097
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:1822
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:1902
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:2036
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:2293
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