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: EvtDecayBase.cc
12 : //
13 : // Description: Store decay parameters for one decay.
14 : //
15 : // Modification history:
16 : //
17 : // RYD September 30, 1997 Module created
18 : //
19 : //------------------------------------------------------------------------
20 : //
21 : #include "EvtGenBase/EvtPatches.hh"
22 : #include <iostream>
23 : #include <fstream>
24 : #include <stdlib.h>
25 : #include <ctype.h>
26 : #include "EvtGenBase/EvtStatus.hh"
27 : #include "EvtGenBase/EvtDecayBase.hh"
28 : #include "EvtGenBase/EvtParticle.hh"
29 : #include "EvtGenBase/EvtPDL.hh"
30 : #include "EvtGenBase/EvtReport.hh"
31 : #include "EvtGenBase/EvtSpinType.hh"
32 : #include <vector>
33 : using std::endl;
34 : using std::fstream;
35 : void EvtDecayBase::checkQ() {
36 : int i;
37 : int q=0;
38 : int qpar;
39 :
40 : //If there are no daughters (jetset etc) then we do not
41 : //want to do this test. Why? Because sometimes the parent
42 : //will have a nonzero charge.
43 :
44 0 : if ( _ndaug != 0) {
45 0 : for(i=0; i<_ndaug; i++ ) {
46 0 : q += EvtPDL::chg3(_daug[i]);
47 : }
48 0 : qpar = EvtPDL::chg3(_parent);
49 :
50 0 : if ( q != qpar ) {
51 0 : report(Severity::Error,"EvtGen") <<_modelname.c_str()<< " generator expected "
52 0 : << " charge to be conserved, found:"<<endl;
53 0 : report(Severity::Error,"EvtGen") << "Parent charge of "<<(qpar/3)<<endl;
54 0 : report(Severity::Error,"EvtGen") << "Sum of daughter charge of "<<(q/3)<<endl;
55 0 : report(Severity::Error,"EvtGen") << "The parent is "<< EvtPDL::name(_parent).c_str()<<endl;
56 0 : for(i=0; i<_ndaug; i++ ) {
57 0 : report(Severity::Error,"EvtGen") << "Daughter "<< EvtPDL::name(_daug[i]).c_str()<<endl;
58 : }
59 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
60 :
61 0 : ::abort();
62 : }
63 : }
64 0 : }
65 :
66 :
67 : double EvtDecayBase::getProbMax( double prob ) {
68 :
69 : int i;
70 :
71 : //diagnostics
72 0 : sum_prob+=prob;
73 0 : if (prob>max_prob) max_prob=prob;
74 :
75 :
76 0 : if ( defaultprobmax && ntimes_prob<=500 ) {
77 : //We are building up probmax with this iteration
78 0 : ntimes_prob += 1;
79 0 : if ( prob > probmax ) { probmax = prob;}
80 0 : if (ntimes_prob==500) {
81 0 : probmax*=1.2;
82 0 : }
83 0 : return 1000000.0*prob;
84 : }
85 :
86 0 : if ( prob> probmax*1.0001) {
87 :
88 0 : report(Severity::Info,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
89 0 : report(Severity::Info,"") << "("<<_modelname.c_str()<<") ";
90 0 : report(Severity::Info,"") << EvtPDL::name(_parent).c_str()<<" -> ";
91 0 : for(i=0;i<_ndaug;i++){
92 0 : report(Severity::Info,"") << EvtPDL::name(_daug[i]).c_str() << " ";
93 : }
94 0 : report(Severity::Info,"") << endl;
95 :
96 0 : if (defaultprobmax) probmax = prob;
97 :
98 : }
99 :
100 0 : ntimes_prob += 1;
101 :
102 :
103 0 : return probmax;
104 :
105 0 : } //getProbMax
106 :
107 :
108 : double EvtDecayBase::resetProbMax(double prob) {
109 :
110 0 : report(Severity::Info,"EvtGen") << "Reseting prob max\n";
111 0 : report(Severity::Info,"EvtGen") << "prob > probmax:("<<prob<<">"<<probmax<<")";
112 0 : report(Severity::Info,"") << "("<<_modelname.c_str()<<")";
113 0 : report(Severity::Info,"") << EvtPDL::getStdHep(_parent)<<"->";
114 :
115 0 : for( int i=0;i<_ndaug;i++){
116 0 : report(Severity::Info,"") << EvtPDL::getStdHep(_daug[i]) << " ";
117 : }
118 0 : report(Severity::Info,"") << endl;
119 :
120 0 : probmax = 0.0;
121 0 : defaultprobmax = 0;
122 0 : ntimes_prob = 0;
123 :
124 0 : return prob;
125 :
126 : }
127 :
128 :
129 : std::string EvtDecayBase::commandName(){
130 0 : return std::string("");
131 : }
132 :
133 : void EvtDecayBase::command(std::string){
134 0 : report(Severity::Error,"EvtGen") << "Should never call EvtDecayBase::command"<<endl;
135 0 : ::abort();
136 : }
137 :
138 : std::string EvtDecayBase::getParamName(int i) {
139 0 : switch(i) {
140 : case 0:
141 0 : return "param00";
142 : case 1:
143 0 : return "param01";
144 : case 2:
145 0 : return "param02";
146 : case 3:
147 0 : return "param03";
148 : case 4:
149 0 : return "param04";
150 : case 5:
151 0 : return "param05";
152 : case 6:
153 0 : return "param06";
154 : case 7:
155 0 : return "param07";
156 : case 8:
157 0 : return "param08";
158 : case 9:
159 0 : return "param09";
160 : default:
161 0 : return "";
162 : }
163 0 : }
164 :
165 : std::string EvtDecayBase::getParamDefault(int /*i*/) {
166 0 : return "";
167 : }
168 :
169 : void EvtDecayBase::init() {
170 :
171 : //This default version of init does nothing;
172 : //A specialized version of this function can be
173 : //supplied for each decay model to do initialization.
174 :
175 0 : return;
176 :
177 : }
178 :
179 : void EvtDecayBase::initProbMax() {
180 :
181 : //This function is called if the decay does not have a
182 : //specialized initialization.
183 : //The default is to set the maximum
184 : //probability to 0 and the number of times called to 0
185 : //and defaultprobmax to 1 such that the decay will be
186 : //generated many many times
187 : //in order to generate a reasonable maximum probability
188 : //for the decay.
189 :
190 0 : defaultprobmax=1;
191 0 : ntimes_prob = 0;
192 0 : probmax = 0.0;
193 :
194 0 : } //initProbMax
195 :
196 :
197 : void EvtDecayBase::saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug,
198 : int narg,std::vector<std::string>& args,
199 : std::string name,
200 : double brfr) {
201 :
202 : int i;
203 :
204 0 : _brfr=brfr;
205 0 : _ndaug=ndaug;
206 0 : _narg=narg;
207 0 : _parent=ipar;
208 :
209 0 : _dsum=0;
210 :
211 0 : if (_ndaug>0) {
212 0 : _daug=new EvtId [_ndaug];
213 0 : for(i=0;i<_ndaug;i++){
214 0 : _daug[i]=daug[i];
215 0 : _dsum+=daug[i].getAlias();
216 : }
217 : }
218 : else{
219 0 : _daug=0;
220 : }
221 :
222 0 : if (_narg>0) {
223 0 : _args=new std::string[_narg+1];
224 0 : for(i=0;i<_narg;i++){
225 0 : _args[i]=args[i];
226 : }
227 : }
228 : else{
229 0 : _args = 0;
230 : }
231 :
232 0 : _modelname=name;
233 :
234 0 : this->init();
235 0 : this->initProbMax();
236 :
237 0 : if (_chkCharge){
238 0 : this->checkQ();
239 0 : }
240 :
241 :
242 0 : if (defaultprobmax){
243 0 : report(Severity::Info,"EvtGen") << "No default probmax for ";
244 0 : report(Severity::Info,"") << "("<<_modelname.c_str()<<") ";
245 0 : report(Severity::Info,"") << EvtPDL::name(_parent).c_str()<<" -> ";
246 0 : for(i=0;i<_ndaug;i++){
247 0 : report(Severity::Info,"") << EvtPDL::name(_daug[i]).c_str() << " ";
248 : }
249 0 : report(Severity::Info,"") << endl;
250 0 : report(Severity::Info,"") << "This is fine for development, but must be provided for production."<<endl;
251 0 : report(Severity::Info,"EvtGen") << "Never fear though - the decay will use the \n";
252 0 : report(Severity::Info,"EvtGen") << "500 iterations to build up a good probmax \n";
253 0 : report(Severity::Info,"EvtGen") << "before accepting a decay. "<<endl;
254 0 : }
255 :
256 0 : }
257 :
258 0 : EvtDecayBase::EvtDecayBase() {
259 :
260 : //the default is that the user module does _not_ set
261 : // any probmax.
262 0 : defaultprobmax=1;
263 0 : ntimes_prob = 0;
264 0 : probmax = 0.0;
265 :
266 0 : _photos=0;
267 0 : _verbose=0;
268 0 : _summary=0;
269 0 : _parent=EvtId(-1,-1);
270 0 : _ndaug=0;
271 0 : _narg=0;
272 0 : _daug=0;
273 0 : _args=0;
274 0 : _argsD=0;
275 0 : _modelname="**********";
276 :
277 : //Default is to check that charge is conserved
278 0 : _chkCharge=1;
279 :
280 : //statistics collection!
281 :
282 0 : max_prob=0.0;
283 0 : sum_prob=0.0;
284 :
285 0 : }
286 :
287 :
288 :
289 : void EvtDecayBase::printSummary() const {
290 0 : if (ntimes_prob>0) {
291 :
292 0 : report(Severity::Info,"EvtGen") << "Calls = "<<ntimes_prob<<" eff: "<<
293 0 : sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax;
294 0 : report(Severity::Info,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : ";
295 0 : }
296 :
297 0 : printInfo();
298 0 : }
299 :
300 :
301 : void EvtDecayBase::printInfo() const {
302 0 : report(Severity::Info,"") << EvtPDL::name(_parent).c_str()<<" -> ";
303 0 : for(int i=0;i<_ndaug;i++){
304 0 : report(Severity::Info,"") << EvtPDL::name(_daug[i]).c_str() << " ";
305 : }
306 0 : report(Severity::Info,"") << " ("<<_modelname.c_str()<<")"<< endl;
307 0 : }
308 :
309 :
310 0 : EvtDecayBase::~EvtDecayBase() {
311 :
312 0 : if (_daug!=0){
313 0 : delete [] _daug;
314 : }
315 :
316 0 : if (_args!=0){
317 0 : delete [] _args;
318 : }
319 :
320 0 : if (_argsD!=0){
321 0 : delete [] _argsD;
322 : }
323 :
324 :
325 0 : }
326 :
327 : void EvtDecayBase::setProbMax(double prbmx){
328 :
329 0 : defaultprobmax=0;
330 0 : probmax=prbmx;
331 :
332 0 : }
333 :
334 : void EvtDecayBase::noProbMax(){
335 :
336 0 : defaultprobmax=0;
337 :
338 0 : }
339 :
340 :
341 : double EvtDecayBase::findMaxMass(EvtParticle *p) {
342 :
343 :
344 0 : double maxOkMass=EvtPDL::getMaxMass(p->getId());
345 :
346 : //protect against vphotons
347 0 : if ( maxOkMass < 0.0000000001 ) return 10000000.;
348 : //and against already determined masses
349 0 : if ( p->hasValidP4() ) maxOkMass=p->mass();
350 :
351 0 : EvtParticle *par=p->getParent();
352 0 : if ( par ) {
353 0 : double maxParMass=findMaxMass(par);
354 : size_t i;
355 : double minDaugMass=0.;
356 0 : for(i=0;i<par->getNDaug();i++){
357 0 : EvtParticle *dau=par->getDaug(i);
358 0 : if ( dau!=p) {
359 : // it might already have a mass
360 0 : if ( dau->isInitialized() || dau->hasValidP4() )
361 0 : minDaugMass+=dau->mass();
362 : else
363 : //give it a bit of phase space
364 0 : minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
365 : }
366 : }
367 0 : if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
368 0 : }
369 : return maxOkMass;
370 0 : }
371 :
372 :
373 : // given list of daughters ( by number ) returns a
374 : // list of viable masses.
375 :
376 : void EvtDecayBase::findMass(EvtParticle *p) {
377 :
378 : //Need to also check that this mass does not screw
379 : //up the parent
380 : //This code assumes that for the ith daughter, 0..i-1
381 : //already have a mass
382 0 : double maxOkMass=findMaxMass(p);
383 :
384 : int count=0;
385 : double mass;
386 : bool massOk=false;
387 : size_t i;
388 0 : while (!massOk) {
389 0 : count++;
390 0 : if ( count > 10000 ) {
391 0 : report(Severity::Info,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
392 0 : report(Severity::Info,"EvtGen") << "Now printing parent and/or grandparent tree\n";
393 0 : if ( p->getParent() ) {
394 0 : if ( p->getParent()->getParent() ) {
395 0 : p->getParent()->getParent()->printTree();
396 0 : report(Severity::Info,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
397 0 : report(Severity::Info,"EvtGen") << p->getParent()->mass() <<endl;
398 0 : }
399 : else{
400 0 : p->getParent()->printTree();
401 0 : report(Severity::Info,"EvtGen") << p->getParent()->mass() <<endl;
402 : }
403 : }
404 0 : else p->printTree();
405 0 : report(Severity::Info,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
406 0 : if ( p->getNDaug() ) {
407 0 : for (i=0; i<p->getNDaug(); i++) {
408 0 : report(Severity::Info,"EvtGen") << p->getDaug(i)->mass()<<" ";
409 : }
410 0 : report(Severity::Info,"EvtGen") << endl;
411 0 : }
412 0 : if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
413 0 : report(Severity::Info,"EvtGen") << "taking a default value\n";
414 0 : p->setMass(maxOkMass);
415 0 : return;
416 : }
417 0 : assert(0);
418 : }
419 0 : mass = EvtPDL::getMass(p->getId());
420 : //Just need to check that this mass is > than
421 : //the mass of all daughters
422 : double massSum=0.;
423 0 : if ( p->getNDaug() ) {
424 0 : for (i=0; i<p->getNDaug(); i++) {
425 0 : massSum+= p->getDaug(i)->mass();
426 : }
427 : }
428 : //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters
429 0 : if (p->getNDaug()<2) massOk=true;
430 0 : if ( p->getParent() ) {
431 0 : if ( p->getParent()->getNDaug()==1 ) massOk=true;
432 : }
433 0 : if ( !massOk ) {
434 0 : if (massSum < mass) massOk=true;
435 0 : if ( mass> maxOkMass) massOk=false;
436 : }
437 : }
438 :
439 0 : p->setMass(mass);
440 :
441 0 : }
442 :
443 :
444 : void EvtDecayBase::findMasses(EvtParticle *p, int ndaugs,
445 : EvtId daugs[10], double masses[10]) {
446 :
447 : int i;
448 : double mass_sum;
449 :
450 : int count=0;
451 :
452 0 : if (!( p->firstornot() )) {
453 0 : for (i = 0; i < ndaugs; i++ ) {
454 0 : masses[i] = p->getDaug(i)->mass();
455 : } //for
456 : } //if
457 : else {
458 0 : p->setFirstOrNot();
459 : // if only one daughter do it
460 :
461 0 : if (ndaugs==1) {
462 0 : masses[0]=p->mass();
463 0 : return;
464 : }
465 :
466 : //until we get a combo whose masses are less than _parent mass.
467 : do {
468 : mass_sum = 0.0;
469 :
470 0 : for (i = 0; i < ndaugs; i++ ) {
471 0 : masses[i] = EvtPDL::getMass(daugs[i]);
472 0 : mass_sum = mass_sum + masses[i];
473 : }
474 :
475 0 : count++;
476 :
477 :
478 0 : if(count==10000) {
479 0 : report(Severity::Error,"EvtGen") <<"Decaying particle:"<<
480 0 : EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
481 0 : report(Severity::Error,"EvtGen") <<"To the following daugthers"<<endl;
482 0 : for (i = 0; i < ndaugs; i++ ) {
483 0 : report(Severity::Error,"EvtGen") <<
484 0 : EvtPDL::name(daugs[i]).c_str() << endl;
485 : }
486 0 : report(Severity::Error,"EvtGen") << "Has been rejected "<<count
487 0 : << " times, will now take minimal masses "
488 0 : << " of daugthers"<<endl;
489 :
490 : mass_sum=0.;
491 0 : for (i = 0; i < ndaugs; i++ ) {
492 0 : masses[i] = EvtPDL::getMinMass(daugs[i]);
493 0 : mass_sum = mass_sum + masses[i];
494 : }
495 0 : if (mass_sum > p->mass()){
496 0 : report(Severity::Error,"EvtGen") << "Parent mass="<<p->mass()
497 0 : << "to light for daugthers."<<endl
498 0 : << "Will throw the event away."<<endl;
499 : //dont terminate - start over on the event.
500 0 : EvtStatus::setRejectFlag();
501 : mass_sum=0.;
502 : // ::abort();
503 0 : }
504 :
505 : }
506 0 : } while ( mass_sum > p->mass());
507 : } //else
508 :
509 0 : return;
510 0 : }
511 :
512 : void EvtDecayBase::checkNArg(int a1, int a2, int a3, int a4) {
513 :
514 0 : if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
515 0 : report(Severity::Error,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
516 0 : report(Severity::Error,"EvtGen") << a1<<endl;;
517 0 : if ( a2>-1) {
518 0 : report(Severity::Error,"EvtGen") << " or " << a2<<endl;
519 0 : }
520 0 : if ( a3>-1) {
521 0 : report(Severity::Error,"EvtGen") << " or " << a3<<endl;
522 0 : }
523 0 : if ( a4>-1) {
524 0 : report(Severity::Error,"EvtGen") << " or " << a4<<endl;
525 0 : }
526 0 : report(Severity::Error,"EvtGen") << " arguments but found:"<< _narg << endl;
527 0 : printSummary();
528 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
529 0 : ::abort();
530 :
531 : }
532 :
533 0 : }
534 : void EvtDecayBase::checkNDaug(int d1, int d2){
535 :
536 0 : if ( _ndaug != d1 && _ndaug != d2 ) {
537 0 : report(Severity::Error,"EvtGen") << _modelname.c_str() << " generator expected ";
538 0 : report(Severity::Error,"EvtGen") << d1;
539 0 : if ( d2>-1) {
540 0 : report(Severity::Error,"EvtGen") << " or " << d2;
541 0 : }
542 0 : report(Severity::Error,"EvtGen") << " daughters but found:"<< _ndaug << endl;
543 0 : printSummary();
544 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
545 0 : ::abort();
546 : }
547 :
548 0 : }
549 :
550 : void EvtDecayBase::checkSpinParent(EvtSpinType::spintype sp) {
551 :
552 0 : EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
553 0 : if ( parenttype != sp ) {
554 0 : report(Severity::Error,"EvtGen") << _modelname.c_str()
555 0 : << " did not get the correct parent spin\n";
556 0 : printSummary();
557 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
558 0 : ::abort();
559 : }
560 :
561 0 : }
562 :
563 : void EvtDecayBase::checkSpinDaughter(int d1, EvtSpinType::spintype sp) {
564 :
565 0 : EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
566 0 : if ( parenttype != sp ) {
567 0 : report(Severity::Error,"EvtGen") << _modelname.c_str()
568 0 : << " did not get the correct daughter spin d="
569 0 : << d1 << endl;
570 0 : printSummary();
571 0 : report(Severity::Error,"EvtGen") << "Will terminate execution!"<<endl;
572 0 : ::abort();
573 : }
574 :
575 0 : }
576 :
577 : double* EvtDecayBase::getArgs() {
578 :
579 0 : if ( _argsD ) return _argsD;
580 : //The user has asked for a list of doubles - the arguments
581 : //better all be doubles...
582 0 : if ( _narg==0 ) return _argsD;
583 :
584 0 : _argsD = new double[_narg];
585 :
586 : int i;
587 0 : char * tc;
588 0 : for(i=0;i<_narg;i++) {
589 0 : _argsD[i] = strtod(_args[i].c_str(),&tc);
590 : }
591 0 : return _argsD;
592 0 : }
593 :
594 : double EvtDecayBase::getArg(unsigned int j) {
595 :
596 : // Verify string
597 :
598 0 : if (getParentId().getId() == 25) {
599 : int i = 0 ;
600 : ++i;
601 0 : }
602 :
603 0 : const char* str = _args[j].c_str();
604 : int i = 0;
605 0 : while(str[i]!=0){
606 0 : if (isalpha(str[i]) && str[i]!='e') {
607 :
608 0 : report(Severity::Info,"EvtGen") << "String " << str << " is not a number" << endl;
609 0 : assert(0);
610 : }
611 0 : i++;
612 : }
613 :
614 : char** tc=0;
615 0 : double result = strtod(_args[j].c_str(),tc);
616 :
617 0 : if (_storedArgs.size() < j+1 ){ // then store the argument's value
618 0 : _storedArgs.push_back(result);
619 0 : }
620 :
621 0 : return result;
622 0 : }
623 :
624 :
625 :
626 :
627 :
628 : bool EvtDecayBase::matchingDecay(const EvtDecayBase &other) const {
629 :
630 0 : if ( _ndaug != other._ndaug) return false;
631 0 : if ( _parent != other._parent) return false;
632 :
633 0 : std::vector<int> useDs;
634 0 : for ( int i=0; i<_ndaug; i++) useDs.push_back(0);
635 :
636 0 : for ( int i=0; i<_ndaug; i++) {
637 : bool foundIt=false;
638 0 : for ( int j=0; j<_ndaug; j++) {
639 0 : if ( useDs[j] == 1 ) continue;
640 0 : if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
641 : foundIt=true;
642 0 : useDs[j]=1;
643 0 : break;
644 : }
645 : }
646 0 : if ( foundIt==false) return false;
647 0 : }
648 0 : for ( int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
649 :
650 0 : return true;
651 :
652 0 : }
|
