ROOTPWA
waveSetGenerator.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2010
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-------------------------------------------------------------------------
22 // File and Version Information:
23 // $Rev:: $: revision of last commit
24 // $Author:: $: author of last commit
25 // $Date:: $: date of last commit
26 //
27 // Description:
28 // class that generates all possible waves given the final state
29 // particles, a list of isobars and constraints on I, J, L, and S
30 //
31 //
32 // Author List:
33 // Boris Grube TUM (original author)
34 //
35 //
36 //-------------------------------------------------------------------------
37 
38 
39 #include <boost/assign.hpp>
40 
41 #include "libconfig.h++"
42 
43 #include "physUtils.hpp"
44 #include "spinUtils.hpp"
45 #include "conversionUtils.hpp"
46 #include "libConfigUtils.hpp"
47 #include "arrayUtils.hpp"
48 #include "particleDataTable.h"
49 #include "waveDescription.h"
50 #include "waveSetGenerator.h"
51 
52 
53 using namespace std;
54 using namespace boost;
55 using namespace boost::tuples;
56 using namespace libconfig;
57 using namespace rpwa;
58 
59 
60 bool waveSetGenerator::_debug = false;
61 
62 
63 waveSetGenerator::waveSetGenerator()
64 {
65  reset();
66 }
67 
68 
69 waveSetGenerator::~waveSetGenerator()
70 { }
71 
72 
73 bool
74 waveSetGenerator::setWaveSetParameters(const string& templateKeyFileName)
75 {
76  // construct template decay topology
77  _templateTopo.reset();
78  waveDescription waveDesc;
79  if ( not waveDesc.parseKeyFile(templateKeyFileName)
80  or not waveDesc.constructDecayTopology(_templateTopo, true)) {
81  printWarn << "problems constructing template decay topology from key file "
82  << "'" << templateKeyFileName << "'. cannot generate wave set." << endl;
83  return false;
84  }
85  // read wave set parameters from template key file
86  printInfo << "reading wave set parameters from key file '" << templateKeyFileName << "'" << endl;
87  Config key;
88  if (not parseLibConfigFile(templateKeyFileName, key, _debug)) {
89  printWarn << "problems reading wave set parameters from '" << templateKeyFileName << "'. "
90  << "cannot generate wave set." << endl;
91  return false;
92  }
93  // find and parse group with wave set parameters
94  const Setting* waveSetParKey = findLibConfigGroup(key.getRoot(), "waveSetParameters", false);
95  if (waveSetParKey) {
96  const Setting* isoSpinRangeKey = findLibConfigArray(*waveSetParKey, "isospinRange", false);
97  if (isoSpinRangeKey)
98  _isospinRange = make_tuple<int, int>((*isoSpinRangeKey)[0], (*isoSpinRangeKey)[1]);
99  const Setting* JRangeKey = findLibConfigArray(*waveSetParKey, "JRange", false);
100  if (JRangeKey)
101  _JRange = make_tuple<int, int>((*JRangeKey)[0], (*JRangeKey)[1]);
102  const Setting* MRangeKey = findLibConfigArray(*waveSetParKey, "MRange", false);
103  waveSetParKey->lookupValue("reflectivity", _reflectivity );
104  waveSetParKey->lookupValue("useReflectivity", _useReflectivity );
105  waveSetParKey->lookupValue("allowSpinExotics", _allowSpinExotics);
106  if (MRangeKey)
107  _spinProjRange = make_tuple<int, int>((*MRangeKey)[0], (*MRangeKey)[1]);
108  const Setting* LRangeKey = findLibConfigArray(*waveSetParKey, "LRange", false);
109  if (LRangeKey)
110  _LRange = make_tuple<int, int>((*LRangeKey)[0], (*LRangeKey)[1]);
111  const Setting* SRangeKey = findLibConfigArray(*waveSetParKey, "SRange", false);
112  if (SRangeKey)
113  _SRange = make_tuple<int, int>((*SRangeKey)[0], (*SRangeKey)[1]);
114  const Setting* isobarBlackListKey = findLibConfigArray(*waveSetParKey,
115  "isobarBlackList", false);
116  if (isobarBlackListKey) {
117  _isobarBlackList.clear();
118  for (int i = 0; i < isobarBlackListKey->getLength(); ++i)
119  _isobarBlackList.push_back((*isobarBlackListKey)[i]);
120  }
121  const Setting* isobarWhiteListKey = findLibConfigArray(*waveSetParKey,
122  "isobarWhiteList", false);
123  if (isobarWhiteListKey) {
124  _isobarWhiteList.clear();
125  for (int i = 0; i < isobarWhiteListKey->getLength(); ++i)
126  _isobarWhiteList.push_back((*isobarWhiteListKey)[i]);
127  }
128  waveSetParKey->lookupValue("requireMinIsobarMass", _requireMinIsobarMass );
129  waveSetParKey->lookupValue("isobarMassWindowSigma", _isobarMassWindowSigma);
130  }
131  if (_debug)
132  printDebug << "parameters of " << *this;
133  return true;
134 }
135 
136 
137 size_t
138 waveSetGenerator::generateWaveSet()
139 {
140  if (not _templateTopo) {
141  printWarn << "template decay topology was not yet constructed. "
142  << "did you call waveSetGenerator::setWaveSetParameters()?" << endl;
143  return 0;
144  }
145  if (not _templateTopo->checkTopology()) {
146  printWarn << "ill-formed template decay topology. cannot generate wave set." << endl;
147  return 0;
148  }
149  // order nodes depth-first
150  vector<nodeDesc> startNds
151  = _templateTopo->sortNodesDfs(_templateTopo->XIsobarDecayVertex());
152  // create decay topologies of all subdecays
153  vector<isobarDecayTopology> subDecays(startNds.size());
154  for (size_t i = 0; i < startNds.size(); ++i)
155  subDecays[i] = _templateTopo->subDecay(startNds[i]);
156  // reverse order, because decay trees are build starting from the final state nodes
157  reverse(startNds.begin (), startNds.end ());
158  reverse(subDecays.begin(), subDecays.end());
159 
160  // create all possible subdecays
161  map<nodeDesc, vector<isobarDecayTopology> > decayPossibilities; // all possible decay graphs starting at certain node
162  for (size_t iStart = 0; iStart < startNds.size(); ++iStart) {
163 
164  // special case for final state nodes
165  if (_templateTopo->isFsVertex(_templateTopo->vertex(startNds[iStart]))) {
166  // final state vertices have no daughter tracks; subdecay
167  // topology consist of just the final state vertex
168  decayPossibilities[startNds[iStart]].push_back(subDecays[iStart]);
169  continue;
170  }
171 
172  if (_debug)
173  printDebug << "generating decay topologies for subdecay[" << iStart << "]: "
174  << subDecays[iStart];
175 
176  // get daughter decay topologies
177  vector<vector<isobarDecayTopology>* > daughterDecays;
178  adjIterator iNd, iNdEnd;
179  for (tie(iNd, iNdEnd) = _templateTopo->adjacentVertices(startNds[iStart]); iNd != iNdEnd; ++iNd)
180  daughterDecays.push_back(&decayPossibilities[*iNd]);
181 
182  // loop over all combinations of daughter decays
183  size_t iDaughter[2];
184  for (iDaughter[0] = 0; iDaughter[0] < daughterDecays[0]->size(); ++iDaughter[0])
185  for (iDaughter[1] = 0; iDaughter[1] < daughterDecays[1]->size(); ++iDaughter[1]) {
186 
187  // get daughter particles from the respective decay topologies
188  const nodeDesc topNodes[2] =
189  {(*daughterDecays[0])[iDaughter[0]].topNode(),
190  (*daughterDecays[1])[iDaughter[1]].topNode()};
191  const particlePtr daughters[2] =
192  {(*daughterDecays[0])[iDaughter[0]].vertex(topNodes[0])->inParticles()[0],
193  (*daughterDecays[1])[iDaughter[1]].vertex(topNodes[1])->inParticles()[0]};
194  if (_debug)
195  printDebug << "combining decay[" << iDaughter[0] << "] of first daughter "
196  << "'" << daughters[0]->name() << "' "
197  << "with decay[" << iDaughter[1] << "] of second daughter "
198  << "'" << daughters[1]->name() << "'" << endl
199  << (*daughterDecays[0])[iDaughter[0]]
200  << (*daughterDecays[1])[iDaughter[1]];
201 
202  // copy parent vertex
203  isobarDecayVertexPtr parentVertex
204  (new isobarDecayVertex(*static_pointer_cast<isobarDecayVertex>
205  (_templateTopo->vertex(startNds[iStart]))));
206  particlePtr parent = parentVertex->parent();
207  // set daughters of parent vertex
208  parentVertex->daughter1() = daughters[0];
209  parentVertex->daughter2() = daughters[1];
210  // join daughter subdecays and parent vertex
211  isobarDecayTopology parentDecay
212  = _templateTopo->joinDaughterDecays(parentVertex,
213  (*daughterDecays[0])[iDaughter[0]],
214  (*daughterDecays[1])[iDaughter[1]]);
215 
216  // calculate parent quantum numbers fixed by daughter quantum numbers
217  const int parentBaryonNmb = daughters[0]->baryonNmb() + daughters[1]->baryonNmb();
218  const int parentCharge = daughters[0]->charge() + daughters[1]->charge();
219  const int parentStrangeness = daughters[0]->strangeness() + daughters[1]->strangeness();
220  const int parentCharm = daughters[0]->charm() + daughters[1]->charm();
221  const int parentBeauty = daughters[0]->beauty() + daughters[1]->beauty();
222  const int parentG = daughters[0]->G() * daughters[1]->G();
223  // daughter quantum numbers that define ranges of possible parent quantum numbers
224  const int daughterI [2] = {daughters[0]->isospin(), daughters[1]->isospin() };
225  const int daughterI_z[2] = {daughters[0]->isospinProj(), daughters[1]->isospinProj()};
226  const int daughterJ [2] = {daughters[0]->J(), daughters[1]->J() };
227 
228  // loop over all allowed combinations of daughter quantum numbers
229  // loop over allowed total spins
230  int S, SMax;
231  for (tie(S, SMax) = getSpinRange(daughterJ[0], daughterJ[1], _SRange); S <= SMax; S += 2) {
232  if (_debug)
233  printDebug << "setting total spin = " << spinQn(S) << " "
234  << "(max spin = " << spinQn(SMax) << ")" << endl;
235 
236  // loop over allowed relative orbital angular momenta
237  for (int L = max(0, get<0>(_LRange)); L <= get<1>(_LRange); L += 2) {
238  if (_debug)
239  printDebug << "setting relative orbital angular momentum = " << spinQn(L) << " "
240  << "(max L = " << spinQn(get<1>(_LRange)) << ")" << endl;
241  const int parentP = daughters[0]->P() * daughters[1]->P() * powMinusOne(L / 2);
242 
243  // loop over allowed total angular momenta
244  int parentJ, parentJMax;
245  for (tie(parentJ, parentJMax) = getSpinRange(L, S, _JRange);
246  parentJ <= parentJMax; parentJ += 2) {
247  if (_debug)
248  printDebug << "setting parent spin = " << spinQn(parentJ) << " "
249  << "(max spin = " << spinQn(parentJMax) << ")" << endl;
250 
251  // loop over allowed isospins
252  int parentI, parentIMax;
253  for (tie(parentI, parentIMax) = getSpinRange(daughterI[0], daughterI[1], _isospinRange);
254  parentI <= parentIMax; parentI += 2) {
255  if (_debug)
256  printDebug << "setting parent isospin = " << spinQn(parentI) << " "
257  << "(max isospin = " << spinQn(parentIMax) << ")" << endl;
258 
259  // get I_z from Gell-Mann-Nishijima formula (see PDG 2008 eq. 14.1)
260  const int parentI_z = 2 * parentCharge - (parentBaryonNmb + parentStrangeness
261  + parentCharm + parentBeauty);
262 
263  // C-parity
264  // defined only for mesons with zero flavor quantum numbers
265  int parentC = 0;
266  if ((parentBaryonNmb == 0) and (parentI_z == 0) and (parentStrangeness == 0)
267  and (parentCharm == 0) and (parentBeauty == 0))
268  parentC = parentG * powMinusOne(parentI / 2);
269  if (_debug)
270  printDebug << "trying isobar quantum numbers IG(JPC) = "
271  << spinQn(parentI) << sign(parentG)
272  << "(" << spinQn(parentJ) << sign(parentP) << sign(parentC) << ")"
273  << ", with L = " << spinQn(L) << ", S = " << spinQn(S) << endl;
274 
275  // check whether charge state is allowed using
276  if (not spinStatesCanCouple(daughterI, daughterI_z, parentI, parentI_z)) {
277  if (_debug)
278  cout << " rejected, because (I, I_z)[0] = (" << spinQn(daughterI[0])
279  << ", " << spinQn(daughterI_z[0]) << ") and (I, I_z)[1] = ("
280  << spinQn(daughterI[1]) << ", " << spinQn(daughterI_z[1])
281  << ") cannot couple to parent (I, I_z) = ("
282  << spinQn(parentI) << ", " << spinQn(parentI_z) << ")" << endl;
283  continue;
284  }
285  if (not _allowSpinExotics)
286  // check whether quantum number combination is spin-exotic
287  if (igjpIsExotic(parentI, parentG, parentJ, parentP)) {
288  if (_debug)
289  cout << " rejected, because quantum number combination IG(JPC) = "
290  << spinQn(parentI) << sign(parentG) << "(" << spinQn(parentJ)
291  << sign(parentP) << sign(parentC) << ") is spin-exotic" << endl;
292  continue;
293  }
294 
295  // enforce rules for particle-antiparticle systems
296  if ( daughters[1]->antiPartProperties()
297  == *(static_pointer_cast<particleProperties>(daughters[0]))) {
298  // !NOTE! since equality test function is virtual,
299  // order in the above comparison is essential
300  // (see S.U. Chung: "C- and G-Parity: a New Definition and Applications, Version V"
301  // eq. 32)
302  if (parentG != powMinusOne((parentI + L + S) / 2)) {
303  if (_debug)
304  cout << " rejected, because I[= " << spinQn(parentI) << "], "
305  << "L[= " << spinQn(L) << "], and S[= " << spinQn(S) << "] "
306  << "are not consistent with G-parity[= " << sign(parentG) << "]. "
307  << "for particle-antiparticle system ["
308  << daughters[0]->name() << ", " << daughters[1]->name() << "] "
309  << "(-)^(I + L + S) = " << sign(powMinusOne((parentI + L + S) / 2))
310  << " must be equal to G-parity." << endl;
311  continue;
312  }
313  // I_z of the parent and C = (-)^(L + S) should be fulfilled automatically
314  // (see S.U. Chung: "C- and G-Parity: a New Definition and Applications, Version V"
315  // eq. 31)
316  assert(parentI_z == 0);
317  assert(parentC == powMinusOne((L + S) / 2));
318  }
319 
320  // enforce rules for particle-particle systems
321  // (see S.U. Chung: "C- and G-Parity: a New Definition and Applications, Version V"
322  // eq. 36)
323  if (*(daughters[0]) == *(daughters[1]))
324  if (powMinusOne((parentI + L + S) / 2) != powMinusOne(daughterI[0])) {
325  if (_debug)
326  cout << " rejected, because system with two identical particles ["
327  << daughters[0]->name() << ", " << daughters[1]->name() << "] must have "
328  << ((powMinusOne(daughterI[0]) == +1) ? "even" : "odd")
329  << " (I[= " << spinQn(parentI) << "] + L[= " << spinQn(L) << "] "
330  << "+ S[= " << spinQn(S) << "]) = "
331  << spinQn(parentI) + spinQn(L) + spinQn(S) << endl;
332  continue;
333  }
334 
335  // set isobar properties
336  particleProperties isobarProp(*parent);
337  isobarProp.setBaryonNmb(parentBaryonNmb);
338  isobarProp.setSCB (parentStrangeness, parentCharm, parentBeauty);
339  isobarProp.setIGJPC (parentI, parentG, parentJ, parentP, parentC);
340 
341  // create all allowed decay topologies
342  if (parent->isXParticle()) {
343  // for X loop over allowed M and reflectivity states
344  tuple<int, int> parentMRange =
345  make_tuple(spinQnLargerEqual((_useReflectivity) ?
346  spinQnSmallerEqual(parentJ, 0) : -parentJ,
347  get<0>(_spinProjRange)),
348  spinQnSmallerEqual(parentJ, get<1>(_spinProjRange)));
349  tuple<int, int> parentReflRange(0, 0);
350  if (_useReflectivity) {
351  if (_reflectivity == 0)
352  parentReflRange = make_tuple(-1, +1);
353  else
354  parentReflRange = make_tuple(signum(_reflectivity), signum(_reflectivity));
355  }
356  for (int parentM = get<0>(parentMRange); parentM <= get<1>(parentMRange);
357  parentM += 2) {
358  if (_debug)
359  printDebug << "setting X spin projection to " << spinQn(parentM) << " "
360  << "(max spin projection = " << spinQn(get<1>(parentMRange))
361  << ")" << endl;
362  for (int parentRefl = get<0>(parentReflRange);
363  parentRefl <= get<1>(parentReflRange); parentRefl += 2) {
364  if (_debug and (parentRefl != 0))
365  printDebug << "setting X reflectivity to "
366  << parityQn(parentRefl) << endl;
367 
368  // check whether amplitude would be zero for given reflectivity
369  if ( _useReflectivity and (parentM == 0)
370  and (reflectivityFactor(parentJ, parentP, parentM, parentRefl) == +1)) {
371  if (_debug)
372  printDebug << "rejected, because for quantum number combination "
373  << "JP(M refl) = " << spinQn(parentJ) << sign(parentP) << "("
374  << spinQn(parentM) << sign(parentRefl) << ") "
375  << "amplitude is always zero"<< endl;
376  continue;
377  }
378 
379  // clone topology
380  const isobarDecayTopologyPtr& decayCopy
381  = createNewDecayTopology(parentDecay, parentVertex, L, S,
382  isobarProp, parentCharge);
383  decayCopy->XDecayVertex()->inParticles()[0]->setSpinProj(parentM);
384  decayCopy->XDecayVertex()->inParticles()[0]->setReflectivity(parentRefl);
385  decayPossibilities[startNds[iStart]].push_back(*decayCopy);
386  }
387  }
388  } else {
389  // for intermediate decay vertices find isobar candidates in particle data table
390  const double minIsobarMass = (_requireMinIsobarMass) ?
391  (daughters[0]->mass() - _isobarMassWindowSigma * daughters[0]->width())
392  + (daughters[1]->mass() - _isobarMassWindowSigma * daughters[1]->width())
393  : 0;
394 
395  const particleProperties::decayMode
396  decay(assign::list_of(daughters[0]->name())(daughters[1]->name()), L, S);
397  vector<const particleProperties*> possibleIsobars
398  = particleDataTable::entriesMatching(isobarProp, "allQn",
399  minIsobarMass, _isobarMassWindowSigma,
400  _isobarWhiteList, _isobarBlackList,
401  decay, _forceDecayCheck);
402  if (_debug) {
403  if (possibleIsobars.empty())
404  printDebug << "rejected, because no matching isobar could be found for "
405  << isobarProp.qnSummary() << " in particle data table" << endl;
406  else
407  printDebug << "found " << possibleIsobars.size() << " isobar candidate(s) for "
408  << isobarProp.qnSummary() << " in particle data table" << endl;
409  }
410 
411  // loop over isobar candidates
412  for (size_t iIsobar = 0; iIsobar < possibleIsobars.size(); ++iIsobar) {
413  if (_debug)
414  printDebug << "setting isobar = '"
415  << possibleIsobars[iIsobar]->name() << "'" << endl;
416  // clone topology
417  const isobarDecayTopologyPtr& decayCopy
418  = createNewDecayTopology(parentDecay, parentVertex, L, S,
419  *possibleIsobars[iIsobar], parentCharge);
420  decayPossibilities[startNds[iStart]].push_back(*decayCopy);
421  }
422  }
423  } // isospin loop
424  } // L-S coupling loop
425  } // L loop
426  } // S loop
427  } // loop over daughter decays
428  } // loop over all start nodes
429 
430  // extract decays for X-decay vertex and add production vertex
431  _waveSet = decayPossibilities[startNds.back()];
432  for (size_t i = 0; i < _waveSet.size(); ++i) {
433  // clone production vertex and set X-particle
434  const productionVertexPtr newProdVert(static_pointer_cast<rpwa::productionVertex>
435  (_templateTopo->productionVertex()->clone(false, false)));
436  const particlePtr& newX = _waveSet[i].XIsobarDecayVertex()->parent();
437  newProdVert->outParticles()[0] = newX;
438  // add production vertex
439  _waveSet[i].setProductionVertex(newProdVert);
440  // correct quantum numbers
441  _waveSet[i].calcIsobarCharges();
442  _waveSet[i].productionVertex()->setXFlavorQN(); // sets baryon nmb, S, C, and B of X
443  }
444 
445  // post process wave set
446  set<size_t> boseSymWaveIndices = findBoseSymDecays();
447  if (boseSymWaveIndices.size() > 0) {
448  _waveSet.erase(remove_if(_waveSet.begin(), _waveSet.end(),
449  rpwa::isInListOfIndices<isobarDecayTopology>(_waveSet.begin(),
450  boseSymWaveIndices)),
451  _waveSet.end());
452  printSucc << "removed " << boseSymWaveIndices.size() << " Bose duplicates "
453  << "from wave set" << endl;
454  }
455  return _waveSet.size();
456 }
457 
458 
459 bool
460 waveSetGenerator::writeKeyFiles(const string& dirName,
461  const bool newKeyFileNameConvention)
462 {
463  size_t countSuccess = 0;
464  for (size_t i = 0; i < _waveSet.size(); ++i) {
465  const string keyFileName = dirName + "/"
466  + waveDescription::waveNameFromTopology(_waveSet[i], newKeyFileNameConvention) + ".key";
467  if (waveDescription::writeKeyFile(keyFileName, _waveSet[i]))
468  ++countSuccess;
469  }
470  printInfo << "wrote " << countSuccess << " out of " << _waveSet.size() << " key files" << endl;
471  if (countSuccess != _waveSet.size()) {
472  printWarn << "writing of " << _waveSet.size() - countSuccess << " key files failed" << endl;
473  return false;
474  }
475  return true;
476 }
477 
478 
479 void
480 waveSetGenerator::reset()
481 {
482  _isospinRange = make_tuple(0, 2);
483  _JRange = make_tuple(0, 0);
484  _spinProjRange = make_tuple(0, 0);
485  _reflectivity = 0;
486  _useReflectivity = false;
487  _allowSpinExotics = false;
488  _LRange = make_tuple(0, 0);
489  _SRange = make_tuple(0, 0);
490  _requireMinIsobarMass = false;
491  _forceDecayCheck = false;
492  _isobarMassWindowSigma = 0;
493  _isobarBlackList.clear();
494  _isobarWhiteList.clear();
495  _templateTopo.reset();
496  _waveSet.clear();
497 }
498 
499 
500 ostream&
501 waveSetGenerator::print(ostream& out) const
502 {
503  out << "wave set generator:" << endl
504  << " isospin range .............. [" << spinQn(get<0>(_isospinRange)) << ", "
505  << spinQn(get<1>(_isospinRange)) << "]" << endl
506  << " J range .................... [" << spinQn(get<0>(_JRange)) << ", "
507  << spinQn(get<1>(_JRange)) << "]" << endl
508  << " M range .................... [" << spinQn(get<0>(_spinProjRange)) << ", "
509  << spinQn(get<1>(_spinProjRange)) << "]" << endl
510  << " allow spin-exotics ......... " << yesNo(_allowSpinExotics) << endl;
511  if (_useReflectivity)
512  out << " generate reflectivity ...... " << ((_reflectivity == 0) ? "both"
513  : sign(_reflectivity)) << endl;
514  out << " L range .................... [" << spinQn(get<0>(_LRange)) << ", "
515  << spinQn(get<1>(_LRange)) << "]" << endl
516  << " S range .................... [" << spinQn(get<0>(_SRange)) << ", "
517  << spinQn(get<1>(_SRange)) << "]" << endl
518  << " require min. isobar mass ... " << yesNo(_requireMinIsobarMass) << endl
519  << " isobar mass window par. .... " << _isobarMassWindowSigma << " [Gamma]" << endl
520  << " force decay checks ......... " << yesNo(_forceDecayCheck) << endl;
521  out << " isobar black list:";
522  if (_isobarBlackList.size() == 0)
523  out << " empty" << endl;
524  else {
525  out << endl;
526  for (size_t i = 0; i < _isobarBlackList.size(); ++i)
527  out << " " << _isobarBlackList[i] << endl;
528  }
529  out << " isobar white list:";
530  if (_isobarWhiteList.size() == 0)
531  out << " all isobars" << endl;
532  else {
533  out << endl;
534  for (size_t i = 0; i < _isobarWhiteList.size(); ++i)
535  out << " " << _isobarWhiteList[i] << endl;
536  }
537  if (_templateTopo)
538  out << "template " << *_templateTopo;
539  return out;
540 }
541 
542 
543 const isobarDecayTopologyPtr
544 waveSetGenerator::createNewDecayTopology(const isobarDecayTopology& parentDecay,
545  const isobarDecayVertexPtr& parentVertex,
546  const int L,
547  const int S,
548  const particleProperties& isobar,
549  const int parentCharge)
550 {
551  // clone topology
552  const isobarDecayTopologyPtr& decayCopy = parentDecay.clone(false, false);
553  // set parent vertex quantum numbers
554  const isobarDecayVertexPtr& parentVertexCopy
555  = static_pointer_cast<isobarDecayVertex>(decayCopy->vertex(parentDecay.node(parentVertex)));
556  parentVertexCopy->setL(L);
557  parentVertexCopy->setS(S);
558  // set parent particle
559  const particlePtr& parentCopy = parentVertexCopy->parent();
560  parentCopy->setProperties(isobar);
561  parentCopy->setCharge(parentCharge);
562  if (_debug)
563  printDebug << "created decay topology for " << *parentVertexCopy << endl;
564  return decayCopy;
565 }
566 
567 
568 set<size_t>
569 waveSetGenerator::findBoseSymDecays() const
570 {
571  printInfo << "looking for Bose duplicates in wave set..." << endl;
572  set<size_t> waveIndicesToRemove;
573  set<size_t> waveIndicesToKeep;
574  for (size_t waveIndex = 0; waveIndex < _waveSet.size(); ++waveIndex) {
575 
576  // don't consider entries that are already tagged
577  set<size_t>::const_iterator indexEntry = waveIndicesToRemove.find(waveIndex);
578  if (indexEntry != waveIndicesToRemove.end())
579  continue;
580  indexEntry = waveIndicesToKeep.find(waveIndex);
581  if (indexEntry != waveIndicesToKeep.end())
582  continue;
583 
584  // get Bose-symmetric vertices
585  isobarDecayTopologyPtr topo = _waveSet[waveIndex].clone();
586  const vector<unsigned int> boseSymVertIds = topo->findIsobarBoseSymVertices();
587  const size_t nmbBoseSymVertIds = boseSymVertIds.size();
588  if (nmbBoseSymVertIds > 1) {
589  printErr << "the code has not yet been tested for decay topologies with "
590  << nmbBoseSymVertIds << " Bose-symmetric decay vertices. aborting." << endl;
591  throw;
592  }
593  if (nmbBoseSymVertIds != 1)
594  continue;
595 
596  // make sure the two daughters are not the same particle
597  isobarDecayVertexPtr vert = topo->isobarDecayVertices()[boseSymVertIds[0]];
598  const particlePtr daughters[2] = {vert->daughter1(), vert->daughter2()};
599  if (*(daughters[0]) == *(daughters[1]))
600  continue;
601  const decayTopology daughterTopos[2] = {topo->subDecay(topo->toNode(daughters[0])),
602  topo->subDecay(topo->toNode(daughters[1]))};
603  if (_debug)
604  printDebug << "wave[" << waveIndex << "] "
605  << waveDescription::waveNameFromTopology(*topo, true)
606  << " has Bose-symmetric vertex " << *vert << ". "
607  << "looking for partner waves in wave set." << endl;
608 
609  // find indices of all decay vertices that are not below the Bose-symmetric vertex
610  vector<unsigned int> nonSymVertIds;
611  {
612  // get graph nodes below Bose-symmetric vertex
613  const vector<nodeDesc> nodes = topo->sortNodesDfs(vert);
614  for (unsigned int i = 0; i < topo->nmbDecayVertices(); ++i) {
615  const nodeDesc nd = topo->node(topo->isobarDecayVertices()[i]);
616  bool isSymNode = false;
617  for (unsigned int j = 0; j < nodes.size(); ++j)
618  if (nd == nodes[j]) {
619  isSymNode = true;
620  break;
621  }
622  if (not isSymNode)
623  nonSymVertIds.push_back(i);
624  }
625  }
626 
627  // look for Bose-symmetric topology
628  // this uses the fact that all waves have the same decay toplogy,
629  // because they were generated from the same template
630  for (size_t i = 0; i < _waveSet.size(); ++i) {
631  if (i == waveIndex)
632  continue;
633  isobarDecayTopologyPtr symTopo = _waveSet[i].clone();
634 
635  // compare vertices that are not below the Bose-symmetric vertex
636  bool foundSymTopo = true;
637  for (size_t j = 0; j < nonSymVertIds.size(); ++j)
638  if ( *(topo->isobarDecayVertices ()[nonSymVertIds[j]])
639  != *(symTopo->isobarDecayVertices()[nonSymVertIds[j]])) {
640  foundSymTopo = false;
641  break;
642  }
643  if (not foundSymTopo)
644  continue;
645  // compare Bose-symmetric vertex
646  isobarDecayVertexPtr symVert = symTopo->isobarDecayVertices()[boseSymVertIds[0]];
647  // swap daughters 1 and 2
648  swap(symVert->outParticles()[0], symVert->outParticles()[1]);
649  if (*symVert != *vert)
650  continue;
651  swap(symVert->outParticles()[1], symVert->outParticles()[0]); // swap back, because otherwise printout is confusing
652  // compare daughter topologies
653  const particlePtr symDaughters [2] = {symVert->daughter2(), symVert->daughter1()}; // !NOTE! reversed order
654  const decayTopology symDaughterTopos[2]
655  = {symTopo->subDecay(symTopo->toNode(symDaughters[0])),
656  symTopo->subDecay(symTopo->toNode(symDaughters[1]))};
657  for (unsigned int k = 0; k < 2; ++k)
658  for (unsigned int l = 0; l < daughterTopos[k].nmbDecayVertices(); ++l)
659  if ( *(daughterTopos [k].decayVertices()[l])
660  != *(symDaughterTopos[k].decayVertices()[l])) {
661  foundSymTopo = false;
662  break;
663  }
664  if (not foundSymTopo)
665  continue;
666  if (_debug)
667  printDebug << "found Bose-partner wave[" << i << "] "
668  << waveDescription::waveNameFromTopology(*symTopo, true) << endl;
669  // keep topology with the heavier first particle
670  if (daughters[0]->mass() > daughters[1]->mass()) {
671  if ( not waveIndicesToRemove.insert(i).second
672  or not waveIndicesToKeep.insert(waveIndex).second) {
673  printErr << "wave[" << i << "]([" << waveIndex << "]) was tagged for "
674  << "removal(for keeping) twice. this should never happen. aborting." << endl;
675  throw;
676  }
677  printInfo << "tagged wave " << waveDescription::waveNameFromTopology(*symTopo, true)
678  << " as superfluous" << endl;
679  } else {
680  if ( not waveIndicesToRemove.insert(waveIndex).second
681  or not waveIndicesToKeep.insert(i).second) {
682  printErr << "wave[" << waveIndex << "]([" << i << "]) was tagged for "
683  << "removal(for keeping) twice. this should never happen. aborting." << endl;
684  throw;
685  }
686  printInfo << "tagged wave " << waveDescription::waveNameFromTopology(*topo, true)
687  << " as superfluous" << endl;
688  }
689  } // inner loop over wave set
690 
691  } // outer loop over wave set
692 
693  printSucc << "found " << waveIndicesToRemove.size() << " Bose duplicates "
694  << "in wave set" << endl;
695  return waveIndicesToRemove;
696 }
Generated by   doxygen 1.8.1.2