ROOTPWA
testDecayTopology.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 // basic test program for vertex and decay topology
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <fstream>
39 #include <sstream>
40 
41 #include "TVector3.h"
42 #include "TSystem.h"
43 
44 #include "../particleDataTable.h"
45 #include "../particle.h"
46 #include "decayGraph.hpp"
47 #include "decayTopology.h"
48 #include "diffractiveDissVertex.h"
49 #include "isobarDecayVertex.h"
50 #include "isobarDecayTopology.h"
51 #include "waveDescription.h"
52 #include "isobarHelicityAmplitude.h"
53 
54 
55 using namespace std;
56 using namespace rpwa;
57 using namespace boost;
58 
59 
60 typedef decayGraph<interactionVertex, particle> graphType;
61 
62 
63 int
64 main(int argc, char** argv)
65 {
66  printSvnVersion();
67 
68  // switch on debug output
69  // particle::setDebug(true);
70  // graphType::setDebug(true);
71  // decayTopologyGraphType::setDebug(true);
72  // decayTopology::setDebug(true);
73  // interactionVertex::setDebug(true);
74  // diffractiveDissVertex::setDebug(true);
75  // isobarDecayVertex::setDebug(true);
76  // isobarDecayTopology::setDebug(true);
77 
78  particleDataTable& pdt = particleDataTable::instance();
79  pdt.readFile();
80 
81  // test construction of vertices
82  if (0) {
83  TVector3 mom;
84  mom = TVector3(1, 2, 3);
85  particlePtr beam = createParticle("pi-");
86  beam->setMomentum(mom);
87  particlePtr target = createParticle("p+");
88  particlePtr X = createParticle("X-");
89  printInfo << "created particles: " << endl
90  << *beam << endl
91  << *target << endl
92  << *X << endl;
93  diffractiveDissVertexPtr vert1 = createDiffractiveDissVertex(beam, target, X);
94  printInfo << "created vertex: " << endl
95  << *vert1 << endl;
96  diffractiveDissVertexPtr vert2(vert1);
97  printInfo << "copied vertex: " << endl
98  << *vert2 << endl;
99 
100  mom = TVector3(3, 4, 5);
101  particlePtr daughter1 = createParticle("pi-");
102  daughter1->setMomentum(mom);
103  mom = TVector3(4, 5, 6);
104  particlePtr daughter2 = createParticle("pi0");
105  daughter2->setMomentum(mom);
106  isobarDecayVertexPtr vert3 = createIsobarDecayVertex(X, daughter1, daughter2, 1, 2);
107  printInfo << "created vertex: " << endl
108  << *vert3 << endl;
109  isobarDecayVertexPtr vert4(vert3);
110  printInfo << "copied vertex: " << endl
111  << *vert4 << endl;
112  }
113 
114  if (0) {
115  particlePtr pi0 = createParticle("pi-");
116  particlePtr pi1 = createParticle("pi+");
117  particlePtr pi2 = createParticle("pi-");
118  particlePtr pi3 = createParticle("pi+");
119  particlePtr pi4 = createParticle("pi-");
120  // define isobars
121  particlePtr sigma = createParticle("sigma");
122  particlePtr a1 = createParticle("a1(1260)+");
123  particlePtr f1 = createParticle("f1(1285)");
124  // define X-system
125  // 2I G 2J P C 2M
126  particlePtr X = createParticle("X+", 2, +1, 4, +1, -1, 2);
127  // define production vertex
128  particlePtr beam = createParticle("pi-");
129  particlePtr target = createParticle("p+");
130  diffractiveDissVertexPtr prodVert = createDiffractiveDissVertex(beam, target, X);
131  // define vertices
132  isobarDecayVertexPtr vert0 = createIsobarDecayVertex(X, pi4, f1, 0, 3);
133  isobarDecayVertexPtr vert1 = createIsobarDecayVertex(f1, pi2, a1, 2, 2);
134  isobarDecayVertexPtr vert2 = createIsobarDecayVertex(a1, pi3, sigma, 2, 0);
135  isobarDecayVertexPtr vert3 = createIsobarDecayVertex(sigma, pi0, pi1, 0, 0);
136 
137  // construct graph
138  graphType g;
139  g.name() = "test graph";
140  g.addVertex(vert0);
141  g.addVertex(vert1);
142  g.addVertex(vert2);
143  g.addVertex(vert3);
144  printInfo << g;
145 
146  for (graphType::nodeIterator i = g.nodes().first; i != g.nodes().second; ++i) {
147  const isobarDecayVertexPtr v = static_pointer_cast<isobarDecayVertex>(g.vertex(*i));
148  g.name(v) = v->parent()->name();
149  }
150  for (graphType::edgeIterator i = g.edges().first; i != g.edges().second; ++i) {
151  const particlePtr& p = g.particle(*i);
152  stringstream n;
153  n << "edge [" << g.index(*i) << "] " << p->name();
154  g.name(*i) = n.str();
155  }
156  g.print(cout, g.nodeNameMap(), g.edgeNameMap());
157  {
158  ofstream graphVizFile("decay.dot");
159  g.writeGraphViz(graphVizFile, g.nodeNameMap(), g.edgeNameMap());
160  gSystem->Exec("dot -Tps -o decay.ps decay.dot");
161  }
162 
163  graphType g2 = g;
164  g2.name() = "graph copy";
165  for (graphType::adjIterator i = g2.adjacentVertices(vert1).first;
166  i != g2.adjacentVertices(vert1).second; ++i) {
167  g2.name (*i) += " !bar!";
168  g2.color(*i) = white_color;
169  printInfo << "vert1 adjacent vertex[" << *i << "]: " << *g2[*i] << endl;
170  }
171  printInfo << "nmbInParticles(vert1): " << g2.nmbInEdges(vert1) << endl;
172  for (graphType::inEdgeIterator i = g2.incomingEdges(vert1).first;
173  i != g2.incomingEdges(vert1).second; ++i) {
174  g2.name (*i) += " !blah!";
175  g2.color(*i) = gray_color;
176  printInfo << "vert1 in edge[" << *i << "]: " << *g2[*i] << endl;
177  }
178  printInfo << "nmbOutParticles(vert1): " << g2.nmbOutEdges(vert1) << endl;
179  for (graphType::outEdgeIterator i = g2.outgoingEdges(vert1).first;
180  i != g2.outgoingEdges(vert1).second; ++i) {
181  g2.name (*i) += " !blah2!";
182  g2.color(*i) = black_color;
183  printInfo << "vert1 out edge[" << *i << "]: " << *g2[*i] << endl;
184  }
185  isobarDecayVertexPtr dummyVert = createIsobarDecayVertex(X, pi4, f1, 0, 3);
186  printInfo << "isNode(vert0) = " << g2.isNode(vert0) << ", "
187  << "isNode(dummyVert) = " << g2.isNode(dummyVert) << endl;
188  particlePtr dummyPart = createParticle("X+", 1, +1, 4, +1, -1, 2);
189  printInfo << "isEdge(f1) = " << g2.isEdge(f1) << ", "
190  << "isEdge(dummyPart) = " << g2.isEdge(dummyPart) << endl;
191  cout << "fromVertex(f1): " << *g2.fromVertex(f1) << endl;
192  cout << "toVertex(f1): " << *g2.toVertex (f1) << endl;
193  cout << "particleExists(vert0, vert1) = " << g2.particleExists(vert0, vert1) << ", "
194  << "particleExists(vert1, vert0) = " << g2.particleExists(vert1, vert0) << endl;
195  cout << "particleConnects(f1, vert0, vert1) = " << g2.particleConnects(f1, vert0, vert1) << endl;
196  cout << "particleConnects(a1, vert0, vert1) = " << g2.particleConnects(a1, vert0, vert1) << endl;
197  printInfo << endl;
198  for (graphType::nodeIterator i = g2.nodes().first; i != g2.nodes().second; ++i)
199  cout << "node " << *i << ": name = '" << g2.name(*i) << "', "
200  << "index = " << g2.index(*i) << ", "
201  << "color = " << g2.color(*i) << ", "
202  << endl;
203  for (graphType::edgeIterator i = g2.edges().first; i != g2.edges().second; ++i)
204  cout << "edge " << *i << ": name = '" << g2.name(*i) << "', "
205  << "index = " << g2.index(*i) << ", "
206  << "color = " << g2.color(*i) << ", "
207  << endl;
208  g2.print(cout, g2.nodeNameMap());
209  g2.clear();
210  printInfo << "after clear:" << endl << g2;
211 
212  decayTopologyGraphType g3;
213  g3.name() = "test graph 3";
214  g3.addVertex(vert0);
215  g3.addVertex(vert1);
216  g3.addVertex(vert2);
217  g3.addVertex(vert3);
218  g3.addVertex(prodVert);
219  g3.addVertex(createFsVertex(pi0));
220  g3.addVertex(createFsVertex(pi1));
221  g3.addVertex(createFsVertex(pi2));
222  g3.addVertex(createFsVertex(pi3));
223  g3.addVertex(createFsVertex(pi4));
224  printInfo << g3;
225  decayTopology topo(g3);
226  topo.name() = "topo graph copy";
227  printInfo << topo;
228 
229  vector<isobarDecayVertexPtr> decayVertices;
230  decayVertices.push_back(vert3);
231  decayVertices.push_back(vert1);
232  decayVertices.push_back(vert2);
233  decayVertices.push_back(vert0);
234  vector<particlePtr> fsParticles;
235  fsParticles.push_back(pi0);
236  fsParticles.push_back(pi1);
237  fsParticles.push_back(pi2);
238  fsParticles.push_back(pi3);
239  fsParticles.push_back(pi4);
240  isobarDecayTopology topo2(prodVert, decayVertices, fsParticles);
241  topo2.name() = "topo";
242  printInfo << topo2;
243  const bool topologyOkay = topo2.checkTopology();
244  printInfo << "topology okay = " << topologyOkay << endl;
245  printInfo << "performing consistency checks on isobar decay topology" << endl;
246  const bool decayConsistent = topo2.checkConsistency();
247  printInfo << "decay consistent = " << decayConsistent << endl;
248  topo2.calcIsobarLzVec();
249 
250  {
251  cout << "-------------------------------------------------------------------------------"
252  << endl;
253  printInfo << "testing cloning" << endl;
254  isobarDecayTopologyPtr topo3 = topo2.clone();
255  isobarDecayVertexPtr v = topo3->isobarDecayVertices()[0];
256  particlePtr p = v->inParticles()[0];
257  v->setL(2);
258  v->setS(2);
259  cout << *p << endl;
260  p->setG(-1);
261  p->setCharge(-1);
262  p->setC(1);
263  //topo3->checkTopology();
264  //topo3->checkConsistency();
265  cout << *topo3;
266  decayTopologyGraphType::nodeIterator iNode, iNodeEnd;
267  cout << "!!! topo2 vertex pointers: ";
268  for (tie(iNode, iNodeEnd) = topo2.nodes(); iNode != iNodeEnd; ++iNode)
269  cout << topo2.vertex(*iNode) << " ";
270  cout << endl;
271  cout << "!!! topo3 vertex pointers: ";
272  for (tie(iNode, iNodeEnd) = topo3->nodes(); iNode != iNodeEnd; ++iNode)
273  cout << topo3->vertex(*iNode) << " ";
274  cout << endl;
275  cout << "!!! topo3 interaction vertex pointers: ";
276  for (unsigned int i = 0; i < topo3->nmbDecayVertices(); ++i)
277  cout << topo3->decayVertices()[i] << " ";
278  cout << endl;
279  cout << "!!! topo3 isobar decay vertex pointers: ";
280  for (unsigned int i = 0; i < topo3->nmbDecayVertices(); ++i)
281  cout << topo3->isobarDecayVertices()[i] << " ";
282  cout << endl;
283  decayTopologyGraphType::edgeIterator iEdge, iEdgeEnd;
284  cout << "!!! topo2 particle pointers: ";
285  for (tie(iEdge, iEdgeEnd) = topo2.edges(); iEdge != iEdgeEnd; ++iEdge)
286  cout << topo2.particle(*iEdge) << " ";
287  cout << endl;
288  cout << "!!! topo3 particle pointers: ";
289  for (tie(iEdge, iEdgeEnd) = topo3->edges(); iEdge != iEdgeEnd; ++iEdge)
290  cout << topo3->particle(*iEdge) << " ";
291  cout << endl;
292  cout << "!!! topo3 FS particle pointers: ";
293  for (unsigned int i = 0; i < topo3->nmbFsParticles(); ++i)
294  cout << topo3->fsParticles()[i] << " ";
295  cout << endl;
296  }
297 
298  {
299  cout << "-------------------------------------------------------------------------------"
300  << endl;
301  printInfo << "testing subdecays and merge" << endl;
302  isobarDecayTopology subDecay1 = topo2.subDecay(2);
303  subDecay1.decayTopologyGraphType::printPointers(cout);
304  subDecay1.decayTopology::print(cout);
305  subDecay1.checkTopology();
306  cout << "subdecay from node[2]: " << subDecay1;
307  const isobarDecayTopology& subDecay2 = topo2.subDecay(1);
308  cout << "subdecay from node[1]: " << subDecay2;
309  const isobarDecayTopology& subDecay3 = topo2.subDecay(9);
310  cout << "subdecay from node[9]: " << subDecay3;
311  subDecay1.addDecay(subDecay3);
312  cout << "subdecay from node[2] + subdecay from node[9]: " << subDecay1;
313  isobarDecayTopology decay;
314  // 2I G 2J P C 2M
315  particlePtr X2 = createParticle("X-", 2, -1, 4, +1, 0, 2);
316  isobarDecayVertexPtr vertX= createIsobarDecayVertex(X2, pi4, f1, 2, 2);
317  decay.addVertex(vertX);
318  subDecay1.addDecay(decay);
319  cout << "subdecay from node[2] + subdecay from node[9] + vert0: " << subDecay1;
320  const isobarDecayTopology& subDecayX = topo2.subDecay(2);
321  isobarDecayTopology newDecay = isobarDecayTopology::joinDaughterDecays(vertX, subDecayX, subDecay3);
322  diffractiveDissVertexPtr prodVertX = createDiffractiveDissVertex(beam, target, X2);
323  newDecay.setProductionVertex(prodVertX);
324  cout << "joined graph: " << newDecay;
325  newDecay.checkTopology();
326  newDecay.checkConsistency();
327  }
328  }
329 
330  if (1) {
331  const string keyFileName = "testFindIsobarBoseSymVertices.key";
332 
333  waveDescription waveDesc;
334  isobarAmplitudePtr amp;
335  if (not waveDesc.parseKeyFile(keyFileName) or not waveDesc.constructAmplitude(amp)) {
336  printErr << "problems constructing amplitude. exiting." << endl;
337  exit(1);
338  }
339  isobarDecayTopologyPtr topo = amp->decayTopology();
340  printInfo << *topo;
341 
342  topo->findIsobarBoseSymVertices();
343  }
344 }
Generated by   doxygen 1.8.1.2