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