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