Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtParticleDecayList.cpp
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: EvtDecayParm.cc
12//
13// Description: Store decay parameters for one decay.
14//
15// Modification history:
16//
17// RYD April 5, 1997 Module created
18//
19//------------------------------------------------------------------------
20//
21#include "EvtGenBase/EvtPatches.hh"
22#include "EvtGenBase/EvtPatches.hh"
23#include <iostream>
24#include <fstream>
25#include <stdlib.h>
26#include <ctype.h>
27#include "EvtGenBase/EvtParticleDecayList.hh"
28#include "EvtGenBase/EvtParticle.hh"
29#include "EvtGenBase/EvtRandom.hh"
30#include "EvtGenBase/EvtReport.hh"
31#include "EvtGenBase/EvtPDL.hh"
32#include "EvtGenBase/EvtStatus.hh"
33using std::endl;
34using std::fstream;
35
36EvtParticleDecayList::EvtParticleDecayList(const EvtParticleDecayList &o) {
37 _nmode=o._nmode;
38 _rawbrfrsum=o._rawbrfrsum;
39 _decaylist=new EvtParticleDecayPtr[_nmode];
40
41 int i;
42 for(i=0;i<_nmode;i++){
43 _decaylist[i]=new EvtParticleDecay;
44
45 EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
46
47 EvtDecayBase *tModelNew=tModel->clone();
48 if (tModel->getPHOTOS()){
49 tModelNew->setPHOTOS();
50 }
51 if (tModel->verbose()){
52 tModelNew->setVerbose();
53 }
54 if (tModel->summary()){
55 tModelNew->setSummary();
56 }
57 std::vector<std::string> args;
58 int j;
59 for(j=0;j<tModel->getNArg();j++){
60 args.push_back(tModel->getArgStr(j));
61 }
62 tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
63 tModel->getDaugs(),
64 tModel->getNArg(),
65 args,
66 tModel->getModelName(),
67 tModel->getBranchingFraction());
68 _decaylist[i]->setDecayModel(tModelNew);
69
70 _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
71 _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
72 }
73
74
75}
76
77
78EvtParticleDecayList::~EvtParticleDecayList(){
79
80 int i;
81 for(i=0;i<_nmode;i++){
82 delete _decaylist[i];
83 }
84
85 if (_decaylist!=0) delete [] _decaylist;
86}
87
88void EvtParticleDecayList::printSummary(){
89
90 int i;
91 for(i=0;i<_nmode;i++){
92 _decaylist[i]->printSummary();
93 }
94
95}
96
97void EvtParticleDecayList::removeDecay(){
98
99 int i;
100 for(i=0;i<_nmode;i++){
101 delete _decaylist[i];
102 }
103
104 delete [] _decaylist;
105 _decaylist=0;
106 _nmode=0;
107 _rawbrfrsum=0.0;
108
109}
110
0ca57c2f 111EvtDecayBase* EvtParticleDecayList::getDecayModel(int imode) {
112
113 EvtDecayBase* theModel(0);
114 if (imode >= 0 && imode < _nmode) {
115 EvtParticleDecay* theDecay = _decaylist[imode];
116 if (theDecay != 0) {
117 theModel = theDecay->getDecayModel();
118 }
119 }
120
121 return theModel;
122
123}
da0e9ce3 124
125EvtDecayBase* EvtParticleDecayList::getDecayModel(EvtParticle *p){
126
127 if (p->getNDaug()!=0) {
128 assert(p->getChannel()>=0);
129 return getDecay(p->getChannel()).getDecayModel();
130 }
131 if (p->getChannel() >(-1) ) {
132 return getDecay(p->getChannel()).getDecayModel();
133 }
134
135
136 if (getNMode()==0 ) {
137 return 0;
138 }
139 if (getRawBrfrSum()<0.00000001 ) {
140 return 0;
141 }
142
143 if (getNMode()==1) {
144 p->setChannel(0);
145 return getDecay(0).getDecayModel();
146 }
147
148 if (p->getChannel() >(-1)) {
149 report(ERROR,"EvtGen") << "Internal error!!!"<<endl;
150 ::abort();
151 }
152
153 int j;
154
155 for (j=0;j<10000000;j++){
156
157 double u=EvtRandom::Flat();
158
159 int i;
160 bool breakL=false;
161 for (i=0;i<getNMode();i++) {
162
163 if ( breakL ) continue;
164 if (u<getDecay(i).getBrfrSum()) {
165 breakL=true;
166 //special case for decay of on particel to another
167 // e.g. K0->K0S
168
169 if (getDecay(i).getDecayModel()->getNDaug()==1 ) {
170 p->setChannel(i);
171 return getDecay(i).getDecayModel();
172 }
173
174 if ( p->hasValidP4() ) {
175 if (getDecay(i).getMassMin() < p->mass() ) {
176 p->setChannel(i);
177 return getDecay(i).getDecayModel();
178 }
179 }
180 else{
181 //Lange apr29-2002 - dont know the mass yet
182 p->setChannel(i);
183 return getDecay(i).getDecayModel();
184 }
185 }
186 }
187 }
188
189 //Ok, we tried 10000000 times above to pick a decay channel that is
190 //kinematically allowed! Now we give up and search all channels!
191 //if that fails, the particle will not be decayed!
192
193 report(ERROR,"EvtGen") << "Tried 10000000 times to generate decay of "
194 << EvtPDL::name(p->getId())<< " with mass="<<p->mass()<<endl;
195 report(ERROR,"EvtGen") << "Will take first kinematically allowed decay in the decay table"
196 <<endl;
197
198 int i;
199
200 //Need to check that we don't use modes with 0 branching fractions.
201 double previousBrSum=0.0;
202 for (i=0;i<getNMode();i++) {
203 if(getDecay(i).getBrfrSum()!=previousBrSum){
204 if ( getDecay(i).getMassMin() < p->mass() ) {
205 p->setChannel(i);
206 return getDecay(i).getDecayModel();
207 }
208 }
209 previousBrSum=getDecay(i).getBrfrSum();
210 }
211
212 report(ERROR,"EvtGen") << "Could not decay:"
213 <<EvtPDL::name(p->getId()).c_str()
214 <<" with mass:"<<p->mass()
215 <<" will throw event away! "<<endl;
216
217 EvtStatus::setRejectFlag();
218 return 0;
219
220}
221
222
223void EvtParticleDecayList::setNMode(int nmode){
224
225 EvtParticleDecayPtr* _decaylist_new= new EvtParticleDecayPtr[nmode];
226
227 if (_nmode!=0){
228 report(ERROR,"EvtGen") << "Error _nmode not equal to zero!!!"<<endl;
229 ::abort();
230 }
231 if (_decaylist!=0) {
232 delete [] _decaylist;
233 }
234 _decaylist=_decaylist_new;
235 _nmode=nmode;
236
237}
238
239
0ca57c2f 240EvtParticleDecay& EvtParticleDecayList::getDecay(int nchannel) const {
da0e9ce3 241 if (nchannel>=_nmode) {
242 report(ERROR,"EvtGen") <<"Error getting channel:"
243 <<nchannel<<" with only "<<_nmode
244 <<" stored!"<<endl;
245 ::abort();
246 }
247 return *(_decaylist[nchannel]);
248}
249
250void EvtParticleDecayList::makeChargeConj(EvtParticleDecayList* conjDecayList){
251
252 _rawbrfrsum=conjDecayList->_rawbrfrsum;
253
254 setNMode(conjDecayList->_nmode);
255
256 int i;
257
258 for(i=0;i<_nmode;i++){
259 _decaylist[i]=new EvtParticleDecay;
260 _decaylist[i]->chargeConj(conjDecayList->_decaylist[i]);
261 }
262
263}
264
265void EvtParticleDecayList::addMode(EvtDecayBase* decay, double brfrsum,
266 double massmin){
267
268 EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode+1];
269
270 int i;
271 for(i=0;i<_nmode;i++){
272 newlist[i]=_decaylist[i];
273 }
274
275 _rawbrfrsum=brfrsum;
276
277 newlist[_nmode]=new EvtParticleDecay;
278
279 newlist[_nmode]->setDecayModel(decay);
280 newlist[_nmode]->setBrfrSum(brfrsum);
281 newlist[_nmode]->setMassMin(massmin);
282
283 EvtDecayBase *newDec=newlist[_nmode]->getDecayModel();
284 for(i=0;i<_nmode;i++){
285 if ( newDec->matchingDecay(*(newlist[i]->getDecayModel())) ) {
286
287 //sometimes its ok..
288 if ( newDec->getModelName() == "JETSET" || newDec->getModelName() == "PYTHIA" ) continue;
289 if ( newDec->getModelName() == "JSCONT" || newDec->getModelName() == "PYCONT" ) continue;
290 if ( newDec->getModelName() == "PYGAGA" ) continue;
291 if ( newDec->getModelName() == "LUNDAREALAW" ) continue;
0ca57c2f 292 if ( newDec->getModelName() == "TAUOLA") continue;
da0e9ce3 293 report(ERROR,"EvtGen") << "Two matching decays with same parent in decay table\n";
294 report(ERROR,"EvtGen") << "Please fix that\n";
295 report(ERROR,"EvtGen") << "Parent " << EvtPDL::name(newDec->getParentId()).c_str() << endl;
296 for (int j=0; j<newDec->getNDaug(); j++)
297 report(ERROR,"EvtGen") << "Daughter " << EvtPDL::name(newDec->getDaug(j)).c_str() << endl;
298 assert(0);
299 }
300 }
301
302 if (_nmode!=0){
303 delete [] _decaylist;
304 }
305
306 if ( ( _nmode == 0 ) && ( _decaylist != 0 ) ) delete [] _decaylist ;
307
308 _nmode++;
309
310 _decaylist=newlist;
311
312}
313
314
315void EvtParticleDecayList::finalize(){
316
317 if (_nmode>0) {
318 if ( _rawbrfrsum< 0.000001 ) {
319 report(ERROR,"EvtGen") << "Please give me a "
320 << "branching fraction sum greater than 0\n";
321 assert(0);
322 }
323 if (fabs(_rawbrfrsum-1.0)>0.0001) {
324 report(INFO,"EvtGen") <<"Warning, sum of branching fractions for "
325 <<EvtPDL::name(_decaylist[0]->getDecayModel()->getParentId()).c_str()
326 <<" is "<<_rawbrfrsum<<endl;
327 report(INFO,"EvtGen") << "rescaled to one! "<<endl;
328
329 }
330
331 int i;
332
333 for(i=0;i<_nmode;i++){
334 double brfrsum=_decaylist[i]->getBrfrSum()/_rawbrfrsum;
335 _decaylist[i]->setBrfrSum(brfrsum);
336 }
337
338 }
339
340}
341
342
343EvtParticleDecayList& EvtParticleDecayList::operator=(const EvtParticleDecayList &o) {
344 if (this != &o) {
345 removeDecay();
346 _nmode=o._nmode;
347 _rawbrfrsum=o._rawbrfrsum;
348 _decaylist=new EvtParticleDecayPtr[_nmode];
349
350 int i;
351 for(i=0;i<_nmode;i++){
352 _decaylist[i]=new EvtParticleDecay;
353
354 EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
355
356 EvtDecayBase *tModelNew=tModel->clone();
357 if (tModel->getPHOTOS()){
358 tModelNew->setPHOTOS();
359 }
360 if (tModel->verbose()){
361 tModelNew->setVerbose();
362 }
363 if (tModel->summary()){
364 tModelNew->setSummary();
365 }
366 std::vector<std::string> args;
367 int j;
368 for(j=0;j<tModel->getNArg();j++){
369 args.push_back(tModel->getArgStr(j));
370 }
371 tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
372 tModel->getDaugs(),
373 tModel->getNArg(),
374 args,
375 tModel->getModelName(),
376 tModel->getBranchingFraction());
377 _decaylist[i]->setDecayModel(tModelNew);
378
379 //_decaylist[i]->setDecayModel(tModel);
380 _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
381 _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
382 }
383 }
384 return *this;
385
386
387}
388
389
390void EvtParticleDecayList::removeMode(EvtDecayBase* decay) {
391 // here we will delete a decay with the same final state particles
392 // and recalculate the branching fractions for the remaining modes
393 int match = -1;
394 int i;
395 double match_bf;
396
397 for(i=0;i<_nmode;i++){
398 if ( decay->matchingDecay(*(_decaylist[i]->getDecayModel())) ) {
399 match = i;
400 }
401 }
402
403 if (match < 0) {
404 report(ERROR,"EvtGen") << " Attempt to remove undefined mode for" << endl
405 << "Parent " << EvtPDL::name(decay->getParentId()).c_str() << endl
406 << "Daughters: ";
407 for (int j=0; j<decay->getNDaug(); j++)
408 report(ERROR,"") << EvtPDL::name(decay->getDaug(j)).c_str() << " ";
409 report(ERROR,"") << endl;
410 ::abort();
411 }
412
413 if (match == 0) {
414 match_bf = _decaylist[match]->getBrfrSum();
415 } else {
416 match_bf = (_decaylist[match]->getBrfrSum()
417 -_decaylist[match-1]->getBrfrSum());
418 }
419
420 double divisor = 1-match_bf;
421 if (divisor < 0.000001 && _nmode > 1) {
422 report(ERROR,"EvtGen") << "Removing requested mode leaves "
423 << EvtPDL::name(decay->getParentId()).c_str()
424 << " with zero sum branching fraction," << endl
425 << "but more than one decay mode remains. Aborting."
426 << endl;
427 ::abort();
428 }
429
430 EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode-1];
431
432 for(i=0;i<match;i++){
433 newlist[i]=_decaylist[i];
434 newlist[i]->setBrfrSum(newlist[i]->getBrfrSum()/divisor);
435 }
436 for(i=match+1; i<_nmode; i++) {
437 newlist[i-1]=_decaylist[i];
438 newlist[i-1]->setBrfrSum((newlist[i-1]->getBrfrSum()-match_bf)/divisor);
439 }
440
441
442 delete [] _decaylist;
443
444 _nmode--;
445
446 _decaylist=newlist;
447
448 if (_nmode == 0) {
449 delete [] _decaylist;
450 }
451
452}
453
454
0ca57c2f 455bool EvtParticleDecayList::isJetSet() const {
da0e9ce3 456 int i ;
457 EvtDecayBase * decayer ;
458
459 for ( i = 0 ;
460 i < getNMode() ;
461 i++ ) {
462 decayer = getDecay( i ).getDecayModel ( ) ;
463 if ( decayer -> getModelName() == "PYTHIA" ) return true ;
464 }
465
466 return false ;
467}