evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtDecayTable.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen 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 * EvtGen 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 EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
22 
25 #include "EvtGenBase/EvtModel.hh"
27 #include "EvtGenBase/EvtPDL.hh"
28 #include "EvtGenBase/EvtParser.hh"
31 #include "EvtGenBase/EvtPatches.hh"
32 #include "EvtGenBase/EvtRadCorr.hh"
33 #include "EvtGenBase/EvtRandom.hh"
34 #include "EvtGenBase/EvtReport.hh"
36 
37 #include <ctype.h>
38 #include <fstream>
39 #include <iomanip>
40 #include <iostream>
41 #include <sstream>
42 #include <stdlib.h>
43 #include <string.h>
44 
45 using std::endl;
46 using std::fstream;
47 using std::ifstream;
48 
50 {
51  _decaytable.clear();
52 }
53 
55 {
56  _decaytable.clear();
57 }
58 
60 {
61  static EvtDecayTable* theDecayTable = 0;
62 
63  if ( theDecayTable == 0 ) {
64  theDecayTable = new EvtDecayTable();
65  }
66 
67  return theDecayTable;
68 }
69 
70 int EvtDecayTable::getNMode( int ipar )
71 {
72  return _decaytable[ipar].getNMode();
73 }
74 
75 EvtDecayBase* EvtDecayTable::getDecay( int ipar, int imode )
76 {
77  return _decaytable[ipar].getDecayModel( imode );
78 }
79 
81 {
82  for ( size_t i = 0; i < EvtPDL::entries(); i++ ) {
83  _decaytable[i].printSummary();
84  }
85 }
86 
88 {
89  int partnum;
90 
91  partnum = p->getId().getAlias();
92 
93  if ( _decaytable[partnum].getNMode() == 0 )
94  return 0;
95  return _decaytable[partnum].getDecayModel( p );
96 }
97 
98 void EvtDecayTable::readDecayFile( const std::string dec_name, bool verbose )
99 {
100  if ( _decaytable.size() < EvtPDL::entries() )
101  _decaytable.resize( EvtPDL::entries() );
102  EvtModel& modelist = EvtModel::instance();
103  EvtExtGeneratorCommandsTable* extGenCommands =
105  std::string colon( ":" ), equals( "=" );
106 
107  int i;
108 
109  EvtGenReport( EVTGEN_INFO, "EvtGen" )
110  << "In readDecayFile, reading:" << dec_name.c_str() << endl;
111 
112  ifstream fin;
113 
114  fin.open( dec_name.c_str() );
115  if ( !fin ) {
116  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
117  << "Could not open " << dec_name.c_str() << endl;
118  }
119  fin.close();
120 
121  EvtParser parser;
122  parser.read( dec_name );
123 
124  int itok;
125 
126  int hasend = 0;
127 
128  std::string token;
129 
130  for ( itok = 0; itok < parser.getNToken(); itok++ ) {
131  token = parser.getToken( itok );
132 
133  if ( token == "End" )
134  hasend = 1;
135  }
136 
137  if ( !hasend ) {
138  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
139  << "Could not find an 'End' in " << dec_name.c_str() << endl;
140  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
141  << "Will terminate execution." << endl;
142  ::abort();
143  }
144 
145  std::string model, parent, sdaug;
146 
147  EvtId ipar;
148 
149  int n_daugh;
150  EvtId daught[MAX_DAUG];
151  double brfr;
152 
153  int itoken = 0;
154 
155  std::vector<EvtModelAlias> modelAliasList;
156 
157  do {
158  token = parser.getToken( itoken++ );
159 
160  //Easy way to turn off photos... Lange September 5, 2000
161  if ( token == "noPhotos" ) {
163  if ( verbose )
164  EvtGenReport( EVTGEN_INFO, "EvtGen" )
165  << "As requested, PHOTOS will be turned off." << endl;
166  } else if ( token == "yesPhotos" ) {
168  if ( verbose )
169  EvtGenReport( EVTGEN_INFO, "EvtGen" )
170  << "As requested, PHOTOS will be turned on for all decays."
171  << endl;
172  } else if ( token == "normalPhotos" ) {
174  if ( verbose )
175  EvtGenReport( EVTGEN_INFO, "EvtGen" )
176  << "As requested, PHOTOS will be turned on only when requested."
177  << endl;
178  } else if ( token == "Alias" ) {
179  std::string newname;
180  std::string oldname;
181 
182  newname = parser.getToken( itoken++ );
183  oldname = parser.getToken( itoken++ );
184 
185  EvtId id = EvtPDL::getId( oldname );
186 
187  if ( id == EvtId( -1, -1 ) ) {
188  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
189  << "Unknown particle name:" << oldname.c_str()
190  << " on line " << parser.getLineofToken( itoken ) << endl;
191  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
192  << "Will terminate execution!" << endl;
193  ::abort();
194  }
195 
196  EvtPDL::alias( id, newname );
197  if ( _decaytable.size() < EvtPDL::entries() )
198  _decaytable.resize( EvtPDL::entries() );
199 
200  } else if ( token == "ModelAlias" ) {
201  std::vector<std::string> modelArgList;
202 
203  std::string aliasName = parser.getToken( itoken++ );
204  std::string modelName = parser.getToken( itoken++ );
205 
206  std::string nameTemp;
207  do {
208  nameTemp = parser.getToken( itoken++ );
209  if ( nameTemp != ";" ) {
210  modelArgList.push_back( nameTemp );
211  }
212  } while ( nameTemp != ";" );
213  EvtModelAlias newAlias( aliasName, modelName, modelArgList );
214  modelAliasList.push_back( newAlias );
215  } else if ( token == "ChargeConj" ) {
216  std::string aname;
217  std::string abarname;
218 
219  aname = parser.getToken( itoken++ );
220  abarname = parser.getToken( itoken++ );
221 
222  EvtId a = EvtPDL::getId( aname );
223  EvtId abar = EvtPDL::getId( abarname );
224 
225  if ( a == EvtId( -1, -1 ) ) {
226  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
227  << "Unknown particle name:" << aname.c_str() << " on line "
228  << parser.getLineofToken( itoken ) << endl;
229  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
230  << "Will terminate execution!" << endl;
231  ::abort();
232  }
233 
234  if ( abar == EvtId( -1, -1 ) ) {
235  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
236  << "Unknown particle name:" << abarname.c_str()
237  << " on line " << parser.getLineofToken( itoken ) << endl;
238  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
239  << "Will terminate execution!" << endl;
240  ::abort();
241  }
242 
243  EvtPDL::aliasChgConj( a, abar );
244 
245  } else if ( token == "JetSetPar" ) {
246  // Check if any old Pythia 6 commands are present
247  std::string pythiaCommand = parser.getToken( itoken++ );
248 
249  Command command;
250 
251  // The old command format is NAME(INT)=VALUE
252  int i1 = pythiaCommand.find_first_of( "(" );
253  int i2 = pythiaCommand.find_first_of( ")" );
254  int i3 = pythiaCommand.find_first_of( "=" );
255 
256  std::string pythiaModule = pythiaCommand.substr( 0, i1 );
257  std::string pythiaParam = pythiaCommand.substr( i1 + 1, i2 - i1 - 1 );
258  std::string pythiaValue = pythiaCommand.substr( i3 + 1 );
259 
260  command["MODULE"] = pythiaModule;
261  command["PARAM"] = pythiaParam;
262  command["VALUE"] = pythiaValue;
263 
264  command["GENERATOR"] = "Both";
265  command["VERSION"] = "PYTHIA6";
266 
267  extGenCommands->addCommand( "PYTHIA", command );
268 
269  } else if ( modelist.isCommand( token ) ) {
270  std::string cnfgstr;
271 
272  cnfgstr = parser.getToken( itoken++ );
273 
274  modelist.storeCommand( token, cnfgstr );
275 
276  } else if ( token == "PythiaGenericParam" ||
277  token == "PythiaAliasParam" || token == "PythiaBothParam" ) {
278  // Read in any Pythia 8 commands, which will be of the form
279  // pythia<type>Param module:param=value, with no spaces in the parameter
280  // string! Here, <type> specifies whether the command is for generic
281  // decays, alias decays, or both.
282 
283  // Pythia 6 commands will be defined by the old JetSetPar command
284  // name, which is handled by the modelist.isCommand() statement above.
285 
286  std::string pythiaCommand = parser.getToken( itoken++ );
287  std::string pythiaModule( "" ), pythiaParam( "" ), pythiaValue( "" );
288 
289  // Separate out the string into the 3 sections using the delimiters
290  // ":" and "=".
291 
292  std::vector<std::string> pComVect1 =
293  this->splitString( pythiaCommand, colon );
294 
295  if ( pComVect1.size() == 2 ) {
296  pythiaModule = pComVect1[0];
297 
298  std::string pCom2 = pComVect1[1];
299 
300  std::vector<std::string> pComVect2 = this->splitString( pCom2,
301  equals );
302 
303  if ( pComVect2.size() == 2 ) {
304  pythiaParam = pComVect2[0];
305  pythiaValue = pComVect2[1];
306  }
307  }
308 
309  // Define the Pythia 8 command and pass it to the external generator
310  // command list.
311  Command command;
312  if ( token == "PythiaGenericParam" ) {
313  command["GENERATOR"] = "Generic";
314  } else if ( token == "PythiaAliasParam" ) {
315  command["GENERATOR"] = "Alias";
316  } else {
317  command["GENERATOR"] = "Both";
318  }
319 
320  command["MODULE"] = pythiaModule;
321  command["PARAM"] = pythiaParam;
322  command["VALUE"] = pythiaValue;
323 
324  command["VERSION"] = "PYTHIA8";
325  extGenCommands->addCommand( "PYTHIA", command );
326 
327  } else if ( token == "CDecay" ) {
328  std::string name;
329 
330  name = parser.getToken( itoken++ );
331  ipar = EvtPDL::getId( name );
332 
333  if ( ipar == EvtId( -1, -1 ) ) {
334  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
335  << "Unknown particle name:" << name.c_str() << " on line "
336  << parser.getLineofToken( itoken - 1 ) << endl;
337  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
338  << "Will terminate execution!" << endl;
339  ::abort();
340  }
341 
342  EvtId cipar = EvtPDL::chargeConj( ipar );
343 
344  if ( _decaytable[ipar.getAlias()].getNMode() != 0 ) {
345  if ( verbose )
346  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
347  << "Redefined decay of " << name.c_str() << " in CDecay"
348  << endl;
349 
350  _decaytable[ipar.getAlias()].removeDecay();
351  }
352 
353  //take contents of cipar and conjugate and store in ipar
354  _decaytable[ipar.getAlias()].makeChargeConj(
355  &_decaytable[cipar.getAlias()] );
356 
357  } else if ( token == "Define" ) {
358  std::string name;
359 
360  name = parser.getToken( itoken++ );
361  // value=atof(parser.getToken(itoken++).c_str());
362 
363  EvtSymTable::define( name, parser.getToken( itoken++ ) );
364 
365  //New code Lange April 10, 2001 - allow the user
366  //to change particle definitions of EXISTING
367  //particles on the fly
368  } else if ( token == "Particle" ) {
369  std::string pname;
370  pname = parser.getToken( itoken++ );
371  if ( verbose )
372  EvtGenReport( EVTGEN_INFO, "EvtGen" ) << pname.c_str() << endl;
373  //There should be at least the mass
374  double newMass = atof( parser.getToken( itoken++ ).c_str() );
375  EvtId thisPart = EvtPDL::getId( pname );
376  double newWidth = EvtPDL::getMeanMass( thisPart );
377  if ( parser.getNToken() > 3 )
378  newWidth = atof( parser.getToken( itoken++ ).c_str() );
379 
380  //Now make the change!
381  EvtPDL::reSetMass( thisPart, newMass );
382  EvtPDL::reSetWidth( thisPart, newWidth );
383 
384  if ( verbose )
385  EvtGenReport( EVTGEN_INFO, "EvtGen" )
386  << "Changing particle properties of " << pname.c_str()
387  << " Mass=" << newMass << " Width=" << newWidth << endl;
388 
389  } else if ( token == "ChangeMassMin" ) {
390  std::string pname;
391  pname = parser.getToken( itoken++ );
392  double tmass = atof( parser.getToken( itoken++ ).c_str() );
393 
394  EvtId thisPart = EvtPDL::getId( pname );
395  EvtPDL::reSetMassMin( thisPart, tmass );
396  if ( verbose )
397  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
398  << "Refined minimum mass for "
399  << EvtPDL::name( thisPart ).c_str() << " to be " << tmass
400  << endl;
401 
402  } else if ( token == "ChangeMassMax" ) {
403  std::string pname;
404  pname = parser.getToken( itoken++ );
405  double tmass = atof( parser.getToken( itoken++ ).c_str() );
406  EvtId thisPart = EvtPDL::getId( pname );
407  EvtPDL::reSetMassMax( thisPart, tmass );
408  if ( verbose )
409  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
410  << "Refined maximum mass for "
411  << EvtPDL::name( thisPart ).c_str() << " to be " << tmass
412  << endl;
413 
414  } else if ( token == "IncludeBirthFactor" ) {
415  std::string pname;
416  pname = parser.getToken( itoken++ );
417  bool yesno = false;
418  if ( parser.getToken( itoken++ ) == "yes" )
419  yesno = true;
420  EvtId thisPart = EvtPDL::getId( pname );
421  EvtPDL::includeBirthFactor( thisPart, yesno );
422  if ( verbose ) {
423  if ( yesno )
424  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
425  << "Include birth factor for "
426  << EvtPDL::name( thisPart ).c_str() << endl;
427  if ( !yesno )
428  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
429  << "No longer include birth factor for "
430  << EvtPDL::name( thisPart ).c_str() << endl;
431  }
432 
433  } else if ( token == "IncludeDecayFactor" ) {
434  std::string pname;
435  pname = parser.getToken( itoken++ );
436  bool yesno = false;
437  if ( parser.getToken( itoken++ ) == "yes" )
438  yesno = true;
439  EvtId thisPart = EvtPDL::getId( pname );
440  EvtPDL::includeDecayFactor( thisPart, yesno );
441  if ( verbose ) {
442  if ( yesno )
443  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
444  << "Include decay factor for "
445  << EvtPDL::name( thisPart ).c_str() << endl;
446  if ( !yesno )
447  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
448  << "No longer include decay factor for "
449  << EvtPDL::name( thisPart ).c_str() << endl;
450  }
451  } else if ( token == "LSNONRELBW" ) {
452  std::string pname;
453  pname = parser.getToken( itoken++ );
454  EvtId thisPart = EvtPDL::getId( pname );
455  std::string tstr = "NONRELBW";
456  EvtPDL::changeLS( thisPart, tstr );
457  if ( verbose )
458  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
459  << "Change lineshape to non-rel BW for "
460  << EvtPDL::name( thisPart ).c_str() << endl;
461  } else if ( token == "LSFLAT" ) {
462  std::string pname;
463  pname = parser.getToken( itoken++ );
464  EvtId thisPart = EvtPDL::getId( pname );
465  std::string tstr = "FLAT";
466  EvtPDL::changeLS( thisPart, tstr );
467  if ( verbose )
468  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
469  << "Change lineshape to flat for "
470  << EvtPDL::name( thisPart ).c_str() << endl;
471  } else if ( token == "LSMANYDELTAFUNC" ) {
472  std::string pname;
473  pname = parser.getToken( itoken++ );
474  EvtId thisPart = EvtPDL::getId( pname );
475  std::string tstr = "MANYDELTAFUNC";
476  EvtPDL::changeLS( thisPart, tstr );
477  if ( verbose )
478  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
479  << "Change lineshape to spikes for "
480  << EvtPDL::name( thisPart ).c_str() << endl;
481 
482  } else if ( token == "BlattWeisskopf" ) {
483  std::string pname;
484  pname = parser.getToken( itoken++ );
485  double tnum = atof( parser.getToken( itoken++ ).c_str() );
486  EvtId thisPart = EvtPDL::getId( pname );
487  EvtPDL::reSetBlatt( thisPart, tnum );
488  if ( verbose )
489  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
490  << "Redefined Blatt-Weisskopf factor "
491  << EvtPDL::name( thisPart ).c_str() << " to be " << tnum
492  << endl;
493  } else if ( token == "BlattWeisskopfBirth" ) {
494  std::string pname;
495  pname = parser.getToken( itoken++ );
496  double tnum = atof( parser.getToken( itoken++ ).c_str() );
497  EvtId thisPart = EvtPDL::getId( pname );
498  EvtPDL::reSetBlattBirth( thisPart, tnum );
499  if ( verbose )
500  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
501  << "Redefined Blatt-Weisskopf birth factor "
502  << EvtPDL::name( thisPart ).c_str() << " to be " << tnum
503  << endl;
504  } else if ( token == "SetLineshapePW" ) {
505  std::string pname;
506  pname = parser.getToken( itoken++ );
507  EvtId thisPart = EvtPDL::getId( pname );
508  std::string pnameD1 = parser.getToken( itoken++ );
509  EvtId thisD1 = EvtPDL::getId( pnameD1 );
510  std::string pnameD2 = parser.getToken( itoken++ );
511  EvtId thisD2 = EvtPDL::getId( pnameD2 );
512  int pw = atoi( parser.getToken( itoken++ ).c_str() );
513  if ( verbose )
514  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
515  << "Redefined Partial wave for " << pname.c_str() << " to "
516  << pnameD1.c_str() << " " << pnameD2.c_str() << " (" << pw
517  << ")" << endl;
518  EvtPDL::setPWForDecay( thisPart, pw, thisD1, thisD2 );
519  EvtPDL::setPWForBirthL( thisD1, pw, thisPart, thisD2 );
520  EvtPDL::setPWForBirthL( thisD2, pw, thisPart, thisD1 );
521 
522  } else if ( token == "Decay" ) {
523  std::string temp_fcn_new_model;
524 
525  EvtDecayBase* temp_fcn_new;
526 
527  double brfrsum = 0.0;
528 
529  parent = parser.getToken( itoken++ );
530  ipar = EvtPDL::getId( parent );
531 
532  if ( ipar == EvtId( -1, -1 ) ) {
533  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
534  << "Unknown particle name:" << parent.c_str() << " on line "
535  << parser.getLineofToken( itoken - 1 ) << endl;
536  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
537  << "Will terminate execution!" << endl;
538  ::abort();
539  }
540 
541  if ( _decaytable[ipar.getAlias()].getNMode() != 0 ) {
542  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
543  << "Redefined decay of " << parent.c_str() << endl;
544  _decaytable[ipar.getAlias()].removeDecay();
545  }
546 
547  do {
548  token = parser.getToken( itoken++ );
549 
550  if ( token != "Enddecay" ) {
551  i = 0;
552  while ( token.c_str()[i++] != 0 ) {
553  if ( isalpha( token.c_str()[i] ) ) {
554  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
555  << "Expected to find a branching fraction or Enddecay "
556  << "but found:" << token.c_str() << " on line "
557  << parser.getLineofToken( itoken - 1 ) << endl;
558  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
559  << "Possibly to few arguments to model "
560  << "on previous line!" << endl;
561  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
562  << "Will terminate execution!" << endl;
563  ::abort();
564  }
565  }
566 
567  brfr = atof( token.c_str() );
568 
569  int isname =
570  EvtPDL::getId( parser.getToken( itoken ) ).getId() >= 0;
571  int ismodel = modelist.isModel( parser.getToken( itoken ) );
572 
573  if ( !( isname || ismodel ) ) {
574  //see if this is an aliased model
575  for ( size_t iAlias = 0; iAlias < modelAliasList.size();
576  iAlias++ ) {
577  if ( modelAliasList[iAlias].matchAlias(
578  parser.getToken( itoken ) ) ) {
579  ismodel = 2;
580  break;
581  }
582  }
583  }
584 
585  if ( !( isname || ismodel ) ) {
586  EvtGenReport( EVTGEN_INFO, "EvtGen" )
587  << parser.getToken( itoken ).c_str()
588  << " is neither a particle name nor "
589  << "the name of a model. " << endl;
590  EvtGenReport( EVTGEN_INFO, "EvtGen" )
591  << "It was encountered on line "
592  << parser.getLineofToken( itoken )
593  << " of the decay file." << endl;
594  EvtGenReport( EVTGEN_INFO, "EvtGen" )
595  << "Please fix it. Thank you." << endl;
596  EvtGenReport( EVTGEN_INFO, "EvtGen" )
597  << "Be sure to check that the "
598  << "correct case has been used. \n";
599  EvtGenReport( EVTGEN_INFO, "EvtGen" )
600  << "Terminating execution. \n";
601  ::abort();
602 
603  itoken++;
604  }
605 
606  n_daugh = 0;
607 
608  while ( EvtPDL::getId( parser.getToken( itoken ) ).getId() >=
609  0 ) {
610  sdaug = parser.getToken( itoken++ );
611  daught[n_daugh++] = EvtPDL::getId( sdaug );
612  if ( daught[n_daugh - 1] == EvtId( -1, -1 ) ) {
613  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
614  << "Unknown particle name:" << sdaug.c_str()
615  << " on line "
616  << parser.getLineofToken( itoken ) << endl;
617  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
618  << "Will terminate execution!" << endl;
619  ::abort();
620  }
621  }
622 
623  model = parser.getToken( itoken++ );
624 
625  int photos = 0;
626  int verbose = 0;
627  int summary = 0;
628 
629  do {
630  if ( model == "PHOTOS" ) {
631  photos = 1;
632  model = parser.getToken( itoken++ );
633  }
634  if ( model == "VERBOSE" ) {
635  verbose = 1;
636  model = parser.getToken( itoken++ );
637  }
638  if ( model == "SUMMARY" ) {
639  summary = 1;
640  model = parser.getToken( itoken++ );
641  }
642  } while ( model == "PHOTOS" || model == "VERBOSE" ||
643  model == "SUMMARY" );
644 
645  //see if this is an aliased model
646  int foundAnAlias = -1;
647  for ( size_t iAlias = 0; iAlias < modelAliasList.size();
648  iAlias++ ) {
649  if ( modelAliasList[iAlias].matchAlias( model ) ) {
650  foundAnAlias = iAlias;
651  break;
652  }
653  }
654 
655  if ( foundAnAlias == -1 ) {
656  if ( !modelist.isModel( model ) ) {
657  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
658  << "Expected to find a model name,"
659  << "found:" << model.c_str() << " on line "
660  << parser.getLineofToken( itoken ) << endl;
661  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
662  << "Will terminate execution!" << endl;
663  ::abort();
664  }
665  } else {
666  model = modelAliasList[foundAnAlias].getName();
667  }
668 
669  temp_fcn_new_model = model;
670  temp_fcn_new = modelist.getFcn( model );
671 
672  if ( photos ) {
673  temp_fcn_new->setPHOTOS();
674  }
675  if ( verbose ) {
676  temp_fcn_new->setVerbose();
677  }
678  if ( summary ) {
679  temp_fcn_new->setSummary();
680  }
681 
682  std::vector<std::string> temp_fcn_new_args;
683 
684  std::string name;
685  int ierr;
686 
687  if ( foundAnAlias == -1 ) {
688  do {
689  name = parser.getToken( itoken++ );
690  if ( name != ";" ) {
691  temp_fcn_new_args.push_back(
692  EvtSymTable::get( name, ierr ) );
693  if ( ierr ) {
694  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
695  << "Reading arguments and found:"
696  << name.c_str() << " on line:"
697  << parser.getLineofToken( itoken - 1 )
698  << endl;
699  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
700  << "Will terminate execution!" << endl;
701  ::abort();
702  }
703  }
704  //int isname=EvtPDL::getId(name).getId()>=0;
705  int ismodel = modelist.isModel( name );
706  if ( ismodel ) {
707  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
708  << "Expected ';' but found:" << name.c_str()
709  << " on line:"
710  << parser.getLineofToken( itoken - 1 )
711  << endl;
712  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
713  << "Most probable error is omitted ';'."
714  << endl;
715  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
716  << "Will terminate execution!" << endl;
717  ::abort();
718  }
719  } while ( name != ";" );
720  } else {
721  std::vector<std::string> copyMe =
722  modelAliasList[foundAnAlias].getArgList();
723  temp_fcn_new_args = copyMe;
724  itoken++;
725  }
726  //Found one decay.
727 
728  brfrsum += brfr;
729 
730  temp_fcn_new->saveDecayInfo( ipar, n_daugh, daught,
731  temp_fcn_new_args.size(),
732  temp_fcn_new_args,
733  temp_fcn_new_model, brfr );
734 
735  double massmin = 0.0;
736 
737  // for (i=0;i<n_daugh;i++){
738  for ( i = 0; i < temp_fcn_new->nRealDaughters(); i++ ) {
739  if ( EvtPDL::getMinMass( daught[i] ) > 0.0001 ) {
740  massmin += EvtPDL::getMinMass( daught[i] );
741  } else {
742  massmin += EvtPDL::getMeanMass( daught[i] );
743  }
744  }
745 
746  _decaytable[ipar.getAlias()].addMode( temp_fcn_new, brfrsum,
747  massmin );
748  }
749  } while ( token != "Enddecay" );
750 
751  _decaytable[ipar.getAlias()].finalize();
752 
753  }
754  // Allow copying of decays from one particle to another; useful
755  // in combination with RemoveDecay
756  else if ( token == "CopyDecay" ) {
757  std::string newname;
758  std::string oldname;
759 
760  newname = parser.getToken( itoken++ );
761  oldname = parser.getToken( itoken++ );
762 
763  EvtId newipar = EvtPDL::getId( newname );
764  EvtId oldipar = EvtPDL::getId( oldname );
765 
766  if ( oldipar == EvtId( -1, -1 ) ) {
767  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
768  << "Unknown particle name:" << oldname.c_str()
769  << " on line " << parser.getLineofToken( itoken ) << endl;
770  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
771  << "Will terminate execution!" << endl;
772  ::abort();
773  }
774  if ( newipar == EvtId( -1, -1 ) ) {
775  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
776  << "Unknown particle name:" << newname.c_str()
777  << " on line " << parser.getLineofToken( itoken ) << endl;
778  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
779  << "Will terminate execution!" << endl;
780  ::abort();
781  }
782  if ( _decaytable[newipar.getAlias()].getNMode() != 0 ) {
783  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
784  << "Redefining decay of " << newname << endl;
785  _decaytable[newipar.getAlias()].removeDecay();
786  }
787  _decaytable[newipar.getAlias()] = _decaytable[oldipar.getAlias()];
788  }
789  // Enable decay deletion; intended primarily for aliases
790  // Peter Onyisi, March 2008
791  else if ( token == "RemoveDecay" ) {
792  parent = parser.getToken( itoken++ );
793  ipar = EvtPDL::getId( parent );
794 
795  if ( ipar == EvtId( -1, -1 ) ) {
796  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
797  << "Unknown particle name:" << parent.c_str() << " on line "
798  << parser.getLineofToken( itoken - 1 ) << endl;
799  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
800  << "Will terminate execution!" << endl;
801  ::abort();
802  }
803 
804  if ( _decaytable[ipar.getAlias()].getNMode() == 0 ) {
805  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
806  << "No decays to delete for " << parent.c_str() << endl;
807  } else {
808  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
809  << "Deleting selected decays of " << parent.c_str() << endl;
810  }
811 
812  do {
813  token = parser.getToken( itoken );
814 
815  if ( token != "Enddecay" ) {
816  n_daugh = 0;
817  while ( EvtPDL::getId( parser.getToken( itoken ) ).getId() >=
818  0 ) {
819  sdaug = parser.getToken( itoken++ );
820  daught[n_daugh++] = EvtPDL::getId( sdaug );
821  if ( daught[n_daugh - 1] == EvtId( -1, -1 ) ) {
822  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
823  << "Unknown particle name:" << sdaug.c_str()
824  << " on line "
825  << parser.getLineofToken( itoken ) << endl;
826  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
827  << "Will terminate execution!" << endl;
828  ::abort();
829  }
830  }
831  token = parser.getToken( itoken );
832  if ( token != ";" ) {
833  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
834  << "Expected ';' but found:" << token << " on line:"
835  << parser.getLineofToken( itoken - 1 ) << endl;
836  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
837  << "Most probable error is omitted ';'." << endl;
838  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
839  << "Will terminate execution!" << endl;
840  ::abort();
841  }
842  token = parser.getToken( itoken++ );
843  EvtDecayBase* temp_fcn_new = modelist.getFcn( "PHSP" );
844  std::vector<std::string> temp_fcn_new_args;
845  std::string temp_fcn_new_model( "PHSP" );
846  temp_fcn_new->saveDecayInfo( ipar, n_daugh, daught, 0,
847  temp_fcn_new_args,
848  temp_fcn_new_model, 0. );
849  _decaytable[ipar.getAlias()].removeMode( temp_fcn_new );
850  }
851  } while ( token != "Enddecay" );
852  itoken++;
853  } else if ( token != "End" ) {
854  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
855  << "Found unknown command:'" << token.c_str() << "' on line "
856  << parser.getLineofToken( itoken ) << endl;
857  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
858  << "Will terminate execution!" << endl;
859  ::abort();
860  }
861 
862  } while ( ( token != "End" ) && itoken != parser.getNToken() );
863 
864  //Now we may need to reset the minimum mass for some particles????
865 
866  for ( size_t ii = 0; ii < EvtPDL::entries(); ii++ ) {
867  EvtId temp( ii, ii );
868  int nModTot = getNMode( ii );
869  //no decay modes
870  if ( nModTot == 0 )
871  continue;
872  //0 width?
873  if ( EvtPDL::getWidth( temp ) < 0.0000001 )
874  continue;
875  int jj;
876  double minMass = EvtPDL::getMaxMass( temp );
877  for ( jj = 0; jj < nModTot; jj++ ) {
878  double tmass = _decaytable[ii].getDecay( jj ).getMassMin();
879  if ( tmass < minMass )
880  minMass = tmass;
881  }
882  if ( minMass > EvtPDL::getMinMass( temp ) ) {
883  if ( verbose )
884  EvtGenReport( EVTGEN_INFO, "EvtGen" )
885  << "Given allowed decays, resetting minMass "
886  << EvtPDL::name( temp ).c_str() << " "
887  << EvtPDL::getMinMass( temp ) << " to " << minMass << endl;
888  EvtPDL::reSetMassMin( temp, minMass );
889  }
890  }
891 }
892 
893 void EvtDecayTable::readXMLDecayFile( const std::string dec_name, bool verbose )
894 {
895  if ( _decaytable.size() < EvtPDL::entries() )
896  _decaytable.resize( EvtPDL::entries() );
897  EvtModel& modelist = EvtModel::instance();
898  EvtExtGeneratorCommandsTable* extGenCommands =
900 
901  EvtParserXml parser;
902  parser.open( dec_name );
903 
904  EvtId ipar;
905  std::string decayParent = "";
906  double brfrSum = 0.;
907  std::vector<EvtModelAlias> modelAliasList;
908  bool endReached = false;
909 
910  while ( parser.readNextTag() ) {
911  //TAGS FOUND UNDER DATA
912  if ( parser.getParentTagTitle() == "data" ) {
913  if ( parser.getTagTitle() == "photos" ) {
914  std::string usage = parser.readAttribute( "usage" );
915  if ( usage == "always" ) {
917  if ( verbose )
918  EvtGenReport( EVTGEN_INFO, "EvtGen" )
919  << "As requested, PHOTOS will be turned on for all decays."
920  << endl;
921  } else if ( usage == "never" ) {
923  if ( verbose )
924  EvtGenReport( EVTGEN_INFO, "EvtGen" )
925  << "As requested, PHOTOS will be turned off."
926  << endl;
927  } else {
929  if ( verbose )
930  EvtGenReport( EVTGEN_INFO, "EvtGen" )
931  << "As requested, PHOTOS will be turned on only when requested."
932  << endl;
933  }
934 
935  } else if ( parser.getTagTitle() == "alias" ) {
936  std::string alias = parser.readAttribute( "name" );
937  std::string particle = parser.readAttribute( "particle" );
938  checkParticle( particle );
939  EvtId id = EvtPDL::getId( particle );
940 
941  EvtPDL::alias( id, alias );
942  if ( _decaytable.size() < EvtPDL::entries() )
943  _decaytable.resize( EvtPDL::entries() );
944 
945  } else if ( parser.getTagTitle() == "modelAlias" ) {
946  std::vector<std::string> modelArgList;
947 
948  std::string alias = parser.readAttribute( "name" );
949  std::string model = parser.readAttribute( "model" );
950  std::string paramStr = parser.readAttribute( "params" );
951  std::istringstream paramStream( paramStr );
952 
953  std::string param;
954 
955  if ( paramStr == "" ) {
956  EvtDecayBase* fcn = modelist.getFcn( model );
957  int i( 0 );
958  std::string paramName = fcn->getParamName( 0 );
959  while ( paramName != "" ) {
960  param = parser.readAttribute( paramName,
961  fcn->getParamDefault( i ) );
962  if ( param == "" )
963  break;
964  modelArgList.push_back( param );
965  ++i;
966  paramName = fcn->getParamName( i );
967  }
968  } else {
969  while ( std::getline( paramStream, param, ' ' ) ) {
970  modelArgList.push_back( param );
971  }
972  }
973  EvtModelAlias newAlias( alias, model, modelArgList );
974  modelAliasList.push_back( newAlias );
975 
976  } else if ( parser.getTagTitle() == "chargeConj" ) {
977  std::string particle = parser.readAttribute( "particle" );
978  std::string conjugate = parser.readAttribute( "conjugate" );
979 
980  EvtId a = EvtPDL::getId( particle );
981  EvtId abar = EvtPDL::getId( conjugate );
982 
983  checkParticle( particle );
984  checkParticle( conjugate );
985 
986  EvtPDL::aliasChgConj( a, abar );
987 
988  } else if ( parser.getTagTitle() == "conjDecay" ) {
989  std::string particle = parser.readAttribute( "particle" );
990 
991  EvtId a = EvtPDL::getId( particle );
992  EvtId abar = EvtPDL::chargeConj( a );
993 
994  checkParticle( particle );
995  checkParticle( abar.getName() );
996 
997  if ( _decaytable[a.getAlias()].getNMode() != 0 ) {
998  if ( verbose )
999  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1000  << "Redefined decay of " << particle.c_str()
1001  << " in ConjDecay" << endl;
1002 
1003  _decaytable[a.getAlias()].removeDecay();
1004  }
1005 
1006  //take contents of abar and conjugate and store in a
1007  _decaytable[a.getAlias()].makeChargeConj(
1008  &_decaytable[abar.getAlias()] );
1009 
1010  } else if ( parser.getTagTitle() == "define" ) {
1011  std::string name = parser.readAttribute( "name" );
1012  std::string value = parser.readAttribute( "value" );
1013  EvtSymTable::define( name, value );
1014 
1015  } else if ( parser.getTagTitle() == "particle" ) {
1016  std::string name = parser.readAttribute( "name" );
1017  double mass = parser.readAttributeDouble( "mass" );
1018  double width = parser.readAttributeDouble( "width" );
1019  double minMass = parser.readAttributeDouble( "massMin" );
1020  double maxMass = parser.readAttributeDouble( "massMax" );
1021  std::string birthFactor = parser.readAttribute(
1022  "includeBirthFactor" );
1023  std::string decayFactor = parser.readAttribute(
1024  "includeDecayFactor" );
1025  std::string lineShape = parser.readAttribute( "lineShape" );
1026  double blattWeisskopfD = parser.readAttributeDouble(
1027  "blattWeisskopfFactor" );
1028  double blattWeisskopfB = parser.readAttributeDouble(
1029  "blattWeisskopfBirth" );
1030 
1031  EvtId thisPart = EvtPDL::getId( name );
1032  checkParticle( name );
1033 
1034  if ( mass != -1 ) {
1035  EvtPDL::reSetMass( thisPart, mass );
1036  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1037  << "Refined mass for "
1038  << EvtPDL::name( thisPart ).c_str() << " to be " << mass
1039  << endl;
1040  }
1041  if ( width != -1 ) {
1042  EvtPDL::reSetWidth( thisPart, width );
1043  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1044  << "Refined width for "
1045  << EvtPDL::name( thisPart ).c_str() << " to be "
1046  << width << endl;
1047  }
1048  if ( minMass != -1 ) {
1049  EvtPDL::reSetMassMin( thisPart, minMass );
1050  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1051  << "Refined minimum mass for "
1052  << EvtPDL::name( thisPart ).c_str() << " to be "
1053  << minMass << endl;
1054  }
1055  if ( maxMass != -1 ) {
1056  EvtPDL::reSetMassMax( thisPart, maxMass );
1057  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1058  << "Refined maximum mass for "
1059  << EvtPDL::name( thisPart ).c_str() << " to be "
1060  << maxMass << endl;
1061  }
1062  if ( !birthFactor.empty() ) {
1063  EvtPDL::includeBirthFactor( thisPart,
1064  stringToBoolean( birthFactor ) );
1065  if ( verbose ) {
1066  if ( stringToBoolean( birthFactor ) ) {
1067  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1068  << "Include birth factor for "
1069  << EvtPDL::name( thisPart ).c_str() << endl;
1070  } else {
1071  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1072  << "No longer include birth factor for "
1073  << EvtPDL::name( thisPart ).c_str() << endl;
1074  }
1075  }
1076  }
1077  if ( !decayFactor.empty() ) {
1078  EvtPDL::includeDecayFactor( thisPart,
1079  stringToBoolean( decayFactor ) );
1080  if ( verbose ) {
1081  if ( stringToBoolean( decayFactor ) ) {
1082  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1083  << "Include decay factor for "
1084  << EvtPDL::name( thisPart ).c_str() << endl;
1085  } else {
1086  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1087  << "No longer include decay factor for "
1088  << EvtPDL::name( thisPart ).c_str() << endl;
1089  }
1090  }
1091  }
1092  if ( !lineShape.empty() ) {
1093  EvtPDL::changeLS( thisPart, lineShape );
1094  if ( verbose )
1095  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1096  << "Change lineshape to " << lineShape << " for "
1097  << EvtPDL::name( thisPart ).c_str() << endl;
1098  }
1099  if ( blattWeisskopfD != -1 ) {
1100  EvtPDL::reSetBlatt( thisPart, blattWeisskopfD );
1101  if ( verbose )
1102  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1103  << "Redefined Blatt-Weisskopf factor "
1104  << EvtPDL::name( thisPart ).c_str() << " to be "
1105  << blattWeisskopfD << endl;
1106  }
1107  if ( blattWeisskopfB != -1 ) {
1108  EvtPDL::reSetBlattBirth( thisPart, blattWeisskopfB );
1109  if ( verbose )
1110  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1111  << "Redefined Blatt-Weisskopf birth factor "
1112  << EvtPDL::name( thisPart ).c_str() << " to be "
1113  << blattWeisskopfB << endl;
1114  }
1115  } else if ( parser.getTagTitle() == "lineShapePW" ) {
1116  std::string parent = parser.readAttribute( "parent" );
1117  std::string daug1 = parser.readAttribute( "daug1" );
1118  std::string daug2 = parser.readAttribute( "daug2" );
1119  int pw = parser.readAttributeInt( "pw" );
1120 
1121  checkParticle( parent );
1122  checkParticle( daug1 );
1123  checkParticle( daug2 );
1124 
1125  EvtId thisPart = EvtPDL::getId( parent );
1126  EvtId thisD1 = EvtPDL::getId( daug1 );
1127  EvtId thisD2 = EvtPDL::getId( daug2 );
1128 
1129  EvtPDL::setPWForDecay( thisPart, pw, thisD1, thisD2 );
1130  EvtPDL::setPWForBirthL( thisD1, pw, thisPart, thisD2 );
1131  EvtPDL::setPWForBirthL( thisD2, pw, thisPart, thisD1 );
1132  if ( verbose )
1133  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1134  << "Redefined Partial wave for " << parent.c_str()
1135  << " to " << daug1.c_str() << " " << daug2.c_str()
1136  << " (" << pw << ")" << endl;
1137 
1138  } else if ( parser.getTagTitle() == "decay" ) { //start of a particle
1139  brfrSum = 0.;
1140  decayParent = parser.readAttribute( "name" );
1141  checkParticle( decayParent );
1142  ipar = EvtPDL::getId( decayParent );
1143 
1144  if ( _decaytable[ipar.getAlias()].getNMode() != 0 ) {
1145  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1146  << "Redefined decay of " << decayParent.c_str() << endl;
1147  _decaytable[ipar.getAlias()].removeDecay();
1148  }
1149 
1150  } else if ( parser.getTagTitle() == "copyDecay" ) {
1151  std::string particle = parser.readAttribute( "particle" );
1152  std::string copy = parser.readAttribute( "copy" );
1153 
1154  EvtId newipar = EvtPDL::getId( particle );
1155  EvtId oldipar = EvtPDL::getId( copy );
1156 
1157  checkParticle( particle );
1158  checkParticle( copy );
1159 
1160  if ( _decaytable[newipar.getAlias()].getNMode() != 0 ) {
1161  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1162  << "Redefining decay of " << particle << endl;
1163  _decaytable[newipar.getAlias()].removeDecay();
1164  }
1165  _decaytable[newipar.getAlias()] = _decaytable[oldipar.getAlias()];
1166 
1167  } else if ( parser.getTagTitle() == "removeDecay" ) {
1168  decayParent = parser.readAttribute( "particle" );
1169  checkParticle( decayParent );
1170  ipar = EvtPDL::getId( decayParent );
1171 
1172  if ( _decaytable[ipar.getAlias()].getNMode() == 0 ) {
1173  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1174  << "No decays to delete for " << decayParent.c_str()
1175  << endl;
1176  } else {
1177  EvtGenReport( EVTGEN_DEBUG, "EvtGen" )
1178  << "Deleting selected decays of " << decayParent.c_str()
1179  << endl;
1180  }
1181 
1182  } else if ( parser.getTagTitle() == "pythiaParam" ) {
1183  Command command;
1184  command["GENERATOR"] = parser.readAttribute( "generator" );
1185  command["MODULE"] = parser.readAttribute( "module" );
1186  command["PARAM"] = parser.readAttribute( "param" );
1187  command["VALUE"] = parser.readAttribute( "value" );
1188  command["VERSION"] = "PYTHIA8";
1189  extGenCommands->addCommand( "PYTHIA", command );
1190 
1191  } else if ( parser.getTagTitle() == "pythia6Param" ) {
1192  Command command;
1193  command["GENERATOR"] = parser.readAttribute( "generator" );
1194  command["MODULE"] = parser.readAttribute( "module" );
1195  command["PARAM"] = parser.readAttribute( "param" );
1196  command["VALUE"] = parser.readAttribute( "value" );
1197  command["VERSION"] = "PYTHIA6";
1198  extGenCommands->addCommand( "PYTHIA", command );
1199 
1200  } else if ( parser.getTagTitle() == "/data" ) { //end of data
1201  endReached = true;
1202  parser.close();
1203  break;
1204  } else if ( parser.getTagTitle() == "Title" ||
1205  parser.getTagTitle() == "Details" ||
1206  parser.getTagTitle() == "Author" ||
1207  parser.getTagTitle() == "Version"
1208  //the above tags are expected to be in the XML decay file but are not used by EvtGen
1209  || parser.getTagTitle() == "dalitzDecay" ||
1210  parser.getTagTitle() == "copyDalitz" ) {
1211  //the above tags are only used by EvtGenModels/EvtDalitzTable
1212  } else {
1213  EvtGenReport( EVTGEN_INFO, "EvtGen" )
1214  << "Unknown tag " << parser.getTagTitle()
1215  << " found in XML decay file near line "
1216  << parser.getLineNumber() << ". Tag will be ignored."
1217  << endl;
1218  }
1219  //TAGS FOUND UNDER DECAY
1220  } else if ( parser.getParentTagTitle() == "decay" ) {
1221  if ( parser.getTagTitle() == "channel" ) { //start of a channel
1222  int nDaughters = 0;
1223  EvtId daughter[MAX_DAUG];
1224 
1225  EvtDecayBase* temp_fcn_new;
1226  std::string temp_fcn_new_model;
1227  std::vector<std::string> temp_fcn_new_args;
1228 
1229  double brfr = parser.readAttributeDouble( "br" );
1230  std::string daugStr = parser.readAttribute( "daughters" );
1231  std::istringstream daugStream( daugStr );
1232  std::string model = parser.readAttribute( "model" );
1233  std::string paramStr = parser.readAttribute( "params" );
1234  std::istringstream paramStream( paramStr );
1235  bool decVerbose = parser.readAttributeBool( "verbose" );
1236  bool decPhotos = parser.readAttributeBool( "photos" );
1237  bool decSummary = parser.readAttributeBool( "summary" );
1238 
1239  std::string daugh;
1240  while ( std::getline( daugStream, daugh, ' ' ) ) {
1241  checkParticle( daugh );
1242  daughter[nDaughters++] = EvtPDL::getId( daugh );
1243  }
1244 
1245  int modelAlias = -1;
1246  for ( size_t iAlias = 0; iAlias < modelAliasList.size();
1247  iAlias++ ) {
1248  if ( modelAliasList[iAlias].matchAlias( model ) ) {
1249  modelAlias = iAlias;
1250  break;
1251  }
1252  }
1253 
1254  if ( modelAlias == -1 ) {
1255  if ( !modelist.isModel( model ) ) {
1256  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1257  << "Expected to find a model name near line "
1258  << parser.getLineNumber() << ","
1259  << "found:" << model.c_str() << endl;
1260  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1261  << "Will terminate execution!" << endl;
1262  ::abort();
1263  }
1264  } else {
1265  model = modelAliasList[modelAlias].getName();
1266  }
1267 
1268  temp_fcn_new_model = model;
1269  temp_fcn_new = modelist.getFcn( model );
1270 
1271  if ( decPhotos )
1272  temp_fcn_new->setPHOTOS();
1273  if ( decVerbose )
1274  temp_fcn_new->setVerbose();
1275  if ( decSummary )
1276  temp_fcn_new->setSummary();
1277 
1278  int ierr;
1279  if ( modelAlias == -1 ) {
1280  std::string param;
1281  if ( paramStr == "" ) {
1282  int i( 0 );
1283  std::string paramName = temp_fcn_new->getParamName( 0 );
1284  while ( paramName != "" ) {
1285  param = parser.readAttribute(
1286  paramName, temp_fcn_new->getParamDefault( i ) );
1287  if ( param == "" )
1288  break; //params must be added in order so we can't just skip the missing ones
1289  temp_fcn_new_args.push_back(
1290  EvtSymTable::get( param, ierr ) );
1291  if ( ierr ) {
1292  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1293  << "Reading arguments near line "
1294  << parser.getLineNumber()
1295  << " and found:" << param.c_str() << endl;
1296  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1297  << "Will terminate execution!" << endl;
1298  ::abort();
1299  }
1300  ++i;
1301  paramName = temp_fcn_new->getParamName( i );
1302  }
1303 
1304  } else { //if the params are not set seperately
1305  while ( std::getline( paramStream, param, ' ' ) ) {
1306  temp_fcn_new_args.push_back(
1307  EvtSymTable::get( param, ierr ) );
1308  if ( ierr ) {
1309  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1310  << "Reading arguments near line "
1311  << parser.getLineNumber()
1312  << " and found:" << param.c_str() << endl;
1313  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1314  << "Will terminate execution!" << endl;
1315  ::abort();
1316  }
1317  }
1318  }
1319  } else {
1320  std::vector<std::string> copyMe =
1321  modelAliasList[modelAlias].getArgList();
1322  temp_fcn_new_args = copyMe;
1323  }
1324 
1325  brfrSum += brfr;
1326 
1327  temp_fcn_new->saveDecayInfo( ipar, nDaughters, daughter,
1328  temp_fcn_new_args.size(),
1329  temp_fcn_new_args,
1330  temp_fcn_new_model, brfr );
1331 
1332  double massMin = 0.0;
1333 
1334  for ( int i = 0; i < temp_fcn_new->nRealDaughters(); i++ ) {
1335  if ( EvtPDL::getMinMass( daughter[i] ) > 0.0001 ) {
1336  massMin += EvtPDL::getMinMass( daughter[i] );
1337  } else {
1338  massMin += EvtPDL::getMeanMass( daughter[i] );
1339  }
1340  }
1341 
1342  _decaytable[ipar.getAlias()].addMode( temp_fcn_new, brfrSum,
1343  massMin );
1344 
1345  } else if ( parser.getTagTitle() == "/decay" ) { //end of a particle
1346  _decaytable[ipar.getAlias()].finalize();
1347  } else
1348  EvtGenReport( EVTGEN_INFO, "EvtGen" )
1349  << "Unexpected tag " << parser.getTagTitle()
1350  << " found in XML decay file near line "
1351  << parser.getLineNumber() << ". Tag will be ignored."
1352  << endl;
1353  //TAGS FOUND UNDER REMOVEDECAY
1354  } else if ( parser.getParentTagTitle() == "removeDecay" ) {
1355  if ( parser.getTagTitle() == "channel" ) { //start of a channel
1356  int nDaughters = 0;
1357  EvtId daughter[MAX_DAUG];
1358 
1359  std::string daugStr = parser.readAttribute( "daughters" );
1360  std::istringstream daugStream( daugStr );
1361 
1362  std::string daugh;
1363  while ( std::getline( daugStream, daugh, ' ' ) ) {
1364  checkParticle( daugh );
1365  daughter[nDaughters++] = EvtPDL::getId( daugh );
1366  }
1367 
1368  EvtDecayBase* temp_fcn_new = modelist.getFcn( "PHSP" );
1369  std::vector<std::string> temp_fcn_new_args;
1370  std::string temp_fcn_new_model( "PHSP" );
1371  temp_fcn_new->saveDecayInfo( ipar, nDaughters, daughter, 0,
1372  temp_fcn_new_args,
1373  temp_fcn_new_model, 0. );
1374  _decaytable[ipar.getAlias()].removeMode( temp_fcn_new );
1375  } else if ( parser.getTagTitle() != "/removeDecay" ) {
1376  EvtGenReport( EVTGEN_INFO, "EvtGen" )
1377  << "Unexpected tag " << parser.getTagTitle()
1378  << " found in XML decay file near line "
1379  << parser.getLineNumber() << ". Tag will be ignored."
1380  << endl;
1381  }
1382  }
1383  } //while lines in file
1384 
1385  if ( !endReached ) {
1386  EvtGenReport( EVTGEN_INFO, "EvtGen" )
1387  << "Either the decay file ended prematurely or the file is badly formed.\n"
1388  << "Error occured near line" << parser.getLineNumber() << endl;
1389  ::abort();
1390  }
1391 
1392  //Now we may need to reset the minimum mass for some particles????
1393  for ( size_t ii = 0; ii < EvtPDL::entries(); ii++ ) {
1394  EvtId temp( ii, ii );
1395  int nModTot = getNMode( ii );
1396  //no decay modes
1397  if ( nModTot == 0 )
1398  continue;
1399  //0 width?
1400  if ( EvtPDL::getWidth( temp ) < 0.0000001 )
1401  continue;
1402  int jj;
1403  double minMass = EvtPDL::getMaxMass( temp );
1404  for ( jj = 0; jj < nModTot; jj++ ) {
1405  double tmass = _decaytable[ii].getDecay( jj ).getMassMin();
1406  if ( tmass < minMass )
1407  minMass = tmass;
1408  }
1409  if ( minMass > EvtPDL::getMinMass( temp ) ) {
1410  if ( verbose )
1411  EvtGenReport( EVTGEN_INFO, "EvtGen" )
1412  << "Given allowed decays, resetting minMass "
1413  << EvtPDL::name( temp ).c_str() << " "
1414  << EvtPDL::getMinMass( temp ) << " to " << minMass << endl;
1415  EvtPDL::reSetMassMin( temp, minMass );
1416  }
1417  }
1418 }
1419 
1420 bool EvtDecayTable::stringToBoolean( std::string valStr )
1421 {
1422  return ( valStr == "true" || valStr == "1" || valStr == "on" ||
1423  valStr == "yes" );
1424 }
1425 
1426 void EvtDecayTable::checkParticle( std::string particle )
1427 {
1428  if ( EvtPDL::getId( particle ) == EvtId( -1, -1 ) ) {
1429  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1430  << "Unknown particle name:" << particle.c_str() << endl;
1431  EvtGenReport( EVTGEN_ERROR, "EvtGen" )
1432  << "Will terminate execution!" << endl;
1433  ::abort();
1434  }
1435 }
1436 
1438 {
1439  int aliasInt = id.getAlias();
1440 
1441  EvtDecayBase* theModel = this->findDecayModel( aliasInt, modeInt );
1442 
1443  return theModel;
1444 }
1445 
1446 EvtDecayBase* EvtDecayTable::findDecayModel( int aliasInt, int modeInt )
1447 {
1448  EvtDecayBase* theModel( 0 );
1449 
1450  if ( aliasInt >= 0 && aliasInt < (int)EvtPDL::entries() ) {
1451  theModel = _decaytable[aliasInt].getDecayModel( modeInt );
1452  }
1453 
1454  return theModel;
1455 }
1456 
1458 {
1459  bool hasPythia = this->hasPythia( id.getAlias() );
1460  return hasPythia;
1461 }
1462 
1463 bool EvtDecayTable::hasPythia( int aliasInt )
1464 {
1465  bool hasPythia( false );
1466  if ( aliasInt >= 0 && aliasInt < (int)EvtPDL::entries() ) {
1467  hasPythia = _decaytable[aliasInt].isJetSet();
1468  }
1469 
1470  return hasPythia;
1471 }
1472 
1474 {
1475  int nModes = this->getNModes( id.getAlias() );
1476  return nModes;
1477 }
1478 
1479 int EvtDecayTable::getNModes( int aliasInt )
1480 {
1481  int nModes( 0 );
1482 
1483  if ( aliasInt >= 0 && aliasInt < (int)EvtPDL::entries() ) {
1484  nModes = _decaytable[aliasInt].getNMode();
1485  }
1486 
1487  return nModes;
1488 }
1489 
1490 int EvtDecayTable::findChannel( EvtId parent, std::string model, int ndaug,
1491  EvtId* daugs, int narg, std::string* args )
1492 {
1493  int i, j, right;
1494  EvtId daugs_scratch[50];
1495  int nmatch, k;
1496 
1497  for ( i = 0; i < _decaytable[parent.getAlias()].getNMode(); i++ ) {
1498  right = 1;
1499 
1500  right = right && model == _decaytable[parent.getAlias()]
1501  .getDecay( i )
1502  .getDecayModel()
1503  ->getModelName();
1504  right = right && ( ndaug == _decaytable[parent.getAlias()]
1505  .getDecay( i )
1506  .getDecayModel()
1507  ->getNDaug() );
1508  right = right && ( narg == _decaytable[parent.getAlias()]
1509  .getDecay( i )
1510  .getDecayModel()
1511  ->getNArg() );
1512 
1513  if ( right ) {
1514  for ( j = 0; j < ndaug; j++ ) {
1515  daugs_scratch[j] = daugs[j];
1516  }
1517 
1518  nmatch = 0;
1519 
1520  for ( j = 0; j < _decaytable[parent.getAlias()]
1521  .getDecay( i )
1522  .getDecayModel()
1523  ->getNDaug();
1524  j++ ) {
1525  for ( k = 0; k < ndaug; k++ ) {
1526  if ( daugs_scratch[k] == _decaytable[parent.getAlias()]
1527  .getDecay( i )
1528  .getDecayModel()
1529  ->getDaug( j ) ) {
1530  daugs_scratch[k] = EvtId( -1, -1 );
1531  nmatch++;
1532  break;
1533  }
1534  }
1535  }
1536 
1537  right = right && ( nmatch == ndaug );
1538 
1539  for ( j = 0; j < _decaytable[parent.getAlias()]
1540  .getDecay( i )
1541  .getDecayModel()
1542  ->getNArg();
1543  j++ ) {
1544  right = right && ( args[j] == _decaytable[parent.getAlias()]
1545  .getDecay( i )
1546  .getDecayModel()
1547  ->getArgStr( j ) );
1548  }
1549  }
1550  if ( right )
1551  return i;
1552  }
1553  return -1;
1554 }
1555 
1556 int EvtDecayTable::inChannelList( EvtId parent, int ndaug, EvtId* daugs )
1557 {
1558  int i, j, k;
1559  EvtId daugs_scratch[MAX_DAUG];
1560 
1561  int dsum = 0;
1562  for ( i = 0; i < ndaug; i++ ) {
1563  dsum += daugs[i].getAlias();
1564  }
1565 
1566  int nmatch;
1567 
1568  int ipar = parent.getAlias();
1569 
1570  int nmode = _decaytable[ipar].getNMode();
1571 
1572  for ( i = 0; i < nmode; i++ ) {
1573  EvtDecayBase* thedecaymodel =
1574  _decaytable[ipar].getDecay( i ).getDecayModel();
1575 
1576  if ( thedecaymodel->getDSum() == dsum ) {
1577  int nd = thedecaymodel->getNDaug();
1578 
1579  if ( ndaug == nd ) {
1580  for ( j = 0; j < ndaug; j++ ) {
1581  daugs_scratch[j] = daugs[j];
1582  }
1583  nmatch = 0;
1584  for ( j = 0; j < nd; j++ ) {
1585  for ( k = 0; k < ndaug; k++ ) {
1586  if ( EvtId( daugs_scratch[k] ) ==
1587  thedecaymodel->getDaug( j ) ) {
1588  daugs_scratch[k] = EvtId( -1, -1 );
1589  nmatch++;
1590  break;
1591  }
1592  }
1593  }
1594  if ( ( nmatch == ndaug ) &&
1595  ( !( ( thedecaymodel->getModelName() == "JETSET" ) ||
1596  ( thedecaymodel->getModelName() == "PYTHIA" ) ) ) ) {
1597  return i;
1598  }
1599  }
1600  }
1601  }
1602 
1603  return -1;
1604 }
1605 
1606 std::vector<std::string> EvtDecayTable::splitString( std::string& theString,
1607  std::string& splitter )
1608 {
1609  // Code from STLplus
1610  std::vector<std::string> result;
1611 
1612  if ( !theString.empty() && !splitter.empty() ) {
1613  for ( std::string::size_type offset = 0;; ) {
1614  std::string::size_type found = theString.find( splitter, offset );
1615 
1616  if ( found != std::string::npos ) {
1617  std::string tmpString = theString.substr( offset, found - offset );
1618  if ( tmpString.size() > 0 ) {
1619  result.push_back( tmpString );
1620  }
1621  offset = found + splitter.size();
1622  } else {
1623  std::string tmpString =
1624  theString.substr( offset, theString.size() - offset );
1625  if ( tmpString.size() > 0 ) {
1626  result.push_back( tmpString );
1627  }
1628  break;
1629  }
1630  }
1631  }
1632 
1633  return result;
1634 }
static EvtExtGeneratorCommandsTable * getInstance()
std::vector< std::string > splitString(std::string &theString, std::string &splitter)
static void setAlwaysRadCorr()
Definition: EvtRadCorr.cpp:77
static void reSetMass(EvtId i, double mass)
Definition: EvtPDL.cpp:397
static std::string name(EvtId i)
Definition: EvtPDL.cpp:382
int getLineNumber()
Definition: EvtParserXml.hh:37
int getNMode(int ipar)
static void define(const std::string &name, std::string d)
Definition: EvtSymTable.cpp:40
virtual int nRealDaughters()
static void reSetMassMin(EvtId i, double mass)
Definition: EvtPDL.cpp:407
static EvtDecayTable * getInstance()
std::string getName() const
Definition: EvtId.cpp:41
virtual std::string getParamDefault(int i)
void setPHOTOS()
Definition: EvtDecayBase.hh:70
static void reSetWidth(EvtId i, double width)
Definition: EvtPDL.cpp:402
virtual std::string getParamName(int i)
EvtDecayBase * getFcn(std::string model_name)
Definition: EvtModel.cpp:47
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
double readAttributeDouble(std::string attribute, double defaultValue=-1.)
int read(const std::string filename)
Definition: EvtParser.cpp:62
static void setPWForDecay(EvtId i, int spin, EvtId d1, EvtId d2)
Definition: EvtPDL.cpp:442
void addCommand(std::string extGenerator, Command command)
EvtDecayBase * getDecay(int ipar, int imode)
static void reSetBlatt(EvtId i, double blatt)
Definition: EvtPDL.cpp:417
static void includeDecayFactor(EvtId i, bool yesno)
Definition: EvtPDL.cpp:432
static double getMeanMass(EvtId i)
Definition: EvtPDL.cpp:314
EvtId getId() const
bool open(std::string filename)
std::string getTagTitle()
Definition: EvtParserXml.hh:35
static void alias(EvtId num, const std::string &newname)
Definition: EvtPDL.cpp:256
bool readAttributeBool(std::string attribute, bool defaultValue=false)
static std::string get(const std::string &name, int &ierr)
Definition: EvtSymTable.cpp:55
static void reSetBlattBirth(EvtId i, double blatt)
Definition: EvtPDL.cpp:422
Definition: EvtId.hh:27
bool stringToBoolean(std::string valStr)
EvtDecayBase * getDecayFunc(EvtParticle *p)
int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
bool readNextTag()
int getDSum() const
Definition: EvtDecayBase.hh:80
std::string getParentTagTitle()
void saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
int getNModes(int aliasInt)
static double getMinMass(EvtId i)
Definition: EvtPDL.cpp:342
EvtDecayBase * findDecayModel(int aliasInt, int modeInt)
static size_t entries()
Definition: EvtPDL.cpp:387
std::map< std::string, std::string > Command
static EvtModel & instance()
Definition: EvtModel.hh:56
static void changeLS(EvtId i, std::string &newLS)
Definition: EvtPDL.cpp:437
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
std::vector< EvtParticleDecayList > _decaytable
int findChannel(EvtId parent, std::string model, int ndaug, EvtId *daugs, int narg, std::string *args)
std::string getModelName() const
Definition: EvtDecayBase.hh:79
static void setNormalRadCorr()
Definition: EvtRadCorr.cpp:87
static double getWidth(EvtId i)
Definition: EvtPDL.cpp:352
void readXMLDecayFile(const std::string dec_name, bool verbose=true)
int readAttributeInt(std::string attribute, int defaultValue=-1)
int getNDaug() const
Definition: EvtDecayBase.hh:65
static void aliasChgConj(EvtId a, EvtId abar)
Definition: EvtPDL.cpp:191
static void reSetMassMax(EvtId i, double mass)
Definition: EvtPDL.cpp:412
const int MAX_DAUG
Definition: EvtParticle.hh:42
void readDecayFile(const std::string dec_name, bool verbose=true)
void storeCommand(std::string cmd, std::string cnfgstr)
Definition: EvtModel.cpp:92
static void setPWForBirthL(EvtId i, int spin, EvtId par, EvtId othD)
Definition: EvtPDL.cpp:447
void setVerbose()
Definition: EvtDecayBase.hh:71
int getAlias() const
Definition: EvtId.hh:44
int isCommand(std::string cmd)
Definition: EvtModel.cpp:84
int getLineofToken(int i)
Definition: EvtParser.cpp:57
void setSummary()
Definition: EvtDecayBase.hh:72
int getNToken()
Definition: EvtParser.cpp:47
void checkParticle(std::string particle)
static EvtId chargeConj(EvtId id)
Definition: EvtPDL.cpp:207
int isModel(std::string name)
Definition: EvtModel.cpp:76
bool hasPythia(int aliasInt)
static double getMaxMass(EvtId i)
Definition: EvtPDL.cpp:337
static void includeBirthFactor(EvtId i, bool yesno)
Definition: EvtPDL.cpp:427
static void setNeverRadCorr()
Definition: EvtRadCorr.cpp:82
const std::string & getToken(int i)
Definition: EvtParser.cpp:52
std::string readAttribute(std::string attribute, std::string defaultValue="")
EvtId getDaug(int i) const
Definition: EvtDecayBase.hh:67