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