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