9 #include <b2bii/modules/B2BIIMdstInput/B2BIIConvertMdstModule.h>
11 #include <framework/datastore/StoreObjPtr.h>
12 #include <framework/datastore/RelationArray.h>
13 #include <framework/database/Database.h>
14 #include <framework/pcore/ProcHandler.h>
16 #include <mdst/dataobjects/HitPatternVXD.h>
17 #include <mdst/dataobjects/HitPatternCDC.h>
20 #include <framework/gearbox/Unit.h>
21 #include <framework/gearbox/Const.h>
22 #include <framework/geometry/B2Vector3.h>
23 #include <analysis/dataobjects/ParticleExtraInfoMap.h>
26 #include <framework/dataobjects/EventMetaData.h>
27 #include <framework/dataobjects/Helix.h>
28 #include <framework/dataobjects/UncertainHelix.h>
31 #include <b2bii/utility/BelleMdstToGenHepevt.h>
34 #include <Math/RotationY.h>
35 #include <Math/Vector3D.h>
36 #include <Math/Vector4D.h>
44 #include "belle_legacy/eid/eid.h"
48 #include "belle_legacy/kid/kid_acc.h"
49 #include "belle_legacy/kid/kid_cdc.h"
53 #include "belle_legacy/findKs/findKs.h"
56 #ifdef HAVE_NISKSFINDER
57 #include "belle_legacy/nisKsFinder/nisKsFinder.h"
60 #ifdef HAVE_GOODLAMBDA
61 #include "belle_legacy/findLambda/findLambda.h"
64 #include "belle_legacy/benergy/BeamEnergy.h"
65 #include "belle_legacy/ip/IpProfile.h"
66 #include "belle_legacy/tables/evtcls.h"
67 #include "belle_legacy/tables/trg.h"
77 bool approximatelyEqual(
float a,
float b,
float epsilon)
79 return fabs(a - b) <= ((fabs(a) < fabs(b) ? fabs(b) : fabs(a)) * epsilon);
82 double adjustAngleRange(
double phi)
84 phi = phi - int(phi / TMath::TwoPi()) * TMath::TwoPi();
85 return phi - int(phi / TMath::Pi()) * TMath::TwoPi();
88 void fill7x7ErrorMatrix(
const TrackFitResult* tfr, TMatrixDSym& error7x7,
const double mass,
const double bField)
92 double d0 = tfr->
getD0();
100 double cosPhi0 = TMath::Cos(phi0);
101 double sinPhi0 = TMath::Sin(phi0);
105 rho = 1.0 / alpha / omega;
109 double energy = TMath::Sqrt(mass * mass + (1.0 + tanl * tanl) * rho * rho);
121 const int iOmega = 2;
125 TMatrixD jacobian(7, 5);
128 jacobian(iPx, iPhi0) = - fabs(rho) * sinPhi0;
129 jacobian(iPx, iOmega) = -
charge * rho * rho * cosPhi0 * alpha;
130 jacobian(iPy, iPhi0) = fabs(rho) * cosPhi0;
131 jacobian(iPy, iOmega) = -
charge * rho * rho * sinPhi0 * alpha;
132 jacobian(iPz, iOmega) = -
charge * rho * rho * tanl * alpha;
133 jacobian(iPz, iTanl) = fabs(rho);
134 if (omega != 0 && energy != 0) {
135 jacobian(iE, iOmega) = - (1.0 + tanl * tanl) * rho * rho / omega / energy;
136 jacobian(iE, iTanl) = tanl * rho * rho / energy;
138 jacobian(iE, iOmega) = (DBL_MAX);
139 jacobian(iE, iTanl) = (DBL_MAX);
141 jacobian(iX, iD0) = sinPhi0;
142 jacobian(iX, iPhi0) = d0 * cosPhi0;
143 jacobian(iY, iD0) = - cosPhi0;
144 jacobian(iY, iPhi0) = d0 * sinPhi0;
145 jacobian(iZ, iZ0) = 1.0;
149 error7x7 = error5x5.Similarity(jacobian);
161 m_mcMatchingMode(c_Direct)
164 setDescription(
"Converts Belle mDST objects (Panther tables and records) to Belle II mDST objects.");
167 "Convert beam parameters or use information stored in "
168 "Belle II database.",
true);
170 "Use 6x6 (position, momentum) covariance matrix for charged tracks instead of 5x5 (helix parameters) covariance matrix",
false);
172 "MC matching mode: 'Direct', or 'GeneratorLevel'",
173 std::string(
"Direct"));
175 "clusters with a E9/E25 value above this threshold are classified as neutral even if tracks are matched to their connected region (matchType == 2)",
179 addParam(
"nisKsInfo",
m_nisEnable,
"Flag to switch on conversion of nisKsFinder info",
true);
181 addParam(
"TrkExtra",
m_convertTrkExtra,
" Flag to switch on conversion of first_x,y,z and last_x,y,z from Mdst_trk_fit",
true);
185 B2DEBUG(1,
"B2BIIConvertMdst: Constructor done.");
203 B2INFO(
"B2BIIConvertMdst: initialized.");
205 B2WARNING(
"nisKsFinder output has been disabled. ksnbVLike, ksnbNoLam, ksnbStandard will not be converted.");
210 B2DEBUG(99,
"[B2BIIConvertMdstModule::initializeDataStore] initialization of DataStore started");
217 m_v0s.registerInDataStore();
258 B2DEBUG(99,
"[B2BIIConvertMdstModule::initializeDataStore] initialization of DataStore ended");
264 B2DEBUG(99,
"B2BIIConvertMdst: beginRun called.");
268 Belle::BeamEnergy::begin_run();
270 Belle::BeamEnergy::dump();
273 Belle::IpProfile::begin_run();
275 Belle::IpProfile::dump();
276 bool usableIP = Belle::IpProfile::usable();
277 B2DEBUG(99,
"B2BIIConvertMdst: IpProfile is usable = " << usableIP);
282 Belle::eid::init_data();
283 Belle::eid::show_use(
"ALL");
292 Belle::Belle_event_Manager& evman = Belle::Belle_event_Manager::get_manager();
293 Belle::Belle_event& evt = evman[0];
295 if (evt.ExpMC() == 2)
308 B2INFO(
"No database entry for this run yet, create one");
316 B2ERROR(
"BeamParameters from condition database are different from converted "
317 "ones, overriding database. Did you make sure the globaltag B2BII is used?");
319 B2INFO(
"BeamSpot, BoostVector, and InvariantMass from condition database are different from converted "
320 "ones, overriding database");
323 B2FATAL(
"Cannot reliably override the Database content in parallel processing "
324 "mode, please run the conversion in single processing mode");
376 const double Eher = Belle::BeamEnergy::E_HER();
377 const double Eler = Belle::BeamEnergy::E_LER();
378 const double crossingAngle = Belle::BeamEnergy::Cross_angle();
379 const double angleLer = M_PI;
380 const double angleHer = crossingAngle;
382 TMatrixDSym covariance(0);
383 HepLorentzVector p_beam = Belle::BeamEnergy::p_beam();
386 ROOT::Math::PxPyPzEVector P_her(0.0, 0.0, TMath::Sqrt(Eher * Eher - mass_e * mass_e), Eher);
387 ROOT::Math::RotationY rotateAroundYAxis(angleHer);
388 P_her = rotateAroundYAxis(P_her);
389 ROOT::Math::PxPyPzEVector P_ler(0.0, 0.0, TMath::Sqrt(Eler * Eler - mass_e * mass_e), Eler);
390 rotateAroundYAxis.SetAngle(angleLer);
391 P_ler = rotateAroundYAxis(P_ler);
394 ROOT::Math::PxPyPzEVector P_beam = P_her + P_ler;
399 B2DEBUG(99,
"Beam Energy: E_HER = " << Eher <<
"; E_LER = " << Eler <<
"; angle = " << crossingAngle);
400 B2DEBUG(99,
"Beam Momentum (pre-convert) : P_X = " << p_beam.px() <<
"; P_Y = " << p_beam.py() <<
"; P_Z = " << p_beam.pz());
401 B2DEBUG(99,
"Beam Momentum (post-convert) : P_X = " << P_beam.Px() <<
"; P_Y = " << P_beam.Py() <<
"; P_Z = " << P_beam.Pz());
406 if (!Belle::IpProfile::usable()) {
411 TVector3(std::numeric_limits<double>::quiet_NaN(),
412 std::numeric_limits<double>::quiet_NaN(),
413 std::numeric_limits<double>::quiet_NaN()
414 ), TMatrixTSym<double>()
420 CLHEP::HepSymMatrix ipErr;
423 ip = Belle::IpProfile::position();
424 ipErr = Belle::IpProfile::position_err();
427 Belle::IpProfile::set_evtbin_number();
431 ip = Belle::IpProfile::e_position();
432 ipErr = Belle::IpProfile::e_position_err();
437 TMatrixDSym cov(ipErr.num_col());
438 for (
int i = 0; i < ipErr.num_row(); ++i) {
439 for (
int j = 0; j < ipErr.num_col(); ++j) {
440 cov(i, j) = ipErr(i + 1, j + 1);
452 Belle::Mdst_charged_Manager& m = Belle::Mdst_charged_Manager::get_manager();
453 for (Belle::Mdst_charged_Manager::iterator chargedIterator = m.begin(); chargedIterator != m.end(); ++chargedIterator) {
454 Belle::Mdst_charged belleTrack = *chargedIterator;
468 const Belle::Gen_hepevt& hep0 = get_hepevt(belleTrack);
471 const Belle::Gen_hepevt* hep =
nullptr;
477 hep = &gen_level(hep0);
485 tracksToMCParticles.
add(track->getArrayIndex(), matchedMCParticle);
489 B2DEBUG(99,
"Can not find MCParticle corresponding to this gen_hepevt (Panther ID = " << hep->get_ID() <<
")");
490 B2DEBUG(99,
"Gen_hepevt: Panther ID = " << hep->get_ID() <<
"; idhep = " << hep->idhep() <<
"; isthep = " << hep->isthep());
500 ksPList->initialize(310, ksPList.
getName());
505 lambda0PList->initialize(3122, lambda0PList.
getName());
508 antiLambda0PList.
create();
509 antiLambda0PList->initialize(-3122, antiLambda0PList.
getName());
511 antiLambda0PList->bindAntiParticleList(*lambda0PList);
516 convGammaPList->initialize(22, convGammaPList.
getName());
519 Belle::Mdst_vee2_Manager& m = Belle::Mdst_vee2_Manager::get_manager();
520 for (Belle::Mdst_vee2_Manager::iterator vee2Iterator = m.begin(); vee2Iterator != m.end(); ++vee2Iterator) {
521 Belle::Mdst_vee2 belleV0 = *vee2Iterator;
524 Belle::Mdst_charged belleTrackP = belleV0.chgd(0);
526 Belle::Mdst_charged belleTrackM = belleV0.chgd(1);
534 switch (belleV0.kind()) {
560 B2WARNING(
"Conversion of vee2 candidate of unknown kind! kind = " << belleV0.kind());
564 int trackID[2] = {0, 0};
566 Belle::Mdst_charged_Manager& charged_mag = Belle::Mdst_charged_Manager::get_manager();
567 for (std::vector<Belle::Mdst_charged>::iterator chgIterator = charged_mag.begin(); chgIterator != charged_mag.end();
569 if (belleV0.chgd(0).get_ID() >= 1 && trackID[0] == 0 && belleV0.chgd(0).get_ID() == chgIterator->get_ID()) {
570 trackID[0] = (int)(chgIterator->get_ID());
573 if (belleV0.chgd(1).get_ID() >= 1 && trackID[1] == 0 && belleV0.chgd(1).get_ID() == chgIterator->get_ID()) {
574 trackID[1] = (int)(chgIterator->get_ID());
581 HepPoint3D dauPivot(belleV0.vx(), belleV0.vy(), belleV0.vz());
582 int trackFitPIndex = -1;
583 int trackFitMIndex = -1;
585 CLHEP::HepLorentzVector momentumP;
586 CLHEP::HepSymMatrix error7x7P(7, 0);
587 HepPoint3D positionP;
588 TMatrixFSym errMatrixP(7);
589 CLHEP::HepLorentzVector momentumM;
590 CLHEP::HepSymMatrix error7x7M(7, 0);
591 HepPoint3D positionM;
592 TMatrixFSym errMatrixM(7);
593 CLHEP::HepSymMatrix error5x5(5, 0);
594 if (trackID[0] >= 1) {
595 if (belleV0.daut()) {
596 std::vector<float> helixParam(5);
597 std::vector<float> helixError(15);
600 auto trackFitP =
m_trackFitResults.appendNew(helixParam, helixError, pTypeP, 0.5, -1, -1, 0);
601 trackFitPIndex = trackFitP->getArrayIndex();
608 for (
unsigned i = 0; i < 7; i++)
609 for (
unsigned j = 0; j < 7; j++)
610 errMatrixP(i, j) = error7x7P[i][j];
612 daughterP =
Particle(trackID[0] - 1, tmpTFR, pTypeP);
613 daughterP.
updateMomentum(ROOT::Math::PxPyPzEVector(momentumP.px(), momentumP.py(), momentumP.pz(), momentumP.e()),
614 ROOT::Math::XYZVector(positionP.x(), positionP.y(), positionP.z()),
618 Belle::Mdst_trk_fit& trk_fit = charged_mag[trackID[0] - 1].trk().mhyp(belleHypP);
619 double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
621 std::vector<float> helixParam(5);
622 std::vector<float> helixError(15);
623 convertHelix(trk_fit, HepPoint3D(0., 0., 0.), helixParam, helixError);
626 if (helixParam[2] == 0) {
627 B2WARNING(
"Helix parameter for curvature == 0. Skipping Track! The parameter is: " << helixParam[2] <<
"...");
631 auto trackFitP =
m_trackFitResults.appendNew(helixParam, helixError, pTypeP, pValue, -1, -1, 0);
633 trackFitPIndex = trackFitP->getArrayIndex();
635 daughterP =
Particle(trackID[0] - 1, trackFitP, pTypeP);
638 helixParam, error5x5,
639 momentumP, positionP, error7x7P);
641 for (
unsigned i = 0; i < 7; i++)
642 for (
unsigned j = 0; j < 7; j++)
643 errMatrixP(i, j) = error7x7P[i][j];
645 daughterP.
updateMomentum(ROOT::Math::PxPyPzEVector(momentumP.px(), momentumP.py(), momentumP.pz(), momentumP.e()),
646 ROOT::Math::XYZVector(positionP.x(), positionP.y(), positionP.z()),
650 if (trackID[1] >= 1) {
651 if (belleV0.daut()) {
652 std::vector<float> helixParam(5);
653 std::vector<float> helixError(15);
656 auto trackFitM =
m_trackFitResults.appendNew(helixParam, helixError, pTypeM, 0.5, -1, -1, 0);
657 trackFitMIndex = trackFitM->getArrayIndex();
663 for (
unsigned i = 0; i < 7; i++)
664 for (
unsigned j = 0; j < 7; j++)
665 errMatrixM(i, j) = error7x7M[i][j];
667 daughterM =
Particle(trackID[1] - 1, tmpTFR, pTypeM);
668 daughterM.
updateMomentum(ROOT::Math::PxPyPzEVector(momentumM.px(), momentumM.py(), momentumM.pz(), momentumM.e()),
669 ROOT::Math::XYZVector(positionM.x(), positionM.y(), positionM.z()),
673 Belle::Mdst_trk_fit& trk_fit = charged_mag[trackID[1] - 1].trk().mhyp(belleHypM);
674 double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
676 std::vector<float> helixParam(5);
677 std::vector<float> helixError(15);
678 convertHelix(trk_fit, HepPoint3D(0., 0., 0.), helixParam, helixError);
681 if (helixParam[2] == 0) {
682 B2WARNING(
"Helix parameter for curvature == 0. Skipping Track! The parameter is: " << helixParam[2] <<
"...");
686 auto trackFitM =
m_trackFitResults.appendNew(helixParam, helixError, pTypeM, pValue, -1, -1, 0);
688 trackFitMIndex = trackFitM->getArrayIndex();
690 daughterM =
Particle(trackID[1] - 1, trackFitM, pTypeM);
693 helixParam, error5x5,
694 momentumM, positionM, error7x7M);
696 for (
unsigned i = 0; i < 7; i++)
697 for (
unsigned j = 0; j < 7; j++)
698 errMatrixM(i, j) = error7x7M[i][j];
700 daughterM.
updateMomentum(ROOT::Math::PxPyPzEVector(momentumM.px(), momentumM.py(), momentumM.pz(), momentumM.e()),
701 ROOT::Math::XYZVector(positionM.x(), positionM.y(), positionM.z()),
712 m_v0s.appendNew(std::make_pair(trackP, trackFitP), std::make_pair(trackM, trackFitM));
722 newDaugP->addRelationTo(pidP);
724 newDaugP->addRelationTo(mcParticleP);
731 ROOT::Math::PxPyPzEVector v0Momentum(belleV0.px(), belleV0.py(), belleV0.pz(), belleV0.energy());
732 ROOT::Math::XYZVector v0Vertex(belleV0.vx(), belleV0.vy(), belleV0.vz());
738 auto appendVertexFitInfo = [](Belle::Mdst_vee2 & _belle_V0,
Particle & _belle2_V0) {
740 _belle2_V0.addExtraInfo(
"chiSquared", _belle_V0.chisq());
742 _belle2_V0.addExtraInfo(
"ndf", 1);
744 double prob = TMath::Prob(_belle_V0.chisq(), 1);
745 _belle2_V0.setPValue(prob);
749 if (belleV0.kind() == 1) {
754 appendVertexFitInfo(belleV0, KS);
756 ksPList->addParticle(newV0);
759 Belle::FindKs belleKSFinder;
760 belleKSFinder.candidates(belleV0, Belle::IpProfile::position(1));
797 }
else if (belleV0.kind() == 2) {
802 appendVertexFitInfo(belleV0, Lambda0);
804 lambda0PList->addParticle(newV0);
807 Belle::FindLambda lambdaFinder;
808 lambdaFinder.candidates(belleV0, Belle::IpProfile::position(1));
809 newV0->
addExtraInfo(
"goodLambda", lambdaFinder.goodLambda());
810 }
else if (belleV0.kind() == 3) {
811 Particle antiLambda0(v0Momentum, -3122);
815 appendVertexFitInfo(belleV0, antiLambda0);
817 antiLambda0PList->addParticle(newV0);
820 Belle::FindLambda lambdaFinder;
821 lambdaFinder.candidates(belleV0, Belle::IpProfile::position(1));
822 newV0->
addExtraInfo(
"goodLambda", lambdaFinder.goodLambda());
823 }
else if (belleV0.kind() == 4) {
828 appendVertexFitInfo(belleV0, gamma);
830 convGammaPList->addParticle(newV0);
834 if (belleV0.kind() <= 3) {
835 Belle::nisKsFinder ksnb;
836 double protIDP =
atcPID(pidP, 2, 4);
837 double protIDM =
atcPID(pidM, 2, 4);
838 ksnb.candidates(belleV0, Belle::IpProfile::position(1), momentumP, protIDP, protIDM);
843 if (belleV0.kind() == 1)
859 Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
860 if (genMgr.count() == 0)
863 typedef std::pair<MCParticleGraph::GraphParticle*, Belle::Gen_hepevt> halfFamily;
864 halfFamily currFamily;
866 std::queue < halfFamily > heritancesQueue;
872 for (Belle::Gen_hepevt_Manager::iterator genIterator = genMgr.begin();
873 genIterator != genMgr.end(); ++genIterator) {
874 Belle::Gen_hepevt hep = *genIterator;
876 if (!(hep.moFirst() == 0 && hep.moLast() == 0))
884 for (
int iDaughter = hep.daFirst(); iDaughter <= hep.daLast();
886 if (iDaughter == 0) {
887 B2DEBUG(95,
"Trying to access generated daughter with Panther ID == 0");
890 currFamily.first = graphParticle;
891 currFamily.second = genMgr(Belle::Panther_ID(iDaughter));
892 heritancesQueue.push(currFamily);
897 while (!heritancesQueue.empty()) {
898 currFamily = heritancesQueue.front();
899 heritancesQueue.pop();
902 Belle::Gen_hepevt& currDaughter = currFamily.second;
905 if (currDaughter.idhep() == 0)
919 int nGrandChildren = currDaughter.daLast() - currDaughter.daFirst() + 1;
921 if (nGrandChildren > 0 && currDaughter.daFirst() != 0) {
922 for (
int igrandchild = currDaughter.daFirst(); igrandchild <= currDaughter.daLast(); ++igrandchild) {
923 if (igrandchild == 0) {
924 B2DEBUG(95,
"Trying to access generated daughter with Panther ID == 0");
928 family.first = graphDaughter;
929 family.second = genMgr(Belle::Panther_ID(igrandchild));
930 heritancesQueue.push(family);
947 Belle::Mdst_ecl_Manager& ecl_manager = Belle::Mdst_ecl_Manager::get_manager();
948 Belle::Mdst_ecl_aux_Manager& ecl_aux_manager = Belle::Mdst_ecl_aux_Manager::get_manager();
950 for (Belle::Mdst_ecl_Manager::iterator eclIterator = ecl_manager.begin(); eclIterator != ecl_manager.end(); ++eclIterator) {
953 Belle::Mdst_ecl mdstEcl = *eclIterator;
954 Belle::Mdst_ecl_aux mdstEclAux(ecl_aux_manager(mdstEcl.get_ID()));
965 B2EclCluster->setConnectedRegionId(B2EclCluster->getArrayIndex() + 1);
966 B2EclCluster->setClusterId(1);
973 const Belle::Gen_hepevt& hep0 = get_hepevt(mdstEcl);
976 const Belle::Gen_hepevt* hep =
nullptr;
982 hep = &gen_level(hep0);
989 eclClustersToMCParticles.
add(B2EclCluster->getArrayIndex(), matchedMCParticleID);
992 B2DEBUG(79,
"Cannot find MCParticle corresponding to this gen_hepevt (Panther ID = " << hep->get_ID() <<
")");
993 B2DEBUG(79,
"Gen_hepevt: Panther ID = " << hep->get_ID() <<
"; idhep = " << hep->idhep() <<
"; isthep = " << hep->isthep());
1006 Belle::Mdst_klm_cluster_Manager& klm_cluster_manager = Belle::Mdst_klm_cluster_Manager::get_manager();
1008 for (Belle::Mdst_klm_cluster_Manager::iterator klmC_Ite = klm_cluster_manager.begin(); klmC_Ite != klm_cluster_manager.end();
1012 Belle::Mdst_klm_cluster mdstKlm_cluster = *klmC_Ite;
1035 plist->initialize(22,
"gamma:mdst");
1038 Belle::Mdst_gamma_Manager& gamma_manager = Belle::Mdst_gamma_Manager::get_manager();
1040 for (Belle::Mdst_gamma_Manager::iterator gammaIterator = gamma_manager.begin(); gammaIterator != gamma_manager.end();
1044 Belle::Mdst_gamma mdstGamma = *gammaIterator;
1045 Belle::Mdst_ecl mdstEcl = mdstGamma.ecl();
1059 plist->addParticle(B2Gamma);
1066 if (matchedMCParticle)
1076 plist->initialize(111,
"pi0:mdst");
1079 Belle::Mdst_pi0_Manager& pi0_manager = Belle::Mdst_pi0_Manager::get_manager();
1080 for (Belle::Mdst_pi0_Manager::iterator pi0Iterator = pi0_manager.begin(); pi0Iterator != pi0_manager.end(); ++pi0Iterator) {
1083 Belle::Mdst_pi0 mdstPi0 = *pi0Iterator;
1084 Belle::Mdst_gamma mdstGamma1 = mdstPi0.gamma(0);
1085 Belle::Mdst_gamma mdstGamma2 = mdstPi0.gamma(1);
1086 if (!mdstGamma1 || !mdstGamma2)
1089 ROOT::Math::PxPyPzEVector p4(mdstPi0.px(), mdstPi0.py(), mdstPi0.pz(), mdstPi0.energy());
1097 if (!B2Gamma1 || !B2Gamma2)
1111 double prob = TMath::Prob(mdstPi0.chisq(), 1);
1115 plist->addParticle(B2Pi0);
1128 plist->initialize(
Const::Klong.getPDGCode(),
"K_L0:mdst");
1130 Belle::Mdst_klong_Manager& klong_manager = Belle::Mdst_klong_Manager::get_manager();
1131 for (Belle::Mdst_klong_Manager::iterator klong_Ite = klong_manager.begin(); klong_Ite != klong_manager.end(); ++klong_Ite) {
1134 Belle::Mdst_klong mdstKlong = *klong_Ite;
1135 Belle::Mdst_klm_cluster mdstKlm = mdstKlong.klmc();
1147 B2KlmCluster->
setClusterPosition(mdstKlong.cos_x(), mdstKlong.cos_y(), mdstKlong.cos_z());
1154 plist->addParticle(B2Klong);
1165 Belle::Gen_hepevt_Manager& GenMgr = Belle::Gen_hepevt_Manager::get_manager();
1166 const double dang(15. / 180.*M_PI);
1168 for (Belle::Gen_hepevt_Manager::iterator klong_hep_it = GenMgr.begin(); klong_hep_it != GenMgr.end(); ++klong_hep_it) {
1172 CLHEP::HepLorentzVector gp4(klong_hep_it->PX(), klong_hep_it->PY(), klong_hep_it->PZ(), klong_hep_it->E());
1174 int bestRecKlongID(0);
1176 for (Belle::Mdst_klong_Manager::iterator klong_rec_it = klong_manager.begin(); klong_rec_it != klong_manager.end();
1180 if ((*klong_rec_it).ecl())
1182 CLHEP::Hep3Vector klp3(klong_rec_it->cos_x(), klong_rec_it->cos_y(), klong_rec_it->cos_z());
1184 if (cos(gp4.theta() - klp3.theta()) > cos(dang) && cos(gp4.phi() - klp3.phi()) > cos(dang)) {
1186 double tmp_sum = cos(gp4.theta() - klp3.theta()) + cos(gp4.phi() - klp3.phi());
1187 if (tmp_sum > sum) {
1196 particlesToMCParticles.
add(bestRecKlongID, matchedMCParticleID);
1211 Belle::Evtcls_flag_Manager& EvtFlagMgr = Belle::Evtcls_flag_Manager::get_manager();
1212 Belle::Evtcls_flag2_Manager& EvtFlag2Mgr = Belle::Evtcls_flag2_Manager::get_manager();
1215 Belle::Evtcls_hadronic_flag_Manager& EvtHadFlagMgr = Belle::Evtcls_hadronic_flag_Manager::get_manager();
1217 std::string name =
"evtcls_flag";
1218 std::string name_had =
"evtcls_hadronic_flag";
1220 std::vector<Belle::Evtcls_flag>::iterator eflagIterator = EvtFlagMgr.begin();
1221 std::vector<Belle::Evtcls_flag2>::iterator eflag2Iterator = EvtFlag2Mgr.begin();
1222 std::vector<Belle::Evtcls_hadronic_flag>::iterator ehadflagIterator = EvtHadFlagMgr.begin();
1225 std::vector<int> flag(20);
1226 for (
int index = 0; index < 20; ++index) {
1228 if (index == 14 || index == 16)
continue;
1229 std::string iVar = name + std::to_string(index);
1232 m_evtInfo->addExtraInfo(iVar, (*eflagIterator).flag(index));
1235 m_evtInfo->addExtraInfo(iVar, (*eflag2Iterator).flag(index - 10));
1237 B2DEBUG(99,
"evtcls_flag(" << index <<
") = " <<
m_evtInfo->getExtraInfo(iVar));
1241 for (
int index = 0; index < 6; ++index) {
1242 std::string iVar = name_had + std::to_string(index);
1243 m_evtInfo->addExtraInfo(iVar, (*ehadflagIterator).hadronic_flag(index));
1244 B2DEBUG(99,
"evtcls_hadronic_flag(" << index <<
") = " <<
m_evtInfo->getExtraInfo(iVar));
1259 if (
event->getExperiment() <= 27) {
1261 Belle::Rectrg_summary_Manager& RecTrgSummaryMgr = Belle::Rectrg_summary_Manager::get_manager();
1262 std::vector<Belle::Rectrg_summary>::iterator eflagIterator = RecTrgSummaryMgr.begin();
1263 std::string name_summary =
"rectrg_summary_m_final";
1266 for (
int index = 0; index < 2; ++index) {
1267 std::string iVar = name_summary + std::to_string(index);
1268 m_evtInfo->addExtraInfo(iVar, (*eflagIterator).final(index));
1269 B2DEBUG(99,
"m_final(" << index <<
") = " <<
m_evtInfo->getExtraInfo(iVar));
1273 Belle::Rectrg_summary3_Manager& RecTrgSummary3Mgr = Belle::Rectrg_summary3_Manager::get_manager();
1275 std::string name_summary3 =
"rectrg_summary3_m_final";
1277 std::vector<Belle::Rectrg_summary3>::iterator eflagIterator3 = RecTrgSummary3Mgr.begin();
1280 for (
int index = 0; index < 3; ++index) {
1281 std::string iVar = name_summary3 + std::to_string(index);
1282 m_evtInfo->addExtraInfo(iVar, (*eflagIterator3).final(index));
1283 B2DEBUG(99,
"m_final(" << index <<
") = " <<
m_evtInfo->getExtraInfo(iVar));
1297 static Belle::kid_acc acc_pdf(0);
1300 const double pmass[5] = { 0.00051099907, 0.105658389, 0.13956995, 0.493677, 0.93827231 };
1302 CLHEP::Hep3Vector mom(chg.px(), chg.py(), chg.pz());
1303 double cos_theta = mom.cosTheta();
1304 double pval = mom.mag();
1306 double npe = chg.acc().photo_electron();
1307 double beta = pval / sqrt(pval * pval + pmass[idp] * pmass[idp]);
1308 double pdfval = acc_pdf.npe2pdf(cos_theta, beta, npe);
1316 CLHEP::Hep3Vector mom(chg.px(), chg.py(), chg.pz());
1317 double pval = mom.mag();
1319 Belle::kid_cdc kidCdc(5);
1320 float factor0 = kidCdc.factor0();
1321 float factor1 = kidCdc.factor1(idp, pval);
1323 if (factor0 == 1.0 && factor1 == 1.0)
return chg.trk().pid(idp);
1325 double m = chg.trk().dEdx() / factor0;
1326 double e = chg.trk().dEdx_exp(idp) * factor1;
1327 double s = chg.trk().sigma_dEdx(idp);
1328 double val = 1. / sqrt(2.*M_PI) / s * exp(-0.5 * (m - e) * (m - e) / s / s);
1335 bool discard_allzero)
1337 if (discard_allzero) {
1338 const double max_l = *std::max_element(likelihoods, likelihoods +
c_nHyp);
1344 for (
int i = 0; i <
c_nHyp; i++) {
1345 float logl = log(likelihoods[i]);
1355 track->addRelationTo(pid);
1361 double likelihoods[
c_nHyp];
1365 for (
int i = 0; i <
c_nHyp; i++) {
1366 accL[i] = tofL[i] = cdcL[i] = 1.0;
1370 const auto& acc = belleTrack.acc();
1371 if (acc and acc.quality() == 0) {
1372 for (
int i = 0; i <
c_nHyp; i++)
1373 accL[i] = likelihoods[i] =
acc_pid(belleTrack, i);
1380 const Belle::Mdst_tof& tof = belleTrack.tof();
1381 if (tof and tof.quality() == 0) {
1382 for (
int i = 0; i <
c_nHyp; i++)
1383 tofL[i] = likelihoods[i] = tof.pid(i);
1390 const Belle::Mdst_trk& trk = belleTrack.trk();
1391 if (trk.dEdx() > 0) {
1392 for (
int i = 0; i <
c_nHyp; i++) {
1393 likelihoods[i] = trk.pid(i);
1394 cdcL[i] =
cdc_pid(belleTrack, i);
1408 Belle::eid electronID(belleTrack);
1409 float eclID_e_pdf = electronID.pdf_e_ecl();
1410 float eclID_h_pdf = electronID.pdf_h_ecl();
1411 float atcID_e_pdf = electronID.atc_pid_pdf(
true, accL, tofL, cdcL);
1412 float atcID_h_pdf = electronID.atc_pid_pdf(
false, accL, tofL, cdcL);
1415 float eclProb = eclID_e_pdf / (eclID_e_pdf + eclID_h_pdf);
1416 float atcProb = atcID_e_pdf / (atcID_e_pdf + atcID_h_pdf);
1418 if (atcProb > 0.999999) atcProb = 0.999999;
1420 double eidCombinedSig = eclProb * atcProb;
1421 double eidCombinedBkg = (1. - eclProb) * (1. - atcProb);
1423 likelihoods[0] = eidCombinedSig;
1425 likelihoods[2] = eidCombinedBkg;
1439 int muid_trackid = belleTrack.muid_ID();
1443 Belle::Mdst_klm_mu_ex_Manager& ex_mgr = Belle::Mdst_klm_mu_ex_Manager::get_manager();
1444 Belle::Mdst_klm_mu_ex& ex = ex_mgr(Belle::Panther_ID(muid_trackid));
1447 if (ex.Chi_2() > 0) {
1449 likelihoods[1] = ex.Muon_likelihood();
1450 likelihoods[2] = ex.Pion_likelihood();
1451 likelihoods[3] = ex.Kaon_likelihood();
1456 for (
int i = 0; i < 5; i++)
1457 if (likelihoods[i] < 0)
1482 const HepPoint3D& newPivot,
1483 std::vector<float>& helixParams,
1484 CLHEP::HepSymMatrix& error5x5,
1485 CLHEP::HepLorentzVector& momentum,
1486 HepPoint3D& position,
1487 CLHEP::HepSymMatrix& error7x7,
const double dPhi)
1489 const HepPoint3D pivot(trk_fit.pivot_x(),
1493 CLHEP::HepVector a(5);
1494 a[0] = trk_fit.helix(0);
1495 a[1] = trk_fit.helix(1);
1496 a[2] = trk_fit.helix(2);
1497 a[3] = trk_fit.helix(3);
1498 a[4] = trk_fit.helix(4);
1499 CLHEP::HepSymMatrix Ea(5, 0);
1500 Ea[0][0] = trk_fit.error(0);
1501 Ea[1][0] = trk_fit.error(1);
1502 Ea[1][1] = trk_fit.error(2);
1503 Ea[2][0] = trk_fit.error(3);
1504 Ea[2][1] = trk_fit.error(4);
1505 Ea[2][2] = trk_fit.error(5);
1506 Ea[3][0] = trk_fit.error(6);
1507 Ea[3][1] = trk_fit.error(7);
1508 Ea[3][2] = trk_fit.error(8);
1509 Ea[3][3] = trk_fit.error(9);
1510 Ea[4][0] = trk_fit.error(10);
1511 Ea[4][1] = trk_fit.error(11);
1512 Ea[4][2] = trk_fit.error(12);
1513 Ea[4][3] = trk_fit.error(13);
1514 Ea[4][4] = trk_fit.error(14);
1516 Belle::Helix helix(pivot, a, Ea);
1519 if (helix.kappa() > 0)
1524 if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
1525 helix.pivot(newPivot);
1526 momentum = helix.momentum(dPhi, mass, position, error7x7);
1528 if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
1529 helix.pivot(HepPoint3D(0., 0., 0.));
1530 momentum = helix.momentum(dPhi, mass, position, error7x7);
1532 momentum = helix.momentum(dPhi, mass, position, error7x7);
1542 const HepPoint3D& newPivot,
1543 std::vector<float>& helixParams, std::vector<float>& helixError)
1545 const HepPoint3D pivot(trk_fit.pivot_x(),
1549 CLHEP::HepVector a(5);
1550 a[0] = trk_fit.helix(0);
1551 a[1] = trk_fit.helix(1);
1552 a[2] = trk_fit.helix(2);
1553 a[3] = trk_fit.helix(3);
1554 a[4] = trk_fit.helix(4);
1555 CLHEP::HepSymMatrix Ea(5, 0);
1556 Ea[0][0] = trk_fit.error(0);
1557 Ea[1][0] = trk_fit.error(1);
1558 Ea[1][1] = trk_fit.error(2);
1559 Ea[2][0] = trk_fit.error(3);
1560 Ea[2][1] = trk_fit.error(4);
1561 Ea[2][2] = trk_fit.error(5);
1562 Ea[3][0] = trk_fit.error(6);
1563 Ea[3][1] = trk_fit.error(7);
1564 Ea[3][2] = trk_fit.error(8);
1565 Ea[3][3] = trk_fit.error(9);
1566 Ea[4][0] = trk_fit.error(10);
1567 Ea[4][1] = trk_fit.error(11);
1568 Ea[4][2] = trk_fit.error(12);
1569 Ea[4][3] = trk_fit.error(13);
1570 Ea[4][4] = trk_fit.error(14);
1572 Belle::Helix helix(pivot, a, Ea);
1574 if (newPivot.x() != 0. || newPivot.y() != 0. || newPivot.z() != 0.) {
1575 helix.pivot(newPivot);
1577 if (pivot.x() != 0. || pivot.y() != 0. || pivot.z() != 0.) {
1578 helix.pivot(HepPoint3D(0., 0., 0.));
1582 CLHEP::HepSymMatrix error5x5(5, 0);
1585 unsigned int size = 5;
1586 unsigned int counter = 0;
1587 for (
unsigned int i = 0; i < size; i++)
1588 for (
unsigned int j = i; j < size; j++)
1589 helixError[counter++] = error5x5[i][j];
1594 CLHEP::HepVector a(5);
1595 CLHEP::HepSymMatrix Ea(5, 0);
1601 helixParams[0] = a[0];
1604 helixParams[1] = adjustAngleRange(a[1] + TMath::Pi() / 2.0);
1610 helixParams[3] = a[3];
1613 helixParams[4] = a[4];
1615 unsigned int size = 5;
1616 for (
unsigned int i = 0; i < size; i++) {
1617 for (
unsigned int j = 0; j < size; j++) {
1618 error5x5[i][j] = Ea[i][j];
1624 if (std::isinf(error5x5[i][j])) {
1625 B2DEBUG(99,
"Helix covariance matrix element found to be infinite. Setting value to DBL_MAX/2.0.");
1626 error5x5[i][j] = DBL_MAX / 2.0;
1634 Belle::Mdst_trk& trk = belleTrack.trk();
1636 for (
int mhyp = 0 ; mhyp <
c_nHyp; ++mhyp) {
1638 double thisMass = pType.
getMass();
1640 Belle::Mdst_trk_fit& trk_fit = trk.mhyp(mhyp);
1643 std::vector<float> helixParam(5);
1645 CLHEP::HepSymMatrix error5x5(5, 0);
1647 CLHEP::HepLorentzVector momentum;
1649 CLHEP::HepSymMatrix error7x7(7, 0);
1651 HepPoint3D position;
1654 helixParam, error5x5,
1655 momentum, position, error7x7, 0.0);
1657 std::vector<float> helixError(15);
1658 unsigned int size = 5;
1659 unsigned int counter = 0;
1660 for (
unsigned int i = 0; i < size; i++)
1661 for (
unsigned int j = i; j < size; j++)
1662 helixError[counter++] = error5x5[i][j];
1664 double pValue = TMath::Prob(trk_fit.chisq(), trk_fit.ndf());
1671 for (
unsigned int i = 0; i < 3; i++)
1672 cdcNHits += trk_fit.nhits(i);
1680 double path_length = 0;
1681 double tof_sigma = 0;
1685 double dedx = trk.dEdx();
1686 short dedx_qual = trk.quality_dedx();
1688 const Belle::Mdst_tof& tof_obj = belleTrack.tof();
1690 tof = tof_obj.tof();
1691 path_length = tof_obj.path_length();
1692 tof_qual = tof_obj.quality();
1693 tof_sigma = tof_obj.sigma_tof();
1696 const Belle::Mdst_acc& acc_obj = belleTrack.acc();
1698 acc_ph = acc_obj.photo_electron();
1699 acc_qual = acc_obj.quality();
1703 auto cdcExtraInfo =
m_belleTrkExtra.appendNew(trk_fit.first_x(), trk_fit.first_y(), trk_fit.first_z(),
1704 trk_fit.last_x(), trk_fit.last_y(), trk_fit.last_z(),
1705 tof, path_length, tof_qual, tof_sigma,
1706 acc_ph, acc_qual, dedx, dedx_qual);
1707 track->addRelationTo(cdcExtraInfo);
1710 int svdHitPattern = trk_fit.hit_svd();
1714 std::bitset<32> svdBitSet(svdHitPattern);
1718 unsigned short svdLayers;
1722 std::bitset<32> svdUMask(
static_cast<std::string
>(
"00000000000000000000000000000011"));
1724 std::bitset<32> svdVMask;
1727 if (
event->getExperiment() <= 27) {
1728 svdVMask = svdUMask << 6;
1731 svdVMask = svdUMask << 8;
1736 for (
unsigned short layerId = 0; layerId < svdLayers; layerId++) {
1737 unsigned short uHits = (svdBitSet & svdUMask).count();
1738 unsigned short vHits = (svdBitSet & svdVMask).count();
1739 patternVxd.
setSVDLayer(layerId + 3, uHits, vHits);
1748 TMatrixDSym cartesianCovariance(6);
1749 for (
unsigned i = 0; i < 7; i++) {
1752 for (
unsigned j = 0; j < 7; j++) {
1760 BFIELD, cartesianCovariance, pValue);
1762 TMatrixDSym helixCovariance = helixFromCartesian.
getCovariance();
1765 for (
unsigned int i = 0; i < 5; ++i)
1766 for (
unsigned int j = i; j < 5; ++j)
1767 helixError[counter++] = helixCovariance(i, j);
1772 track->setTrackFitResultIndex(pType, trackFit->getArrayIndex());
1800 B2WARNING(
"Trying to convert Gen_hepevt with idhep = " << idHep <<
1801 ". This should never happen.");
1804 mcParticle->
setPDG(idHep);
1807 if (genHepevt.isthep() > 0) {
1811 mcParticle->
setMass(genHepevt.M());
1813 ROOT::Math::PxPyPzEVector p4(genHepevt.PX(), genHepevt.PY(), genHepevt.PZ(), genHepevt.E());
1820 if (genHepevt.daFirst() > 0) {
1821 Belle::Gen_hepevt_Manager& genMgr = Belle::Gen_hepevt_Manager::get_manager();
1822 Belle::Gen_hepevt daughterParticle = genMgr(Belle::Panther_ID(genHepevt.daFirst()));
1827 mcParticle->
setDecayTime(std::numeric_limits<float>::infinity());
1842 eclCluster->
setPhi(ecl.phi());
1844 eclCluster->
setR(ecl.r());
1847 double covarianceMatrix[6];
1848 covarianceMatrix[0] = ecl.error(0);
1849 covarianceMatrix[1] = ecl.error(1);
1850 covarianceMatrix[2] = ecl.error(2);
1851 covarianceMatrix[3] = ecl.error(3);
1852 covarianceMatrix[4] = ecl.error(4);
1853 covarianceMatrix[5] = ecl.error(5);
1856 eclCluster->
setLAT(eclAux.width());
1860 eclCluster->
setTime(eclAux.property(0));
1867 klmCluster->
setLayers(klm_cluster.layers());
1880 Belle::Mdst_ecl_trk_Manager& m = Belle::Mdst_ecl_trk_Manager::get_manager();
1881 Belle::Mdst_charged_Manager& chgMg = Belle::Mdst_charged_Manager::get_manager();
1886 std::vector<int> insert_order_types = {1, 2, 0};
1887 for (
auto& insert_type : insert_order_types) {
1888 for (Belle::Mdst_ecl_trk_Manager::iterator ecltrkIterator = m.begin(); ecltrkIterator != m.end(); ++ecltrkIterator) {
1889 Belle::Mdst_ecl_trk mECLTRK = *ecltrkIterator;
1891 if (mECLTRK.type() != insert_type)
1894 Belle::Mdst_ecl mdstEcl = mECLTRK.ecl();
1895 Belle::Mdst_trk mTRK = mECLTRK.trk();
1903 for (Belle::Mdst_charged_Manager::iterator chgIterator = chgMg.begin(); chgIterator != chgMg.end(); ++chgIterator) {
1904 Belle::Mdst_charged mChar = *chgIterator;
1905 Belle::Mdst_trk mTRK_in_charged = mChar.trk();
1907 if (mTRK_in_charged.get_ID() == mTRK.get_ID()) {
1909 tracksToECLClusters.
add(mChar.get_ID() - 1, mdstEcl.get_ID() - 1, 1.0);
1924 Belle::Mdst_klm_cluster_Manager& klm_cluster_manager = Belle::Mdst_klm_cluster_Manager::get_manager();
1927 for (Belle::Mdst_klm_cluster_Manager::iterator klmC_Ite = klm_cluster_manager.begin(); klmC_Ite != klm_cluster_manager.end();
1930 Belle::Mdst_klm_cluster mdstKlm_cluster = *klmC_Ite;
1931 Belle::Mdst_trk mTRK = mdstKlm_cluster.trk();
1932 Belle::Mdst_ecl mECL = mdstKlm_cluster.ecl();
1934 if (mTRK) klmClustersToTracks.
add(mdstKlm_cluster.get_ID() - 1, mTRK.get_ID() - 1);
1935 if (mECL) klmClustersToEclClusters.
add(mdstKlm_cluster.get_ID() - 1, mECL.get_ID() - 1);
1952 const int mask = 0x00f00000;
1953 int high_bits =
id & mask;
1954 if (high_bits == 0 || high_bits == mask)
return id;
2014 int bellePDGCode = belleMC.idhep();
2015 int belleIIPDGCode = mcP->
getPDG();
2017 if (bellePDGCode == 0)
2018 B2WARNING(
"[B2BIIConvertMdstModule] " << objectName <<
" matched to Gen_hepevt with idhep = 0.");
2020 if (bellePDGCode != belleIIPDGCode)
2021 B2WARNING(
"[B2BIIConvertMdstModule] " << objectName <<
" matched to different MCParticle! " << bellePDGCode <<
" vs. " <<
2024 const double belleMomentum[] = { belleMC.PX(), belleMC.PY(), belleMC.PZ() };
2027 for (
unsigned i = 0; i < 3; i++) {
2028 double relDev = (belle2Momentum[i] - belleMomentum[i]) / belleMomentum[i];
2030 if (relDev > 1e-3) {
2031 B2WARNING(
"[B2BIIConvertMdstModule] " << objectName <<
" matched to different MCParticle!");
2032 B2INFO(
" - Gen_hepevt [" << bellePDGCode <<
"] px/py/pz = " << belleMC.PX() <<
"/" << belleMC.PY() <<
"/" << belleMC.PZ());
2033 B2INFO(
" - TrackFitResult [" << belleIIPDGCode <<
"] px/py/pz = " << mcP->
get4Vector().Px() <<
"/" << mcP->
get4Vector().Py() <<
"/"
2041 CLHEP::HepLorentzVector& momentum, HepPoint3D& position, CLHEP::HepSymMatrix& error)
2043 const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2044 CLHEP::HepVector a(5);
2045 CLHEP::HepSymMatrix Ea(5, 0);
2047 a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2048 a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2049 a[4] = vee.daut().helix_p(4);
2050 Ea[0][0] = vee.daut().error_p(0); Ea[1][0] = vee.daut().error_p(1);
2051 Ea[1][1] = vee.daut().error_p(2); Ea[2][0] = vee.daut().error_p(3);
2052 Ea[2][1] = vee.daut().error_p(4); Ea[2][2] = vee.daut().error_p(5);
2053 Ea[3][0] = vee.daut().error_p(6); Ea[3][1] = vee.daut().error_p(7);
2054 Ea[3][2] = vee.daut().error_p(8); Ea[3][3] = vee.daut().error_p(9);
2055 Ea[4][0] = vee.daut().error_p(10); Ea[4][1] = vee.daut().error_p(11);
2056 Ea[4][2] = vee.daut().error_p(12); Ea[4][3] = vee.daut().error_p(13);
2057 Ea[4][4] = vee.daut().error_p(14);
2059 a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2060 a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2061 a[4] = vee.daut().helix_m(4);
2062 Ea[0][0] = vee.daut().error_m(0); Ea[1][0] = vee.daut().error_m(1);
2063 Ea[1][1] = vee.daut().error_m(2); Ea[2][0] = vee.daut().error_m(3);
2064 Ea[2][1] = vee.daut().error_m(4); Ea[2][2] = vee.daut().error_m(5);
2065 Ea[3][0] = vee.daut().error_m(6); Ea[3][1] = vee.daut().error_m(7);
2066 Ea[3][2] = vee.daut().error_m(8); Ea[3][3] = vee.daut().error_m(9);
2067 Ea[4][0] = vee.daut().error_m(10); Ea[4][1] = vee.daut().error_m(11);
2068 Ea[4][2] = vee.daut().error_m(12); Ea[4][3] = vee.daut().error_m(13);
2069 Ea[4][4] = vee.daut().error_m(14);
2072 Belle::Helix helix(pivot, a, Ea);
2075 momentum = helix.momentum(0., pType.
getMass(), position, error);
2079 std::vector<float>& helixError)
2081 const HepPoint3D pivot(vee.vx(), vee.vy(), vee.vz());
2082 CLHEP::HepVector a(5);
2083 CLHEP::HepSymMatrix Ea(5, 0);
2085 a[0] = vee.daut().helix_p(0); a[1] = vee.daut().helix_p(1);
2086 a[2] = vee.daut().helix_p(2); a[3] = vee.daut().helix_p(3);
2087 a[4] = vee.daut().helix_p(4);
2088 Ea[0][0] = vee.daut().error_p(0);
2089 Ea[1][0] = vee.daut().error_p(1);
2090 Ea[1][1] = vee.daut().error_p(2);
2091 Ea[2][0] = vee.daut().error_p(3);
2092 Ea[2][1] = vee.daut().error_p(4);
2093 Ea[2][2] = vee.daut().error_p(5);
2094 Ea[3][0] = vee.daut().error_p(6);
2095 Ea[3][1] = vee.daut().error_p(7);
2096 Ea[3][2] = vee.daut().error_p(8);
2097 Ea[3][3] = vee.daut().error_p(9);
2098 Ea[4][0] = vee.daut().error_p(10);
2099 Ea[4][1] = vee.daut().error_p(11);
2100 Ea[4][2] = vee.daut().error_p(12);
2101 Ea[4][3] = vee.daut().error_p(13);
2102 Ea[4][4] = vee.daut().error_p(14);
2104 a[0] = vee.daut().helix_m(0); a[1] = vee.daut().helix_m(1);
2105 a[2] = vee.daut().helix_m(2); a[3] = vee.daut().helix_m(3);
2106 a[4] = vee.daut().helix_m(4);
2107 Ea[0][0] = vee.daut().error_m(0);
2108 Ea[1][0] = vee.daut().error_m(1);
2109 Ea[1][1] = vee.daut().error_m(2);
2110 Ea[2][0] = vee.daut().error_m(3);
2111 Ea[2][1] = vee.daut().error_m(4);
2112 Ea[2][2] = vee.daut().error_m(5);
2113 Ea[3][0] = vee.daut().error_m(6);
2114 Ea[3][1] = vee.daut().error_m(7);
2115 Ea[3][2] = vee.daut().error_m(8);
2116 Ea[3][3] = vee.daut().error_m(9);
2117 Ea[4][0] = vee.daut().error_m(10);
2118 Ea[4][1] = vee.daut().error_m(11);
2119 Ea[4][2] = vee.daut().error_m(12);
2120 Ea[4][3] = vee.daut().error_m(13);
2121 Ea[4][4] = vee.daut().error_m(14);
2124 Belle::Helix helix(pivot, a, Ea);
2127 helix.pivot(HepPoint3D(0., 0., 0.));
2129 CLHEP::HepSymMatrix error5x5(5, 0);
2132 unsigned int size = 5;
2133 unsigned int counter = 0;
2134 for (
unsigned int i = 0; i < size; i++)
2135 for (
unsigned int j = i; j < size; j++)
2136 helixError[counter++] = error5x5[i][j];
2140 const HepPoint3D& position,
2141 const CLHEP::HepSymMatrix& error,
2142 const short int charge,
2145 const uint64_t hitPatternCDCInitializer,
2146 const uint32_t hitPatternVXDInitializer,
2149 ROOT::Math::XYZVector pos(position.x(), position.y(), position.z());
2150 ROOT::Math::XYZVector mom(momentum.px(), momentum.py(), momentum.pz());
2152 TMatrixDSym errMatrix(6);
2153 for (
unsigned i = 0; i < 7; i++) {
2156 for (
unsigned j = 0; j < 7; j++) {
2167 return TrackFitResult(pos, mom, errMatrix, charge, pType, pValue,
BFIELD, hitPatternCDCInitializer, hitPatternVXDInitializer, ndf);
2172 if (!pid)
return 0.5;
2179 if (acc_sig + acc_bkg > 0.0)
2180 acc = acc_sig / (acc_sig + acc_bkg);
2187 double tof_all = tof_sig + tof_bkg;
2189 tof = tof_sig / tof_all;
2190 if (tof < 0.001) tof = 0.001;
2191 if (tof > 0.999) tof = 0.999;
2199 double cdc_all = cdc_sig + cdc_bkg;
2201 cdc = cdc_sig / cdc_all;
2202 if (cdc < 0.001) cdc = 0.001;
2203 if (cdc > 0.999) cdc = 0.999;
2207 double pid_sig = acc * tof * cdc;
2208 double pid_bkg = (1. - acc) * (1. - tof) * (1. - cdc);
2210 return pid_sig / (pid_sig + pid_bkg);
2216 B2DEBUG(99,
"B2BIIConvertMdst: endRun done.");
2222 B2DEBUG(99,
"B2BIIConvertMdst: terminate called");
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 ...
B2BIIConvertMdstModule()
Constructor.
std::map< int, int > mdstKlmToKLMCluster
map of Mdst_klm Panther IDs and corresponding KLMCluster StoreArray indices
void convertMdstChargedTable()
Reads and converts all entries of Mdst_charged (Mdst_trk and Mdst_trk_fit) Panther table to Track (Tr...
void setTracksToECLClustersRelations()
Sets Track -> ECLCluster relations.
OptionalDBObjPtr< CollisionBoostVector > m_collisionBoostVectorDB
CollisionBoostVector for boost vector.
StoreObjPtr< EventExtraInfo > m_evtInfo
Event Extra Info.
const int ERRMCONV[7]
CONVERSION OF TRACK ERROR MATRIX ELEMENTS.
StoreArray< KLMCluster > m_klmClusters
KLM clusters.
void belleVeeDaughterToCartesian(const Belle::Mdst_vee2 &vee, const int charge, const Const::ParticleType &pType, CLHEP::HepLorentzVector &momentum, HepPoint3D &position, CLHEP::HepSymMatrix &error)
Fills 4-momentum, position and 7x7 error matrix from Belle Vee daughter.
@ c_Direct
Direct matching.
@ c_GeneratorLevel
Match to generator-level particles.
virtual void initialize() override
Initialize the module.
void setLikelihoods(PIDLikelihood *pid, Const::EDetector det, double likelihoods[c_nHyp], bool discard_allzero=false)
Add given Belle likelihoods (not log-likelihoods, in Belle hypothesis order) for given detector to pi...
StoreArray< V0 > m_v0s
V0-particles.
virtual void event() override
Called for each event.
void convertGenHepEvtTable()
Reads and converts all entries of Gen_hepevt Panther table to MCParticle dataobjects and adds them to...
double cdc_pid(const Belle::Mdst_charged &chg, int idp)
Returns CDC likelihood for given hypothesis idp.
void convertRecTrgTable()
Reads and converts m_final from rectrg_summary3.
void setKLMClustersRelations()
Sets KLMCluster -> Track and ECLCluster relations.
bool m_convertEvtcls
Flag to switch on conversion of Evtcls table.
const double KAPPA2OMEGA
Conversion factor for Kappa -> Omega helix parameters.
virtual void endRun() override
Called when the current run is finished.
bool m_nisEnable
Flag to switch on conversion of nisKsFinder info.
StoreArray< TrackFitResult > m_trackFitResults
Track fir results.
BeamSpot m_beamSpot
Interaction Point of the beam.
StoreArray< Particle > m_particles
Particles.
virtual void terminate() override
Terminates the module.
OptionalDBObjPtr< CollisionInvariantMass > m_collisionInvMDB
CollisionInvariantMass for Invariant Mass of Beam.
void convertHelix(Belle::Helix &helix, std::vector< float > &helixParams, CLHEP::HepSymMatrix &error5x5)
Converts Belle's Helix parameters and it's covariance matrix to Belle II's version.
void convertMdstPi0Table()
Reads all entries of Mdst_Pi0 Panther table, creates a particle list 'pi0:mdst' and adds them to Stor...
std::map< int, int > genHepevtToMCParticle
map of Gen_hepevt Panther IDs and corresponding MCParticle StoreArray indices
void testMCRelation(const Belle::Gen_hepevt &belleMC, const MCParticle *mcP, const std::string &objectName)
Checks if the reconstructed object (Track, ECLCluster, ...) was matched to the same MCParticle.
void convertIPProfile(bool beginRun=false)
Stores the IPProfiles in BeamSpot (currently in DataStore)
void convertMdstECLObject(const Belle::Mdst_ecl &ecl, const Belle::Mdst_ecl_aux &eclAux, ECLCluster *eclCluster)
Converts Mdst_ecl(_aux) record to ECLCluster object.
CollisionInvariantMass m_collisionInvM
CollisionInvariantMass for the invariant mass of the beam.
bool m_convertBeamParameters
Convert beam parameters or use information stored in Belle II database.
static double acc_pid(const Belle::Mdst_charged &chg, int idp)
Returns ACC likelihood for given hypothesis idp.
virtual void beginRun() override
Module functions to be called from event process.
StoreArray< Track > m_tracks
Tracks.
void convertBeamEnergy()
Stores beam parameters (energy, angles) in CollisionInvariantMass and CollisionBoostVector (currently...
void convertPIDData(const Belle::Mdst_charged &belleTrack, const Track *track)
Get PID information for belleTrack and add it to PIDLikelihood (with relation from track).
void convertMdstKLMObject(const Belle::Mdst_klm_cluster &klm, KLMCluster *klmCluster)
Converts Mdst_klm_cluster record to KLMCluster object.
Belle2::MCParticleGraph m_particleGraph
MCParticle Graph to build Belle2 MC Particles.
void convertMdstKLongTable()
Reads all entries of Mdst_Klong Panther table, creates a particle list 'K_L0:mdst' and adds them to S...
OptionalDBObjPtr< BeamSpot > m_beamSpotDB
BeamSpot for IP.
bool m_realData
flag that tells whether given data sample is for real data or MC
StoreArray< ECLCluster > m_eclClusters
ECL clusters.
void convertGenHepevtObject(const Belle::Gen_hepevt &genHepevt, MCParticleGraph::GraphParticle *mcParticle)
Converts Gen_hepevt record to MCParticleGraph object.
void initializeDataStore()
Initializes Belle II DataStore.
void convertMdstECLTable()
Reads and converts all entries of Mdst_ecl(_aux) Panther table to ECLCluster dataobjects and adds the...
CollisionBoostVector m_collisionBoostVector
CollisionBoostVector for bosst vector of the beam.
std::string m_mcMatchingModeString
MC matching mode.
int m_lastIPProfileBin
variable to tell us which IPProfile bin was active last time we looked
void convertMdstKLMTable()
Reads and converts all entries of Mdst_klm_cluster Panther table to KLMCluster dataobjects and adds t...
bool m_convertRecTrg
Flag to switch on conversion of rectrg_summary3.
std::map< int, int > mdstGammaToParticle
map of gamma Panther IDs and corresponding Particle StoreArray indices
static const Const::ChargedStable c_belleHyp_to_chargedStable[c_nHyp]
maps Belle hypotheses to Const::ChargedStable (from http://belle.kek.jp/secured/wiki/doku....
std::map< int, int > mdstKlongToParticle
map of Klong Panther IDs and corresponding Particle StoreArray indices
double m_matchType2E9oE25Threshold
E9/E25 threshold value clusters with a value above this threshold are classified as neutral even if t...
void convertMdstChargedObject(const Belle::Mdst_charged &belleTrack, Track *track)
Converts Mdst_charged (Mdst_trk(_fit)) record to Track (TrackFitResult) object.
StoreArray< MCParticle > m_mcParticles
MC particles.
void convertMdstVee2Table()
Reads and converts all entries of Mdst_vee2 Panther table to V0 dataobjects and adds them to StoreArr...
StoreArray< PIDLikelihood > m_pidLikelihoods
output PIDLikelihood array.
std::map< int, int > mdstEclToECLCluster
map of Mdst_ecl Panther IDs and corresponding ECLCluster StoreArray indices
TrackFitResult createTrackFitResult(const CLHEP::HepLorentzVector &momentum, const HepPoint3D &position, const CLHEP::HepSymMatrix &error, const short int charge, const Const::ParticleType &pType, const float pValue, const uint64_t hitPatternCDCInitializer, const uint32_t hitPatternVXDInitializer, const uint16_t ndf)
Creates TrackFitResult and fills it.
virtual ~B2BIIConvertMdstModule() override
Destructor.
const double BFIELD
B filed in TESLA.
MCMatchingMode m_mcMatchingMode
C matching mode.
StoreArray< BelleTrkExtra > m_belleTrkExtra
Belle CDC extra information.
bool m_use6x6CovarianceMatrix4Tracks
flag that tells which form of covariance matrix should be used in the conversion of charged tracks
int recoverMoreThan24bitIDHEP(int id)
Helper function to recover falsely set idhep info in GenHepEvt list.
bool m_convertTrkExtra
Flag to switch on conversion of first(last)_{x,y,z} of mdst_trk_fit.
int getHelixParameters(const Belle::Mdst_trk_fit &trk_fit, const double mass, const HepPoint3D &newPivot, std::vector< float > &helixParams, CLHEP::HepSymMatrix &error5x5, CLHEP::HepLorentzVector &momentum, HepPoint3D &position, CLHEP::HepSymMatrix &error7x7, const double dPhi=0.0)
Fills Helix parameters (converted to Belle II version), 5x5 error matrix, 4-momentum,...
static const int c_nHyp
Number of Belle track hypotheses (see c_belleHyp_to_chargedStable).
void convertEvtclsTable()
Reads and converts all entries of evtcls Panther table.
This class contains the beam spot position and size modeled as a gaussian distribution in space.
void setIP(const TVector3 &ipPosition, const TMatrixDSym &covariance)
Set the IP position and its error matrix.
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.
The ParticleType class for identifying different particle types.
int getPDGCode() const
PDG code.
double getMass() const
Particle mass.
A class for sets of detector IDs whose content is limited to restricted set of valid detector IDs.
static const ChargedStable muon
muon particle
EDetector
Enum for identifying the detector components (detector and subdetector).
static const ChargedStable pion
charged pion particle
static const ParticleType Klong
K^0_L particle.
static const double speedOfLight
[cm/ns]
static const double electronMass
electron mass
static const ChargedStable proton
proton particle
static const ChargedStable kaon
charged kaon particle
static const ParticleType photon
photon particle
static const ChargedStable electron
electron particle
static const ChargedStable deuteron
deuteron particle
void setE9oE21(double E9oE21)
Set E9/E21 energy ratio.
void setTheta(double theta)
Set Theta of Shower (radian).
void setPhi(double phi)
Set Phi of Shower (radian).
void setdeltaL(double deltaL)
Set deltaL for shower shape.
void setEnergyHighestCrystal(double energyhighestcrystal)
Set energy of highest energetic crystal (GeV).
void setEnergyRaw(double energyraw)
Set Uncorrect Energy deposited (GeV).
void setTime(double time)
Set time information.
void setCovarianceMatrix(double covArray[6])
Set covariance matrix (3x3), i.e.
void setLAT(double LAT)
Set Lateral distribution parameter.
void setNumberOfCrystals(double noc)
Set number of crystals (sum of weights).
void setEnergy(double energy)
Set Corrected Energy (GeV).
void setIsTrack(bool istrack)
Set m_isTrack true if the cluster matches with a track.
void setR(double r)
Set R (in cm).
static double getAlpha(const double bZ)
Calculates the alpha value for a given magnetic field in Tesla.
Hit pattern of CDC hits within a track.
void setNHits(unsigned short nHits)
Sets the 8 MSBs to the total number of hits in the CDC.
ULong64_t getInteger() const
Getter for underlying integer type.
Hit pattern of the VXD within a track.
unsigned int getInteger() const
Getter for the underlying integer.
void setSVDLayer(const unsigned short layerId, unsigned short uHits, unsigned short vHits)
Set the number of hits in a specific layer of the SVD.
A class that describes the interval of experiments/runs for which an object in the database is valid.
void setLayers(int layers)
Set number of layers with hits.
void setClusterPosition(float globalX, float globalY, float globalZ)
Set global position.
void setInnermostLayer(int innermostLayer)
Set number of the innermost layer with hits.
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.
@ c_PrimaryParticle
bit 0: Particle is primary particle.
void setDecayTime(float time)
Set decay time.
void setMass(float mass)
Set particle mass.
void setDecayVertex(const ROOT::Math::XYZVector &vertex)
Set decay vertex.
void setValidVertex(bool valid)
Set indication wether vertex and time information is valid or just default.
void setProductionVertex(const ROOT::Math::XYZVector &vertex)
Set production vertex position.
ROOT::Math::PxPyPzEVector get4Vector() const
Return 4Vector of particle.
void setPDG(int pdg)
Set PDG code of the particle.
int getPDG() const
Return PDG code of particle.
void set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets the 4Vector of particle.
void setStatus(unsigned short int status)
Set Status code for the particle.
void setProductionTime(float time)
Set production time.
void setDescription(const std::string &description)
Sets the description of the module.
Class to collect log likelihoods from TOP, ARICH, dEdx, ECL and KLM aimed for output to mdst includes...
Class to store reconstructed particles.
void appendDaughter(const Particle *daughter, const bool updateType=true, const int daughterProperty=c_Ordinary)
Appends index of daughter to daughters index array.
void setVertex(const ROOT::Math::XYZVector &vertex)
Sets position (decay vertex)
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
void setPValue(double pValue)
Sets chi^2 probability of fit.
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.
static bool parallelProcessingUsed()
Returns true if multiple processes have been spawned, false in single-core mode.
Low-level class to create/modify relations between StoreArrays.
void add(index_type from, index_type to, weight_type weight=1.0)
Add a new element to the relation.
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
const std::string & getName() const
Return name under which the object is saved in the DataStore.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
bool create(bool replace=false)
Create a default object in the data store.
Type-safe access to single objects in the data store.
Values of the result of a track fit with a given particle hypothesis.
Helix getHelix() const
Conversion to framework Helix (without covariance).
TMatrixDSym getCovariance5() const
Getter for covariance matrix of perigee parameters in matrix form.
short getChargeSign() const
Return track charge (1 or -1).
double getOmega() const
Getter for omega.
double getD0() const
Getter for d0.
double getTanLambda() const
Getter for tanLambda.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
ROOT::Math::XYZVector getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
double getPhi0() const
Getter for phi0.
Class that bundles various TrackFitResults.
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]
REG_MODULE(B2BIIConvertBeamParams)
Register the module.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
static Database & Instance()
Instance of a singleton Database.
static DBStore & Instance()
Instance of a singleton DBStore.
bool storeData(const std::string &name, TObject *object, const IntervalOfValidity &iov)
Store an object in the database.
void addConstantOverride(const std::string &name, TObject *obj, bool oneRun=false)
Add constant override payload.
B2Vector3< double > B2Vector3D
typedef for common usage with double
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.
Abstract base class for different kinds of events.