BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDecayTable.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtDecayTable.cc
12//
13// Description:
14//
15// Modification history:
16//
17// DJL/RYD September 25, 1996 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtPatches.hh"
22
23#include "EvtDecayBase.hh"
24#include "EvtDecayParm.hh"
25#include "EvtDecayTable.hh"
26#include "EvtModel.hh"
27#include "EvtModelAlias.hh"
28#include "EvtPDL.hh"
29#include "EvtParser.hh"
30#include "EvtParticle.hh"
31#include "EvtParticleDecayList.hh"
32#include "EvtRadCorr.hh"
33#include "EvtRandom.hh"
34#include "EvtReport.hh"
35#include "EvtSymTable.hh"
36#include <ctype.h>
37#include <fstream>
38#include <iomanip>
39#include <iostream>
40#include <stdlib.h>
41#include <string.h>
42#include <vector>
43using std::endl;
44using std::fstream;
45using std::ifstream;
46
47static std::vector<EvtParticleDecayList> decaytable;
48
49int EvtDecayTable::getNMode( int ipar ) { return decaytable[ipar].getNMode(); }
50
51EvtDecayBase* getDecay( int ipar, int imode ) {
52 return decaytable[ipar].getDecayModel( imode );
53}
54
55EvtDecayBase* EvtDecayTable::gettheDecay( EvtId parent, int imode ) {
56 int ipar = parent.getAlias();
57 EvtDecayBase* thedecaymodel = decaytable[ipar].getDecay( imode ).getDecayModel();
58 return thedecaymodel;
59}
60
61void EvtDecayTable::printSummary() {
62
63 int i;
64
65 for ( i = 0; i < EvtPDL::entries(); i++ ) { decaytable[i].printSummary(); }
66}
67
68EvtDecayBase* EvtDecayTable::getDecayFunc( EvtParticle* p ) {
69 int partnum;
70
71 partnum = p->getId().getAlias();
72
73 if ( decaytable[partnum].getNMode() == 0 ) return 0;
74 return decaytable[partnum].getDecayModel( p );
75}
76
77void EvtDecayTable::readDecayFile( const std::string dec_name ) {
78
79 if ( decaytable.size() < EvtPDL::entries() ) decaytable.resize( EvtPDL::entries() );
80 EvtModel& modelist = EvtModel::instance();
81 int i;
82
83 report( INFO, "EvtGen" ) << "In readDecayFile, reading:" << dec_name.c_str() << endl;
84
85 ifstream fin;
86
87 fin.open( dec_name.c_str() );
88 if ( !fin ) { report( ERROR, "EvtGen" ) << "Could not open " << dec_name.c_str() << endl; }
89 fin.close();
90
91 EvtParser parser;
92 parser.Read( dec_name );
93
94 int itok;
95
96 int hasend = 0;
97
98 std::string token;
99
100 for ( itok = 0; itok < parser.getNToken(); itok++ )
101 {
102
103 token = parser.getToken( itok );
104
105 if ( token == "End" ) hasend = 1;
106 }
107
108 if ( !hasend )
109 {
110 report( ERROR, "EvtGen" ) << "Could not find an 'End' in " << dec_name.c_str() << endl;
111 report( ERROR, "EvtGen" ) << "Will terminate execution." << endl;
112 ::abort();
113 }
114
115 std::string model, parent, sdaug;
116
117 EvtId ipar;
118
119 int n_daugh;
120 EvtId daught[MAX_DAUG];
121 double brfr;
122
123 int itoken = 0;
124
125 std::vector<EvtModelAlias> modelAliasList;
126
127 do {
128
129 token = parser.getToken( itoken++ );
130
131 // Easy way to turn off photos... Lange September 5, 2000
132 if ( token == "noPhotos" )
133 {
134 EvtRadCorr::setNeverRadCorr();
135 report( INFO, "EvtGen" ) << "As requested, PHOTOS will be turned off." << endl;
136 }
137 else if ( token == "yesPhotos" )
138 {
139 EvtRadCorr::setAlwaysRadCorr();
140 report( INFO, "EvtGen" ) << "As requested, PHOTOS will be turned on." << endl;
141 }
142 else if ( token == "normalPhotos" )
143 {
144 EvtRadCorr::setNormalRadCorr();
145 report( INFO, "EvtGen" ) << "As requested, PHOTOS will be turned on only when requested."
146 << endl;
147 }
148 else if ( token == "Alias" )
149 {
150
151 std::string newname;
152 std::string oldname;
153
154 newname = parser.getToken( itoken++ );
155 oldname = parser.getToken( itoken++ );
156
157 EvtId id = EvtPDL::getId( oldname );
158
159 if ( id == EvtId( -1, -1 ) )
160 {
161 report( ERROR, "EvtGen" ) << "Unknown particle name:" << oldname.c_str() << " on line "
162 << parser.getLineofToken( itoken ) << endl;
163 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
164 ::abort();
165 }
166
167 EvtPDL::alias( id, newname );
168 if ( decaytable.size() < EvtPDL::entries() ) decaytable.resize( EvtPDL::entries() );
169 }
170 else if ( token == "ModelAlias" )
171 {
172 std::vector<std::string> modelArgList;
173
174 std::string aliasName = parser.getToken( itoken++ );
175 std::string modelName = parser.getToken( itoken++ );
176
177 std::string nameTemp;
178 do {
179 nameTemp = parser.getToken( itoken++ );
180 if ( nameTemp != ";" ) { modelArgList.push_back( nameTemp ); }
181 } while ( nameTemp != ";" );
182 EvtModelAlias newAlias( aliasName, modelName, modelArgList );
183 modelAliasList.push_back( newAlias );
184 }
185 else if ( token == "ChargeConj" )
186 {
187
188 std::string aname;
189 std::string abarname;
190
191 aname = parser.getToken( itoken++ );
192 abarname = parser.getToken( itoken++ );
193
194 EvtId a = EvtPDL::getId( aname );
195 EvtId abar = EvtPDL::getId( abarname );
196
197 if ( a == EvtId( -1, -1 ) )
198 {
199 report( ERROR, "EvtGen" ) << "Unknown particle name:" << aname.c_str() << " on line "
200 << parser.getLineofToken( itoken ) << endl;
201 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
202 ::abort();
203 }
204
205 if ( abar == EvtId( -1, -1 ) )
206 {
207 report( ERROR, "EvtGen" ) << "Unknown particle name:" << abarname.c_str()
208 << " on line " << parser.getLineofToken( itoken ) << endl;
209 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
210 ::abort();
211 }
212
213 EvtPDL::aliasChgConj( a, abar );
214 }
215 else if ( modelist.isCommand( token ) )
216 {
217
218 std::string cnfgstr;
219
220 cnfgstr = parser.getToken( itoken++ );
221
222 modelist.storeCommand( token, cnfgstr );
223 }
224 else if ( token == "CDecay" )
225 {
226
227 std::string name;
228
229 name = parser.getToken( itoken++ );
230 ipar = EvtPDL::getId( name );
231
232 if ( ipar == EvtId( -1, -1 ) )
233 {
234 report( ERROR, "EvtGen" ) << "Unknown particle name:" << name.c_str() << " on line "
235 << parser.getLineofToken( itoken - 1 ) << endl;
236 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
237 ::abort();
238 }
239
240 EvtId cipar = EvtPDL::chargeConj( ipar );
241
242 if ( decaytable[ipar.getAlias()].getNMode() != 0 )
243 {
244
245 report( DEBUG, "EvtGen" )
246 << "Redefined decay of " << name.c_str() << " in CDecay" << endl;
247
248 decaytable[ipar.getAlias()].removeDecay();
249 }
250
251 // take contents of cipar and conjugate and store in ipar
252 decaytable[ipar.getAlias()].makeChargeConj( &decaytable[cipar.getAlias()] );
253 }
254 else if ( token == "Define" )
255 {
256
257 std::string name;
258
259 name = parser.getToken( itoken++ );
260 // value=atof(parser.getToken(itoken++).c_str());
261
262 EvtSymTable::Define( name, parser.getToken( itoken++ ) );
263
264 // New code Lange April 10, 2001 - allow the user
265 // to change particle definitions of EXISTING
266 // particles on the fly
267 }
268 else if ( token == "Particle" )
269 {
270
271 std::string pname;
272 pname = parser.getToken( itoken++ );
273 report( INFO, "EvtGen" ) << pname.c_str() << endl;
274 // There should be at least the mass
275 double newMass = atof( parser.getToken( itoken++ ).c_str() );
276 EvtId thisPart = EvtPDL::getId( pname );
277 double newWidth = EvtPDL::getMeanMass( thisPart );
278 if ( parser.getNToken() > 3 ) newWidth = atof( parser.getToken( itoken++ ).c_str() );
279
280 // Now make the change!
281 EvtPDL::reSetMass( thisPart, newMass );
282 EvtPDL::reSetWidth( thisPart, newWidth );
283
284 report( INFO, "EvtGen" ) << "Changing particle properties of " << pname.c_str()
285 << " Mass=" << newMass << " Width=" << newWidth << endl;
286 }
287 else if ( token == "ChangeMassMin" )
288 {
289 std::string pname;
290 pname = parser.getToken( itoken++ );
291 double tmass = atof( parser.getToken( itoken++ ).c_str() );
292
293 EvtId thisPart = EvtPDL::getId( pname );
294 EvtPDL::reSetMassMin( thisPart, tmass );
295 report( DEBUG, "EvtGen" ) << "Refined minimum mass for "
296 << EvtPDL::name( thisPart ).c_str() << " to be " << tmass
297 << endl;
298 }
299 else if ( token == "ChangeMassMax" )
300 {
301 std::string pname;
302 pname = parser.getToken( itoken++ );
303 double tmass = atof( parser.getToken( itoken++ ).c_str() );
304 EvtId thisPart = EvtPDL::getId( pname );
305 EvtPDL::reSetMassMax( thisPart, tmass );
306 report( DEBUG, "EvtGen" ) << "Refined maximum mass for "
307 << EvtPDL::name( thisPart ).c_str() << " to be " << tmass
308 << endl;
309 }
310 else if ( token == "IncludeBirthFactor" )
311 {
312 std::string pname;
313 pname = parser.getToken( itoken++ );
314 bool yesno = false;
315 if ( parser.getToken( itoken++ ).c_str() == "yes" ) yesno = true;
316 EvtId thisPart = EvtPDL::getId( pname );
317 EvtPDL::includeBirthFactor( thisPart, yesno );
318 if ( yesno )
319 report( DEBUG, "EvtGen" )
320 << "Include birth factor for " << EvtPDL::name( thisPart ).c_str() << endl;
321 if ( !yesno )
322 report( DEBUG, "EvtGen" ) << "No longer include birth factor for "
323 << EvtPDL::name( thisPart ).c_str() << endl;
324 }
325 else if ( token == "IncludeDecayFactor" )
326 {
327 std::string pname;
328 pname = parser.getToken( itoken++ );
329 bool yesno = false;
330 if ( parser.getToken( itoken++ ).c_str() == "yes" ) yesno = true;
331 EvtId thisPart = EvtPDL::getId( pname );
332 EvtPDL::includeDecayFactor( thisPart, yesno );
333 if ( yesno )
334 report( DEBUG, "EvtGen" )
335 << "Include decay factor for " << EvtPDL::name( thisPart ).c_str() << endl;
336 if ( !yesno )
337 report( DEBUG, "EvtGen" ) << "No longer include decay factor for "
338 << EvtPDL::name( thisPart ).c_str() << endl;
339 }
340 else if ( token == "AddFactorPn" )
341 {
342 std::string pname;
343 pname = parser.getToken( itoken++ );
344 double factor = atof( parser.getToken( itoken++ ).c_str() );
345 EvtId thisPart = EvtPDL::getId( pname );
346 EvtPDL::addFactorPn( thisPart, factor );
347 report( DEBUG, "EvtGen" ) << "Include momentum factor Pn= " << factor << " for "
348 << EvtPDL::name( thisPart ).c_str() << endl;
349 }
350 else if ( token == "LSNONRELBW" )
351 {
352 std::string pname;
353 pname = parser.getToken( itoken++ );
354 EvtId thisPart = EvtPDL::getId( pname );
355 std::string tstr = "NONRELBW";
356 EvtPDL::changeLS( thisPart, tstr );
357 report( DEBUG, "EvtGen" ) << "Change lineshape to non-rel BW for "
358 << EvtPDL::name( thisPart ).c_str() << endl;
359 }
360 else if ( token == "SP6LSFIX" )
361 {
362 std::string pname;
363 pname = parser.getToken( itoken++ );
364 EvtId thisPart = EvtPDL::getId( pname );
365 EvtPDL::fixLSForSP8( thisPart );
366 report( DEBUG, "EvtGen" ) << "Fixed lineshape for SP6 --from M.Baak "
367 << EvtPDL::name( thisPart ).c_str() << endl;
368 }
369 else if ( token == "LSFLAT" )
370 {
371 std::string pname;
372 pname = parser.getToken( itoken++ );
373 EvtId thisPart = EvtPDL::getId( pname );
374 std::string tstr = "FLAT";
375 EvtPDL::changeLS( thisPart, tstr );
376 report( DEBUG, "EvtGen" ) << "Change lineshape to flat for "
377 << EvtPDL::name( thisPart ).c_str() << endl;
378 }
379 else if ( token == "LSMANYDELTAFUNC" )
380 {
381 std::string pname;
382 pname = parser.getToken( itoken++ );
383 EvtId thisPart = EvtPDL::getId( pname );
384 std::string tstr = "MANYDELTAFUNC";
385 EvtPDL::changeLS( thisPart, tstr );
386 report( DEBUG, "EvtGen" ) << "Change lineshape to spikes for "
387 << EvtPDL::name( thisPart ).c_str() << endl;
388 }
389 else if ( token == "BlattWeisskopf" )
390 {
391 std::string pname;
392 pname = parser.getToken( itoken++ );
393 double tnum = atof( parser.getToken( itoken++ ).c_str() );
394 EvtId thisPart = EvtPDL::getId( pname );
395 EvtPDL::reSetBlatt( thisPart, tnum );
396 report( DEBUG, "EvtGen" ) << "Redefined Blatt-Weisskopf factor "
397 << EvtPDL::name( thisPart ).c_str() << " to be " << tnum
398 << endl;
399 }
400 else if ( token == "SetLineshapePW" )
401 {
402 std::string pname;
403 pname = parser.getToken( itoken++ );
404 EvtId thisPart = EvtPDL::getId( pname );
405 std::string pnameD1 = parser.getToken( itoken++ );
406 EvtId thisD1 = EvtPDL::getId( pnameD1 );
407 std::string pnameD2 = parser.getToken( itoken++ );
408 EvtId thisD2 = EvtPDL::getId( pnameD2 );
409 int pw = atoi( parser.getToken( itoken++ ).c_str() );
410 report( DEBUG, "EvtGen" ) << "Redefined Partial wave for " << pname.c_str() << " to "
411 << pnameD1.c_str() << " " << pnameD2.c_str() << " (" << pw
412 << ")" << endl;
413 EvtPDL::setPWForDecay( thisPart, pw, thisD1, thisD2 );
414 }
415 else if ( token == "Decay" )
416 {
417
418 std::string temp_fcn_new_model;
419
420 EvtDecayBase* temp_fcn_new;
421
422 double brfrsum = 0.0;
423
424 parent = parser.getToken( itoken++ );
425 ipar = EvtPDL::getId( parent );
426
427 if ( ipar == EvtId( -1, -1 ) )
428 {
429 report( ERROR, "EvtGen" ) << "Unknown particle name:" << parent.c_str() << " on line "
430 << parser.getLineofToken( itoken - 1 ) << endl;
431 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
432 ::abort();
433 }
434
435 if ( decaytable[ipar.getAlias()].getNMode() != 0 )
436 {
437 report( DEBUG, "EvtGen" ) << "Redefined decay of " << parent.c_str() << endl;
438 decaytable[ipar.getAlias()].removeDecay();
439 }
440
441 do {
442
443 token = parser.getToken( itoken++ );
444
445 if ( token != "Enddecay" )
446 {
447
448 i = 0;
449 while ( token.c_str()[i++] != 0 )
450 {
451 if ( isalpha( token.c_str()[i] ) )
452 {
453 report( ERROR, "EvtGen" ) << "Expected to find a branching fraction or Enddecay "
454 << "but found:" << token.c_str() << " on line "
455 << parser.getLineofToken( itoken - 1 ) << endl;
456 report( ERROR, "EvtGen" ) << "Possibly to few arguments to model "
457 << "on previous line!" << endl;
458 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
459 ::abort();
460 }
461 }
462
463 brfr = atof( token.c_str() );
464
465 int isname = EvtPDL::getId( parser.getToken( itoken ) ).getId() >= 0;
466 int ismodel = modelist.isModel( parser.getToken( itoken ) );
467
468 if ( !( isname || ismodel ) )
469 {
470 // see if this is an aliased model
471 int iAlias;
472 for ( iAlias = 0; iAlias < modelAliasList.size(); iAlias++ )
473 {
474 if ( modelAliasList[iAlias].matchAlias( parser.getToken( itoken ) ) )
475 {
476 ismodel = 2;
477 break;
478 }
479 }
480 }
481
482 if ( !( isname || ismodel ) )
483 {
484
485 report( INFO, "EvtGen" )
486 << parser.getToken( itoken ).c_str() << " is neither a particle name nor "
487 << "the name of a model. " << endl;
488 report( INFO, "EvtGen" )
489 << "It was encountered on line " << parser.getLineofToken( itoken )
490 << " of the decay file." << endl;
491 report( INFO, "EvtGen" ) << "Please fix it. Thank you." << endl;
492 report( INFO, "EvtGen" ) << "Be sure to check that the "
493 << "correct case has been used. \n";
494 report( INFO, "EvtGen" ) << "Terminating execution. \n";
495 ::abort();
496
497 itoken++;
498 }
499
500 n_daugh = 0;
501
502 while ( EvtPDL::getId( parser.getToken( itoken ) ).getId() >= 0 )
503 {
504 sdaug = parser.getToken( itoken++ );
505 daught[n_daugh++] = EvtPDL::getId( sdaug );
506 if ( daught[n_daugh - 1] == EvtId( -1, -1 ) )
507 {
508 report( ERROR, "EvtGen" )
509 << "Unknown particle name:" << sdaug.c_str() << " on line "
510 << parser.getLineofToken( itoken ) << endl;
511 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
512 ::abort();
513 }
514 }
515
516 model = parser.getToken( itoken++ );
517
518 int photos = 0;
519 int verbose = 0;
520 int summary = 0;
521
522 do {
523 if ( model == "PHOTOS" )
524 {
525 photos = 1;
526 model = parser.getToken( itoken++ );
527 }
528 if ( model == "VERBOSE" )
529 {
530 verbose = 1;
531 model = parser.getToken( itoken++ );
532 }
533 if ( model == "SUMMARY" )
534 {
535 summary = 1;
536 model = parser.getToken( itoken++ );
537 }
538 } while ( model == "PHOTOS" || model == "VERBOSE" || model == "SUMMARY" );
539
540 // see if this is an aliased model
541 int iAlias;
542 int foundAnAlias = -1;
543 for ( iAlias = 0; iAlias < modelAliasList.size(); iAlias++ )
544 {
545 if ( modelAliasList[iAlias].matchAlias( model ) )
546 {
547 foundAnAlias = iAlias;
548 break;
549 }
550 }
551
552 if ( foundAnAlias == -1 )
553 {
554 if ( !modelist.isModel( model ) )
555 {
556 report( ERROR, "EvtGen" ) << "Expected to find a model name,"
557 << "found:" << model.c_str() << " on line "
558 << parser.getLineofToken( itoken ) << endl;
559 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
560 ::abort();
561 }
562 }
563 else { model = modelAliasList[foundAnAlias].getName(); }
564
565 temp_fcn_new_model = model;
566 temp_fcn_new = modelist.getFcn( model );
567
568 if ( photos ) { temp_fcn_new->setPHOTOS(); }
569 if ( verbose ) { temp_fcn_new->setVerbose(); }
570 if ( summary ) { temp_fcn_new->setSummary(); }
571
572 std::vector<std::string> temp_fcn_new_args;
573
574 std::string name;
575 int ierr;
576
577 if ( foundAnAlias == -1 )
578 {
579 do {
580 name = parser.getToken( itoken++ );
581 if ( name != ";" )
582 {
583 temp_fcn_new_args.push_back( EvtSymTable::Get( name, ierr ) );
584 if ( ierr )
585 {
586 report( ERROR, "EvtGen" )
587 << "Reading arguments and found:" << name.c_str()
588 << " on line:" << parser.getLineofToken( itoken - 1 ) << endl;
589 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
590 ::abort();
591 }
592 }
593 // int isname=EvtPDL::getId(name).getId()>=0;
594 int ismodel = modelist.isModel( name );
595 if ( ismodel )
596 {
597 report( ERROR, "EvtGen" )
598 << "Expected ';' but found:" << name.c_str()
599 << " on line:" << parser.getLineofToken( itoken - 1 ) << endl;
600 report( ERROR, "EvtGen" ) << "Most probable error is omitted ';'." << endl;
601 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
602 ::abort();
603 }
604 } while ( name != ";" );
605 }
606 else
607 {
608 std::vector<std::string> copyMe = modelAliasList[foundAnAlias].getArgList();
609 temp_fcn_new_args = copyMe;
610 itoken++;
611 }
612 // Found one decay.
613
614 brfrsum += brfr;
615
616 temp_fcn_new->saveDecayInfo( ipar, n_daugh, daught, temp_fcn_new_args.size(),
617 temp_fcn_new_args, temp_fcn_new_model, brfr );
618
619 double massmin = 0.0;
620
621 // for (i=0;i<n_daugh;i++){
622 for ( i = 0; i < temp_fcn_new->nRealDaughters(); i++ )
623 {
624 if ( EvtPDL::getMinMass( daught[i] ) > 0.0001 )
625 { massmin += EvtPDL::getMinMass( daught[i] ); }
626 else { massmin += EvtPDL::getMeanMass( daught[i] ); }
627 }
628
629 decaytable[ipar.getAlias()].addMode( temp_fcn_new, brfrsum, massmin );
630 }
631 } while ( token != "Enddecay" );
632
633 decaytable[ipar.getAlias()].finalize();
634 }
635 // Allow copying of decays from one particle to another; useful
636 // in combination with RemoveDecay
637 else if ( token == "CopyDecay" )
638 {
639 std::string newname;
640 std::string oldname;
641
642 newname = parser.getToken( itoken++ );
643 oldname = parser.getToken( itoken++ );
644
645 EvtId newipar = EvtPDL::getId( newname );
646 EvtId oldipar = EvtPDL::getId( oldname );
647
648 if ( oldipar == EvtId( -1, -1 ) )
649 {
650 report( ERROR, "EvtGen" ) << "Unknown particle name:" << oldname.c_str() << " on line "
651 << parser.getLineofToken( itoken ) << endl;
652 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
653 ::abort();
654 }
655 if ( newipar == EvtId( -1, -1 ) )
656 {
657 report( ERROR, "EvtGen" ) << "Unknown particle name:" << newname.c_str() << " on line "
658 << parser.getLineofToken( itoken ) << endl;
659 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
660 ::abort();
661 }
662 if ( decaytable[newipar.getAlias()].getNMode() != 0 )
663 {
664 report( DEBUG, "EvtGen" ) << "Redefining decay of " << newname << endl;
665 decaytable[newipar.getAlias()].removeDecay();
666 }
667 decaytable[newipar.getAlias()] = decaytable[oldipar.getAlias()];
668 }
669 // Enable decay deletion; intended primarily for aliases
670 // Peter Onyisi, March 2008
671 else if ( token == "RemoveDecay" )
672 {
673 parent = parser.getToken( itoken++ );
674 ipar = EvtPDL::getId( parent );
675
676 if ( ipar == EvtId( -1, -1 ) )
677 {
678 report( ERROR, "EvtGen" ) << "Unknown particle name:" << parent.c_str() << " on line "
679 << parser.getLineofToken( itoken - 1 ) << endl;
680 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
681 ::abort();
682 }
683
684 if ( decaytable[ipar.getAlias()].getNMode() == 0 )
685 { report( DEBUG, "EvtGen" ) << "No decays to delete for " << parent.c_str() << endl; }
686 else
687 {
688 report( DEBUG, "EvtGen" ) << "Deleting selected decays of " << parent.c_str() << endl;
689 }
690
691 do {
692 token = parser.getToken( itoken );
693
694 if ( token != "Enddecay" )
695 {
696 n_daugh = 0;
697 while ( EvtPDL::getId( parser.getToken( itoken ) ).getId() >= 0 )
698 {
699 sdaug = parser.getToken( itoken++ );
700 daught[n_daugh++] = EvtPDL::getId( sdaug );
701 if ( daught[n_daugh - 1] == EvtId( -1, -1 ) )
702 {
703 report( ERROR, "EvtGen" )
704 << "Unknown particle name:" << sdaug.c_str() << " on line "
705 << parser.getLineofToken( itoken ) << endl;
706 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
707 ::abort();
708 }
709 }
710 token = parser.getToken( itoken );
711 if ( token != ";" )
712 {
713 report( ERROR, "EvtGen" )
714 << "Expected ';' but found:" << token
715 << " on line:" << parser.getLineofToken( itoken - 1 ) << endl;
716 report( ERROR, "EvtGen" ) << "Most probable error is omitted ';'." << endl;
717 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
718 ::abort();
719 }
720 token = parser.getToken( itoken++ );
721 EvtDecayBase* temp_fcn_new = modelist.getFcn( "PHSP" );
722 std::vector<std::string> temp_fcn_new_args;
723 std::string temp_fcn_new_model( "PHSP" );
724 temp_fcn_new->saveDecayInfo( ipar, n_daugh, daught, 0, temp_fcn_new_args,
725 temp_fcn_new_model, 0. );
726 decaytable[ipar.getAlias()].removeMode( temp_fcn_new );
727 }
728 } while ( token != "Enddecay" );
729 itoken++;
730 }
731 else if ( token != "End" )
732 {
733
734 report( ERROR, "EvtGen" ) << "Found unknown command:'" << token.c_str() << "' on line "
735 << parser.getLineofToken( itoken ) << endl;
736 report( ERROR, "EvtGen" ) << "Will terminate execution!" << endl;
737 ::abort();
738 }
739
740 } while ( ( token != "End" ) && itoken != parser.getNToken() );
741
742 // Now we may need to reset the minimum mass for some particles????
743
744 int ii;
745 for ( ii = 0; ii < EvtPDL::entries(); ii++ )
746 {
747 EvtId temp( ii, ii );
748 int nModTot = getNMode( ii );
749 // no decay modes
750 if ( nModTot == 0 ) continue;
751 // 0 width?
752 if ( EvtPDL::getWidth( temp ) < 0.0000001 ) continue;
753 int jj;
754 double minMass = EvtPDL::getMaxMass( temp );
755 for ( jj = 0; jj < nModTot; jj++ )
756 {
757 double tmass = decaytable[ii].getDecay( jj ).getMassMin();
758 if ( tmass < minMass ) minMass = tmass;
759 }
760 if ( minMass > EvtPDL::getMinMass( temp ) )
761 {
762
763 report( INFO, "EvtGen" ) << "Given allowed decays, resetting minMass "
764 << EvtPDL::name( temp ).c_str() << " "
765 << EvtPDL::getMinMass( temp ) << " to " << minMass << endl;
766 EvtPDL::reSetMassMin( temp, minMass );
767 }
768 }
769}
770
771int EvtDecayTable::findChannel( EvtId parent, std::string model, int ndaug, EvtId* daugs,
772 int narg, std::string* args ) {
773
774 int i, j, right;
775 EvtId daugs_scratch[50];
776 int nmatch, k;
777
778 for ( i = 0; i < decaytable[parent.getAlias()].getNMode(); i++ )
779 {
780
781 right = 1;
782
783 right =
784 right &&
785 model == decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getModelName();
786 right =
787 right &&
788 ( ndaug == decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getNDaug() );
789 right = right &&
790 ( narg == decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getNArg() );
791
792 if ( right )
793 {
794
795 for ( j = 0; j < ndaug; j++ ) { daugs_scratch[j] = daugs[j]; }
796
797 nmatch = 0;
798
799 for ( j = 0; j < decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getNDaug();
800 j++ )
801 {
802
803 for ( k = 0; k < ndaug; k++ )
804 {
805 if ( daugs_scratch[k] ==
806 decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getDaug( j ) )
807 {
808 daugs_scratch[k] = EvtId( -1, -1 );
809 nmatch++;
810 break;
811 }
812 }
813 }
814
815 right = right && ( nmatch == ndaug );
816
817 for ( j = 0; j < decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getNArg();
818 j++ )
819 {
820 right =
821 right &&
822 ( args[j] ==
823 decaytable[parent.getAlias()].getDecay( i ).getDecayModel()->getArgStr( j ) );
824 }
825 }
826 if ( right ) return i;
827 }
828 return -1;
829}
830
831int EvtDecayTable::inChannelList( EvtId parent, int ndaug, EvtId* daugs ) {
832
833 int i, j, k;
834 // std::cout<<"=============MAX_DAUG = "<<MAX_DAUG<<endl;
835 EvtId daugs_scratch[MAX_DAUG];
836
837 int dsum = 0;
838 for ( i = 0; i < ndaug; i++ ) { dsum += daugs[i].getAlias(); }
839
840 int nmatch;
841
842 int ipar = parent.getAlias();
843
844 int nmode = decaytable[ipar].getNMode();
845
846 for ( i = 0; i < nmode; i++ )
847 {
848
849 EvtDecayBase* thedecaymodel = decaytable[ipar].getDecay( i ).getDecayModel();
850
851 if ( thedecaymodel->getDSum() == dsum )
852 {
853
854 int nd = thedecaymodel->getNDaug();
855
856 if ( ndaug == nd )
857 {
858 for ( j = 0; j < ndaug; j++ ) { daugs_scratch[j] = daugs[j]; }
859 nmatch = 0;
860 for ( j = 0; j < nd; j++ )
861 {
862 for ( k = 0; k < ndaug; k++ )
863 {
864 if ( EvtId( daugs_scratch[k] ) == thedecaymodel->getDaug( j ) )
865 {
866 daugs_scratch[k] = EvtId( -1, -1 );
867 nmatch++;
868 break;
869 }
870 }
871 }
872 if ( ( nmatch == ndaug ) && ( !( ( thedecaymodel->getModelName() == "JETSET" ) ||
873 ( thedecaymodel->getModelName() == "LUNDCHARM" ) ||
874 ( thedecaymodel->getModelName() == "PYTHIA" ) ) ) )
875 { return i; }
876 }
877 }
878 }
879
880 return -1;
881}
getDecay
EvtDecayBase * getDecay(int ipar, int imode)
Definition EvtDecayTable.cc:51
MAX_DAUG
const int MAX_DAUG
Definition EvtParticle.hh:38
report
ostream & report(Severity severity, const char *facility)
Definition EvtReport.cc:34
DEBUG
@ DEBUG
Definition EvtReport.hh:53
ERROR
@ ERROR
Definition EvtReport.hh:49
INFO
@ INFO
Definition EvtReport.hh:52
EvtDecayBase::getModelName
std::string getModelName()
Definition EvtDecayBase.hh:79
EvtDecayBase::saveDecayInfo
void saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
Definition EvtDecayBase.cc:151
EvtDecayBase::nRealDaughters
virtual int nRealDaughters()
Definition EvtDecayBase.hh:105
EvtDecayBase::setSummary
void setSummary()
Definition EvtDecayBase.hh:74
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition EvtDecayBase.hh:69
EvtDecayBase::setVerbose
void setVerbose()
Definition EvtDecayBase.hh:73
EvtDecayTable::gettheDecay
static EvtDecayBase * gettheDecay(EvtId parent, int imode)
Definition EvtDecayTable.cc:55
EvtDecayTable::findChannel
static int findChannel(EvtId parent, std::string model, int ndaug, EvtId *daugs, int narg, std::string *args)
Definition EvtDecayTable.cc:771
EvtDecayTable::getNMode
static int getNMode(int ipar)
Definition EvtDecayTable.cc:49
EvtDecayTable::inChannelList
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition EvtDecayTable.cc:831
EvtDecayTable::getDecayFunc
static EvtDecayBase * getDecayFunc(EvtParticle *)
Definition EvtDecayTable.cc:68
EvtDecayTable::getDecay
static EvtDecayBase * getDecay(int ipar, int imode)
EvtDecayTable::readDecayFile
static void readDecayFile(const std::string dec_name)
Definition EvtDecayTable.cc:77
EvtDecayTable::printSummary
static void printSummary()
Definition EvtDecayTable.cc:61
EvtId
Definition EvtId.hh:27
EvtId::getAlias
int getAlias() const
Definition EvtId.hh:42
EvtId::getId
int getId() const
Definition EvtId.hh:40
EvtModel::getFcn
EvtDecayBase * getFcn(std::string model_name)
Definition EvtModel.cc:44
EvtModel::storeCommand
void storeCommand(std::string cmd, std::string cnfgstr)
Definition EvtModel.cc:89
EvtModel::isModel
int isModel(std::string name)
Definition EvtModel.cc:77
EvtModel::instance
static EvtModel & instance()
Definition EvtModel.hh:56
EvtModel::isCommand
int isCommand(std::string cmd)
Definition EvtModel.cc:83
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition EvtPDL.hh:59
EvtPDL::reSetBlatt
static void reSetBlatt(EvtId i, double blatt)
Definition EvtPDL.hh:84
EvtPDL::addFactorPn
static void addFactorPn(EvtId i, double factor)
Definition EvtPDL.hh:93
EvtPDL::fixLSForSP8
static void fixLSForSP8(EvtId i)
Definition EvtPDL.hh:103
EvtPDL::changeLS
static void changeLS(EvtId i, std::string &newLS)
Definition EvtPDL.hh:96
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition EvtPDL.hh:43
EvtPDL::reSetWidth
static void reSetWidth(EvtId i, double width)
Definition EvtPDL.hh:75
EvtPDL::entries
static int entries()
Definition EvtPDL.hh:73
EvtPDL::reSetMass
static void reSetMass(EvtId i, double mass)
Definition EvtPDL.hh:74
EvtPDL::alias
static void alias(EvtId num, const std::string &newname)
Definition EvtPDL.cc:244
EvtPDL::includeDecayFactor
static void includeDecayFactor(EvtId i, bool yesno)
Definition EvtPDL.hh:90
EvtPDL::reSetMassMax
static void reSetMassMax(EvtId i, double mass)
Definition EvtPDL.hh:81
EvtPDL::name
static std::string name(EvtId i)
Definition EvtPDL.hh:70
EvtPDL::chargeConj
static EvtId chargeConj(EvtId id)
Definition EvtPDL.cc:193
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition EvtPDL.hh:56
EvtPDL::setPWForDecay
static void setPWForDecay(EvtId i, int spin, EvtId d1, EvtId d2)
Definition EvtPDL.hh:99
EvtPDL::aliasChgConj
static void aliasChgConj(EvtId a, EvtId abar)
Definition EvtPDL.cc:176
EvtPDL::includeBirthFactor
static void includeBirthFactor(EvtId i, bool yesno)
Definition EvtPDL.hh:87
EvtPDL::getMaxMass
static double getMaxMass(EvtId i)
Definition EvtPDL.hh:55
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:272
EvtPDL::reSetMassMin
static void reSetMassMin(EvtId i, double mass)
Definition EvtPDL.hh:78
EvtParser::getLineofToken
int getLineofToken(int i)
Definition EvtParser.cc:52
EvtParser::getNToken
int getNToken()
Definition EvtParser.cc:48
EvtParser::getToken
const std::string & getToken(int i)
Definition EvtParser.cc:50
EvtParser::Read
int Read(const std::string filename)
Definition EvtParser.cc:54
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:100
EvtRadCorr::setAlwaysRadCorr
static void setAlwaysRadCorr()
Definition EvtRadCorr.cc:64
EvtRadCorr::setNormalRadCorr
static void setNormalRadCorr()
Definition EvtRadCorr.cc:72
EvtRadCorr::setNeverRadCorr
static void setNeverRadCorr()
Definition EvtRadCorr.cc:68
EvtSymTable::Define
static void Define(const std::string &name, std::string d)
Definition EvtSymTable.cc:38
EvtSymTable::Get
static std::string Get(const std::string &name, int &ierr)
Definition EvtSymTable.cc:54