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