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