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