]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TEvtGen/EvtGenBase/EvtDecayBase.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
L1phase shift corrected
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtDecayBase.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: 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>
33using std::endl;
34using std::fstream;
35void 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
67double 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
108double 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
129std::string EvtDecayBase::commandName(){
130 return std::string("");
131}
132
133void EvtDecayBase::command(std::string){
134 report(ERROR,"EvtGen") << "Should never call EvtDecayBase::command"<<endl;
135 ::abort();
136}
137
138
139
140void EvtDecayBase::init() {
141
142 //This default version of init does nothing;
143 //A specialized version of this function can be
144 //supplied for each decay model to do initialization.
145
146 return;
147
148}
149
150void EvtDecayBase::initProbMax() {
151
152 //This function is called if the decay does not have a
153 //specialized initialization.
154 //The default is to set the maximum
155 //probability to 0 and the number of times called to 0
156 //and defaultprobmax to 1 such that the decay will be
157 //generated many many times
158 //in order to generate a reasonable maximum probability
159 //for the decay.
160
161 defaultprobmax=1;
162 ntimes_prob = 0;
163 probmax = 0.0;
164
165} //initProbMax
166
167
168void EvtDecayBase::saveDecayInfo(EvtId ipar, int ndaug, EvtId *daug,
169 int narg,std::vector<std::string>& args,
170 std::string name,
171 double brfr) {
172
173 int i;
174
175 _brfr=brfr;
176 _ndaug=ndaug;
177 _narg=narg;
178 _parent=ipar;
179
180 _dsum=0;
181
182 if (_ndaug>0) {
183 _daug=new EvtId [_ndaug];
184 for(i=0;i<_ndaug;i++){
185 _daug[i]=daug[i];
186 _dsum+=daug[i].getAlias();
187 }
188 }
189 else{
190 _daug=0;
191 }
192
193 if (_narg>0) {
194 _args=new std::string[_narg+1];
195 for(i=0;i<_narg;i++){
196 _args[i]=args[i];
197 }
198 }
199 else{
200 _args = 0;
201 }
202
203 _modelname=name;
204
205 this->init();
206 this->initProbMax();
207
208 if (_chkCharge){
209 this->checkQ();
210 }
211
212
213 if (defaultprobmax){
214 report(INFO,"EvtGen") << "No default probmax for ";
215 report(INFO,"") << "("<<_modelname.c_str()<<") ";
216 report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
217 for(i=0;i<_ndaug;i++){
218 report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
219 }
220 report(INFO,"") << endl;
221 report(INFO,"") << "This is fine for development, but must be provided for production."<<endl;
222 report(INFO,"EvtGen") << "Never fear though - the decay will use the \n";
223 report(INFO,"EvtGen") << "500 iterations to build up a good probmax \n";
224 report(INFO,"EvtGen") << "before accepting a decay. "<<endl;
225 }
226
227}
228
229EvtDecayBase::EvtDecayBase() {
230
231 //the default is that the user module does _not_ set
232 // any probmax.
233 defaultprobmax=1;
234 ntimes_prob = 0;
235 probmax = 0.0;
236
237 _photos=0;
238 _verbose=0;
239 _summary=0;
240 _parent=EvtId(-1,-1);
241 _ndaug=0;
242 _narg=0;
243 _daug=0;
244 _args=0;
245 _argsD=0;
246 _modelname="**********";
247
248 //Default is to check that charge is conserved
249 _chkCharge=1;
250
251 //statistics collection!
252
253 max_prob=0.0;
254 sum_prob=0.0;
255
256}
257
258
259
260void EvtDecayBase::printSummary() const {
261 if (ntimes_prob>0) {
262
263 report(INFO,"EvtGen") << "Calls = "<<ntimes_prob<<" eff: "<<
264 sum_prob/(probmax*ntimes_prob)<<" frac. max:"<<max_prob/probmax;
265 report(INFO,"") <<" probmax:"<<probmax<<" max:"<<max_prob<<" : ";
266 }
267
268 printInfo();
269}
270
271
272void EvtDecayBase::printInfo() const {
273 report(INFO,"") << EvtPDL::name(_parent).c_str()<<" -> ";
274 for(int i=0;i<_ndaug;i++){
275 report(INFO,"") << EvtPDL::name(_daug[i]).c_str() << " ";
276 }
277 report(INFO,"") << " ("<<_modelname.c_str()<<")"<< endl;
278}
279
280
281EvtDecayBase::~EvtDecayBase() {
282
283 if (_daug!=0){
284 delete [] _daug;
285 }
286
287 if (_args!=0){
288 delete [] _args;
289 }
290
291 if (_argsD!=0){
292 delete [] _argsD;
293 }
294
295
296}
297
298void EvtDecayBase::setProbMax(double prbmx){
299
300 defaultprobmax=0;
301 probmax=prbmx;
302
303}
304
305void EvtDecayBase::noProbMax(){
306
307 defaultprobmax=0;
308
309}
310
311
312double EvtDecayBase::findMaxMass(EvtParticle *p) {
313
314
315 double maxOkMass=EvtPDL::getMaxMass(p->getId());
316
317 //protect against vphotons
318 if ( maxOkMass < 0.0000000001 ) return 10000000.;
319 //and against already determined masses
320 if ( p->hasValidP4() ) maxOkMass=p->mass();
321
322 EvtParticle *par=p->getParent();
323 if ( par ) {
324 double maxParMass=findMaxMass(par);
325 size_t i;
326 double minDaugMass=0.;
327 for(i=0;i<par->getNDaug();i++){
328 EvtParticle *dau=par->getDaug(i);
329 if ( dau!=p) {
330 // it might already have a mass
331 if ( dau->isInitialized() || dau->hasValidP4() )
332 minDaugMass+=dau->mass();
333 else
334 //give it a bit of phase space
335 minDaugMass+=1.000001*EvtPDL::getMinMass(dau->getId());
336 }
337 }
338 if ( maxOkMass>(maxParMass-minDaugMass)) maxOkMass=maxParMass-minDaugMass;
339 }
340 return maxOkMass;
341}
342
343
344// given list of daughters ( by number ) returns a
345// list of viable masses.
346
347void EvtDecayBase::findMass(EvtParticle *p) {
348
349 //Need to also check that this mass does not screw
350 //up the parent
351 //This code assumes that for the ith daughter, 0..i-1
352 //already have a mass
353 double maxOkMass=findMaxMass(p);
354
355 int count=0;
356 double mass;
357 bool massOk=false;
358 size_t i;
359 while (!massOk) {
360 count++;
361 if ( count > 10000 ) {
362 report(INFO,"EvtGen") << "Can not find a valid mass for: " << EvtPDL::name(p->getId()).c_str() <<endl;
363 report(INFO,"EvtGen") << "Now printing parent and/or grandparent tree\n";
364 if ( p->getParent() ) {
365 if ( p->getParent()->getParent() ) {
366 p->getParent()->getParent()->printTree();
367 report(INFO,"EvtGen") << p->getParent()->getParent()->mass() <<endl;
368 report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
369 }
370 else{
371 p->getParent()->printTree();
372 report(INFO,"EvtGen") << p->getParent()->mass() <<endl;
373 }
374 }
375 else p->printTree();
376 report(INFO,"EvtGen") << "maxokmass=" << maxOkMass << " " << EvtPDL::getMinMass(p->getId()) << " " << EvtPDL::getMaxMass(p->getId())<<endl;
377 if ( p->getNDaug() ) {
378 for (i=0; i<p->getNDaug(); i++) {
379 report(INFO,"EvtGen") << p->getDaug(i)->mass()<<" ";
380 }
381 report(INFO,"EvtGen") << endl;
382 }
383 if ( maxOkMass >= EvtPDL::getMinMass(p->getId()) ) {
384 report(INFO,"EvtGen") << "taking a default value\n";
385 p->setMass(maxOkMass);
386 return;
387 }
388 assert(0);
389 }
390 mass = EvtPDL::getMass(p->getId());
391 //Just need to check that this mass is > than
392 //the mass of all daughters
393 double massSum=0.;
394 if ( p->getNDaug() ) {
395 for (i=0; i<p->getNDaug(); i++) {
396 massSum+= p->getDaug(i)->mass();
397 }
398 }
399 //some special cases are handled with 0 (stable) or 1 (k0->ks/kl) daughters
400 if (p->getNDaug()<2) massOk=true;
401 if ( p->getParent() ) {
402 if ( p->getParent()->getNDaug()==1 ) massOk=true;
403 }
404 if ( !massOk ) {
405 if (massSum < mass) massOk=true;
406 if ( mass> maxOkMass) massOk=false;
407 }
408 }
409
410 p->setMass(mass);
411
412}
413
414
415void EvtDecayBase::findMasses(EvtParticle *p, int ndaugs,
416 EvtId daugs[10], double masses[10]) {
417
418 int i;
419 double mass_sum;
420
421 int count=0;
422
423 if (!( p->firstornot() )) {
424 for (i = 0; i < ndaugs; i++ ) {
425 masses[i] = p->getDaug(i)->mass();
426 } //for
427 } //if
428 else {
429 p->setFirstOrNot();
430 // if only one daughter do it
431
432 if (ndaugs==1) {
433 masses[0]=p->mass();
434 return;
435 }
436
437 //until we get a combo whose masses are less than _parent mass.
438 do {
439 mass_sum = 0.0;
440
441 for (i = 0; i < ndaugs; i++ ) {
442 masses[i] = EvtPDL::getMass(daugs[i]);
443 mass_sum = mass_sum + masses[i];
444 }
445
446 count++;
447
448
449 if(count==10000) {
450 report(ERROR,"EvtGen") <<"Decaying particle:"<<
451 EvtPDL::name(p->getId()).c_str()<<" (m="<<p->mass()<<")"<<endl;
452 report(ERROR,"EvtGen") <<"To the following daugthers"<<endl;
453 for (i = 0; i < ndaugs; i++ ) {
454 report(ERROR,"EvtGen") <<
455 EvtPDL::name(daugs[i]).c_str() << endl;
456 }
457 report(ERROR,"EvtGen") << "Has been rejected "<<count
458 << " times, will now take minimal masses "
459 << " of daugthers"<<endl;
460
461 mass_sum=0.;
462 for (i = 0; i < ndaugs; i++ ) {
463 masses[i] = EvtPDL::getMinMass(daugs[i]);
464 mass_sum = mass_sum + masses[i];
465 }
466 if (mass_sum > p->mass()){
467 report(ERROR,"EvtGen") << "Parent mass="<<p->mass()
468 << "to light for daugthers."<<endl
469 << "Will throw the event away."<<endl;
470 //dont terminate - start over on the event.
471 EvtStatus::setRejectFlag();
472 mass_sum=0.;
473 // ::abort();
474 }
475
476 }
477 } while ( mass_sum > p->mass());
478 } //else
479
480 return;
481}
482
483void EvtDecayBase::checkNArg(int a1, int a2, int a3, int a4) {
484
485 if ( _narg != a1 && _narg != a2 && _narg != a3 && _narg != a4 ) {
486 report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected "<<endl;
487 report(ERROR,"EvtGen") << a1<<endl;;
488 if ( a2>-1) {
489 report(ERROR,"EvtGen") << " or " << a2<<endl;
490 }
491 if ( a3>-1) {
492 report(ERROR,"EvtGen") << " or " << a3<<endl;
493 }
494 if ( a4>-1) {
495 report(ERROR,"EvtGen") << " or " << a4<<endl;
496 }
497 report(ERROR,"EvtGen") << " arguments but found:"<< _narg << endl;
498 printSummary();
499 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
500 ::abort();
501
502 }
503
504}
505void EvtDecayBase::checkNDaug(int d1, int d2){
506
507 if ( _ndaug != d1 && _ndaug != d2 ) {
508 report(ERROR,"EvtGen") << _modelname.c_str() << " generator expected ";
509 report(ERROR,"EvtGen") << d1;
510 if ( d2>-1) {
511 report(ERROR,"EvtGen") << " or " << d2;
512 }
513 report(ERROR,"EvtGen") << " daughters but found:"<< _ndaug << endl;
514 printSummary();
515 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
516 ::abort();
517 }
518
519}
520
521void EvtDecayBase::checkSpinParent(EvtSpinType::spintype sp) {
522
523 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
524 if ( parenttype != sp ) {
525 report(ERROR,"EvtGen") << _modelname.c_str()
526 << " did not get the correct parent spin\n";
527 printSummary();
528 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
529 ::abort();
530 }
531
532}
533
534void EvtDecayBase::checkSpinDaughter(int d1, EvtSpinType::spintype sp) {
535
536 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getDaug(d1));
537 if ( parenttype != sp ) {
538 report(ERROR,"EvtGen") << _modelname.c_str()
539 << " did not get the correct daughter spin d="
540 << d1 << endl;
541 printSummary();
542 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
543 ::abort();
544 }
545
546}
547
548double* EvtDecayBase::getArgs() {
549
550 if ( _argsD ) return _argsD;
551 //The user has asked for a list of doubles - the arguments
552 //better all be doubles...
553 if ( _narg==0 ) return _argsD;
554
555 _argsD = new double[_narg];
556
557 int i;
558 char * tc;
559 for(i=0;i<_narg;i++) {
560 _argsD[i] = strtod(_args[i].c_str(),&tc);
561 }
562 return _argsD;
563}
564
565double EvtDecayBase::getArg(unsigned int j) {
566
567 // Verify string
568
569 if (getParentId().getId() == 25) {
570 int i;
571 ++i;
572 }
573
574 const char* str = _args[j].c_str();
575 int i = 0;
576 while(str[i]!=0){
577 if (isalpha(str[i]) && str[i]!='e') {
578
579 report(INFO,"EvtGen") << "String " << str << " is not a number" << endl;
580 assert(0);
581 }
582 i++;
583 }
584
585 char** tc=0;
586 double result = strtod(_args[j].c_str(),tc);
587
588 if (_storedArgs.size() < j+1 ){ // then store the argument's value
589 _storedArgs.push_back(result);
590 }
591
592 return result;
593}
594
595
596
597
598
599bool EvtDecayBase::matchingDecay(const EvtDecayBase &other) const {
600
601 if ( _ndaug != other._ndaug) return false;
602 if ( _parent != other._parent) return false;
603
604 std::vector<int> useDs;
605 for ( int i=0; i<_ndaug; i++) useDs.push_back(0);
606
607 for ( int i=0; i<_ndaug; i++) {
608 bool foundIt=false;
609 for ( int j=0; j<_ndaug; j++) {
610 if ( useDs[j] == 1 ) continue;
611 if ( _daug[i] == other._daug[j] && _daug[i].getAlias() == other._daug[j].getAlias()) {
612 foundIt=true;
613 useDs[j]=1;
614 break;
615 }
616 }
617 if ( foundIt==false) return false;
618 }
619 for ( int i=0; i<_ndaug; i++) if ( useDs[i]==0) return false;
620
621 return true;
622
623}