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