Belle II Software development
TauDecayModeModule.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 <analysis/modules/TauDecayMode/TauDecayModeModule.h>
10#include <framework/logging/Logger.h>
11#include <framework/particledb/EvtGenDatabasePDG.h>
12#include <vector>
13#include <algorithm>
14#include <TParticlePDG.h>
15#include <map>
16#include <fstream>
17#include <Math/LorentzRotation.h>
18#include <Math/Boost.h>
19using namespace std;
20using namespace Belle2;
21
22std::map<std::string, int> make_map(const std::string& file, int chg)
23{
24 std::string fileName;
25 if (file == "") {
26 if (chg < 0) {
27 B2INFO("Missing input mapping file: use mp_file_minus=basf2.find_file('data/analysis/modules/TauDecayMode/map_tauminus.txt') TauDecayMode.param('file_minus', mp_file_minus) to classify with default mapping.");
28 } else {
29 B2INFO("Missing input mapping file: use mp_file_plus=basf2.find_file('data/analysis/modules/TauDecayMode/map_tauplus.txt') TauDecayMode.param('file_plus', mp_file_plus) to classify with default mapping.");
30 }
31 } else {
32 fileName = file;
33 }
34
35 ifstream f;
36 f.open(fileName);
37 std::string line;
38 std::string key;
39 int value;
40 std::map <std::string, int> map_tau;
41 while (f.good()) {
42 getline(f, line);
43 istringstream ss(line);
44 ss >> key >> value;
45 map_tau[key] = value;
46 }
47 f.close();
48 return map_tau;
49}
50
51//-----------------------------------------------------------------
52// Register the Module
53//-----------------------------------------------------------------
54REG_MODULE(TauDecayMode);
55//-----------------------------------------------------------------
56// Implementation
57//-----------------------------------------------------------------
58
60 m_mprong(0), m_pprong(0), m_megstar(0), m_pegstar(0),
64{
65 // Set module properties
66 setDescription("Module to identify generated tau pair decays, using MCParticle information.\n"
67 "By default, each tau decay is numbered as TauolaBelle2DecayMode [Ref: BELLE2-NOTE-PH-2020-055]");
68 //Parameter definition
69 addParam("printmode", m_printmode, "Printout more information from each event", std::string("default"));
70 addParam("file_minus", m_file_minus, "Path for an alternative mapping for tau- decays", std::string(""));
71 addParam("file_plus", m_file_plus, "Path for an alternative mapping for tau+ decays", std::string(""));
72}
73
74//
76{
77 mode_decay_minus = make_map(m_file_minus, -1);
78 mode_decay_plus = make_map(m_file_plus, 1);
79 m_tauDecay.registerInDataStore();
80 m_event_metadata.isRequired();
81}
82
83//
85{
86
88 if (tauPair) {
89
91
92 if (m_printmode == "missing") {
93 if (m_mmode == -1) {
94 B2INFO("TauDecayMode:: EventNumber = " << m_event_metadata->getEvent()
95 << " Decay: tau- -> " << m_tauminusdecaymode << " Mode = " << m_mmode);
96 }
97 if (m_pmode == -1) {
98 B2INFO("TauDecayMode:: EventNumber = " << m_event_metadata->getEvent()
99 << " Decay: tau+ -> " << m_tauplusdecaymode << " Mode = " << m_pmode);
100 }
101 }
102
103 if (m_printmode == "all") {
104 B2INFO("TauDecayMode:: EventNumber = " << m_event_metadata->getEvent()
105 << " Decay: tau- -> " << m_tauminusdecaymode << " Mode = " << m_mmode);
106 B2INFO("TauDecayMode:: EventNumber = " << m_event_metadata->getEvent()
107 << " Decay: tau+ -> " << m_tauplusdecaymode << " Mode = " << m_pmode);
108 }
109
110 //
111 if (m_mmode == -1) {
112 m_taum_no = m_taum_no - 1;
114 }
115 if (m_pmode == -1) {
116 m_taup_no = m_taup_no - 1;
118 }
119
122
125
126 } else {
127 m_pmode = -1;
128 m_mmode = -1;
129 m_pprong = -1;
130 m_mprong = -1;
131 m_megstar = -1;
132 m_pegstar = -1;
133 }
134
135 if (!m_tauDecay) m_tauDecay.create();
136 m_tauDecay->addTauMinusIdMode(m_mmode);
137 m_tauDecay->addTauPlusIdMode(m_pmode);
138
139 m_tauDecay->addTauPlusMcProng(m_pprong);
140 m_tauDecay->addTauMinusMcProng(m_mprong);
141
142 m_tauDecay->addTauMinusEgstar(m_megstar);
143 m_tauDecay->addTauPlusEgstar(m_pegstar);
144}
145
147{
148 // Clear local vectors
149 vec_nut.clear();
150 vec_anut.clear();
151 vec_numu.clear();
152 vec_anumu.clear();
153 vec_nue.clear();
154 vec_anue.clear();
155 vec_em.clear();
156 vec_ep.clear();
157 vec_mum.clear();
158 vec_mup.clear();
159 vec_pim.clear();
160 vec_pip.clear();
161 vec_km.clear();
162 vec_kp.clear();
163 vec_apro.clear();
164 vec_pro.clear();
165 vec_pi0.clear();
166 vec_k0s.clear();
167 vec_k0l.clear();
168 vec_eta.clear();
169 vec_omega.clear();
170 vec_etapr.clear();
171 vec_phi.clear();
172 vec_rhom.clear();
173 vec_rhop.clear();
174 vec_rho0.clear();
175 vec_kstarm.clear();
176 vec_kstarp.clear();
177 vec_kstar0.clear();
178 vec_kstar0_br.clear();
179 vec_a1m.clear();
180 vec_a1p.clear();
181 vec_a00_980.clear();
182 vec_a0m_980.clear();
183 vec_a0p_980.clear();
184 vec_a00_1450.clear();
185 vec_a0m_1450.clear();
186 vec_a0p_1450.clear();
187 vec_b1m.clear();
188 vec_b1p.clear();
189 vec_f1.clear();
190 vec_f0.clear();
191 vec_lambda.clear();
192 vec_lmb_br.clear();
193 vec_alpha.clear();
194 vec_gam.clear();
195 vec_radgam_taum.clear();
196 vec_radgam_taup.clear();
197 //
198 map<int, std::vector<int>> map_vec;
199 //
200 map_vec[11] = vec_em;
201 map_vec[-11] = vec_ep;
202 map_vec[12] = vec_nue;
203 map_vec[-12] = vec_anue;
204 map_vec[13] = vec_mum;
205 map_vec[-13] = vec_mup;
206 map_vec[14] = vec_numu;
207 map_vec[-14] = vec_anumu;
208 map_vec[16] = vec_nut;
209 map_vec[-16] = vec_anut;
210 map_vec[-211] = vec_pim;
211 map_vec[211] = vec_pip;
212 map_vec[-321] = vec_km;
213 map_vec[321] = vec_kp;
214 map_vec[-2212] = vec_apro;
215 map_vec[2212] = vec_pro;
216 map_vec[111] = vec_pi0;
217 map_vec[310] = vec_k0s;
218 map_vec[130] = vec_k0l;
219 map_vec[221] = vec_eta;
220 map_vec[223] = vec_omega;
221 map_vec[331] = vec_etapr;
222 map_vec[333] = vec_phi;
223 map_vec[-213] = vec_rhom;
224 map_vec[213] = vec_rhop;
225 map_vec[113] = vec_rho0;
226 map_vec[-323] = vec_kstarm;
227 map_vec[323] = vec_kstarp;
228 map_vec[313] = vec_kstar0;
229 map_vec[-313] = vec_kstar0_br;
230 map_vec[-20213] = vec_a1m;
231 map_vec[20213] = vec_a1p;
232 map_vec[-9000211] = vec_a0m_980;
233 map_vec[9000211] = vec_a0p_980;
234 map_vec[9000111] = vec_a00_980;
235 map_vec[-10211] = vec_a0m_1450;
236 map_vec[10211] = vec_a0p_1450;
237 map_vec[10111] = vec_a00_1450;
238 map_vec[-10213] = vec_b1m;
239 map_vec[10213] = vec_b1p;
240 map_vec[20223] = vec_f1;
241 map_vec[9010221] = vec_f0;
242 map_vec[3122] = vec_lambda;
243 map_vec[-3122] = vec_lmb_br;
244 map_vec[94144] = vec_alpha;
245 map_vec[22] = vec_gam;
246
247 bool TauMinusIsLeptonicDecay = false;
248 bool TauPlusIsLeptonicDecay = false;
249
250 bool elecFirst = true;
251 bool muonFirst = true;
252
253 bool isPiPizGamTauMinusFirst = true;
254 bool isPiPizGamTauPlusFirst = true;
255
256 bool isLFVTauMinus2BodyDecayFirst = true;
257 bool isLFVTauPlus2BodyDecayFirst = true;
258
259 bool isChargedRhoFromTauMinusFirst = true;
260 bool isChargedRhoFromTauPlusFirst = true;
261
262 bool isChargedA1FromTauMinusFirst = true;
263 bool isChargedA1FromTauPlusFirst = true;
264
265 bool isEtaPPGFromTauMinusFirst = true;
266 bool isEtaPPGFromTauPlusFirst = true;
267
268 bool isOmegaPizGamFromTauMinusFirst = true;
269 bool isOmegaPizGamFromTauPlusFirst = true;
270
271 bool isEtaPizPizPizFromTauMinusFirst = true;
272 bool isEtaPizPizPizFromTauPlusFirst = true;
273
276
277 bool isOmegaPimPipFromTauMinusFirst = true;
278 bool isOmegaPimPipFromTauPlusFirst = true;
279
282
283 // Loop of MCParticles
284 for (int i = 0; i < MCParticles.getEntries(); i++) {
285
286 MCParticle& p = *MCParticles[i];
287
288 int pdgid = p.getPDG();
289
290 if (pdgid == -Const::Klong.getPDGCode()) pdgid = -pdgid; // Strange feature in TauolaBelle2
291
292 // Special treatment for photons
293 bool accept_photon = false;
294
295 if (!p.hasStatus(MCParticle::c_PrimaryParticle))
296 continue; // only consider particles coming from generator, e.g. discard particles added by Geant4
297 if (p.isInitial()) continue; // pick e- e+, but not from the incoming beams
298
299 // Check IsLeptonicDecay first
300 if (abs(pdgid) == 12 || abs(pdgid) == 14) { // (anti-)nu_e or (anti-)nu_mu
301 int leptonicdecay = 0;
302 const MCParticle* mother = p.getMother(); // tau- or tau+
303 const vector<MCParticle*> daughters = mother->getDaughters();
304 for (MCParticle* d : daughters) {
305 int dauid = abs(d->getPDG());
306 if (dauid == 11 || dauid == 12 || dauid == 13 || dauid == 14 || dauid == 16) {
307 leptonicdecay++;
308 }
309 }
310 if (leptonicdecay == 3) {
311 if (mother->getPDG() == 15) TauMinusIsLeptonicDecay = true;
312 if (mother->getPDG() == -15) TauPlusIsLeptonicDecay = true;
313 }
314 }
315
316 else if (pdgid == Const::electron.getPDGCode() && elecFirst) {
317 elecFirst = false;
318 const MCParticle* mother = p.getMother();
319 if (not mother) continue; // In some low multiplicity generators there may be no mother
320 const vector<MCParticle*> daughters = mother->getDaughters();
321 int nElMinus = 0;
322 int nElPlus = 0;
323 stringstream elec_ss;
324 for (MCParticle* d : daughters) {
325 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
326 elec_ss << d->getPDG() << " ";
327 if (d->getPDG() == Const::electron.getPDGCode()) nElMinus++;
328 if (d->getPDG() == -Const::electron.getPDGCode()) nElPlus++;
329 }
330 if (nElMinus == 1 && nElPlus == 1) { // use this information in getRecursiveMotherCharge
331 B2DEBUG(19, "Mother of elec pair is = " << mother->getPDG() << " which has daughters : " << elec_ss.str());
332 }
333 }
334
335 else if (pdgid == Const::muon.getPDGCode() && muonFirst) {
336 muonFirst = false;
337 const MCParticle* mother = p.getMother();
338 if (not mother) continue; // In some low multiplicity generators there may be no mother
339 const vector<MCParticle*> daughters = mother->getDaughters();
340 int nMuMinus = 0;
341 int nMuPlus = 0;
342 stringstream muon_ss;
343 for (MCParticle* d : daughters) {
344 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
345 muon_ss << d->getPDG() << " ";
346 if (d->getPDG() == Const::muon.getPDGCode()) nMuMinus++;
347 if (d->getPDG() == -Const::muon.getPDGCode()) nMuPlus++;
348 }
349 if (nMuMinus == 1 && nMuPlus == 1) { // use this information in getRecursiveMotherCharge
350 B2DEBUG(19, "Mother of muon pair is = " << mother->getPDG() << " which has daughters : " << muon_ss.str());
351 }
352 }
353
354 // Photons from PHOTOS are radiative photons
355 else if (pdgid == Const::photon.getPDGCode() && p.hasStatus(MCParticle::c_IsPHOTOSPhoton)) {
357 if (chg < 0) vec_radgam_taum.push_back(i);
358 if (chg > 0) vec_radgam_taup.push_back(i);
359 }
360
361 // Photons from PHOTOS do not define tau decay mode
362 else if (pdgid == Const::photon.getPDGCode() && !p.hasStatus(MCParticle::c_IsPHOTOSPhoton)) {
363 const MCParticle* mother = p.getMother();
364 if (not mother) continue; // In some low multiplicity generators there may be no mother
365 int mothid = abs(mother->getPDG());
366
367 // check if the gamma comes from final state charged particles {e, mu, pi, K, p, b_1}
368 bool isRadiationfromFinalStateChargedParticle = false;
369 if (mothid == 11 ||
370 mothid == 13 ||
371 mothid == 211 ||
372 mothid == 321 ||
373 mothid == 2212 ||
374 mothid == 10213) {
375 isRadiationfromFinalStateChargedParticle = true;
376 }
377
378 // check if the gamma from ChargedRho
379 bool isRadiationFromChargedRho = false;
380 if (mothid == 213) {
381 int chg = getRecursiveMotherCharge(mother);
382 if (chg < 0 && isChargedRhoFromTauMinusFirst) {
383 isChargedRhoFromTauMinusFirst = false;
384 isRadiationFromChargedRho = true;
385 }
386 if (chg > 0 && isChargedRhoFromTauPlusFirst) {
387 isChargedRhoFromTauPlusFirst = false;
388 isRadiationFromChargedRho = true;
389 }
390 }
391
392 // check if the gamma from ChargedA1
393 bool isRadiationFromChargedA1 = false;
394 if (mothid == 20213) {
395 int chg = getRecursiveMotherCharge(mother);
396 if (chg < 0 && isChargedA1FromTauMinusFirst) {
397 isChargedA1FromTauMinusFirst = false;
398 isRadiationFromChargedA1 = true;
399 }
400 if (chg > 0 && isChargedA1FromTauPlusFirst) {
401 isChargedA1FromTauPlusFirst = false;
402 isRadiationFromChargedA1 = true;
403 }
404 }
405
406 // check if the gamma comes from intermediate W+- boson
407 bool isRadiationfromIntermediateWBoson = false;
408 if (mothid == 24) isRadiationfromIntermediateWBoson = true;
409
410 // check if it is a tau- -> pi- omega nu, omega -> pi0 gamma decay
411 // Note: TauolaBelle2 generator treats this a coherent production
412 // e.g. includes omega in the form factor but not as an explicit final state particle
413 bool isPiPizGam = false;
414 if (isRadiationfromIntermediateWBoson) {
415 const vector<MCParticle*> daughters = mother->getDaughters();
416 int nPiSisters = 0;
417 int nPizSisters = 0;
418 int nTotSisters = 0;
419 int nOtherSisters = 0;
420 for (MCParticle* d : daughters) {
421 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
422 nTotSisters++;
423 if (abs(d->getPDG()) == Const::pion.getPDGCode()) {
424 nPiSisters++;
425 } else if (abs(d->getPDG()) == Const::pi0.getPDGCode()) {
426 nPizSisters++;
427 } else if (abs(d->getPDG()) != Const::photon.getPDGCode()) {
428 nOtherSisters++;
429 }
430 }
431 // allow final state radiation from the tau lepton, but ignore information on number of photons for decay mode classifiction
432 if (nTotSisters >= 3 && nPiSisters == 1 && nPizSisters == 1 && nOtherSisters == 0) {
433 int chg = getRecursiveMotherCharge(mother);
434 if (chg < 0 && isPiPizGamTauMinusFirst) {
435 isPiPizGamTauMinusFirst = false;
436 isPiPizGam = true;
437 }
438 if (chg > 0 && isPiPizGamTauPlusFirst) {
439 isPiPizGamTauPlusFirst = false;
440 isPiPizGam = true;
441 }
442 }
443 }
444
445 // check if the gamma comes from tau
446 bool isRadiationfromTau = false;
447 if (mothid == 15) isRadiationfromTau = true;
448
449 // Photons from LeptonicTauDecays are radiative photons not added by PHOTOS
450 if (isRadiationfromTau) {
452 if (chg < 0 && TauMinusIsLeptonicDecay) vec_radgam_taum.push_back(i);
453 if (chg > 0 && TauPlusIsLeptonicDecay) vec_radgam_taup.push_back(i);
454 }
455
456 // check if the gamma comes from 2 body LFV decays of tau, e.g. tau- -> e-/mu- gamma with arbitrary number of extra photon radiations from PHOTOS/FSR
457 bool isLFVTau2BodyDecay = false;
458 if (isRadiationfromTau) {
459 bool hasNeutrinoAsSister = false;
460 int numChargedSister = 0;
461 int numNeutralNonNeutrinoNonPhotonSister = 0;
462 const vector<MCParticle*> daughters = mother->getDaughters();
463 for (MCParticle* d : daughters) {
464 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
465 hasNeutrinoAsSister = find(Neutrinos.begin(), Neutrinos.end(), abs(d->getPDG())) != Neutrinos.end();
466 if (hasNeutrinoAsSister) break;
467 }
468 if (!hasNeutrinoAsSister) {
469 for (MCParticle* d : daughters) {
470 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
471 bool isChargedFinalState = find(finalStatePDGs.begin(), finalStatePDGs.end(), abs(d->getPDG())) != finalStatePDGs.end();
472 if (isChargedFinalState) {
473 numChargedSister++;
474 } else if (d->getPDG() != Const::photon.getPDGCode()) {
475 numNeutralNonNeutrinoNonPhotonSister++;
476 }
477 }
478 if (numChargedSister == 1 && numNeutralNonNeutrinoNonPhotonSister == 0) {
479 if (mother->getPDG() == 15 && isLFVTauMinus2BodyDecayFirst) {
480 isLFVTauMinus2BodyDecayFirst = false;
481 isLFVTau2BodyDecay = true;
482 }
483 if (mother->getPDG() == -15 && isLFVTauPlus2BodyDecayFirst) {
484 isLFVTauPlus2BodyDecayFirst = false;
485 isLFVTau2BodyDecay = true;
486 }
487 }
488 }
489 B2DEBUG(19, "hasNeutrinoAsSister = " << hasNeutrinoAsSister
490 << " numChargedSister = " << numChargedSister
491 << " numNeutralNonNeutrinoNonPhotonSister = " << numNeutralNonNeutrinoNonPhotonSister
492 << " isLFVTau2BodyDecay = " << isLFVTau2BodyDecay);
493 }
494
495 bool isPi0GG = false;
496 bool isEtaGG = false;
497 bool isEtpGG = false;
498 bool isPi0GEE = false;
499 bool isEtaPPG = false;
500 bool isOmPizG = false;
501 if (mothid == 111 || mothid == 221 || mothid == 331 || mothid == 223) {
502 const vector<MCParticle*> daughters = mother->getDaughters();
503 int numberofTotalDaughters = 0;
504 int numberOfPhotonDaughters = 0;
505 int numberOfElectronDaughters = 0;
506 int numberOfPionDaughters = 0;
507 int numberOfPizDaughters = 0;
508 for (MCParticle* d : daughters) {
509 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
510 numberofTotalDaughters ++;
511 if (abs(d->getPDG()) == Const::photon.getPDGCode()) numberOfPhotonDaughters++;
512 if (abs(d->getPDG()) == Const::electron.getPDGCode()) numberOfElectronDaughters++;
513 if (abs(d->getPDG()) == Const::pion.getPDGCode()) numberOfPionDaughters++;
514 if (abs(d->getPDG()) == Const::pi0.getPDGCode()) numberOfPizDaughters++;
515 }
516 if (numberofTotalDaughters == 2 && numberOfPhotonDaughters == 2) {
517 if (mothid == Const::pi0.getPDGCode()) isPi0GG = true;
518 if (mothid == 221) isEtaGG = true;
519 if (mothid == 331) isEtpGG = true;
520 }
521 if (numberofTotalDaughters >= 3 && numberOfPhotonDaughters >= 1 && numberOfElectronDaughters == 2 && numberOfPionDaughters == 0
522 && numberOfPizDaughters == 0) {
523 if (mothid == Const::pi0.getPDGCode()) isPi0GEE = true;
524 }
525 if (numberofTotalDaughters >= 3 && numberOfPhotonDaughters >= 1 && numberOfPizDaughters == 0 && numberOfPionDaughters == 2) {
526 if (mothid == 221) {
527 int chg = getRecursiveMotherCharge(mother);
528 if (chg < 0 && isEtaPPGFromTauMinusFirst) {
529 isEtaPPGFromTauMinusFirst = false;
530 isEtaPPG = true;
531 }
532 if (chg > 0 && isEtaPPGFromTauPlusFirst) {
533 isEtaPPGFromTauPlusFirst = false;
534 isEtaPPG = true;
535 }
536 }
537 }
538 if (numberofTotalDaughters >= 2 && numberOfPhotonDaughters >= 1 && numberOfPizDaughters == 1 && numberOfPionDaughters == 0) {
539 if (mothid == 223) {
540 int chg = getRecursiveMotherCharge(mother);
541 if (chg < 0 && isOmegaPizGamFromTauMinusFirst) {
542 isOmegaPizGamFromTauMinusFirst = false;
543 isOmPizG = true;
544 }
545 if (chg > 0 && isOmegaPizGamFromTauPlusFirst) {
546 isOmegaPizGamFromTauPlusFirst = false;
547 isOmPizG = true;
548 }
549 }
550 }
551 }
552
553 B2DEBUG(19, "isRadiationfromFinalStateChargedParticle = " << isRadiationfromFinalStateChargedParticle);
554 B2DEBUG(19, "isRadiationFromChargedRho = " << isRadiationFromChargedRho);
555 B2DEBUG(19, "isRadiationFromChargedA1 = " << isRadiationFromChargedA1);
556 B2DEBUG(19, "isRadiationfromIntermediateWBoson = " << isRadiationfromIntermediateWBoson << " isPiPizGam = " << isPiPizGam);
557 B2DEBUG(19, "isRadiationfromTau = " << isRadiationfromTau << " isLFVTau2BodyDecay = " << isLFVTau2BodyDecay);
558 B2DEBUG(19, "isPi0GG = " << isPi0GG << " isEtaGG = " << isEtaGG << " isEtpGG = " << isEtpGG << " isPi0GEE = " << isPi0GEE);
559 B2DEBUG(19, "isEtaPPG = " << isEtaPPG << " isOmPizG = " << isOmPizG);
560
561 // accept photons if (isRadiationfromFinalStateChargedParticle is false) or
562 // if (isRadiationfromIntermediateWBoson is false) or
563 // if (isRadiationfromTau is true) {gamma is from 2 body LFV decays of tau} or
564 // if the radiation is coming from any other decay [except pi0 -> gamma gamma, eta-> gamma gamma, eta' -> gamma gamma]
565 if (isRadiationfromFinalStateChargedParticle) {
566 } else if (isRadiationFromChargedRho) {
567 accept_photon = true;
568 } else if (isRadiationFromChargedA1) {
569 accept_photon = true;
570 } else if (isRadiationfromIntermediateWBoson) {
571 if (isPiPizGam) {
572 accept_photon = true;
573 }
574 } else if (isRadiationfromTau) { // accept one photon from tau -> (charged particle) + gamma [no neutrinos]
575 if (isLFVTau2BodyDecay) {
576 accept_photon = true;
577 }
578 } else if (isPi0GG) {
579 } else if (isEtaGG) {
580 } else if (isEtpGG) {
581 } else if (isPi0GEE) {
582 } else if (isEtaPPG) {
583 accept_photon = true;
584 } else if (isOmPizG) {
585 accept_photon = true;
586 }
587 }
588
589 // Without further analysis, it is NOT possible to separate
590 // tau- -> pi- 2pi0 eta (-> pi- pi+ pi0) nu decays [Mode 39] from
591 // tau- -> 2pi- pi+ eta (-> pi0 pi0 pi0) nu decays [Mode 42]
592
593 else if (pdgid == Const::pi0.getPDGCode() && (isEtaPizPizPizFromTauMinusFirst || isEtaPizPizPizFromTauPlusFirst)) {
594 const MCParticle* mother = p.getMother();
595 // In some low multiplicity generators there may be no mother
596 if (mother and (mother->getPDG() == 221)) { // eta -> pi0 pi0 pi0
597 const vector<MCParticle*> daughters = mother->getDaughters();
598 int nPizSisters = 0;
599 for (MCParticle* d : daughters) {
600 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
601 if (d->getPDG() == Const::pi0.getPDGCode()) nPizSisters++;
602 }
603 B2DEBUG(19, "nPizSisters = " << nPizSisters);
604 if (nPizSisters == 3) {
605 int chg = getRecursiveMotherCharge(mother);
606 if (chg < 0 && isEtaPizPizPizFromTauMinusFirst) {
607 isEtaPizPizPizFromTauMinusFirst = false;
609 }
610 if (chg > 0 && isEtaPizPizPizFromTauPlusFirst) {
611 isEtaPizPizPizFromTauPlusFirst = false;
613 }
614 }
615 }
616 }
617
618 // Without further analysis, it is NOT possible to separate
619 // tau- -> pi- 2pi0 omega (-> pi- pi+) nu decays [Mode 49] from
620 // tau- -> 2pi- pi+ 2pi0 nu decays [Mode 66]
621 // Note: TauolaBelle2 treats Mode 66 is coherent production,
622 // e.g. includes omega in the form factor but not as an explicit final state particle.
623 // Note2: omega is explicitly included in following mode:
624 // tau- -> pi- pi0 omega (-> pi- pi+ pi0) nu decays [Mode 130]
625 // which is where Mode 49 is mapped to if there is no omega->pi-pi+ decay
626
627 // Same confusion can also happen for
628 // tau- -> pi- pi0 omega (-> pi- pi+) nu decays [Mode 131]
629 // tau- -> pi- omega (-> pi- pi+ pi0) nu decays [Mode 236]
630
631 else if (pdgid == -Const::pion.getPDGCode() && (isOmegaPimPipFromTauMinusFirst || isOmegaPimPipFromTauPlusFirst)) {
632 const MCParticle* mother = p.getMother();
633 // In some low multiplicity generators there may be no mother
634 if (mother and (mother->getPDG() == 223)) { // omega -> pi- pi+
635 const vector<MCParticle*> daughters = mother->getDaughters();
636 int nOmegaDaughters = 0;
637 int nPimSisters = 0;
638 int nPipSisters = 0;
639 int nPizSisters = 0;
640 for (MCParticle* d : daughters) {
641 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
642 nOmegaDaughters++;
643 if (d->getPDG() == -Const::pion.getPDGCode()) nPimSisters++;
644 if (d->getPDG() == Const::pion.getPDGCode()) nPipSisters++;
645 if (d->getPDG() == Const::pi0.getPDGCode()) nPizSisters++;
646 }
647 B2DEBUG(19,
648 "nOmegaDaughters = " << nOmegaDaughters << " "
649 "nPimSisters = " << nPimSisters << " "
650 "nPipSisters = " << nPipSisters << " "
651 "nPizSisters = " << nPizSisters);
652 if (nOmegaDaughters >= 2 && nPimSisters == 1 && nPipSisters == 1 && nPizSisters == 0) { // allow omega->pi-pi+gama
653 int chg = getRecursiveMotherCharge(mother);
654 if (chg < 0 && isOmegaPimPipFromTauMinusFirst) {
655 isOmegaPimPipFromTauMinusFirst = false;
657 }
658 if (chg > 0 && isOmegaPimPipFromTauPlusFirst) {
659 isOmegaPimPipFromTauPlusFirst = false;
661 }
662 }
663 }
664 }
665 B2DEBUG(19,
666 "isEtaPizPizPizFromTauMinusFirst = " << isEtaPizPizPizFromTauMinusFirst << " "
667 "m_isEtaPizPizPizFromTauMinus = " << m_isEtaPizPizPizFromTauMinus << " "
668 "isEtaPizPizPizFromTauPlusFirst = " << isEtaPizPizPizFromTauPlusFirst << " "
669 "m_isEtaPizPizPizFromTauPlus = " << m_isEtaPizPizPizFromTauPlus
670 );
671 B2DEBUG(19,
672 "isOmegaPimPipFromTauMinusFirst = " << isOmegaPimPipFromTauMinusFirst << " "
673 "m_isOmegaPimPipFromTauMinus = " << m_isOmegaPimPipFromTauMinus << " "
674 "isOmegaPimPipFromTauPlusFirst = " << isOmegaPimPipFromTauPlusFirst << " "
675 "m_isOmegaPimPipFromTauPlus = " << m_isEtaPizPizPizFromTauPlus
676 );
677
678 // Fill up all the vector of particles
679 map<int, std::vector<int>>::iterator ite ;
680 for (ite = map_vec.begin(); ite != map_vec.end(); ++ite) {
681 if (pdgid == ite->first) {
682 if (pdgid == Const::photon.getPDGCode()) {
683 if (accept_photon) {
684 B2DEBUG(19, "Photon accepted");
685 ite-> second.push_back(i);
686 }
687 } else {
688 ite-> second.push_back(i);
689 }
690 break;
691 }
692 }
693
694 } // End of loop over MCParticles
695
696 vec_dau_tauminus.clear();
697 vec_dau_tauplus.clear();
698
699 map<int, std::vector<int>>::iterator itr ;
700 for (itr = map_vec.begin(); itr != map_vec.end(); ++itr) {
701 for (unsigned int i = 0; i < itr-> second.size(); i++) {
702 int ii = itr-> second[i];
704 if (chg < 0) vec_dau_tauminus.push_back(ii);
705 if (chg > 0) vec_dau_tauplus.push_back(ii);
706 }
707 }
708
710 std::string pdgname;
711
712 //make decay string for tau-
714 for (unsigned iorder = 0; iorder < OrderedList.size(); ++iorder) {
715 int ii = OrderedList[iorder];
716 //
717 for (unsigned int i = 0; i < vec_dau_tauminus.size(); i++) {
719 int pdg = p->getPDG();
720 if (pdg == -Const::Klong.getPDGCode()) pdg = -pdg; // Strange Feature in TauolaBelle2
721 if (pdg != ii) continue;
722 m_tauminusdecaymode.append(".");
723 if (pdg == 94144) {
724 pdgname = "alpha";
725 } else {
726 pdgname = databasePDG->GetParticle(pdg)->GetName();
727 }
728 m_tauminusdecaymode.append(pdgname);
729 }
730 }
731 //
733
734 //make decay string for tau+
736 for (unsigned iorder = 0; iorder < OrderedList.size(); ++iorder) {
737 int ii = OrderedList[iorder];
738 //
739 for (unsigned int i = 0; i < vec_dau_tauplus.size(); i++) {
741 int pdg = p->getPDG();
742 if (pdg == -Const::Klong.getPDGCode()) pdg = -pdg; // Strange Feature in TauolaBelle2
743 if (pdg != ii) continue;
744 m_tauplusdecaymode.append(".");
745 if (pdg == 94144) {
746 pdgname = "alpha";
747 } else {
748 pdgname = databasePDG->GetParticle(pdg)->GetName();
749 }
750 m_tauplusdecaymode.append(pdgname);
751 }
752 }
753 //
755
756}
757
758//
759int TauDecayModeModule::TauolaBelle2DecayMode(const std::string& state, int chg)
760{
761 std::map<std::string, int> mode_decay = (chg < 0) ? mode_decay_minus : mode_decay_plus;
763 for (itr = mode_decay.begin(); itr != mode_decay.end(); ++itr) {
764 if (state == itr-> first) {
765 int mode = itr-> second;
766 if (mode == 39) {
767 if (chg < 0 && m_isEtaPizPizPizFromTauMinus) mode = 42;
768 if (chg > 0 && m_isEtaPizPizPizFromTauPlus) mode = 42;
769 }
770 if (mode == 49) {
771 if (chg < 0 && !m_isOmegaPimPipFromTauMinus) mode = 130;
772 if (chg > 0 && !m_isOmegaPimPipFromTauPlus) mode = 130;
773 }
774 if (mode == 131) {
775 if (chg < 0 && !m_isOmegaPimPipFromTauMinus) mode = 236;
776 if (chg > 0 && !m_isOmegaPimPipFromTauPlus) mode = 236;
777 }
778 if (mode == 247) {
779 if (chg < 0 && !m_isEtaPizPizPizFromTauMinus) mode = 15;
780 if (chg > 0 && !m_isEtaPizPizPizFromTauPlus) mode = 15;
781 }
782 return mode;
783 }
784 }
785 return -1;
786}
787
788//
790{
791 const MCParticle* mother = p->getMother();
792 if (mother == nullptr) {
793 return 0;
794 } else if (abs(p->getPDG()) == Const::electron.getPDGCode() && mother->getPDG() == Const::pi0.getPDGCode()) { // pi0 -> e-e+ gamma
795 return 0;
796 } else if (abs(p->getPDG()) == Const::electron.getPDGCode() && mother->getPDG() == 221) { // eta -> e-e+ gamma
797 return 0;
798 } else if (abs(p->getPDG()) == Const::electron.getPDGCode() && mother->getPDG() == 113) { // rho0 -> e-e+
799 return 0;
800 } else if (abs(p->getPDG()) == Const::electron.getPDGCode() && mother->getPDG() == 223) { // omega -> e- e+
801 return 0;
802 } else if (abs(p->getPDG()) == Const::electron.getPDGCode() && mother->getPDG() == 333) { // phi -> e- e+
803 return 0;
804 } else if (abs(p->getPDG()) == Const::muon.getPDGCode() && mother->getPDG() == 221) { // eta -> mu-mu+ gamma
805 return 0;
806 } else if (abs(p->getPDG()) == Const::muon.getPDGCode() && mother->getPDG() == 113) { // rho0 -> mu-mu+
807 return 0;
808 } else if (abs(p->getPDG()) == Const::muon.getPDGCode() && mother->getPDG() == 333) { // phi -> mu-mu+
809 return 0;
810 } else if (mother->getPDG() == 15) {
811 return -1;
812 } else if (mother->getPDG() == -15) {
813 return 1;
814 } else {
815 return getRecursiveMotherCharge(mother);
816 }
817}
818
819//
821{
822 numOfTauPlus = 0;
823 numOfTauMinus = 0;
824 idOfTauPlus = 0;
825 idOfTauMinus = 0;
826 for (int i = 0; i < MCParticles.getEntries(); i++) {
827 MCParticle& p = *MCParticles[i];
828
829 if (p.getStatus() == 1 && p.getPDG() == 15) {
831 idOfTauMinus = p.getIndex();
832 }
833 if (p.getStatus() == 1 && p.getPDG() == -15) {
834 numOfTauPlus++;
835 idOfTauPlus = p.getIndex();
836 }
837 }
838 if (numOfTauPlus == 1 && numOfTauMinus == 1) {
839 tauPair = true;
840 } else tauPair = false;
841}
842
843//
845{
846 int ret = 0;
847 const vector<MCParticle*> daughters = p.getDaughters();
848 if (daughters.empty()) return ret;
849 for (MCParticle* d : daughters) {
850 if (!d->hasStatus(MCParticle::c_PrimaryParticle)) continue;
851 // TODO: Improve how to identify a final state particle.
852 bool isChargedFinalState = find(finalStatePDGs.begin(),
853 finalStatePDGs.end(),
854 abs(d->getPDG())) != finalStatePDGs.end();
855 if (isChargedFinalState) ret++;
856 else ret += getProngOfDecay(*d);
857 }
858 return ret;
859}
860
861// Energy of the radiative photon in tau rest frame
862double TauDecayModeModule::getEgstar(const std::vector<int>& vec_radgam, const MCParticle& p)
863{
864 double egstar = -1.;
865 ROOT::Math::PxPyPzEVector p4_tau = p.get4Vector();
866 B2DEBUG(19, "p4_tau: " << p4_tau << " " << p4_tau.P());
867 ROOT::Math::Boost boost_to_mother_rest_frame(p4_tau.BoostToCM());
868 for (auto i : vec_radgam) {
869 ROOT::Math::PxPyPzEVector p4_gamma = MCParticles[i]->get4Vector();
870 B2DEBUG(19, "p4_gamma: " << p4_gamma << " " << p4_gamma.P());
871 ROOT::Math::PxPyPzEVector p4_gamma_rest = boost_to_mother_rest_frame * p4_gamma;
872 B2DEBUG(19, "p4_gamma_rest: " << p4_gamma_rest << " " << p4_gamma_rest.P());
873 double energy_rest = p4_gamma_rest.E();
874 if (energy_rest > egstar) egstar = energy_rest;
875 }
876 B2DEBUG(19, "egstar: " << egstar);
877 return egstar;
878}
static const ParticleType pi0
neutral pion particle
Definition Const.h:674
static const ChargedStable muon
muon particle
Definition Const.h:660
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 ParticleType photon
photon particle
Definition Const.h:673
static const ChargedStable electron
electron particle
Definition Const.h:659
Replacement for TDatabasePDG that is filled from EvtGen's evt.pdl.
static EvtGenDatabasePDG * Instance()
Instance method that loads the EvtGen table.
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
@ c_IsPHOTOSPhoton
bit 8: Particle is an radiative photon from PHOTOS
Definition MCParticle.h:63
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
std::vector< Belle2::MCParticle * > getDaughters() const
Get vector of all daughter particles, empty vector if none.
Definition MCParticle.cc:52
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
std::vector< int > vec_numu
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_em
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_omega
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_k0s
Variable name of the vector where particles identified in the event are stored.
Int_t m_mprong
Prong of the decay channel of tau-.
std::vector< int > vec_b1p
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_f1
Variable name of the vector where particles identified in the event are stored.
int numOfTauPlus
Number of tau+ in the event.
Int_t m_pprong
Prong of the decay channel of tau+.
StoreObjPtr< EventMetaData > m_event_metadata
event number
void IdentifyTauPair()
Identifies if the event is a generated tau pair.
std::vector< int > vec_apro
Variable name of the vector where particles identified in the event are stored.
double getEgstar(const std::vector< int > &vec_radgam, const MCParticle &mc)
Energy of the radiative photon in tau rest frame.
std::vector< int > vec_pi0
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_alpha
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_radgam_taum
Variable name of the vector where particles identified in the event are stored.
StoreArray< MCParticle > MCParticles
StoreArray of MCParticles.
std::vector< int > vec_anue
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_nue
Variable name of the vector where particles identified in the event are stored.
virtual void initialize() override
Initializes the module.
std::vector< int > vec_a0p_980
Variable name of the vector where particles identified in the event are stored.
std::string m_file_plus
Alternative mapping for tau+.
static constexpr std::array< int, 46 > OrderedList
PDG codes of ORDERED particles.
std::vector< int > vec_a0m_1450
Variable name of the vector where particles identified in the event are stored.
virtual void event() override
Method is called for each event.
bool m_isEtaPizPizPizFromTauMinus
Flag for eta->pi0pi0pi0 decays from tau-.
static constexpr std::array< int, 3 > Neutrinos
PDG codes of neutrinos in final state particles in generation: {nu_e, nu_mu, mu_tau}...
double m_pegstar
Energy of radiative photon from tau+.
std::vector< int > vec_lmb_br
Variable name of the vector where particles identified in the event are stored.
Int_t m_mmode
ID of the decay channel of tau-.
void AnalyzeTauPairEvent()
Analyze a generated tau pair event.
bool m_isOmegaPimPipFromTauMinus
Flag for omega->pi-pi+ decays from tau-.
int getRecursiveMotherCharge(const MCParticle *mc)
Identifies particles coming from tau decays.
std::vector< int > vec_phi
Variable name of the vector where particles identified in the event are stored.
int m_taum_no
number of tau- unclassified events
std::string m_tauminusdecaymode
Variable name for the decay mode of the tau-.
std::vector< int > vec_a1p
Variable name of the vector where particles identified in the event are stored.
bool m_isOmegaPimPipFromTauPlus
Flag for omega->pi-pi+ decays from tau+.
bool tauPair
Boolean variable used to identify tau event.
int TauolaBelle2DecayMode(const std::string &s, int chg)
Classifies the decays of the event and assigns a decay mode.
int idOfTauMinus
Index of the generated tau-.
std::vector< int > vec_lambda
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_anut
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_f0
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_mum
Variable name of the vector where particles identified in the event are stored.
Int_t m_pmode
ID of the decay channel of tau+.
int numOfTauMinus
Number of tau- in the event.
std::vector< int > vec_nut
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_kstar0_br
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_radgam_taup
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_a00_1450
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_kstar0
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_a0m_980
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_a00_980
Variable name of the vector where particles identified in the event are stored.
int getProngOfDecay(const MCParticle &mc)
Identifies the number of charged final state particles in the decay.
std::vector< int > vec_a1m
Variable name of the vector where particles identified in the event are stored.
int idOfTauPlus
Index of the generated tau+.
StoreObjPtr< TauPairDecay > m_tauDecay
pointer to tau pair decay objects
std::vector< int > vec_rhop
Variable name of the vector where particles identified in the event are stored.
std::string m_printmode
Parameter passed by the user to indicate the information to be printed.
std::vector< int > vec_pro
Variable name of the vector where particles identified in the event are stored.
bool m_isEtaPizPizPizFromTauPlus
Flag for eta->pi0pi0pi0 decays from tau+.
std::string m_tauplusdecaymode
Variable name for the decay mode of the tau+.
std::vector< int > vec_pim
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_mup
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_pip
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_kstarp
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_kstarm
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_etapr
Variable name of the vector where particles identified in the event are stored.
double m_megstar
Energy of radiative photon from tau-.
std::vector< int > vec_ep
Variable name of the vector where particles identified in the event are stored.
TauDecayModeModule()
Constructor: Sets the description, the properties and the parameters of the module.
std::vector< int > vec_k0l
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_eta
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_km
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_anumu
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_dau_tauplus
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_rhom
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_gam
Variable name of the vector where particles identified in the event are stored.
std::string m_file_minus
Alternative mapping for tau-.
int m_taup_no
number of tau+ unclassified events
std::vector< int > vec_dau_tauminus
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_kp
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_b1m
Variable name of the vector where particles identified in the event are stored.
std::vector< int > vec_rho0
Variable name of the vector where particles identified in the event are stored.
static constexpr std::array< int, 5 > finalStatePDGs
PDG codes accepted as charged final state particles in generation: {e, mu, pi, K, p}...
std::map< std::string, int > mode_decay_minus
Mapping for tau- decays.
std::map< std::string, int > mode_decay_plus
Mapping for tau+ decays.
std::vector< int > vec_a0p_1450
Variable name of the vector where particles identified in the event are stored.
STL iterator class.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
Abstract base class for different kinds of events.
STL namespace.