00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include <iostream>
00024 #include <fstream>
00025 #include <stdlib.h>
00026 #include <ctype.h>
00027 #include "EvtGenBase/EvtParticleDecayList.hh"
00028 #include "EvtGenBase/EvtParticle.hh"
00029 #include "EvtGenBase/EvtRandom.hh"
00030 #include "EvtGenBase/EvtReport.hh"
00031 #include "EvtGenBase/EvtPDL.hh"
00032 #include "EvtGenBase/EvtStatus.hh"
00033 using std::endl;
00034 using std::fstream;
00035
00036 EvtParticleDecayList::EvtParticleDecayList(const EvtParticleDecayList &o) {
00037 _nmode=o._nmode;
00038 _rawbrfrsum=o._rawbrfrsum;
00039 _decaylist=new EvtParticleDecayPtr[_nmode];
00040
00041 int i;
00042 for(i=0;i<_nmode;i++){
00043 _decaylist[i]=new EvtParticleDecay;
00044
00045 EvtDecayBase *tModel=o._decaylist[i]->getDecayModel();
00046
00047 EvtDecayBase *tModelNew=tModel->clone();
00048 if (tModel->getPHOTOS()){
00049 tModelNew->setPHOTOS();
00050 }
00051 if (tModel->verbose()){
00052 tModelNew->setVerbose();
00053 }
00054 if (tModel->summary()){
00055 tModelNew->setSummary();
00056 }
00057 std::vector<std::string> args;
00058 int j;
00059 for(j=0;j<tModel->getNArg();j++){
00060 args.push_back(tModel->getArgStr(j));
00061 }
00062 tModelNew->saveDecayInfo(tModel->getParentId(),tModel->getNDaug(),
00063 tModel->getDaugs(),
00064 tModel->getNArg(),
00065 args,
00066 tModel->getModelName(),
00067 tModel->getBranchingFraction());
00068 _decaylist[i]->setDecayModel(tModelNew);
00069
00070
00071 _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum());
00072 _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin());
00073 }
00074
00075
00076 }
00077
00078
00079 EvtParticleDecayList::~EvtParticleDecayList(){
00080
00081 int i;
00082 for(i=0;i<_nmode;i++){
00083 delete _decaylist[i];
00084 }
00085
00086 if (_decaylist!=0) delete [] _decaylist;
00087 }
00088
00089 void EvtParticleDecayList::printSummary(){
00090
00091 int i;
00092 for(i=0;i<_nmode;i++){
00093 _decaylist[i]->printSummary();
00094 }
00095
00096 }
00097
00098 void EvtParticleDecayList::removeDecay(){
00099
00100 int i;
00101 for(i=0;i<_nmode;i++){
00102 delete _decaylist[i];
00103 }
00104
00105 delete [] _decaylist;
00106 _decaylist=0;
00107 _nmode=0;
00108 _rawbrfrsum=0.0;
00109
00110 }
00111
00112
00113 EvtDecayBase* EvtParticleDecayList::getDecayModel(EvtParticle *p){
00114
00115 if (p->getNDaug()!=0) {
00116 assert(p->getChannel()>=0);
00117 return getDecay(p->getChannel()).getDecayModel();
00118 }
00119 if (p->getChannel() >(-1) ) {
00120 return getDecay(p->getChannel()).getDecayModel();
00121 }
00122
00123
00124 if (getNMode()==0 ) {
00125 return 0;
00126 }
00127 if (getRawBrfrSum()<0.00000001 ) {
00128 return 0;
00129 }
00130
00131 if (getNMode()==1) {
00132 p->setChannel(0);
00133 return getDecay(0).getDecayModel();
00134 }
00135
00136 if (p->getChannel() >(-1)) {
00137 report(ERROR,"EvtGen") << "Internal error!!!"<<endl;
00138 ::abort();
00139 }
00140
00141 int j;
00142
00143 for (j=0;j<1000;j++){
00144
00145 double u=EvtRandom::Flat();
00146
00147 int i;
00148 bool breakL=false;
00149 for (i=0;i<getNMode();i++) {
00150
00151 if ( breakL ) continue;
00152 if (u<getDecay(i).getBrfrSum()) {
00153 breakL=true;
00154
00155
00156
00157 if (getDecay(i).getDecayModel()->getNDaug()==1 ) {
00158 p->setChannel(i);
00159 return getDecay(i).getDecayModel();
00160 }
00161
00162 if ( p->hasValidP4() ) {
00163
00164 if (getDecay(i).getMassMin() < p->mass() ) {
00165 p->setChannel(i);
00166 return getDecay(i).getDecayModel();
00167 }
00168 }
00169 else{
00170
00171 p->setChannel(i);
00172 return getDecay(i).getDecayModel();
00173 }
00174 }
00175 }
00176 }
00177
00178
00179
00180
00181
00182 int i;
00183
00184 for (i=0;i<getNMode();i++) {
00185
00186 if ( getDecay(i).getMassMin() < p->mass() ) {
00187 return getDecay(i).getDecayModel();
00188 }
00189 }
00190
00191 report(ERROR,"EvtGen") << "Could not decay:"
00192 <<EvtPDL::name(p->getId()).c_str()
00193 <<" with mass:"<<p->mass()
00194 <<" will throw event away! "<<endl;
00195
00196
00197
00198
00199 EvtStatus::setRejectFlag();
00200 return 0;
00201
00202 }
00203
00204
00205 void EvtParticleDecayList::setNMode(int nmode){
00206
00207 EvtParticleDecayPtr* _decaylist_new= new EvtParticleDecayPtr[nmode];
00208
00209 if (_nmode!=0){
00210 report(ERROR,"EvtGen") << "Error _nmode not equal to zero!!!"<<endl;
00211 ::abort();
00212 delete [] _decaylist;
00213 }
00214 _decaylist=_decaylist_new;
00215 _nmode=nmode;
00216
00217 }
00218
00219
00220 EvtParticleDecay& EvtParticleDecayList::getDecay(int nchannel) {
00221 if (nchannel>=_nmode) {
00222 report(ERROR,"EvtGen") <<"Error getting channel:"
00223 <<nchannel<<" with only "<<_nmode
00224 <<" stored!"<<endl;
00225 ::abort();
00226 }
00227 return *(_decaylist[nchannel]);
00228 }
00229
00230 void EvtParticleDecayList::makeChargeConj(EvtParticleDecayList* conjDecayList){
00231
00232 _rawbrfrsum=conjDecayList->_rawbrfrsum;
00233
00234 setNMode(conjDecayList->_nmode);
00235
00236 int i;
00237
00238 for(i=0;i<_nmode;i++){
00239 _decaylist[i]=new EvtParticleDecay;
00240 _decaylist[i]->chargeConj(conjDecayList->_decaylist[i]);
00241 }
00242
00243 }
00244
00245 void EvtParticleDecayList::addMode(EvtDecayBase* decay, double brfrsum,
00246 double massmin){
00247
00248 EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode+1];
00249
00250 int i;
00251 for(i=0;i<_nmode;i++){
00252 newlist[i]=_decaylist[i];
00253 }
00254
00255 _rawbrfrsum=brfrsum;
00256
00257 newlist[_nmode]=new EvtParticleDecay;
00258
00259 newlist[_nmode]->setDecayModel(decay);
00260 newlist[_nmode]->setBrfrSum(brfrsum);
00261 newlist[_nmode]->setMassMin(massmin);
00262
00263 EvtDecayBase *newDec=newlist[_nmode]->getDecayModel();
00264 for(i=0;i<_nmode;i++){
00265 if ( newDec->matchingDecay(*(newlist[i]->getDecayModel())) ) {
00266
00267
00268 if ( newDec->getModelName() == "JETSET" || newDec->getModelName() == "PYTHIA" ) continue;
00269 if ( newDec->getModelName() == "JSCONT" || newDec->getModelName() == "PYCONT" ) continue;
00270 if ( newDec->getModelName() == "PYGAGA" ) continue;
00271 report(ERROR,"EvtGen") << "Two matching decays with same parent in decay table\n";
00272 report(ERROR,"EvtGen") << "Please fix that\n";
00273 report(ERROR,"EvtGen") << "Parent " << EvtPDL::name(newDec->getParentId()).c_str() << endl;
00274 for (int j=0; j<newDec->getNDaug(); j++)
00275 report(ERROR,"EvtGen") << "Daughter " << EvtPDL::name(newDec->getDaug(j)).c_str() << endl;
00276 assert(0);
00277 }
00278 }
00279
00280 if (_nmode!=0){
00281 delete [] _decaylist;
00282 }
00283
00284 _nmode++;
00285
00286 _decaylist=newlist;
00287
00288 }
00289
00290
00291 void EvtParticleDecayList::finalize(){
00292
00293 if (_nmode>0) {
00294 if ( _rawbrfrsum< 0.000001 ) {
00295 report(ERROR,"EvtGen") << "Please give me a "
00296 << "branching fraction sum greater than 0\n";
00297 assert(0);
00298 }
00299 if (fabs(_rawbrfrsum-1.0)>0.0001) {
00300 report(INFO,"EvtGen") <<"Warning, sum of branching fractions for "
00301 <<EvtPDL::name(_decaylist[0]->getDecayModel()->getParentId()).c_str()
00302 <<" is "<<_rawbrfrsum<<endl;
00303 report(INFO,"EvtGen") << "rescaled to one! "<<endl;
00304
00305 }
00306
00307 int i;
00308
00309 for(i=0;i<_nmode;i++){
00310 double brfrsum=_decaylist[i]->getBrfrSum()/_rawbrfrsum;
00311 _decaylist[i]->setBrfrSum(brfrsum);
00312 }
00313
00314 }
00315
00316 }
00317 void EvtParticleDecayList::removeMode(EvtDecayBase* decay) {
00318
00319
00320 int match = -1;
00321 int i;
00322 double match_bf;
00323
00324 for(i=0;i<_nmode;i++){
00325 if ( decay->matchingDecay(*(_decaylist[i]->getDecayModel())) ) {
00326 match = i;
00327 }
00328 }
00329
00330 if (match < 0) {
00331 report(ERROR,"EvtGen") << " Attempt to remove undefined mode for" << endl
00332 << "Parent " << EvtPDL::name(decay->getParentId()).c_str() << endl
00333 << "Daughters: ";
00334 for (int j=0; j<decay->getNDaug(); j++)
00335 report(ERROR,"") << EvtPDL::name(decay->getDaug(j)).c_str() << " ";
00336 report(ERROR,"") << endl;
00337 ::abort();
00338 }
00339
00340 if (match == 0) {
00341 match_bf = _decaylist[match]->getBrfrSum();
00342 } else {
00343 match_bf = (_decaylist[match]->getBrfrSum()
00344 -_decaylist[match-1]->getBrfrSum());
00345 }
00346
00347 double divisor = 1-match_bf;
00348 if (divisor < 0.000001 && _nmode > 1) {
00349 report(ERROR,"EvtGen") << "Removing requested mode leaves "
00350 << EvtPDL::name(decay->getParentId()).c_str()
00351 << " with zero sum branching fraction," << endl
00352 << "but more than one decay mode remains. Aborting."
00353 << endl;
00354 ::abort();
00355 }
00356
00357 EvtParticleDecayPtr* newlist=new EvtParticleDecayPtr[_nmode-1];
00358
00359 for(i=0;i<match;i++){
00360 newlist[i]=_decaylist[i];
00361 newlist[i]->setBrfrSum(newlist[i]->getBrfrSum()/divisor);
00362 }
00363 for(i=match+1; i<_nmode; i++) {
00364 newlist[i-1]=_decaylist[i];
00365 newlist[i-1]->setBrfrSum((newlist[i-1]->getBrfrSum()-match_bf)/divisor);
00366 }
00367
00368
00369 delete [] _decaylist;
00370
00371 _nmode--;
00372
00373 _decaylist=newlist;
00374
00375 if (_nmode == 0) {
00376 delete [] _decaylist;
00377 }
00378
00379 }