ROOTPWA
decayTopology.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 // container class that holds all external information for
29 // amplitude calculation
30 // internally the decay process is represented as a graph
31 // the graph is constraint to contain exactly one production
32 // vertex and at least one interaction vertex; in addtion for
33 // each final-state particle a corresponding final-state vertex
34 // is created for internal use
35 //
36 // "final-state" particles are the measured decay daughters;
37 // additional final-state particles that belong to the production
38 // process are handled by the production vertex
39 //
40 //
41 // Author List:
42 // Boris Grube TUM (original author)
43 //
44 //
45 //-------------------------------------------------------------------------
46 
47 
48 #include <vector>
49 #include <map>
50 #include <algorithm>
51 
52 #include "TClonesArray.h"
53 #include "TClass.h"
54 #include "TObjString.h"
55 #include "TVector3.h"
56 
57 #include "reportingUtilsRoot.hpp"
58 #include "conversionUtils.hpp"
59 #include "decayTopology.h"
60 
61 
62 using namespace std;
63 using namespace boost;
64 using namespace rpwa;
65 
66 
67 bool decayTopology::_debug = false;
68 
69 
70 decayTopology::decayTopology()
71  : decayTopologyGraphType(),
72  _prodVertex (),
73  _decayVertices (),
74  _fsParticles (),
75  _fsDataPartIndexMap (),
76  _fsDataPartMomCache ()
77 {
78 }
79 
80 
81 decayTopology::decayTopology(const productionVertexPtr& productionVertex,
82  const vector<interactionVertexPtr>& decayVertices,
83  const vector<particlePtr>& fsParticles,
84  const bool performTopologyCheck)
85 {
86  constructDecay(productionVertex, decayVertices, fsParticles, performTopologyCheck);
87 }
88 
89 
90 decayTopology::decayTopology(const decayTopology& topo)
91 {
92  *this = topo;
93 }
94 
95 
96 decayTopology::decayTopology(const decayTopologyGraphType& graph)
97 {
98  *this = graph;
99 }
100 
101 
102 decayTopology::~decayTopology()
103 { }
104 
105 
106 decayTopology&
107 decayTopology::operator =(const decayTopology& topo)
108 {
109  if (this != &topo) {
110  decayTopologyGraphType::operator =(topo);
111  _prodVertex = topo._prodVertex;
112  _decayVertices = topo._decayVertices;
113  _fsParticles = topo._fsParticles;
114  _fsDataPartIndexMap = topo._fsDataPartIndexMap;
115  _fsDataPartMomCache = topo._fsDataPartMomCache;
116  }
117  return *this;
118 }
119 
120 
121 decayTopology&
122 decayTopology::operator =(const decayTopologyGraphType& graph)
123 {
124  if (this != &graph) {
125  decayTopologyGraphType::operator =(graph);
126  buildInternalData();
127  }
128  return *this;
129 }
130 
131 
132 decayTopology*
133 decayTopology::doClone(const bool cloneFsParticles,
134  const bool cloneProdKinematics) const
135 {
136  if (_debug)
137  printDebug << "cloning decay topology '" << name() << "' "
138  << ((cloneFsParticles ) ? "in" : "ex") << "cluding final-state particles, "
139  << ((cloneProdKinematics) ? "in" : "ex") << "cluding production kinematics particles"
140  << endl;
141  // copy graph data structure
142  decayTopology* topoClone = new decayTopology(*this);
143  // clone nodes
144  nodeIterator iNd, iNdEnd;
145  vector<interactionVertexPtr> newVerts;
146  for (tie(iNd, iNdEnd) = topoClone->nodes(); iNd != iNdEnd; ++iNd) {
147  const interactionVertexPtr& vert = topoClone->vertex(*iNd);
148  // clone vertex
149  if ( isDecayVertex(vert)
150  or isProductionVertex(vert)
151  or (cloneFsParticles and isFsVertex(vert)))
152  topoClone->isDecayVertex(vert);
153  newVerts.push_back(topoClone->cloneNode(*iNd));
154  }
155  // looping over incoming particles of respective vertices and clone them
156  // this also clones dangling particles that have no associated edge
158  for (unsigned int i = 0; i < newVerts.size(); ++i) {
159  if (not topoClone->isProductionVertex(newVerts[i]) or cloneProdKinematics)
160  for (unsigned int j = 0; j < newVerts[i]->nmbInParticles(); ++j) {
161  const particlePtr part = newVerts[i]->inParticles()[j];
162  if (not topoClone->isFsParticle(part) or cloneFsParticles) {
163  if (topoClone->isEdge(part))
164  topoClone->cloneEdge(topoClone->edge(part));
165  else {
166  // dangling particle
167  if (_debug)
168  printDebug << "cloning dangling " << *part << endl;
169  // clone particle and replace it in array of incoming particles in vertex
170  particlePtr newPart(part->clone());
171  newVerts[i]->inParticles()[j] = newPart;
172  }
173  }
174  }
175  }
176  topoClone->name() = "topoClone";
177  return topoClone;
178 }
179 
180 
181 void
182 decayTopology::clear()
183 {
184  decayTopologyGraphType::clear();
185  _prodVertex.reset();
186  _decayVertices.clear();
187  _fsParticles.clear();
188  _fsDataPartIndexMap.clear();
189  _fsDataPartMomCache.clear();
190 }
191 
192 
193 map<string, unsigned int>
194 decayTopology::nmbIndistFsParticles() const
195 {
196  map<string, unsigned int> partMult;
197  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
198  const string partName = fsParticles()[i]->name();
199  map<string, unsigned int>::iterator entry = partMult.find(partName);
200  if (entry != partMult.end())
201  ++(entry->second);
202  else
203  partMult[partName] = 1;
204  }
205  return partMult;
206 }
207 
208 
209 int
210 decayTopology::fsParticlesIntrinsicParity() const
211 {
212  // calculate product of final-state particles
213  int intrinsicParity = 1;
214  for (unsigned int i = 0; i < nmbFsParticles(); ++i)
215  intrinsicParity *= fsParticles()[i]->P();
216  return intrinsicParity;
217 }
218 
219 
220 int
221 decayTopology::spaceInvEigenValue() const
222 {
223  // parity of X is P = P_spatial * P_intrinsic; we want to know P_spatial
224  return XParticle()->P() / fsParticlesIntrinsicParity();
225 }
226 
227 
228 int
229 decayTopology::reflectionEigenValue() const
230 {
231  // eigenvalue of reflection through production plane is r_spatial / P_intrinsic
232  // where r = -1 / refl is the reflectivity eiganvalue
233  return -1 / (XParticle()->reflectivity() * fsParticlesIntrinsicParity());
234 }
235 
236 
237 void
238 decayTopology::transformFsParticles(const TLorentzRotation& L)
239 {
240  for (unsigned int i = 0; i < nmbFsParticles(); ++i)
241  fsParticles()[i]->transform(L);
242 }
243 
244 
245 bool
246 decayTopology::isVertex(const interactionVertexPtr& vert) const
247 {
248  if (not vert)
249  return false;
250  nodeIterator iNd, iNdEnd;
251  for (tie(iNd, iNdEnd) = nodes(); iNd != iNdEnd; ++iNd)
252  if (vert == vertex(*iNd))
253  return true;
254  return false;
255 }
256 
257 
258 bool
259 decayTopology::isParticle(const particlePtr& part) const
260 {
261  if (not part)
262  return false;
263  edgeIterator iEd, iEdEnd;
264  for (tie(iEd, iEdEnd) = edges(); iEd != iEdEnd; ++iEd)
265  if (part == particle(*iEd))
266  return true;
267  return false;
268 }
269 
270 
271 bool
272 decayTopology::isProductionVertex(const interactionVertexPtr& vert) const
273 {
274  if (not vert or not isVertex(vert))
275  return false;
276  return (vert == _prodVertex);
277 }
278 
279 
280 bool
281 decayTopology::isDecayVertex(const interactionVertexPtr& vert) const
282 {
283  if (not vert)
284  return false;
285  if (not isVertex(vert))
286  return false;
287  if (isProductionVertex(vert) or isFsVertex(vert))
288  return false;
289  return true;
290 }
291 
292 
293 int
294 decayTopology::decayVertexIndex(const interactionVertexPtr& vert) const
295 {
296  if (not vert or not isDecayVertex(vert))
297  return -1;
298  for (unsigned int i = 0; i < nmbDecayVertices(); ++i)
299  if (vert == decayVertices()[i])
300  return i;
301  return -1;
302 }
303 
304 
305 bool
306 decayTopology::isFsVertex(const interactionVertexPtr& vert) const
307 {
308  if (not vert)
309  return false;
310  if (dynamic_pointer_cast<fsVertex>(vert))
311  return isVertex(vert);
312  return false;
313 }
314 
315 
316 bool
317 decayTopology::isFsParticle(const particlePtr& part) const
318 {
319  if (not part or not isParticle(part))
320  return false;
321  for (unsigned int i = 0; i < nmbFsParticles(); ++i)
322  if (part == fsParticles()[i])
323  return true;
324  return false;
325 }
326 
327 
328 int
329 decayTopology::fsParticlesIndex(const particlePtr& part) const
330 {
331  if (not part or not isParticle(part))
332  return -1;
333  for (unsigned int i = 0; i < nmbFsParticles(); ++i)
334  if (part == fsParticles()[i])
335  return i;
336  return -1;
337 }
338 
339 
340 bool
341 decayTopology::checkTopology() const
342 {
343  bool topologyIsOkay = true;
344  // make sure there are no dangling particles
345  nodeIterator iNd, iNdEnd;
346  for (tie(iNd, iNdEnd) = nodes(); iNd != iNdEnd; ++iNd) {
347  const interactionVertexPtr& vert = vertex(*iNd);
348  if (not vert) {
349  printWarn << "node[" << *iNd << "] has vertex null pointer" << endl;
350  topologyIsOkay = false;
351  continue;
352  }
353  if (isProductionVertex(vert)) {
354  // for production vertex only X particle has associated edge
355  if (not isEdge(static_pointer_cast<rpwa::productionVertex>(vert)->XParticle())) {
356  printWarn << "X particle of " << *vert << " has no associated edge" << endl;
357  topologyIsOkay = false;
358  } else if (_debug)
359  printDebug << "success: X particle of " << *vert << " has associated edge" << endl;
360  } else {
361  for (unsigned int i = 0; i < vert->nmbInParticles(); ++i)
362  if (not isEdge(vert->inParticles()[i])) {
363  printWarn << "incoming particle[" << i << "] of " << *vert
364  << " has no associated edge" << endl;
365  topologyIsOkay = false;
366  } else if (_debug)
367  printDebug << "success: incoming particle[" << i << "] of " << *vert
368  << " has associated edge" << endl;
369  for (unsigned int i = 0; i < vert->nmbOutParticles(); ++i)
370  if (not isEdge(vert->outParticles()[i])) {
371  printWarn << "outgoing particle[" << i << "] of " << *vert
372  << " has no associated edge" << endl;
373  topologyIsOkay = false;
374  } else if (_debug)
375  printDebug << "success: outgoing particle[" << i << "] of " << *vert
376  << " has associated edge" << endl;
377  }
378  }
379  // check production vertex
380  if (not productionVertex()) {
381  printWarn << "production vertex is null pointer" << endl;
382  topologyIsOkay = false;
383  } else {
384  if (not isNode(_prodVertex) or (vertex(node(_prodVertex)) != _prodVertex)) {
385  printWarn << *_prodVertex << " has no associated node" << endl;
386  topologyIsOkay = false;
387  } else if (_debug)
388  printDebug << "success: " << *_prodVertex << " has associated node" << endl;
389  if (nmbOutEdges(_prodVertex) != 1) {
390  printWarn << *_prodVertex << " has " << nmbOutEdges(_prodVertex) << " != 1 "
391  << "outgoing edges" << endl;
392  topologyIsOkay = false;
393  } else if (nmbInEdges(_prodVertex) != 0) {
394  printWarn << *_prodVertex << " has " << nmbInEdges(_prodVertex) << " != 0 "
395  << "incoming edges" << endl;
396  topologyIsOkay = false;
397  } else if (_debug)
398  printDebug << "success: " << *_prodVertex
399  << " has exactly 1 outgoing and no incoming edges" << endl;
400  }
401  // make sure that all nodes are reachable from top node
402  vector<nodeDesc> sortedNds;
403  unsigned int nmbVertices = nmbDecayVertices() + nmbFsParticles();
404  if (productionVertex()) {
405  sortedNds = sortNodesDfs(node(productionVertex()));
406  ++nmbVertices;
407  } else {
408  nodeDesc topNd = topNode();
409  if (_debug)
410  printDebug << "using node[" << topNd << "] as top node" << endl;
411  sortedNds = sortNodesDfs(topNd);
412  }
413  if (sortedNds.size() != nmbVertices) {
414  printWarn << "number of nodes reachable from top node is "
415  << sortedNds.size() << " != " << nmbVertices << endl;
416  topologyIsOkay = false;
417  } else if (_debug)
418  printDebug << "success: all nodes are reachable from top node" << endl;
419  // make sure all interaction vertices have corresponding nodes
420  for (unsigned int i = 0; i < nmbDecayVertices(); ++i) {
421  const interactionVertexPtr& vert = decayVertices()[i];
422  if (not vert) {
423  printWarn << "interaction vertex[" << i << "] is null pointer" << endl;
424  topologyIsOkay = false;
425  continue;
426  }
427  if (not isNode(vert) or (vertex(node(vert)) != vert)) {
428  printWarn << *vert << " at " << vert << " has no associated node" << endl;
429  topologyIsOkay = false;
430  } else if (_debug)
431  printDebug << "success: " << *vert << " has associated node" << endl;
432  if (nmbOutEdges(vert) < 1) {
433  printWarn << *vert << " has " << nmbOutEdges(vert) << " < 1 "
434  << "outgoing edges" << endl;
435  topologyIsOkay = false;
436  } else if (nmbInEdges(vert) < 1) {
437  printWarn << *vert << " has " << nmbInEdges(vert) << " < 1 "
438  << "incoming edges" << endl;
439  topologyIsOkay = false;
440  } else if (_debug)
441  printDebug << "success: " << *vert
442  << " has at least 1 outgoing and 1 incoming edge" << endl;
443  }
444  // make sure all final-state particles have corresponding edges
445  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
446  const particlePtr& part = fsParticles()[i];
447  if (not part) {
448  printWarn << "final-state particle[" << i << "] is null pointer" << endl;
449  topologyIsOkay = false;
450  continue;
451  }
452  if (not isEdge(part) or (particle(this->edge(part)) != part)) {
453  printWarn << "final state " << *part << " has no associated edge" << endl;
454  topologyIsOkay = false;
455  } else if (_debug)
456  printDebug << "success: final state " << *part << " has associated edge" << endl;
457  if (isEdge(part) and not isFsVertex(toVertex(part))) {
458  printWarn << "vertex associated to final-state particle "
459  << "'" << part->name() << "' is not a final-state vertex" << endl;
460  topologyIsOkay = false;
461  } else if (_debug)
462  printDebug << "success: vertex associated to final-state particle "
463  << "'" << part->name() << "' is a final-state vertex" << endl;
464  }
465  // make sure edges connect the right vertices
466  edgeIterator iEd, iEdEnd;
467  for (tie(iEd, iEdEnd) = edges(); iEd != iEdEnd; ++iEd) {
468  const particlePtr& part = particle(*iEd);
469  // check outgoing particles of the vertex the edge is coming from
470  const interactionVertexPtr& fromVert = fromVertex(*iEd);
471  bool isInFromVert = false;
472  for (unsigned int i = 0; i < fromVert->nmbOutParticles(); ++i)
473  if (part == fromVert->outParticles()[i]) {
474  isInFromVert = true;
475  break;
476  }
477  if (not isInFromVert) {
478  printWarn << "particle '" << part->name() << "' = " << part << " associated to edge " << *iEd
479  << " is not in list of outgoing particles of " << *fromVert << " at " << fromVert
480  << ": ";
481  for (unsigned int i = 0; i < fromVert->nmbOutParticles(); ++i)
482  cout << "[" << i << "] = " << fromVert->outParticles()[i] << " ";
483  cout << endl;
484  topologyIsOkay = false;
485  } else if (_debug)
486  printDebug << "success: found particle '" << part->name() << "' "
487  << "in list of outgoing particles of vertex " << *fromVert << endl;
488  // check outgoing particles of the vertex the edge is coming from
489  const interactionVertexPtr& toVert = toVertex (*iEd);
490  bool isInToVert = false;
491  for (unsigned int i = 0; i < toVert->nmbInParticles(); ++i)
492  if (part == toVert->inParticles()[i]) {
493  isInToVert = true;
494  break;
495  }
496  if (not isInToVert) {
497  printWarn << "particle '" << part->name() << "' = " << part << " associated to edge " << *iEd
498  << " is not in list of incoming particles of " << *toVert << " at " << toVert
499  << ": ";
500  for (unsigned int i = 0; i < toVert->nmbOutParticles(); ++i)
501  cout << "[" << i << "] = " << toVert->outParticles()[i] << " ";
502  cout << endl;
503  topologyIsOkay = false;
504  } else if (_debug)
505  printDebug << "success: found particle '" << part->name() << "' "
506  << "in list of incoming particles of vertex " << *toVert << endl;
507  }
508  // make sure final-state indices make sense
509  // two cases are allowed:
510  // i) all final-state particles have indices at default value -1
511  // ii) all final-state particles have an index != -1; the indices
512  // range from 0 to (number of final-state particles - 1)
513  // sort final-state particles w.r.t. name
514  map<string, vector<particlePtr> > fsPartSorted;
515  bool allIndicesAtDefault = true;
516  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
517  const particlePtr& part = fsParticles()[i];
518  fsPartSorted[part->name()].push_back(part);
519  if (part->index() >= 0)
520  allIndicesAtDefault = false;
521  }
522  if (allIndicesAtDefault) {
523  if (_debug)
524  printDebug << "success: all final-state particles have default indices" << endl;
525  } else
526  for (map<string, vector<particlePtr> >::iterator i = fsPartSorted.begin();
527  i != fsPartSorted.end(); ++i) {
528  // sort particles with same name by their index
529  vector<particlePtr>& parts = i->second;
530  sort(parts.begin(), parts.end(), rpwa::compareIndicesAsc);
531  // check that indices are consecutive and include all unmbers
532  // from 0 to (number of final-state particles - 1)
533  for (unsigned int j = 0; j < parts.size(); ++j)
534  if (parts[j]->index() != (int)j) {
535  printWarn << "indices for final-state particles '" << i->first << "' are not consecutive. "
536  << "expected index [" << j << "], found [" << parts[j]->index() << "]." << endl;
537  topologyIsOkay = false;
538  } else if (_debug)
539  printDebug << "success: final-state particle '" << i->first << "' "
540  << "has expected index [" << j << "]" << endl;
541  }
542  if (_debug)
543  printDebug << "decay topology " << ((topologyIsOkay) ? "passed" : "did not pass")
544  << " all tests" << endl;
545  return topologyIsOkay;
546 }
547 
548 
549 decayTopology
550 decayTopology::subDecay(const nodeDesc& startNd,
551  const bool linkToMotherTopo)
552 {
553  decayTopology subTopo(dfsSubGraph(startNd, linkToMotherTopo));
554  subTopo.name() = "subdecay";
555  return subTopo;
556 }
557 
558 
559 void
560 decayTopology::addDecay(const decayTopology& topo)
561 {
562  decayTopologyGraphType::addGraph(topo);
563  buildInternalData();
564 }
565 
566 
567 void decayTopology::setProductionVertex(const productionVertexPtr& productionVertex)
568 {
569  if (not productionVertex) {
570  printErr << "null pointer for production vertex. aborting." << endl;
571  throw;
572  }
573  if (not productionVertex->XParticle()) {
574  printErr << "null pointer for particle[0] coming out of production vertex. aborting." << endl;
575  throw;
576  }
577  name() = "\"" + productionVertex->XParticle()->qnSummary() + "\"";
578  if (not _prodVertex) {
579  // topology does not have production vertex -> create graph node
580  addVertex(productionVertex);
581  } else {
582  // delete existing production vertex -> update graph node
583  vertex(node(_prodVertex)) = productionVertex;
584  _prodVertex.reset();
585  }
586  _prodVertex = productionVertex;
587 }
588 
589 
590 bool
591 decayTopology::initKinematicsData(const TClonesArray& prodKinPartNames,
592  const TClonesArray& decayKinPartNames)
593 {
594  // two modes are supported:
595  // i) all final-state particles have indices at default value -1: in
596  // case there are more than one final-state particles of the same
597  // type, the momenta are assigned in the same order the particles
598  // are stored in the final-state particle array
599  // ii) all final-state particles have an index != -1: the data are
600  // assigned according to the indices of the particles
601  // in case only part of the final-state particles have non-default
602  // indices, the result is undefined
603 
604  // init production kinematics
605  bool success = productionVertex()->initKinematicsData(prodKinPartNames);
606 
607  // check decay kinematics data
608  const string partClassName = decayKinPartNames.GetClass()->GetName();
609  if (partClassName != "TObjString") {
610  printWarn << "decay kinematics particle names are of type '" << partClassName
611  << "' and not TObjString. cannot read decay kinematics." << endl;
612  success = false;
613  }
614  const int nmbFsPart = decayKinPartNames.GetEntriesFast();
615  if ((nmbFsPart < 0) or ((unsigned int)nmbFsPart != nmbFsParticles())) {
616  printWarn << "array of decay kinematics particle names has wrong size: "
617  << nmbFsPart << " (expected " << nmbFsParticles() << ")" << endl;
618  success = false;
619  }
620  if (not success)
621  return false;
622  // sort decay kinematics particles w.r.t. particle name
623  map<string, vector<unsigned int> > fsDataPartIndices;
624  map<string, unsigned int> fsDataPartCount;
625  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
626  const string name = ((TObjString*)decayKinPartNames[i])->GetString().Data();
627  fsDataPartIndices[name].push_back(i);
628  fsDataPartCount[name] = 0;
629  }
630 
631  // create index map: index of final-state particle in decay graph -> index in data array
632  _fsDataPartIndexMap.clear();
633  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
634  const particlePtr& part = fsParticles()[i];
635  const string partName = part->name();
636  map<string, vector<unsigned int> >::const_iterator entry = fsDataPartIndices.find(partName);
637  if (entry != fsDataPartIndices.end()) {
638  int partIndex = part->index();
639  if (partIndex < 0) {
640  partIndex = fsDataPartCount[partName];
641  ++fsDataPartCount[partName];
642  }
643  if ((unsigned int)partIndex < entry->second.size()) {
644  if (_debug)
645  printDebug << "assigning decay kinematics data index [" << entry->second[partIndex] << "] "
646  << "to final-state particle '" << partName << "'[" << partIndex << "] "
647  << "at index [" << i << "]" << endl;
648  _fsDataPartIndexMap[i] = entry->second[partIndex];
649  } else {
650  printWarn << "index [" << partIndex << "] for final-state particle "
651  << "'" << partName << "' out of range. input data contain only "
652  << entry->second.size() << " entries for this particle type." << endl;
653  success = false;
654  }
655  } else {
656  printWarn << "cannot find entry for final-state particle '" << part->name() << "' "
657  << "in input data." << endl;
658  success = false;
659  }
660  }
661  if (_fsDataPartIndexMap.size() != nmbFsParticles()) {
662  printWarn << "could not find all final-state particles in input data." << endl;
663  success = false;
664  }
665 
666  return success;
667 }
668 
669 
670 bool
671 decayTopology::readKinematicsData(const TClonesArray& prodKinMomenta,
672  const TClonesArray& decayKinMomenta)
673 {
674  // set production kinematics
675  bool success = productionVertex()->readKinematicsData(prodKinMomenta);
676 
677  // check momentum array
678  const int nmbFsPart = decayKinMomenta.GetEntriesFast();
679  if ((nmbFsPart < 0) or ((unsigned int)nmbFsPart != nmbFsParticles())) {
680  printWarn << "array of decay kinematics particle momenta has wrong size: "
681  << nmbFsPart << " (expected " << nmbFsParticles() << "). "
682  << "cannot read decay kinematics." << endl;
683  success = false;
684  }
685  if (not success)
686  return false;
687 
688  // set decay kinematics
689  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
690  const particlePtr& part = fsParticles()[i];
691  const unsigned int partIndex = _fsDataPartIndexMap[i];
692  const TVector3* mom = dynamic_cast<TVector3*>(decayKinMomenta[partIndex]);
693  if (not mom) {
694  printWarn << "decay kinematics data entry [" << partIndex << "] is not of type TVector3. "
695  << "cannot read decay kinematics momentum for particle '" << part->name() << "'. "
696  << "skipping." << endl;
697  success = false;
698  continue;
699  }
700  if (_debug)
701  printDebug << "setting momentum of final-state particle '" << part->name() << "' "
702  << "at index [" << i << "] to " << *mom << " GeV "
703  << "at input data index [" << partIndex << "]" << endl;
704  part->setMomentum(*mom);
705  }
706  fillKinematicsDataCache();
707 
708  return success;
709 }
710 
711 
712 void
713 decayTopology::fillKinematicsDataCache()
714 {
715  // adjust cache size
716  _fsDataPartMomCache.resize(nmbFsParticles(), TVector3());
717  // set decay kinematics
718  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
719  const particlePtr& part = fsParticles()[i];
720  _fsDataPartMomCache[i] = part->momentum();
721  }
722 }
723 
724 
725 bool
726 decayTopology::revertMomenta()
727 {
728  // revert production kinematics
729  bool success = productionVertex()->revertMomenta();
730  // revert decay kinematics
731  if (_fsDataPartMomCache.size() != nmbFsParticles()) {
732  if (_debug)
733  printWarn << "cache for final-state particle momenta has size "
734  << _fsDataPartMomCache.size() << " != # of final-state particles = "
735  << nmbFsParticles() << endl;
736  return false;
737  }
738  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
739  const particlePtr& part = fsParticles()[i];
740  part->setMomentum(_fsDataPartMomCache[i]);
741  if (_debug)
742  printDebug << "resetting momentum of final-state particle '" << part->name() << "'"
743  << "[" << i << "] to " << _fsDataPartMomCache[i] << " GeV" << endl;
744  }
745  return success;
746 }
747 
748 
749 bool
750 decayTopology::revertMomenta(const vector<unsigned int>& fsPartPermMap) // final-state permutation map
751 {
752  // revert production kinematics
753  bool success = productionVertex()->revertMomenta();
754  // revert decay kinematics
755  if ( (_fsDataPartMomCache.size() != nmbFsParticles())
756  or (fsPartPermMap.size () != nmbFsParticles())) {
757  if (_debug)
758  printWarn << "cache size for final-state particle momenta (= "
759  << _fsDataPartMomCache.size() << ") or size of permutation map (= "
760  << fsPartPermMap.size() << ") do not match # of final-state particles (= "
761  << nmbFsParticles() << ")" << endl;
762  return false;
763  }
764  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
765  const unsigned int newIndex = fsPartPermMap[i];
766  const particlePtr& part = fsParticles()[i];
767  part->setMomentum(_fsDataPartMomCache[newIndex]);
768  if (_debug)
769  printDebug << "(re)setting momentum of final-state particle "
770  << "'" << part->name() << "'[" << i << "] "
771  << "to that of '" << fsParticles()[newIndex]->name()
772  << "'[" << newIndex << "] = " << _fsDataPartMomCache[newIndex] << " GeV" << endl;
773  }
774  return success;
775 }
776 
777 
778 ostream&
779 decayTopology::print(ostream& out) const
780 {
781  // print nodes
782  out << "decay topology '" << name() << "' has " << nmbNodes() << " node(s):" << endl;
783  if (productionVertex())
784  out << " production node[" << node(productionVertex()) << "] = "
785  << productionVertex() << ": " << *productionVertex() << endl;
786  else
787  out << " topology has no production node." << endl;
788  for (unsigned int i = 0; i < nmbDecayVertices(); ++i)
789  out << " interaction node[" << node(decayVertices()[i]) << "] = "
790  << decayVertices()[i] << ": " << *decayVertices()[i] << endl;
791  nodeIterator iNd, iNdEnd;
792  for (tie(iNd, iNdEnd) = nodes(); iNd != iNdEnd; ++iNd)
793  if (isFsVertex(vertex(*iNd)))
794  out << " final-state node[" << *iNd << "] = " << vertex(*iNd) << ": "
795  << *vertex(*iNd) << endl;
796  // print edges
797  out << "decay topology '" << name() << "' has " << nmbEdges() << " edge(s):" << endl;
798  edgeIterator iEd, iEdEnd;
799  for (tie(iEd, iEdEnd) = edges(); iEd != iEdEnd; ++iEd) {
800  const particlePtr& part = particle(*iEd);
801  out << " edge[" << *iEd << "] = [" << fromNode(*iEd) << ", "
802  << toNode(*iEd) << "] = " << part << ": '" << part->name() << "'";
803  if (isFsParticle(part))
804  out << " (final state)";
805  out << endl;
806  }
807  return out;
808 }
809 
810 
811 ostream&
812 decayTopology::printProdKinParticles(ostream& out) const
813 {
814  out << "production kinematics particles:" << endl;
815  for (unsigned int i = 0; i < productionVertex()->nmbInParticles(); ++i) {
816  const particlePtr& part = productionVertex()->inParticles()[i];
817  out << " incoming particle[" << i << "]: " << part->qnSummary() << ", "
818  << "index = " << part->index() << ", p = " << part->lzVec() << " GeV" << endl;
819  }
820  for (unsigned int i = 0; i < productionVertex()->nmbOutParticles(); ++i) {
821  const particlePtr& part = productionVertex()->outParticles()[i];
822  out << " outgoing particle[" << i << "]: " << part->qnSummary() << ", "
823  << "index = " << part->index() << ", p = " << part->lzVec() << " GeV" << endl;
824  }
825  return out;
826 }
827 
828 
829 ostream&
830 decayTopology::printDecayKinParticles(ostream& out) const
831 {
832  out << "decay kinematics particles:" << endl;
833  for (unsigned int i = 0; i < nmbFsParticles(); ++i) {
834  const particlePtr& part = fsParticles()[i];
835  out << " particle[" << i << "]: " << part->qnSummary() << ", "
836  << "index = " << part->index() << ", p = " << part->lzVec() << " GeV" << endl;
837  }
838  return out;
839 }
840 
841 
842 decayTopology&
843 decayTopology::constructDecay(const productionVertexPtr& productionVertex,
844  const vector<interactionVertexPtr>& decayVertices,
845  const vector<particlePtr>& fsParticles,
846  const bool performTopologyCheck)
847 {
848  clear();
849  const unsigned int nmbDecayVert = decayVertices.size();
850  const unsigned int nmbFsPart = fsParticles.size();
851  if (_debug)
852  printDebug << "constructing decay topology with " << nmbDecayVert << " decay vertices and "
853  << nmbFsPart << " final-state particles" << endl;
854  if (nmbFsPart < 1) {
855  printErr << "cannot construct decay topology without final-state particles. aborting." << endl;
856  throw;
857  }
858  if (nmbDecayVert < 1) {
859  printWarn << "need at least production and X-decay vertex "
860  << "to construct decay topology. aborting." << endl;
861  throw;
862  }
863 
864  // create graph node for production vertex and store pointer
865  if (not productionVertex) {
866  printErr << "null pointer for production vertex. aborting." << endl;
867  throw;
868  }
869  if (not productionVertex->XParticle()) {
870  printErr << "null pointer for X particle coming out of production vertex. aborting." << endl;
871  throw;
872  }
873  name() = "\"" + productionVertex->XParticle()->qnSummary() + "\"";
874  _prodVertex = productionVertex;
875  addVertex(productionVertex);
876  // create graph nodes for interaction vertices and store pointers
877  for (unsigned int i = 0; i < nmbDecayVert; ++i) {
878  if (not decayVertices[i]) {
879  printErr << "null pointer for decay vertex[" << i << "]. aborting." << endl;
880  throw;
881  }
882  addVertex(decayVertices[i]);
883  }
884  // create final-state nodes and vertices
885  for (unsigned int i = 0; i < nmbFsPart; ++i) {
886  if (not fsParticles[i]) {
887  printErr << "null pointer for final-state particle[" << i << "]. aborting." << endl;
888  throw;
889  }
890  const fsVertexPtr& fsVert = createFsVertex(fsParticles[i]);
891  addVertex(fsVert);
892  }
893  // memorize depth-first sorted interaction vertices and final-state particles
894  vector<nodeDesc> sortedNds = sortNodesDfs(node(_prodVertex));
895  for (unsigned int i = 0; i < sortedNds.size(); ++i) {
896  const interactionVertexPtr& vert = vertex(sortedNds[i]);
897  const fsVertexPtr& fsVert = dynamic_pointer_cast<fsVertex>(vert);
898  if (fsVert) {
899  const particlePtr& part = fsVert->inParticles()[0];
900  for (unsigned int i = 0; i < nmbFsPart; ++i)
901  if (part == fsParticles[i]) {
902  _fsParticles.push_back(part);
903  break;
904  }
905  } else if (vert != _prodVertex) {
906  for (unsigned int i = 0; i < nmbDecayVert; ++i)
907  if (vert == decayVertices[i]) {
908  _decayVertices.push_back(vert);
909  break;
910  }
911  }
912  }
913  // check that topology makes sense
914  if (performTopologyCheck and not checkTopology()) {
915  printErr << "topology has problems that need to be fixed. aborting." << endl;
916  throw;
917  }
918  return *this;
919 }
920 
921 
922 void
923 decayTopology::buildInternalData()
924 {
925  bool success = true;
926  // find production vertex
927  // assumes that production vertex is the only vertex in graph that
928  // has no incoming edges and exactly one outgoing edge
929  _prodVertex.reset();
930  unsigned int nmbProdVertCandidates = 0;
931  nodeIterator iNd, iNdEnd;
932  for (tie(iNd, iNdEnd) = nodes(); iNd != iNdEnd; ++iNd)
933  if ((nmbInEdges(*iNd) == 0) and (nmbOutEdges(*iNd) == 1)) {
934  _prodVertex = dynamic_pointer_cast<rpwa::productionVertex>(vertex(*iNd));
935  if (_prodVertex)
936  ++nmbProdVertCandidates;
937  }
938  // set final-state particles
939  _fsParticles.clear();
940  vector<nodeDesc> sortedNds;
941  if (_prodVertex)
942  sortedNds = sortNodesDfs(node(_prodVertex));
943  else {
944  // take non-FS vertex with no incoming particle as top node
945  nodeDesc topNode;
946  unsigned int nmbTopVertCandidates = 0;
947  for (tie(iNd, iNdEnd) = nodes(); iNd != iNdEnd; ++iNd)
948  if ((nmbInEdges(*iNd) == 0) and not isFsVertex(vertex(*iNd))) {
949  topNode = *iNd;
950  ++nmbTopVertCandidates;
951  }
952  if (nmbTopVertCandidates == 1)
953  sortedNds = sortNodesDfs(topNode);
954  else {
955  if (_debug)
956  printWarn << "ill-formed graph. vertex order will be undefined." << endl;
957  for (tie(iNd, iNdEnd) = nodes(); iNd != iNdEnd; ++iNd)
958  sortedNds.push_back(*iNd);
959  }
960  }
961  for (unsigned int i = 0; i < sortedNds.size(); ++i) {
962  const interactionVertexPtr& vert = vertex(sortedNds[i]);
963  if (isFsVertex(vert))
964  _fsParticles.push_back(vert->inParticles()[0]);
965  }
966  if (nmbProdVertCandidates > 1) {
967  printWarn << "found " << nmbProdVertCandidates << " instead of 0 or 1 candidate "
968  << "for production vertex in graph '" << name() << "'" << endl;
969  success = false;
970  }
971  if (_fsParticles.size() < 1) {
972  printWarn << "cannot find final-state particles in graph '" << name() << "'" << endl;
973  success = false;
974  }
975  if (not success) {
976  printErr << "cannot construct decay topology from graph '" << name() << "'. aborting." << endl;
977  throw;
978  }
979  // find interaction vertices
980  _decayVertices.clear();
981  for (unsigned int i = 0; i < sortedNds.size(); ++i) {
982  const interactionVertexPtr& vert = vertex(sortedNds[i]);
983  if (isDecayVertex(vert))
984  _decayVertices.push_back(vert);
985  }
986 }
987 
988 
989 interactionVertexPtr
990 decayTopology::cloneNode(const nodeDesc& nd,
991  const bool cloneInParticles,
992  const bool cloneOutParticles)
993 {
994  const interactionVertexPtr vert = vertex(nd); // this must not be a reference
995  const bool isProdVert = isProductionVertex(vert);
996  const bool isDecayVert = isDecayVertex(vert);
997  const interactionVertexPtr& newVert = decayTopologyGraphType::cloneNode(nd); // created new vertex for this node
998  // update member variables
999  if (isProdVert)
1000  _prodVertex = static_pointer_cast<rpwa::productionVertex>(newVert);
1001  else if (isDecayVert)
1002  for (unsigned int i = 0; i < nmbDecayVertices(); ++i)
1003  if (vert == _decayVertices[i])
1004  _decayVertices[i] = newVert;
1005  return newVert;
1006 }
1007 
1008 
1009 particlePtr
1010 decayTopology::cloneEdge(const edgeDesc& ed)
1011 {
1012  const particlePtr part = particle(ed); // this must not be a reference
1013  const particlePtr& newPart = decayTopologyGraphType::cloneEdge(ed);
1014  // update member variable
1015  if (isFsParticle(part))
1016  for (unsigned int i = 0; i < nmbFsParticles(); ++i)
1017  if (part == fsParticles()[i])
1018  _fsParticles[i] = newPart;
1019  return newPart;
1020 }
Generated by   doxygen 1.8.1.2