[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtParticleDecayList.cxx

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtDecayParm.cc
//
// Description: Store decay parameters for one decay.
//
// Modification history:
//
//    RYD     April 5, 1997         Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <ctype.h>
#include "EvtGenBase/EvtParticleDecayList.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtStatus.hh"
using std::endl;
using std::fstream;

EvtParticleDecayList::EvtParticleDecayList(const EvtParticleDecayList &o) {
  _nmode=o._nmode;
  _rawbrfrsum=o._rawbrfrsum;
  _decaylist=new EvtParticleDecayPtr[_nmode];

  int i;
  for(i=0;i<_nmode;i++){
    _decaylist[i]=new EvtParticleDecay;  

    EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
    
    EvtDecayBase *tModelNew=tModel->clone();
    if (tModel->getPHOTOS()){
      tModelNew->setPHOTOS();
    }
    if (tModel->verbose()){
      tModelNew->setVerbose();
    }
    if (tModel->summary()){
      tModelNew->setSummary();
    }
    std::vector<std::string> args;
    int j;
    for(j=0;j<tModel->getNArg();j++){
      args.push_back(tModel->getArgStr(j));
    }
    tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
			     tModel->getDaugs(),
			     tModel->getNArg(),
			     args,
			     tModel->getModelName(),
			     tModel->getBranchingFraction());
    _decaylist[i]->setDecayModel(tModelNew);
    
    _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
    _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
  }


}


EvtParticleDecayList::~EvtParticleDecayList(){

  int i;
  for(i=0;i<_nmode;i++){
    delete _decaylist[i];
  }
  
  if (_decaylist!=0) delete [] _decaylist;
}

void EvtParticleDecayList::printSummary(){

  int i;
  for(i=0;i<_nmode;i++){
    _decaylist[i]->printSummary();
  }
  
}

void EvtParticleDecayList::removeDecay(){
  
  int i;
  for(i=0;i<_nmode;i++){
    delete _decaylist[i];
  }
  
  delete [] _decaylist; 
  _decaylist=0; 
  _nmode=0; 
  _rawbrfrsum=0.0;
  
}


EvtDecayBase* EvtParticleDecayList::getDecayModel(EvtParticle *p){

  if (p->getNDaug()!=0) {
    assert(p->getChannel()>=0);
    return getDecay(p->getChannel()).getDecayModel();
  }
  if (p->getChannel() >(-1) ) {
    return getDecay(p->getChannel()).getDecayModel();
  }


  if (getNMode()==0 ) {
    return 0;
  }
  if (getRawBrfrSum()<0.00000001 ) {
    return 0;
  }

  if (getNMode()==1) {
    p->setChannel(0);
    return getDecay(0).getDecayModel();
  }
  
  if (p->getChannel() >(-1)) {
    report(ERROR,"EvtGen") << "Internal error!!!"<<endl;
    ::abort();
  }

  int j;

  for (j=0;j<10000000;j++){

    double u=EvtRandom::Flat();

    int i;
    bool breakL=false;
    for (i=0;i<getNMode();i++) {

      if ( breakL ) continue;
      if (u<getDecay(i).getBrfrSum()) {
	breakL=true;
	//special case for decay of on particel to another
	// e.g. K0->K0S

	if (getDecay(i).getDecayModel()->getNDaug()==1 ) {
	  p->setChannel(i);
	  return getDecay(i).getDecayModel(); 
	} 

	if ( p->hasValidP4() ) {
	  if (getDecay(i).getMassMin() < p->mass() ) {
	    p->setChannel(i);
	    return getDecay(i).getDecayModel(); 
	  }
	}
	else{
	    //Lange apr29-2002 - dont know the mass yet
	    p->setChannel(i);
	    return getDecay(i).getDecayModel(); 
	}
      }
    }
  }

  //Ok, we tried 10000000 times above to pick a decay channel that is
  //kinematically allowed! Now we give up and search all channels!
  //if that fails, the particle will not be decayed!
  
  report(ERROR,"EvtGen") << "Tried 10000000 times to generate decay of "
                         << EvtPDL::name(p->getId())<< " with mass="<<p->mass()<<endl;
  report(ERROR,"EvtGen") << "Will take first kinematically allowed decay in the decay table" 
			 <<endl;

  int i;

  //Need to check that we don't use modes with 0 branching fractions.
  double previousBrSum=0.0;
  for (i=0;i<getNMode();i++) {
      if(getDecay(i).getBrfrSum()!=previousBrSum){
	  if ( getDecay(i).getMassMin() < p->mass() ) {
	      p->setChannel(i);
	      return getDecay(i).getDecayModel(); 
	  }
      }
      previousBrSum=getDecay(i).getBrfrSum();
  }

  report(ERROR,"EvtGen") << "Could not decay:"
			 <<EvtPDL::name(p->getId()).c_str()
			 <<" with mass:"<<p->mass()
			 <<" will throw event away! "<<endl;
  
  EvtStatus::setRejectFlag();
  return 0;

}


void EvtParticleDecayList::setNMode(int nmode){

  EvtParticleDecayPtr* _decaylist_new= new EvtParticleDecayPtr[nmode];

  if (_nmode!=0){
    report(ERROR,"EvtGen") << "Error _nmode not equal to zero!!!"<<endl;
    ::abort();
  }
  if (_decaylist!=0) {
    delete [] _decaylist;
  }
  _decaylist=_decaylist_new;
  _nmode=nmode;

}


EvtParticleDecay& EvtParticleDecayList::getDecay(int nchannel) {
  if (nchannel>=_nmode) {
    report(ERROR,"EvtGen") <<"Error getting channel:"
			   <<nchannel<<" with only "<<_nmode
			   <<" stored!"<<endl;
    ::abort();
  }
  return *(_decaylist[nchannel]);
}

void EvtParticleDecayList::makeChargeConj(EvtParticleDecayList* conjDecayList){

  _rawbrfrsum=conjDecayList->_rawbrfrsum;

  setNMode(conjDecayList->_nmode);
  
  int i;

  for(i=0;i<_nmode;i++){
    _decaylist[i]=new EvtParticleDecay;
    _decaylist[i]->chargeConj(conjDecayList->_decaylist[i]);
  }

}

void EvtParticleDecayList::addMode(EvtDecayBase* decay, double brfrsum,
				   double massmin){

  EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode+1];

  int i;
  for(i=0;i<_nmode;i++){
    newlist[i]=_decaylist[i];
  }

  _rawbrfrsum=brfrsum;

  newlist[_nmode]=new EvtParticleDecay;  

  newlist[_nmode]->setDecayModel(decay);
  newlist[_nmode]->setBrfrSum(brfrsum);
  newlist[_nmode]->setMassMin(massmin);

  EvtDecayBase *newDec=newlist[_nmode]->getDecayModel();
  for(i=0;i<_nmode;i++){
    if ( newDec->matchingDecay(*(newlist[i]->getDecayModel())) ) {

      //sometimes its ok..
      if ( newDec->getModelName() == "JETSET" || newDec->getModelName() == "PYTHIA" ) continue;
      if ( newDec->getModelName() == "JSCONT" || newDec->getModelName() == "PYCONT" ) continue;
      if ( newDec->getModelName() == "PYGAGA"  ) continue;
      if ( newDec->getModelName() == "LUNDAREALAW" ) continue;
      report(ERROR,"EvtGen") << "Two matching decays with same parent in decay table\n";
      report(ERROR,"EvtGen") << "Please fix that\n";
      report(ERROR,"EvtGen") << "Parent " << EvtPDL::name(newDec->getParentId()).c_str() << endl;
      for (int j=0; j<newDec->getNDaug(); j++)
	report(ERROR,"EvtGen") << "Daughter " << EvtPDL::name(newDec->getDaug(j)).c_str() << endl;
      assert(0);
    }
  }

  if (_nmode!=0){
    delete [] _decaylist;
  }

  if ( ( _nmode == 0 ) && ( _decaylist != 0 ) ) delete [] _decaylist ;

  _nmode++;

  _decaylist=newlist;

}


void EvtParticleDecayList::finalize(){

  if (_nmode>0) {
    if ( _rawbrfrsum< 0.000001 ) {
      report(ERROR,"EvtGen") << "Please give me a "
			     <<    "branching fraction sum greater than 0\n";
      assert(0);
    }
    if (fabs(_rawbrfrsum-1.0)>0.0001) {
      report(INFO,"EvtGen") <<"Warning, sum of branching fractions for "
      			    <<EvtPDL::name(_decaylist[0]->getDecayModel()->getParentId()).c_str()
      			    <<" is "<<_rawbrfrsum<<endl;
      report(INFO,"EvtGen") << "rescaled to one! "<<endl;
      
    }

    int i;

    for(i=0;i<_nmode;i++){
      double brfrsum=_decaylist[i]->getBrfrSum()/_rawbrfrsum;
      _decaylist[i]->setBrfrSum(brfrsum);      
    }

  }

}


EvtParticleDecayList& EvtParticleDecayList::operator=(const EvtParticleDecayList &o) {
   if (this != &o) {
      removeDecay();
      _nmode=o._nmode;
      _rawbrfrsum=o._rawbrfrsum;
      _decaylist=new EvtParticleDecayPtr[_nmode];
      
      int i;
      for(i=0;i<_nmode;i++){
       _decaylist[i]=new EvtParticleDecay;  
       
       EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
       
       EvtDecayBase *tModelNew=tModel->clone();
       if (tModel->getPHOTOS()){
          tModelNew->setPHOTOS();
       }
       if (tModel->verbose()){
          tModelNew->setVerbose();
       }
       if (tModel->summary()){
          tModelNew->setSummary();
       }
       std::vector<std::string> args;
       int j;
       for(j=0;j<tModel->getNArg();j++){
          args.push_back(tModel->getArgStr(j));
       }
       tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
                                tModel->getDaugs(),
                                tModel->getNArg(),
                                args,
                                tModel->getModelName(),
                                tModel->getBranchingFraction());
       _decaylist[i]->setDecayModel(tModelNew);
       
       //_decaylist[i]->setDecayModel(tModel);
       _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
       _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
      }
   }
   return *this;


}


void EvtParticleDecayList::removeMode(EvtDecayBase* decay) {
   // here we will delete a decay with the same final state particles
   // and recalculate the branching fractions for the remaining modes
   int match = -1;
   int i;
   double match_bf;

   for(i=0;i<_nmode;i++){
      if ( decay->matchingDecay(*(_decaylist[i]->getDecayModel())) ) {
       match = i;
      }
   }

   if (match < 0) {
      report(ERROR,"EvtGen") << " Attempt to remove undefined mode for" << endl
                           << "Parent " << EvtPDL::name(decay->getParentId()).c_str() << endl
                           << "Daughters: ";
      for (int j=0; j<decay->getNDaug(); j++)
      report(ERROR,"") << EvtPDL::name(decay->getDaug(j)).c_str() << " ";
      report(ERROR,"") << endl;
      ::abort();
   }

   if (match == 0) {
      match_bf = _decaylist[match]->getBrfrSum();
   } else {
      match_bf = (_decaylist[match]->getBrfrSum()
                -_decaylist[match-1]->getBrfrSum());
   }

   double divisor = 1-match_bf;
   if (divisor < 0.000001 && _nmode > 1) {
      report(ERROR,"EvtGen") << "Removing requested mode leaves "
                           <<  EvtPDL::name(decay->getParentId()).c_str()
                           << " with zero sum branching fraction," << endl
                           << "but more than one decay mode remains. Aborting."
                           << endl;
      ::abort();
   }

   EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode-1];

   for(i=0;i<match;i++){
      newlist[i]=_decaylist[i];
      newlist[i]->setBrfrSum(newlist[i]->getBrfrSum()/divisor);
   }
   for(i=match+1; i<_nmode; i++) {
      newlist[i-1]=_decaylist[i];
      newlist[i-1]->setBrfrSum((newlist[i-1]->getBrfrSum()-match_bf)/divisor);
   }


   delete [] _decaylist;

  _nmode--;

  _decaylist=newlist;

  if (_nmode == 0) {
     delete [] _decaylist;
  }

}


bool EvtParticleDecayList::isJetSet() {
  int i ;
  EvtDecayBase * decayer ;
 
  for ( i = 0 ;
        i < getNMode() ;
	i++ ) {
    decayer = getDecay( i ).getDecayModel ( ) ;
    if ( decayer -> getModelName() == "PYTHIA" ) return true ;
  }
  
  return false ;
}
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	}