ROOTPWA
waveDescription.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 reads/writes wave description from/to keyfiles,
29 // constructs decay topologies and amplitudes
30 // can be saved into .root files
31 //
32 //
33 // Author List:
34 // Boris Grube TUM (original author)
35 //
36 //
37 //-------------------------------------------------------------------------
38 
39 
40 #include <fstream>
41 #include <map>
42 #include <cstdio>
43 
44 #include <boost/algorithm/string.hpp>
45 #include <boost/tokenizer.hpp>
46 #include <boost/assign.hpp>
47 
48 #include "TClass.h"
49 
50 #include "libConfigUtils.hpp"
51 #include "conversionUtils.hpp"
52 #include "diffractiveDissVertex.h"
53 #include "leptoProductionVertex.h"
54 #include "isobarHelicityAmplitude.h"
55 #include "isobarCanonicalAmplitude.h"
56 #include "waveDescription.h"
57 
58 
59 using namespace std;
60 using namespace boost;
61 using namespace boost::assign;
62 using namespace libconfig;
63 using namespace rpwa;
64 
65 
66 ClassImp(waveDescription);
67 
68 
69 bool waveDescription::_debug = false;
70 
71 
72 map<string,string> waveDescription::isobars = map_list_of
73  ("pi+", "\\pi^+")
74  ("pi-", "\\pi^-")
75  ("pi+-", "\\pi^\\pm")
76  ("pi-+", "\\pi^\\mp")
77  ("sigma0", "\\sigma")
78  ("rho(770)0", "\\rho^0(770)")
79  ("a1(1260)-", "a_1^-(1260)")
80  ("a2(1320)-", "a_2^-(1320)")
81  ("rho(1450)0", "\\rho^0(1450)")
82  ("rho(1700)0", "\\rho^0(1700)")
83  ("pi(1300)-", "\\pi^-(1300)")
84  ("pi(1800)-", "\\pi^-(1800)")
85  ("pi2(1670)-", "\\pi^-_2(1670)")
86  ("f0(1370)0", "f_0^0(1370)")
87  ("f0(1500)0", "f_0^0(1500)")
88  ("f0(1700)0", "f_0^0(1700)")
89  ("f1(1285)0", "f_1^0(1285)")
90  ("f1(1420)0", "f_1^0(1420)")
91  ("b1(1235)0", "b_1^0(1235)")
92  ("b1(1800)0", "b_1^0(1800)")
93  ("b0(1800)0", "b_0^0(1800)")
94  ("b2(1800)0", "b_2^0(1800)")
95  ("b1(1500)0", "b_1^0(1500)")
96  ("f2(1270)0", "f_2^0(1270)")
97  ("f2(1950)0", "f_2^0(1950)")
98  ("f2(1565)0", "f_2^0(1565)")
99  ("f2(2010)0", "f_2^0(2010)")
100  ("eta(1440)0" , "\\eta^0(1420)")
101  ("eta2(1645)0", "\\eta_2^0(1645)")
102  ("eta1(1600)0", "\\eta_1^0(1600)")
103  ("pi1(1600)-", "\\pi_1^-(1600)")
104  ("rho3(1690)0", "\\rho_3^0(1690)")
105  ("rho(1600)0", "\\rho^0(1600)");
106 
107 
108 waveDescription::waveDescription()
109  : TObject (),
110  _key (0),
111  _keyFileParsed (false),
112  _keyFileLocalCopy("")
113 {
114  //waveDescription::Class()->IgnoreTObjectStreamer(); // don't store TObject's fBits and fUniqueID
115 }
116 
117 
118 waveDescription::~waveDescription()
119 {
120  if (_key)
121  delete _key;
122 }
123 
124 
125 void
126 waveDescription::clear()
127 {
128  _key = 0;
129  _keyFileParsed = false;
130  _keyFileLocalCopy = "";
131 }
132 
133 
134 waveDescription&
135 waveDescription::operator =(const waveDescription& waveDesc)
136 {
137  if (this != &waveDesc) {
138  TObject::operator =(waveDesc);
139  _key = waveDesc._key;
140  _keyFileParsed = waveDesc._keyFileParsed;
141  _keyFileLocalCopy = waveDesc._keyFileLocalCopy;
142  }
143  return *this;
144 }
145 
146 
147 bool
148 waveDescription::parseKeyFile(const string& keyFileName)
149 {
150  _keyFileParsed = false;
151  if (not _key)
152  _key = new Config();
153  if (not parseLibConfigFile(keyFileName, *_key, _debug)) {
154  printWarn << "problems reading key file '" << keyFileName << "'. "
155  << "cannot construct decay topology." << endl;
156  delete _key;
157  _key = 0;
158  return false;
159  }
160  // read key file contents into string
161  {
162  ifstream keyFile(keyFileName.c_str());
163  if (not keyFile or not keyFile.good()) {
164  printWarn << "cannot read from file '" << keyFileName << "'" << endl;
165  return false;
166  }
167  _keyFileLocalCopy = "";
168  string line;
169  while(getline(keyFile, line))
170  _keyFileLocalCopy += line + "\n";
171  }
172  _keyFileParsed = true;
173  return true;
174 }
175 
176 
177 ostream&
178 waveDescription::printKeyFileContents(ostream& out) const
179 {
180  if (_keyFileLocalCopy != "") {
181  typedef tokenizer<char_separator<char> > tokenizer;
182  char_separator<char> separator("\n");
183  tokenizer keyFileLines(_keyFileLocalCopy, separator);
184  unsigned int lineNumber = 0;
185  for (tokenizer::iterator i = keyFileLines.begin(); i != keyFileLines.end(); ++i)
186  out << setw(5) << ++lineNumber << " " << *i << endl;
187  } else
188  out << "key file contents string is empty" << endl;
189  return out;
190 }
191 
192 
193 bool
194 waveDescription::constructDecayTopology(isobarDecayTopologyPtr& topo,
195  const bool fromTemplate) const
196 {
197  if (not _key or not _keyFileParsed) {
198  printWarn << "parsing was not successful. cannot construct decay topology." << endl;
199  return false;
200  }
201 
202  if (topo)
203  topo.reset();
204  topo = isobarDecayTopologyPtr(); // null pointer
205 
206  const Setting& rootKey = _key->getRoot();
207  printInfo << "constructing decay topology from key file" << endl;
208 
209  // find wave group
210  const Setting* decayVertKey = findLibConfigGroup(rootKey, "decayVertex");
211  if (not decayVertKey) {
212  printWarn << "cannot find 'wave' group. cannot construct decay topology." << endl;
213  return false;
214  }
215 
216  // find X quantum numbers group
217  const Setting* XQnKey = findLibConfigGroup(*decayVertKey, "XQuantumNumbers", not fromTemplate);
218  particlePtr X;
219  if (not XQnKey)
220  if (not fromTemplate) {
221  printWarn << "cannot find 'XQuantumNumbers' group. cannot construct decay topology." << endl;
222  return false;
223  } else {
224  // create default X particle
225  if (_debug)
226  printDebug << "cannot find 'XQuantumNumbers' group. "
227  << "creating X particle with default properties." << endl;
228  X = createParticle("X", not fromTemplate);
229  }
230  else
231  // create X particle from key
232  if (not constructXParticle(*XQnKey, X)) {
233  printWarn << "problems constructing X particle. cannot construct decay topology." << endl;
234  return false;
235  }
236 
237  // create production vertex
238  productionVertexPtr prodVert = productionVertexPtr();
239  if (not constructProductionVertex(rootKey, X, prodVert)) {
240  printWarn << "problems constructing production vertex. cannot construct decay topology." << endl;
241  return false;
242  }
243 
244  // find X decay group
245  const Setting* XDecayKey = findLibConfigGroup(*decayVertKey, "XDecay");
246  if (not XDecayKey) {
247  printWarn << "cannot find 'XDecay' group. cannot construct decay topology." << endl;
248  return false;
249  }
250  // traverse decay chain and create final state particles and isobar decay vertices
251  vector<isobarDecayVertexPtr> decayVertices;
252  vector<particlePtr> fsParticles;
253  if (not constructDecayVertex(*XDecayKey, X, decayVertices, fsParticles, fromTemplate)) {
254  printWarn << "problems constructing decay chain. cannot construct decay topology." << endl;
255  return false;
256  }
257 
258  // create decay isobar decay topology
259  topo = createIsobarDecayTopology(prodVert, decayVertices, fsParticles);
260  // backward compatibility: allow sloppy key files, where charges of
261  // isobars are not explicitely defined
262  topo->calcIsobarCharges();
264  //topo->calcIsobarBaryonNmbs();
265  //topo->productionVertex()->setXFlavorQN(); // sets baryon nmb, S, C, and B of X
266 
267  printSucc << "constructed decay topology from key file" << endl;
268  return true;
269 }
270 
271 
272 bool
273 waveDescription::constructAmplitude(isobarAmplitudePtr& amplitude) const
274 {
275  isobarDecayTopologyPtr topo;
276  if (not constructDecayTopology(topo)) {
277  printWarn << "problems constructing decay topology. cannot construct decay amplitude." << endl;
278  return false;
279  }
280  return constructAmplitude(amplitude, topo);
281 }
282 
283 
284 bool
285 waveDescription::constructAmplitude(isobarAmplitudePtr& amplitude,
286  const isobarDecayTopologyPtr& topo) const
287 {
288  if (not topo) {
289  printWarn << "null pointer to decay topology. cannot construct decay amplitude." << endl;
290  return false;
291  }
292  if (not topo->checkTopology() or not topo->checkConsistency()) {
293  printWarn << "decay topology has issues. cannot construct decay amplitude." << endl;
294  return false;
295  }
296  if (amplitude)
297  amplitude.reset();
298  // get amplitude parameters from key file
299  // default values
300  string formalism = "helicity";
301  bool boseSymmetrize = true;
302  bool isospinSymmetrize = true;
303  bool useReflectivityBasis = true;
304  // find amplitude group
305  const Setting& rootKey = _key->getRoot();
306  const libconfig::Setting* amplitudeKey = findLibConfigGroup(rootKey, "amplitude", false);
307  if (amplitudeKey) {
308  if (amplitudeKey->lookupValue("formalism", formalism) and _debug)
309  printDebug << "setting amplitude formalism to '" << formalism << "'" << endl;
310  if (amplitudeKey->lookupValue("boseSymmetrize", boseSymmetrize) and _debug)
311  printDebug << "setting amplitude option 'boseSymmetrize' to "
312  << trueFalse(boseSymmetrize) << endl;
313  if (amplitudeKey->lookupValue("isospinSymmetrize", isospinSymmetrize) and _debug)
314  printDebug << "setting amplitude option 'isospinSymmetrize' to "
315  << trueFalse(isospinSymmetrize) << endl;
316  if (amplitudeKey->lookupValue("useReflectivityBasis", useReflectivityBasis) and _debug)
317  printDebug << "setting amplitude option 'useReflectivityBasis' to "
318  << trueFalse(useReflectivityBasis) << endl;
319  }
320  // construct amplitude
321  amplitude = mapAmplitudeType(formalism, topo);
322  if (_debug)
323  printDebug << "constructed amplitude '"<< amplitude->name() << "': "
324  << enDisabled(boseSymmetrize ) << " Bose symmetrization, "
325  << enDisabled(isospinSymmetrize ) << " isospin symmetrization, "
326  << enDisabled(useReflectivityBasis) << " reflectivity basis" << endl;
327  if (not amplitude) {
328  printWarn << "problems constructing decay amplitude." << endl;
329  return false;
330  }
331  amplitude->enableBoseSymmetrization (boseSymmetrize );
332  amplitude->enableIsospinSymmetrization(isospinSymmetrize );
333  amplitude->enableReflectivityBasis (useReflectivityBasis);
334  return true;
335 }
336 
337 
338 string
339 waveDescription::waveNameFromTopology(isobarDecayTopology topo,
340  const bool newConvention,
341  const isobarDecayVertexPtr& currentVertex)
342 {
343  ostringstream waveName;
344  if (currentVertex == interactionVertexPtr()) {
345  if (not topo.checkTopology() or not topo.checkConsistency()) {
346  printWarn << "decay topology has issues. cannot construct wave name." << endl;
347  return "";
348  }
349  // X quantum numbers
350  const particle& X = *(topo.XParticle());
351  if (newConvention)
352  waveName << "[" << spinQn(X.isospin()) << parityQn(X.G()) << ","
353  << spinQn(X.J()) << parityQn(X.P()) << parityQn(X.C()) << ","
354  << spinQn(X.spinProj()) << parityQn(X.reflectivity()) << "]"
355  << waveNameFromTopology(topo, newConvention, topo.XIsobarDecayVertex());
356  else
357  waveName << spinQn(X.isospin()) << sign(X.G())
358  << spinQn(X.J()) << sign(X.P()) << sign(X.C())
359  << spinQn(X.spinProj()) << sign(X.reflectivity())
360  << waveNameFromTopology(topo, newConvention, static_pointer_cast<isobarDecayVertex>
361  (topo.toVertex(topo.XIsobarDecayVertex()->daughter1())))
362  << "_" << spinQn(topo.XIsobarDecayVertex()->L())
363  << spinQn(topo.XIsobarDecayVertex()->S()) << "_"
364  << waveNameFromTopology(topo, newConvention, static_pointer_cast<isobarDecayVertex>
365  (topo.toVertex(topo.XIsobarDecayVertex()->daughter2())));
366  } else {
367  // recurse down decay chain
368  if (newConvention) {
369  // first daughter
370  waveName << "=[" << currentVertex->daughter1()->name();
371  if (not topo.isFsParticle(currentVertex->daughter1()))
372  waveName << waveNameFromTopology
373  (topo, newConvention,
374  static_pointer_cast<isobarDecayVertex>(topo.toVertex(currentVertex->daughter1())));
375  // L, S
376  waveName << "[" << spinQn(currentVertex->L()) << "," << spinQn(currentVertex->S()) << "]";
377  // second daughter
378  waveName << currentVertex->daughter2()->name();
379  if (not topo.isFsParticle(currentVertex->daughter2()))
380  waveName << waveNameFromTopology
381  (topo, newConvention,
382  static_pointer_cast<isobarDecayVertex>(topo.toVertex(currentVertex->daughter2())));
383  waveName << "]";
384  } else {
385  waveName << ((currentVertex->parent()->charge() != 0) ? currentVertex->parent()->name()
386  : currentVertex->parent()->bareName());
387  isobarDecayTopology subGraph = topo.subDecay(topo.node(currentVertex));
388  if (not topo.isFsVertex(currentVertex) and subGraph.nmbFsParticles() > 2)
389  waveName << "="
390  << waveNameFromTopology(topo, newConvention, static_pointer_cast<isobarDecayVertex>
391  (topo.toVertex(currentVertex->daughter1())))
392  << "_" << spinQn(currentVertex->L()) << spinQn(currentVertex->S()) << "_"
393  << waveNameFromTopology(topo, newConvention, static_pointer_cast<isobarDecayVertex>
394  (topo.toVertex(currentVertex->daughter2())));
395  }
396  }
397  {
398  string name = waveName.str();
399  if (newConvention) {
400  replace_all(name, "(", "_");
401  replace_all(name, ")", "_");
402  } else {
403  replace_all(name, "(", "");
404  replace_all(name, ")", "");
405  }
406  return name;
407  }
408 }
409 
410 
411 string
412 waveDescription::waveLaTeXFromTopology(isobarDecayTopology topo,
413  const isobarDecayVertexPtr& currentVertex)
414 {
415  ostringstream waveLaTeX;
416  if (currentVertex == interactionVertexPtr()) {
417  if (not topo.checkTopology() or not topo.checkConsistency()) {
418  printWarn << "decay topology has issues. cannot construct wave LaTeX." << endl;
419  return "";
420  }
421  // X quantum numbers
422  const particle& X = *(topo.XParticle());
423  waveLaTeX << spinQn(X.isospin()) << "^{"<< parityQn(X.G()) << "}"
424  << spinQn(X.J()) << "^{" << parityQn(X.P()) << parityQn(X.C()) << "}"
425  << spinQn(X.spinProj()) << "^{" << parityQn(X.reflectivity()) << "}\\quad & "
426  << waveLaTeXFromTopology(topo, topo.XIsobarDecayVertex());
427  }
428  else if(not (topo.isFsParticle(currentVertex->daughter1())
429  and topo.isFsParticle(currentVertex->daughter1()))){
430  // recurse down decay chain
431  // do this only if not both daughters are fs partiles
432 
433  bool isXdecay= ( currentVertex == topo.XIsobarDecayVertex() );
434 
435  // first daughter
436  string dau1=isobars[currentVertex->daughter1()->name()];
437  if(dau1.length()<2){
438  dau1="{\\bf ";dau1+=currentVertex->daughter1()->name();dau1+="}";
439  }
440  if(!isXdecay)waveLaTeX << "\\rightarrow\\left\\{ ";
441  waveLaTeX << dau1;
442  if (not topo.isFsParticle(currentVertex->daughter1()))
443  waveLaTeX << waveLaTeXFromTopology
444  (topo,
445  static_pointer_cast<isobarDecayVertex>(topo.toVertex(currentVertex->daughter1())));
446  // L, S
447  waveLaTeX << "\\ells{" << spinQn(currentVertex->L()) << "}{" << spinQn(currentVertex->S()) << "}";
448  // second daughter
449  string dau2=isobars[currentVertex->daughter2()->name()];
450  if(dau2.length()<2){
451  dau2="{\\bf ";dau2+=currentVertex->daughter2()->name();dau2+="}"; }
452  waveLaTeX << dau2;
453  if (not topo.isFsParticle(currentVertex->daughter2()))
454  waveLaTeX << waveLaTeXFromTopology
455  (topo,
456  static_pointer_cast<isobarDecayVertex>(topo.toVertex(currentVertex->daughter2())));
457  if(!isXdecay)waveLaTeX << "\\right\\} ";
458  }
459  return waveLaTeX.str();
460 }
461 
462 
463 bool
464 waveDescription::parseKeyFileLocalCopy()
465 {
466  _keyFileParsed = false;
467  if (not _key)
468  _key = new Config();
469  if (not parseLibConfigString(_keyFileLocalCopy, *_key, _debug)) {
470  printWarn << "problems parsing key file string:" << endl;
471  printKeyFileContents(cout);
472  cout << " cannot construct decay topology." << endl;
473  delete _key;
474  _key = 0;
475  return false;
476  }
477  _keyFileParsed = true;
478  return true;
479 }
480 
481 
482 bool
483 waveDescription::constructXParticle(const Setting& XQnKey,
484  particlePtr& X)
485 {
486  if (_debug)
487  printDebug << "reading X quantum numbers from '" << XQnKey.getPath() << "':" << endl;
488  // get X quantum numbers
489  map<string, int> mandatoryXQn, optionalXQn;
490  mandatoryXQn["isospin" ];
491  optionalXQn ["G" ];
492  mandatoryXQn["J" ];
493  mandatoryXQn["P" ];
494  optionalXQn ["C" ];
495  mandatoryXQn["M" ];
496  optionalXQn ["refl" ];
497  optionalXQn ["baryonNmb" ];
498  optionalXQn ["strangeness"];
499  optionalXQn ["charm" ];
500  optionalXQn ["beauty" ];
501  bool success = true;
502  for (map<string, int>::iterator i = mandatoryXQn.begin(); i != mandatoryXQn.end(); ++i)
503  if (not XQnKey.lookupValue(i->first, i->second)) {
504  printWarn << "cannot find integer field '" << i->first << "in "
505  << "'" << XQnKey.getPath() << "'" << endl;
506  success = false;
507  }
508  for (map<string, int>::iterator i = optionalXQn.begin(); i != optionalXQn.end(); ++i)
509  if (not XQnKey.lookupValue(i->first, i->second))
510  i->second = 0;
511  if (not success) {
512  printWarn << "cannot find X quantum numbers. cannot construct decay topology." << endl;
513  return false;
514  }
515  // create X particle
516  X = createParticle("X",
517  mandatoryXQn["isospin"], optionalXQn["G"],
518  mandatoryXQn["J"], mandatoryXQn["P"], optionalXQn["C"],
519  mandatoryXQn["M"], optionalXQn["refl"]);
520  X->setBaryonNmb(optionalXQn["baryonNmb"]);
521  X->setSCB(optionalXQn["strangeness"], optionalXQn["charm"], optionalXQn["beauty"]);
522  if (_debug)
523  printDebug << "constructed X particle: " << X->qnSummary() << endl;
524  return true;
525 }
526 
527 
528 productionVertexPtr
529 waveDescription::mapProductionVertexType(const Setting& prodVertKey,
530  const string& vertType,
531  const particlePtr& X)
532 {
533  productionVertexPtr prodVert;
534  bool success = true;
535  if ((vertType == "diffractiveDissVertex") or (vertType == "leptoProductionVertex")) {
536  const string prodKinParticleNames[] = {"beam", "target", "recoil"};
537  map<string, particlePtr> prodKinParticles;
538  for (unsigned int i = 0; i < sizeof(prodKinParticleNames) / sizeof(string); ++i)
539  prodKinParticles[prodKinParticleNames[i]] = particlePtr();
540  // construct particles
541  const Setting* beamParticleKey = 0; // needed for leptoproduction vertex
542  for (map<string, particlePtr>::iterator i = prodKinParticles.begin();
543  i != prodKinParticles.end(); ++i) {
544  const Setting* particleKey = findLibConfigGroup(prodVertKey, i->first,
545  (i->first != "recoil") ? true : false);
546  if (particleKey) {
547  if (i->first == "beam")
548  beamParticleKey = particleKey;
549  if (not constructParticle(*particleKey, i->second))
550  success = false;
551  } else if (i->first != "recoil")
552  success = false;
553  }
554  if (success) {
555  if (vertType == "diffractiveDissVertex")
556  prodVert = createDiffractiveDissVertex(prodKinParticles["beam"], prodKinParticles["target"],
557  X, prodKinParticles["recoil"]);
558  else if (vertType == "leptoProductionVertex" ) {
559  prodVert = createLeptoProductionVertex(prodKinParticles["beam"], prodKinParticles["target"],
560  X, prodKinParticles["recoil"]);
561  double beamLongPol = 0;
562  if ((beamParticleKey->lookupValue("longPol", beamLongPol))) {
563  if (_debug)
564  printDebug << "setting polarization of beam " << prodKinParticles["beam"]->qnSummary()
565  << " to " << beamLongPol << endl;
566  } else
567  printWarn << "no polarization is given for beam " << prodKinParticles["beam"]->qnSummary()
568  << ". assuming unpolarized beam." << endl;
569  static_pointer_cast<leptoProductionVertex>(prodVert)->setBeamPol(beamLongPol);
570  }
571  }
572  } else
573  printWarn << "unknown production vertex type '" << vertType << "'" << endl;
574  if (_debug) {
575  if (success and prodVert)
576  printDebug << "constructed " << *prodVert << endl;
577  else
578  printWarn << "problems constructing production vertex of type '" << vertType << "'" << endl;
579  }
580  return prodVert;
581 }
582 
583 
584 bool
585 waveDescription::constructProductionVertex(const Setting& rootKey,
586  const particlePtr& X,
587  productionVertexPtr& prodVert)
588 {
589  // find production vertex group
590  const Setting* prodVertKey = findLibConfigGroup(rootKey, "productionVertex");
591  if (not prodVertKey) {
592  printWarn << "cannot find 'productionVertex' group" << endl;
593  return false;
594  }
595  if (_debug)
596  printDebug << "reading production vertex from '" << prodVertKey->getPath() << "':" << endl;
597  //bool success = true;
598  // get vertex type
599  string vertType;
600  if (not prodVertKey->lookupValue("type", vertType)) {
601  printWarn << "cannot find 'type' entry in '" << prodVertKey->getPath() << "'" << endl;
602  //success = false;
603  }
604  // create production vertex
605  prodVert = mapProductionVertexType(*prodVertKey, vertType, X);
606  return (prodVert) ? true : false;
607 }
608 
609 
610 bool
611 waveDescription::constructParticle(const Setting& particleKey,
612  particlePtr& particle,
613  const bool requirePartInTable)
614 {
615  if (_debug)
616  printDebug << "reading particle information from '" << particleKey.getPath() << "':" << endl;
617  string name;
618  if (particleKey.lookupValue("name", name)) {
619  particle = createParticle(name, requirePartInTable);
620  int spinProj;
621  if (particleKey.lookupValue("spinProj", spinProj))
622  particle->setSpinProj(spinProj);
623  if (_debug)
624  printDebug << "created particle " << particle->qnSummary() << flush;
625  int index;
626  if (particleKey.lookupValue("index", index)) {
627  particle->setIndex(index);
628  if (_debug)
629  cout << " with index [" << index << "]" << flush;
630  }
631  if (_debug)
632  cout << endl;
633  return true;
634  } else {
635  printWarn << "cannot construct particle from entry '" << particleKey.getPath() << "'. "
636  << "name field is missing." << endl;
637  particle = particlePtr();
638  return false;
639  }
640 }
641 
642 
643 massDependencePtr
644 waveDescription::mapMassDependenceType(const string& massDepType)
645 {
646  massDependencePtr massDep;
647  if ( (massDepType == "BreitWigner")
648  or (massDepType == "")) // default mass dependence
649  massDep = createRelativisticBreitWigner();
650  else if (massDepType == "flat")
651  massDep = createFlatMassDependence();
652  else if (massDepType == "f_0(980)")
653  massDep = createF0980BreitWigner();
654  else if (massDepType == "piPiSWaveAuMorganPenningtonM")
655  massDep = createPiPiSWaveAuMorganPenningtonM();
656  else if (massDepType == "piPiSWaveAuMorganPenningtonVes")
657  massDep = createPiPiSWaveAuMorganPenningtonVes();
658  else if (massDepType == "piPiSWaveAuMorganPenningtonKachaev")
659  massDep = createPiPiSWaveAuMorganPenningtonKachaev();
660  else if (massDepType == "rhoPrime")
661  massDep = createRhoPrimeMassDep();
662  else {
663  printWarn << "unknown mass dependence '" << massDepType << "'. using Breit-Wigner." << endl;
664  massDep = createRelativisticBreitWigner();
665  }
666  return massDep;
667 }
668 
669 
670 bool
671 waveDescription::constructDecayVertex(const Setting& parentKey,
672  const particlePtr& parentParticle,
673  vector<isobarDecayVertexPtr>& decayVertices,
674  vector<particlePtr>& fsParticles,
675  const bool fromTemplate)
676 {
677  if (_debug)
678  printDebug << "reading decay vertex information from '" << parentKey.getPath() << "':" << endl;
679  bool success = true;
680 
681  const Setting* fsPartKeys = findLibConfigList(parentKey, "fsParticles", false);
682  vector<particlePtr> fsDaughters;
683  if (fsPartKeys)
684  for (int i = 0; i < fsPartKeys->getLength(); ++i) {
685  particlePtr p;
686  if (constructParticle((*fsPartKeys)[i], p)) {
687  fsDaughters.push_back(p);
688  fsParticles.push_back(p);
689  } else
690  success = false;
691  }
692 
693  const Setting* isobarKeys = findLibConfigList(parentKey, "isobars", false);
694  vector<particlePtr> isobarDaughters;
695  if (isobarKeys)
696  for (int i = 0; i < isobarKeys->getLength(); ++i) {
697  particlePtr p;
698  if (constructParticle((*isobarKeys)[i], p, not fromTemplate))
699  isobarDaughters.push_back(p);
700  else
701  success = false;
702  success &= constructDecayVertex((*isobarKeys)[i], isobarDaughters.back(),
703  decayVertices, fsParticles, fromTemplate);
704  }
705 
706  const unsigned int nmbDaughters = fsDaughters.size() + isobarDaughters.size();
707  if (nmbDaughters != 2) {
708  printWarn << "cannot construct isobar vertex, because number of daughters "
709  << "(" << nmbDaughters << ") in '" << parentKey.getPath() << "' "
710  << "is not 2" << endl;
711  success = false;
712  }
713  if (not success)
714  return false;
715 
716  // get isobar vertex parameters
717  int L = 0, S = 0;
718  if (( not parentKey.lookupValue("L", L)
719  or not parentKey.lookupValue("S", S)) and not fromTemplate)
720  printWarn << "Either L or S are not specified in '" << parentKey.getPath() << "'. "
721  << "using zero." << endl;
722 
723  // get mass dependence
724  massDependencePtr massDep;
725  if (parentParticle->bareName() != "X") {
726  string massDepType = "";
727  const Setting* massDepKey = findLibConfigGroup(parentKey, "massDep", false);
728  if (massDepKey)
729  massDepKey->lookupValue("name", massDepType);
730  massDep = mapMassDependenceType(massDepType);
731  }
732 
733  // if there is 1 final state particle and 1 isobar put them in the
734  // same order as in the key file
735  vector<particlePtr> daughters(2, particlePtr());
736  if ((fsDaughters.size() == 1) and (isobarDaughters.size() == 1)) {
737  if (fsPartKeys->getIndex() < isobarKeys->getIndex()) {
738  daughters[0] = fsDaughters [0];
739  daughters[1] = isobarDaughters[0];
740  } else {
741  daughters[0] = isobarDaughters[0];
742  daughters[1] = fsDaughters [0];
743  }
744  } else if (fsDaughters.size() == 2)
745  daughters = fsDaughters;
746  else if (isobarDaughters.size() == 2)
747  daughters = isobarDaughters;
748 
749  // construct isobar decay vertex
750  decayVertices.push_back(createIsobarDecayVertex(parentParticle, daughters[0],
751  daughters[1], L, S, massDep));
752  if (_debug)
753  printDebug << "constructed " << *decayVertices.back() << endl;
754  return true;
755 }
756 
757 
758 isobarAmplitudePtr
759 waveDescription::mapAmplitudeType(const string& formalismType,
760  const isobarDecayTopologyPtr& topo)
761 {
762  isobarAmplitudePtr amplitude;
763  if (formalismType == "helicity")
764  amplitude = createIsobarHelicityAmplitude(topo);
765  else if (formalismType == "canonical")
766  amplitude = createIsobarCanonicalAmplitude(topo);
767  else {
768  printWarn << "unknown amplitude formalism '" << formalismType << "'. "
769  << "using helicity formalism." << endl;
770  amplitude = createIsobarHelicityAmplitude(topo);
771  }
772  return amplitude;
773 }
774 
775 
776 bool
777 waveDescription::setProductionVertexKeys(Setting& prodVertKey,
778  const productionVertexPtr& prodVert)
779 {
780  const string prodVertName = prodVert->name();
781  if ((prodVertName == "diffractiveDissVertex") or (prodVertName == "leptoProductionVertex")) {
782  if (_debug)
783  printDebug << "setting keys for " << *prodVert << endl;
784  prodVertKey.add("type", Setting::TypeString) = prodVertName;
785  map<string, particlePtr> prodKinParticles;
786  if (prodVertName == "diffractiveDissVertex") {
787  const diffractiveDissVertexPtr& vert = static_pointer_cast<diffractiveDissVertex>(prodVert);
788  prodKinParticles["beam" ] = vert->beam();
789  prodKinParticles["target"] = vert->target();
790  if (vert->recoil()->name() != vert->target()->name())
791  prodKinParticles["recoil"] = vert->recoil();
792  } else if (prodVertName == "leptoProductionVertex") {
793  const leptoProductionVertexPtr& vert = static_pointer_cast<leptoProductionVertex>(prodVert);
794  prodKinParticles["beam" ] = vert->beamLepton();
795  prodKinParticles["target"] = vert->target();
796  if (vert->recoil()->name() != vert->target()->name())
797  prodKinParticles["recoil"] = vert->recoil();
798  }
799  // write particles
800  Setting* beamParticleKey = 0; // needed for leptoproduction vertex
801  for (map<string, particlePtr>::iterator i = prodKinParticles.begin();
802  i != prodKinParticles.end(); ++i) {
803  Setting& particleKey = prodVertKey.add(i->first, Setting::TypeGroup);
804  particleKey.add("name", Setting::TypeString) = i->second->name();
805  if (i->first == "beam")
806  beamParticleKey = &particleKey;
807  }
808  if ((prodVertName == "leptoProductionVertex") and beamParticleKey)
809  beamParticleKey->add("longPol", Setting::TypeFloat)
810  = static_pointer_cast<leptoProductionVertex>(prodVert)->beamPol();
811  return true;
812  } else {
813  printWarn << "setting of keys for production vertex of this type is not yet implemented:"
814  << *prodVert << endl;
815  return false;
816  }
817 }
818 
819 
820 bool
821 waveDescription::setXQuantumNumbersKeys(Setting& XQnKey,
822  const particle& X)
823 {
824  if (_debug)
825  printDebug << "setting quantum number keys for " << X.qnSummary() << endl;
826  XQnKey.add("isospin", Setting::TypeInt) = X.isospin();
827  if (X.G() != 0)
828  XQnKey.add("G", Setting::TypeInt) = X.G();
829  XQnKey.add("J", Setting::TypeInt) = X.J();
830  XQnKey.add("P", Setting::TypeInt) = X.P();
831  if (X.C() != 0)
832  XQnKey.add("C", Setting::TypeInt) = X.C();
833  XQnKey.add("M", Setting::TypeInt) = X.spinProj();
834  if (X.reflectivity() != 0)
835  XQnKey.add("refl", Setting::TypeInt) = X.reflectivity();
836  if (X.baryonNmb() != 0)
837  XQnKey.add("baryonNmb", Setting::TypeInt) = X.baryonNmb();
838  if (X.strangeness() != 0)
839  XQnKey.add("strangeness", Setting::TypeInt) = X.strangeness();
840  if (X.charm() != 0)
841  XQnKey.add("charm", Setting::TypeInt) = X.charm();
842  if (X.beauty() != 0)
843  XQnKey.add("beauty", Setting::TypeInt) = X.beauty();
844  return true;
845 }
846 
847 
848 bool
849 waveDescription::setMassDependence(Setting& isobarDecayKey,
850  const massDependence& massDep)
851 {
852  const string massDepName = massDep.name();
853  if (massDepName == "flatMassDependence")
854  // default for X
855  return true;
856  else if (massDepName == "relativisticBreitWigner")
857  // default mass dependence for isobars
858  return true;
859  else {
860  if (_debug)
861  printDebug << "setting key for '" << massDep << "' mass dependence to "
862  << "'" << massDepName << "'" << endl;
863  Setting& massDepKey = isobarDecayKey.add("massDep", Setting::TypeGroup);
864  massDepKey.add("name", Setting::TypeString) = massDepName;
865  }
866  return true;
867 }
868 
869 
870 bool
871 waveDescription::setXDecayKeys(Setting& parentDecayKey,
872  const isobarDecayTopology& topo,
873  const isobarDecayVertex& vert)
874 {
875  if (_debug)
876  printDebug << "setting keys for decay " << vert << endl;
877  setMassDependence(parentDecayKey, *vert.massDependence());
878  vector<particlePtr> fsParticles;
879  vector<particlePtr> isobars;
880  for (unsigned int i = 0; i < vert.nmbOutParticles(); ++i) {
881  const particlePtr& part = vert.outParticles()[i];
882  if (topo.isFsParticle(part))
883  fsParticles.push_back(part);
884  else
885  isobars.push_back(part);
886  }
887  bool success = true;
888  if (isobars.size() > 0) {
889  Setting& isobarsKey = parentDecayKey.add("isobars", Setting::TypeList);
890  for (unsigned int i = 0; i < isobars.size(); ++i) {
891  Setting& isobarKey = isobarsKey.add(Setting::TypeGroup);
892  isobarKey.add("name", Setting::TypeString) = isobars[i]->name();
893  if (not setXDecayKeys(isobarKey, topo,
894  *static_pointer_cast<isobarDecayVertex>(topo.toVertex(isobars[i]))))
895  success = false;
896  }
897  }
898  parentDecayKey.add("L", Setting::TypeInt) = (int)vert.L();
899  parentDecayKey.add("S", Setting::TypeInt) = (int)vert.S();
900  if (fsParticles.size() > 0) {
901  Setting& fsParticlesKey = parentDecayKey.add("fsParticles", Setting::TypeList);
902  for (unsigned int i = 0; i < fsParticles.size(); ++i) {
903  Setting& fsParticleKey = fsParticlesKey.add(Setting::TypeGroup);
904  fsParticleKey.add("name", Setting::TypeString) = fsParticles[i]->name();
905  }
906  }
907  return success;
908 }
909 
910 
911 bool
912 waveDescription::setKeysFromTopology(Setting& rootKey,
913  const isobarDecayTopology& topo,
914  const bool setProdVert)
915 {
916  if (not topo.checkTopology() or not topo.checkConsistency()) {
917  printWarn << "decay topology has issues. cannot write key file." << endl;
918  return false;
919  }
920  if (_debug)
921  printDebug << "setting keys for " << topo;
922  if (setProdVert) {
923  Setting& prodVertKey = rootKey.add("productionVertex", Setting::TypeGroup);
924  setProductionVertexKeys(prodVertKey, topo.productionVertex());
925  }
926  Setting& decayVertKey = rootKey.add ("decayVertex", Setting::TypeGroup);
927  Setting& XQnKey = decayVertKey.add("XQuantumNumbers", Setting::TypeGroup);
928  Setting& XDecayKey = decayVertKey.add("XDecay", Setting::TypeGroup);
929  if (not setXQuantumNumbersKeys(XQnKey, *(topo.XParticle()))) {
930  printWarn << "problems setting X quantum numbers" << endl;
931  return false;
932  }
933  if (not setXDecayKeys(XDecayKey, topo, *(topo.XIsobarDecayVertex()))) {
934  printWarn << "problems setting X decay" << endl;
935  return false;
936  }
937  return true;
938 }
939 
940 
941 bool
942 waveDescription::setAmplitude(Setting& amplitudeKey,
943  const isobarAmplitude& amplitude)
944 {
945  if (_debug)
946  printDebug << "setting amplitude keys." << endl;
947  string formalism = "";
948  const string amplitudeName = amplitude.name();
949  if (amplitudeName == "isobarCanonicalAmplitude") {
950  formalism = "canonical";
951  } else if (amplitudeName != "isobarHelicityAmplitude") {
952  printWarn << "unknown amplitude type '" << amplitudeName << "'" << endl;
953  return false;
954  }
955  if (formalism != "") {
956  if (_debug)
957  printDebug << "setting 'formalism' key for '" << amplitudeName << "' to "
958  << "'" << formalism << "'" << endl;
959  amplitudeKey.add("formalism", Setting::TypeString) = formalism;
960  }
961  if (not amplitude.boseSymmetrization())
962  amplitudeKey.add("boseSymmetrize", Setting::TypeBoolean) = false;
963  if (not amplitude.isospinSymmetrization())
964  amplitudeKey.add("isospinSymmetrize", Setting::TypeBoolean) = false;
965  if (not amplitude.reflectivityBasis())
966  amplitudeKey.add("useReflectivityBasis", Setting::TypeBoolean) = false;
967  return true;
968 }
969 
970 
971 bool
972 waveDescription::setKeysFromAmplitude(Setting& rootKey,
973  const isobarAmplitude& amplitude,
974  const bool setProdVert)
975 {
976  if (_debug)
977  printInfo << "setting keys for "<< amplitude;
978  if (not setKeysFromTopology(rootKey, *amplitude.decayTopology(), setProdVert)) {
979  printWarn << "problems setting keys for decay topology" << endl;
980  return false;
981  }
982  Setting& amplitudeKey = rootKey.add("amplitude", Setting::TypeGroup);
983  if (not setAmplitude(amplitudeKey, amplitude)) {
984  printWarn << "problems setting amplitude paramaters" << endl;
985  return false;
986  }
987  return true;
988 }
989 
990 
991 bool
992 waveDescription::writeKeyFile(FILE& outStream,
993  const isobarDecayTopology& topo,
994  const bool writeProdVert)
995 {
996  Config key;
997  Setting& rootKey = key.getRoot();
998  if (not setKeysFromTopology(rootKey, topo, writeProdVert)) {
999  printWarn << "problems writing keys for decay topology. cannot write key file." << endl;
1000  return false;
1001  }
1002  try {
1003  key.write(&outStream);
1004  } catch (const FileIOException& ioEx) {
1005  printWarn << "I/O error while writing key file" << endl;
1006  return false;
1007  }
1008  return true;
1009 }
1010 
1011 
1012 bool
1013 waveDescription::writeKeyFile(FILE& outStream,
1014  const isobarAmplitude& amplitude,
1015  const bool writeProdVert)
1016 {
1017  Config key;
1018  Setting& rootKey = key.getRoot();
1019  if (not setKeysFromAmplitude(rootKey, amplitude)) {
1020  printWarn << "problems writing keys for amplitude. cannot write key file." << endl;
1021  return false;
1022  }
1023  try {
1024  key.write(&outStream);
1025  } catch (const FileIOException& ioEx) {
1026  printWarn << "I/O error while writing key file" << endl;
1027  return false;
1028  }
1029  return true;
1030 }
Generated by   doxygen 1.8.1.2