1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtDecayParm.cc
13 // Description: Store decay parameters for one decay.
15 // Modification history:
17 // RYD April 5, 1997 Module created
19 //------------------------------------------------------------------------
21 #include "EvtGenBase/EvtPatches.hh"
22 #include "EvtGenBase/EvtPatches.hh"
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"
36 EvtParticleDecayList::EvtParticleDecayList(const EvtParticleDecayList &o) {
38 _rawbrfrsum=o._rawbrfrsum;
39 _decaylist=new EvtParticleDecayPtr[_nmode];
42 for(i=0;i<_nmode;i++){
43 _decaylist[i]=new EvtParticleDecay;
45 EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
47 EvtDecayBase *tModelNew=tModel->clone();
48 if (tModel->getPHOTOS()){
49 tModelNew->setPHOTOS();
51 if (tModel->verbose()){
52 tModelNew->setVerbose();
54 if (tModel->summary()){
55 tModelNew->setSummary();
57 std::vector<std::string> args;
59 for(j=0;j<tModel->getNArg();j++){
60 args.push_back(tModel->getArgStr(j));
62 tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
66 tModel->getModelName(),
67 tModel->getBranchingFraction());
68 _decaylist[i]->setDecayModel(tModelNew);
70 _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
71 _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
78 EvtParticleDecayList::~EvtParticleDecayList(){
81 for(i=0;i<_nmode;i++){
85 if (_decaylist!=0) delete [] _decaylist;
88 void EvtParticleDecayList::printSummary(){
91 for(i=0;i<_nmode;i++){
92 _decaylist[i]->printSummary();
97 void EvtParticleDecayList::removeDecay(){
100 for(i=0;i<_nmode;i++){
101 delete _decaylist[i];
104 delete [] _decaylist;
112 EvtDecayBase* EvtParticleDecayList::getDecayModel(EvtParticle *p){
114 if (p->getNDaug()!=0) {
115 assert(p->getChannel()>=0);
116 return getDecay(p->getChannel()).getDecayModel();
118 if (p->getChannel() >(-1) ) {
119 return getDecay(p->getChannel()).getDecayModel();
123 if (getNMode()==0 ) {
126 if (getRawBrfrSum()<0.00000001 ) {
132 return getDecay(0).getDecayModel();
135 if (p->getChannel() >(-1)) {
136 report(ERROR,"EvtGen") << "Internal error!!!"<<endl;
142 for (j=0;j<10000000;j++){
144 double u=EvtRandom::Flat();
148 for (i=0;i<getNMode();i++) {
150 if ( breakL ) continue;
151 if (u<getDecay(i).getBrfrSum()) {
153 //special case for decay of on particel to another
156 if (getDecay(i).getDecayModel()->getNDaug()==1 ) {
158 return getDecay(i).getDecayModel();
161 if ( p->hasValidP4() ) {
162 if (getDecay(i).getMassMin() < p->mass() ) {
164 return getDecay(i).getDecayModel();
168 //Lange apr29-2002 - dont know the mass yet
170 return getDecay(i).getDecayModel();
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!
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"
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() ) {
193 return getDecay(i).getDecayModel();
196 previousBrSum=getDecay(i).getBrfrSum();
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;
204 EvtStatus::setRejectFlag();
210 void EvtParticleDecayList::setNMode(int nmode){
212 EvtParticleDecayPtr* _decaylist_new= new EvtParticleDecayPtr[nmode];
215 report(ERROR,"EvtGen") << "Error _nmode not equal to zero!!!"<<endl;
219 delete [] _decaylist;
221 _decaylist=_decaylist_new;
227 EvtParticleDecay& EvtParticleDecayList::getDecay(int nchannel) {
228 if (nchannel>=_nmode) {
229 report(ERROR,"EvtGen") <<"Error getting channel:"
230 <<nchannel<<" with only "<<_nmode
234 return *(_decaylist[nchannel]);
237 void EvtParticleDecayList::makeChargeConj(EvtParticleDecayList* conjDecayList){
239 _rawbrfrsum=conjDecayList->_rawbrfrsum;
241 setNMode(conjDecayList->_nmode);
245 for(i=0;i<_nmode;i++){
246 _decaylist[i]=new EvtParticleDecay;
247 _decaylist[i]->chargeConj(conjDecayList->_decaylist[i]);
252 void EvtParticleDecayList::addMode(EvtDecayBase* decay, double brfrsum,
255 EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode+1];
258 for(i=0;i<_nmode;i++){
259 newlist[i]=_decaylist[i];
264 newlist[_nmode]=new EvtParticleDecay;
266 newlist[_nmode]->setDecayModel(decay);
267 newlist[_nmode]->setBrfrSum(brfrsum);
268 newlist[_nmode]->setMassMin(massmin);
270 EvtDecayBase *newDec=newlist[_nmode]->getDecayModel();
271 for(i=0;i<_nmode;i++){
272 if ( newDec->matchingDecay(*(newlist[i]->getDecayModel())) ) {
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;
289 delete [] _decaylist;
292 if ( ( _nmode == 0 ) && ( _decaylist != 0 ) ) delete [] _decaylist ;
301 void EvtParticleDecayList::finalize(){
304 if ( _rawbrfrsum< 0.000001 ) {
305 report(ERROR,"EvtGen") << "Please give me a "
306 << "branching fraction sum greater than 0\n";
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;
319 for(i=0;i<_nmode;i++){
320 double brfrsum=_decaylist[i]->getBrfrSum()/_rawbrfrsum;
321 _decaylist[i]->setBrfrSum(brfrsum);
329 EvtParticleDecayList& EvtParticleDecayList::operator=(const EvtParticleDecayList &o) {
333 _rawbrfrsum=o._rawbrfrsum;
334 _decaylist=new EvtParticleDecayPtr[_nmode];
337 for(i=0;i<_nmode;i++){
338 _decaylist[i]=new EvtParticleDecay;
340 EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
342 EvtDecayBase *tModelNew=tModel->clone();
343 if (tModel->getPHOTOS()){
344 tModelNew->setPHOTOS();
346 if (tModel->verbose()){
347 tModelNew->setVerbose();
349 if (tModel->summary()){
350 tModelNew->setSummary();
352 std::vector<std::string> args;
354 for(j=0;j<tModel->getNArg();j++){
355 args.push_back(tModel->getArgStr(j));
357 tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
361 tModel->getModelName(),
362 tModel->getBranchingFraction());
363 _decaylist[i]->setDecayModel(tModelNew);
365 //_decaylist[i]->setDecayModel(tModel);
366 _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
367 _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
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
383 for(i=0;i<_nmode;i++){
384 if ( decay->matchingDecay(*(_decaylist[i]->getDecayModel())) ) {
390 report(ERROR,"EvtGen") << " Attempt to remove undefined mode for" << endl
391 << "Parent " << EvtPDL::name(decay->getParentId()).c_str() << endl
393 for (int j=0; j<decay->getNDaug(); j++)
394 report(ERROR,"") << EvtPDL::name(decay->getDaug(j)).c_str() << " ";
395 report(ERROR,"") << endl;
400 match_bf = _decaylist[match]->getBrfrSum();
402 match_bf = (_decaylist[match]->getBrfrSum()
403 -_decaylist[match-1]->getBrfrSum());
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."
416 EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode-1];
418 for(i=0;i<match;i++){
419 newlist[i]=_decaylist[i];
420 newlist[i]->setBrfrSum(newlist[i]->getBrfrSum()/divisor);
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);
428 delete [] _decaylist;
435 delete [] _decaylist;
441 bool EvtParticleDecayList::isJetSet() {
443 EvtDecayBase * decayer ;
448 decayer = getDecay( i ).getDecayModel ( ) ;
449 if ( decayer -> getModelName() == "PYTHIA" ) return true ;