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