#include <EvtParticleDecayList.hh>
Public Member Functions | |
void | addMode (EvtDecayBase *decay, double brfr, double massmin) |
void | alocateDecay (int nmode) |
EvtParticleDecayList (const EvtParticleDecayList &o) | |
EvtParticleDecayList () | |
void | finalize () |
EvtParticleDecay & | getDecay (int nchannel) |
EvtDecayBase * | getDecayModel (int imode) |
EvtDecayBase * | getDecayModel (EvtParticle *p) |
int | getNMode () |
double | getRawBrfrSum () |
void | makeChargeConj (EvtParticleDecayList *conjDecayList) |
void | printSummary () |
void | removeDecay () |
void | removeMode (EvtDecayBase *decay) |
void | setNMode (int nmode) |
void | setRawBrfrSum (double rawbrfrsum) |
~EvtParticleDecayList () | |
Private Attributes | |
EvtParticleDecayPtr * | _decaylist |
int | _nmode |
double | _rawbrfrsum |
|
00032 { 00033 _decaylist=0; 00034 _nmode=0; 00035 _rawbrfrsum=0; 00036 }
|
|
00036 { 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 //_decaylist[i]->setDecayModel(tModel); 00071 _decaylist[i]->setBrfrSum(o._decaylist[i]->getBrfrSum()); 00072 _decaylist[i]->setMassMin(o._decaylist[i]->getMassMin()); 00073 } 00074 00075 00076 }
|
|
00079 { 00080 00081 int i; 00082 for(i=0;i<_nmode;i++){ 00083 delete _decaylist[i]; 00084 } 00085 00086 if (_decaylist!=0) delete [] _decaylist; 00087 }
|
|
00246 { 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 //sometimes its ok.. 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 }
|
|
00059 { 00060 _decaylist= new EvtParticleDecayPtr[nmode]; 00061 }
|
|
00291 { 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 }
|
|
00220 { 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 }
|
|
00047 {return _decaylist[imode]->getDecayModel();}
|
|
00113 { 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 //special case for decay of on particel to another 00155 // e.g. K0->K0S 00156 00157 if (getDecay(i).getDecayModel()->getNDaug()==1 ) { 00158 p->setChannel(i); 00159 return getDecay(i).getDecayModel(); 00160 } 00161 00162 if ( p->hasValidP4() ) { 00163 //report(INFO,"EvtGen") << "amazing " << EvtPDL::name(p->getId()) << " " << getDecay(i).getMassMin() << " "<<p->mass() << " " << i << endl; 00164 if (getDecay(i).getMassMin() < p->mass() ) { 00165 p->setChannel(i); 00166 return getDecay(i).getDecayModel(); 00167 } 00168 } 00169 else{ 00170 //Lange apr29-2002 - dont know the mass yet 00171 p->setChannel(i); 00172 return getDecay(i).getDecayModel(); 00173 } 00174 } 00175 } 00176 } 00177 00178 //Ok, we tried 1000 times above to pick a decay channel that is 00179 //kinematically allowed! Now we give up and search all channels! 00180 //if that fails, the particle will not be decayed! 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 // report(ERROR,"EvtGen") << "Will terminate execution."<<endl; 00197 00198 // ::abort(); 00199 EvtStatus::setRejectFlag(); 00200 return 0; 00201 00202 }
|
|
00042 {return _nmode;}
|
|
00052 {return _rawbrfrsum;}
|
|
00230 { 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 }
|
|
00089 { 00090 00091 int i; 00092 for(i=0;i<_nmode;i++){ 00093 _decaylist[i]->printSummary(); 00094 } 00095 00096 }
|
|
00098 { 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 }
|
|
00317 { 00318 // here we will delete a decay with the same final state particles 00319 // and recalculate the branching fractions for the remaining modes 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 }
|
|
00205 { 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 }
|
|
00053 {_rawbrfrsum=rawbrfrsum;}
|
|
|
|
|
|
|