11 #include <analysis/dataobjects/RestOfEvent.h>
13 #include <framework/datastore/StoreArray.h>
15 #include <analysis/dataobjects/Particle.h>
16 #include <mdst/dataobjects/MCParticle.h>
17 #include <mdst/dataobjects/Track.h>
18 #include <mdst/dataobjects/ECLCluster.h>
19 #include <mdst/dataobjects/KLMCluster.h>
21 #include <analysis/ClusterUtility/ClusterUtils.h>
23 #include <TLorentzVector.h>
30 for (
auto* particleToAdd : particlesToAdd) {
31 std::vector<const Particle*> daughters = particleToAdd->getFinalStateDaughters();
32 for (
auto* daughter : daughters) {
35 if (allParticles[myIndex]->isCopyOf(daughter,
true)) {
41 B2DEBUG(10,
"\t\tAdding particle with PDG " << daughter->getPDGCode());
50 std::vector<const Particle*> result;
54 B2DEBUG(10,
"ROE contains no particles, masks are empty too");
61 bool maskFound =
false;
63 if (mask.getName() == maskName) {
65 source = mask.getParticles();
70 B2FATAL(
"No " << maskName <<
" mask defined in current ROE!");
73 for (
const int index : source) {
74 if ((allParticles[index]->getParticleSource() == Particle::EParticleSourceObject::c_Composite or
75 allParticles[index]->getParticleSource() == Particle::EParticleSourceObject::c_V0) && unpackComposite) {
76 auto fsdaughters = allParticles[index]->getFinalStateDaughters();
77 for (
auto* daughter : fsdaughters) {
78 result.push_back(daughter);
82 result.push_back(allParticles[index]);
88 auto particles =
getParticles(maskName, unpackComposite);
89 std::vector<const Particle*> photons;
90 for (
auto* particle : particles) {
91 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_ECLCluster) {
92 photons.push_back(particle);
99 auto particles =
getParticles(maskName, unpackComposite);
100 std::vector<const Particle*> hadrons;
101 for (
auto* particle : particles) {
102 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_KLMCluster) {
103 hadrons.push_back(particle);
110 bool unpackComposite)
const
112 auto particles =
getParticles(maskName, unpackComposite);
113 std::vector<const Particle*> charged;
114 for (
auto* particle : particles) {
115 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_Track) {
116 if (pdg == 0 || pdg == abs(particle->getPDGCode())) {
117 charged.push_back(particle);
127 if (maskName !=
"" && !
hasMask(maskName)) {
128 B2FATAL(
"No " << maskName <<
" mask defined in current ROE!");
131 std::vector<const Particle*> particlesROE =
getParticles(maskName);
138 B2FATAL(
"Creation of ROE Mask with an empty name is not allowed!");
141 B2FATAL(
"ROE Mask already exists!");
143 Mask elon(name, origin);
152 B2FATAL(
"No " << maskName <<
" mask defined in current ROE!");
154 std::string maskNameToGetParticles = maskName;
155 if (!mask->isValid()) {
156 maskNameToGetParticles =
"";
158 std::vector<const Particle*> allROEParticles =
getParticles(maskNameToGetParticles);
159 std::vector<const Particle*> toKeepinROE;
160 for (
auto* roeParticle : allROEParticles) {
164 toKeepinROE.push_back(roeParticle);
168 if (listType != roeParticle->getParticleSource()) {
169 toKeepinROE.push_back(roeParticle);
170 }
else if (discard) {
172 toKeepinROE.push_back(roeParticle);
176 mask->clearParticles();
177 mask->addParticles(toKeepinROE);
181 const std::shared_ptr<Variable::Cut>& eclCut,
const std::shared_ptr<Variable::Cut>& klmCut,
bool updateExisting)
185 B2FATAL(
"ROE Mask does not exist!");
187 std::string sourceName =
"";
188 if (updateExisting) {
190 sourceName = maskName;
193 std::vector<const Particle*> allROEParticles =
getParticles(sourceName,
false);
194 std::vector<const Particle*> maskedParticles;
196 for (
auto* particle : allROEParticles) {
197 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_Track && (!trackCut || trackCut->
check(particle))) {
198 maskedParticles.push_back(particle);
200 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_ECLCluster && (!eclCut || eclCut->
check(particle))) {
201 maskedParticles.push_back(particle);
203 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_KLMCluster && (!klmCut || klmCut->
check(particle))) {
204 maskedParticles.push_back(particle);
207 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_Composite or
208 particle->getParticleSource() == Particle::EParticleSourceObject::c_V0) {
209 maskedParticles.push_back(particle);
212 mask->clearParticles();
213 mask->addParticles(maskedParticles);
220 B2FATAL(
"ROE Mask does not exist!");
222 std::vector<const Particle*> allROEParticles =
getParticles(name,
false);
223 std::vector<int> indicesToErase;
225 for (
auto* maskParticle : allROEParticles) {
227 for (
auto* daughterV0 : daughtersV0) {
228 if (daughterV0->isCopyOf(maskParticle,
true)) {
233 indicesToErase.push_back(maskParticle->getArrayIndex());
236 if (daughtersV0.size() != indicesToErase.size()) {
237 B2DEBUG(10,
"Only " << indicesToErase.size() <<
" daughters are excluded from mask particles. Abort");
240 std::string toprint =
"We will erase next indices from " + name +
" mask: ";
241 for (
auto& i : indicesToErase) {
242 toprint += std::to_string(i) +
" ";
244 B2DEBUG(10, toprint);
246 mask->addV0(particleV0, indicesToErase);
253 B2FATAL(
"ROE Mask does not exist!");
255 if (!mask->isValid()) {
258 if (particleV0->
getParticleSource() != Particle::EParticleSourceObject::c_Composite and
263 for (
auto* daughter : daughtersV0) {
264 if (daughter->getParticleSource() != Particle::EParticleSourceObject::c_Track) {
268 if (mask->hasV0(particleV0)) {
277 if (mask.getName() == name) {
285 TLorentzVector roe4Vector;
287 for (
const Particle* particle : myParticles) {
289 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_KLMCluster) {
292 roe4Vector += particle->get4Vector();
301 if (mask.getName() == name) {
310 std::vector<const Track*> result;
311 std::vector<const Particle*> allParticles =
getParticles(maskName);
312 for (
auto* particle : allParticles) {
313 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_Track) {
314 result.push_back(particle->getTrack());
321 std::vector<const ECLCluster*> result;
322 std::vector<const Particle*> allParticles =
getParticles(maskName);
323 for (
auto* particle : allParticles) {
325 auto* cluster = particle->getECLCluster();
327 result.push_back(cluster);
334 std::vector<const KLMCluster*> result;
335 std::vector<const Particle*> allParticles =
getParticles(maskName);
336 for (
auto* particle : allParticles) {
337 if (particle->getParticleSource() == Particle::EParticleSourceObject::c_KLMCluster) {
338 result.push_back(particle->getKLMCluster());
345 int nTracks =
getTracks(maskName).size();
351 return nROEECLClusters;
356 return nROEKLMClusters;
367 TLorentzVector roe4VectorTracks;
374 double fractions[n] = {0, 0, 1, 0, 0, 0};
378 for (
auto& roeTrack : roeTracks) {
380 const Track* track = roeTrack;
385 B2ERROR(
"Track with no PID information!");
395 int mcpdg = abs(mcp->getPDG());
406 double muIDBelle = 0.5;
407 if (pid->isAvailable(Const::KLM))
408 muIDBelle = exp(pid->getLogL(
Const::muon, Const::KLM));
413 if (eIDBelle > 0.9 and eIDBelle > muIDBelle)
415 else if (muIDBelle > 0.9 and eIDBelle < muIDBelle)
418 else if (atcPIDBelle_Kpi > 0.6)
430 if (fractions[0] == -1)
431 particlePDG = mcChoice;
433 else if (fractions[0] == -2)
434 particlePDG = pidChoice;
437 particlePDG = fracChoice;
440 const TrackFitResult* tfr = roeTrack->getTrackFitResultWithClosestMass(trackParticle);
443 double tempMass = trackParticle.getMass();
445 TLorentzVector temp4Vector;
446 temp4Vector.SetVect(tempMom);
447 temp4Vector.SetE(TMath::Sqrt(tempMom.Mag2() + tempMass * tempMass));
449 roe4VectorTracks += temp4Vector;
452 return roe4VectorTracks;
458 TLorentzVector roe4VectorECLClusters;
462 for (
auto& roeCluster : roeClusters) {
463 if (roeCluster->isNeutral())
468 return roe4VectorECLClusters;
473 for (
auto* listParticle : particlesToUpdate) {
474 if (roeParticle->
isCopyOf(listParticle,
true)) {
483 std::vector<std::string> maskNames;
486 maskNames.push_back(mask.getName());
500 B2INFO(
"This ROE is nested.");
509 std::string printout =
" -> ";
510 for (
const int index : indices) {
511 printout += std::to_string(index) +
", ";
519 std::set<int> source;
520 if (maskName ==
"") {
524 bool maskFound =
false;
526 if (mask.getName() == maskName) {
528 source = mask.getParticles();
533 B2FATAL(
"No " << maskName <<
" mask defined in current ROE!");
536 int particlePDG = (pdgCode == 0) ?
getPDGCode() : pdgCode;
537 auto isFlavored = (isSelfConjugated) ? Particle::EFlavorType::c_Unflavored : Particle::EFlavorType::c_Flavored;
538 return particles.appendNew(
get4Vector(maskName), particlePDG, isFlavored, std::vector(source.begin(),
539 source.end()), Particle::PropertyFlags::c_IsUnspecified);
546 double acc_sig = exp(pid->getLogL(
Const::kaon, set));
547 double acc_bkg = exp(pid->getLogL(
Const::pion, set));
549 if (acc_sig + acc_bkg > 0.0)
550 acc = acc_sig / (acc_sig + acc_bkg);
554 double tof_sig = exp(pid->getLogL(
Const::kaon, set));
555 double tof_bkg = exp(pid->getLogL(
Const::pion, set));
557 double tof_all = tof_sig + tof_bkg;
559 tof = tof_sig / tof_all;
560 if (tof < 0.001) tof = 0.001;
561 if (tof > 0.999) tof = 0.999;
566 double cdc_sig = exp(pid->getLogL(
Const::kaon, set));
567 double cdc_bkg = exp(pid->getLogL(
Const::pion, set));
569 double cdc_all = cdc_sig + cdc_bkg;
571 cdc = cdc_sig / cdc_all;
572 if (cdc < 0.001) cdc = 0.001;
573 if (cdc > 0.999) cdc = 0.999;
577 double pid_sig = acc * tof * cdc;
578 double pid_bkg = (1. - acc) * (1. - tof) * (1. - cdc);
580 return pid_sig / (pid_sig + pid_bkg);