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