00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #include "MdcRecoUtil/AstStringMap.h"
00029 #include "MdcRecoUtil/Pdt.h"
00030 #include "MdcRecoUtil/PdtEntry.h"
00031
00032 #include <string.h>
00033 #include <stdio.h>
00034 #include <assert.h>
00035
00036 #include <string>
00037 using std::string;
00038 #include <vector>
00039 using std::vector;
00040 #include <map>
00041 using std::map;
00042 #include <iostream>
00043 #include <strstream>
00044 using std::cerr;
00045 using std::cout;
00046 using std::endl;
00047 using std::istrstream;
00048 using std::ostream;
00049 using std::strstream;
00050
00051 #define HBARC (197.327*1.e-3*1.e-13) // GeV*cm
00052
00053 PdtEntry* Pdt::_positiveEntries[5]={0,0,0,0,0};
00054 PdtEntry* Pdt::_negativeEntries[5]={0,0,0,0,0};
00055 PdtEntry* Pdt::_neutralEntries[5]={0,0,0,0,0};
00056
00057
00058 map<int, PdtEntry*>*
00059 Pdt::_entries = new map<int, PdtEntry*>;
00060
00061
00062 AstStringMap< PdtEntry >*
00063 Pdt::_entriesn = new AstStringMap< PdtEntry >;
00064
00065
00066 map<int, PdtEntry*>*
00067 Pdt::_entriesg = new map<int, PdtEntry*>;
00068
00069
00070 map<int, PdtEntry*>*
00071 Pdt::_entriesl = new map<int, PdtEntry*>;
00072
00073 PdtLund::LundType Pdt::antiCode(PdtLund::LundType kf_id){
00074
00075 int kf = (int) kf_id;
00076
00077 int i2 = (kf/10) % 10;
00078 int i3 = (kf/100) % 10;
00079 int i4 = (kf/1000) % 10;
00080
00081 if (i2 == i3 && i4 == 0 && kf > 100)
00082 return lundId(kf);
00083 else if (kf >= 21 && kf <= 23) return lundId(kf);
00084
00085 else if (kf == 310) return lundId(130);
00086 else if (kf == 130) return lundId(310);
00087
00088 else return lundId(-kf);
00089 }
00090
00091 void Pdt::addParticle(const char *pname, PdtLund::LundType id,
00092 float spin, float charge, float mass,
00093 float width, float cut){
00094
00095 PdtEntry* nentry = new PdtEntry(pname, id, spin, charge, mass, width, cut);
00096 (*_entries)[pdgId(id)]=nentry;
00097
00098 if ( (*_entriesn)[pname] != 0 ) {
00099
00100 }
00101
00102 _entriesn->insert( pname )=nentry;
00103 if(PdtGeant::null!=nentry->geantId()){
00104 (*_entriesg)[nentry->geantId()]=nentry;
00105 }
00106
00107 (*_entriesl)[id]=nentry;
00108
00109 if(PdtPid::null!=nentry->pidId()) {
00110 if(nentry->charge()>0) {
00111 _positiveEntries[nentry->pidId()]=nentry;
00112 }else{
00113 _negativeEntries[nentry->pidId()]=nentry;
00114 }
00115 }
00116
00117 if(PdtPid::none!=nentry->pidNeutId()){
00118 if(nentry->charge()==0){
00119 _neutralEntries[nentry->pidNeutId()]=nentry;
00120 }
00121 }
00122
00123 }
00124
00125 void Pdt::addParticle(const char *pname, PdtGeant::GeantType id,
00126 float spin, float charge, float mass,
00127 float width, float cut){
00128
00129 PdtEntry* nentry = new PdtEntry(pname, id, spin, charge, mass, width, cut);
00130
00131 (*_entries)[pdgId(id)]=nentry;
00132 if ( (*_entriesn)[pname] != 0 ) {
00133
00134 }
00135
00136 _entriesn->insert( pname )=nentry;
00137 if(PdtGeant::null!=nentry->geantId()){
00138 (*_entriesg)[id]=nentry;
00139 }
00140
00141 (*_entriesl)[lundId(id)]=nentry;
00142
00143 if(PdtPid::null!=nentry->pidId()) {
00144 if(nentry->charge()>0) {
00145 _positiveEntries[nentry->pidId()]=nentry;
00146 }else{
00147 _negativeEntries[nentry->pidId()]=nentry;
00148 }
00149
00150 }
00151
00152 if(PdtPid::none!=nentry->pidNeutId()){
00153 if(nentry->charge()==0) {
00154 _neutralEntries[nentry->pidNeutId()]=nentry;
00155 }
00156 }
00157
00158 }
00159
00160 void Pdt::addDecay(const char* pname, float bf,
00161 const char* child1, const char* child2, const char* child3,
00162 const char* child4, const char* child5){
00163
00164 PdtEntry *primary = lookup(pname);
00165 if (primary == 0) return;
00166 const char *children[5] = {child1, child2, child3, child4, child5};
00167
00168 vector<PdtEntry*> *kids = new vector<PdtEntry*>;
00169
00170 int nkids;
00171 for(nkids=0; nkids<5 && children[nkids]!=0; nkids++ )
00172 {
00173 PdtEntry* secondary = lookup(children[nkids]);
00174 if( secondary ==0 ) break;
00175 kids->push_back(secondary);
00176 }
00177
00178 primary->addDecay(bf, kids );
00179
00180 }
00181
00182 void Pdt::addDecay(PdtLund::LundType id, float bf,
00183 PdtLund::LundType child1,
00184 PdtLund::LundType child2,
00185 PdtLund::LundType child3,
00186 PdtLund::LundType child4,
00187 PdtLund::LundType child5){
00188
00189 PdtEntry *primary = lookup(id);
00190 if (primary == 0) return;
00191
00192 PdtLund::LundType children[5] = {child1, child2, child3, child4, child5};
00193 vector<PdtEntry*> *kids = new vector<PdtEntry*>;
00194
00195 int nkids;
00196 for(nkids=0; nkids<5 && children[nkids]!=0; nkids++ )
00197 {
00198 PdtEntry* secondary = lookup(children[nkids]);
00199 if( secondary ==0 ) break;
00200 kids->push_back(secondary);
00201 }
00202
00203 primary->addDecay(bf, kids );
00204
00205 }
00206
00207 PdtEntry* Pdt::lookup( const std::string& name ){
00208 return (*_entriesn)[name];
00209 }
00210
00211 PdtEntry* Pdt::lookup(PdtLund::LundType id){
00212
00213 return (*_entriesl)[id];
00214 }
00215
00216 PdtEntry* Pdt::lookup(PdtGeant::GeantType id){
00217
00218 return (*_entriesg)[id];
00219 }
00220
00221 PdtEntry* Pdt::lookup(PdtPdg::PdgType id){
00222
00223 return (*_entries)[id];
00224 }
00225
00226 PdtEntry* Pdt::lookup(PdtPid::PidType id, int charge){
00227
00228
00229 if(id==PdtPid::null) return lookup(PdtLund::null);
00230 if(charge>0) return _positiveEntries[id];
00231 return _negativeEntries[id];
00232
00233 }
00234
00235 PdtEntry* Pdt::lookup(PdtPid::PidNeutralType id, int){
00236
00237
00238
00239 if(id==PdtPid::none) return lookup(PdtLund::null);
00240 return _neutralEntries[id];
00241
00242 }
00243
00244 void Pdt::printOn(ostream& cout){
00245
00246
00247 cout << "Number of particles in tables: "<<endl
00248 << "By pdt number :" << _entries->size() <<endl
00249 << "By name :" << _entriesn->size()<<endl
00250 << "By geant number :" << _entriesg->size()<<endl
00251 << "By lund number :" <<_entriesl->size()<<endl;
00252
00253 cout << '\n';
00254 int i;
00255 for (i=0;i<72;i++) cout << '-';
00256 cout<<'\n';
00257
00258 cout << " Particle Data Table\n\n";
00259 cout << " id name mass width spin charge lifetime\n";
00260 cout << " (GeV) (GeV) (cm) \n";
00261 for (i=0; i<72; i++) cout << '-';
00262 cout << '\n';
00263
00264 PdtEntry *p;
00265
00266 map<int, PdtEntry*>::iterator iter=_entries->begin();
00267
00268 while ( iter != _entries->end()) {
00269 p = iter->second;
00270 assert(0!=p);
00271 p->printOn(cout);
00272 p->printBFOn(cout);
00273 ++iter;
00274 }
00275 cout << '\n';
00276 for (i=0; i<72; i++) cout << '-';
00277 cout << '\n';
00278
00279 map<int, PdtEntry*>::iterator iter1=_entries->begin();
00280
00281 while ( iter1 != _entries->end()) {
00282 p=iter1->second;
00283 cout << "The conjuate of "<< p->name() << " is ";
00284 cout << conjugate(p)->name() <<endl;
00285 iter1++;
00286 }
00287 }
00288
00289
00290
00291 void Pdt::readMCppTable(string filenm){
00292
00293 string path = filenm;
00294 char line[512];
00295 FILE *ifl = fopen(path.c_str(),"r");
00296
00297 if(0 == ifl) cout<<" Pdt::readMCppTable: please copy " <<
00298 filenm << " in to your run directory!"<<
00299 std::endl;
00300 assert(0!=ifl);
00301
00302 while ( fgets(line,512,ifl) )
00303 {
00304 if (strlen(line) >= 511)
00305 {
00306 cerr << "Pdt.read : input line is too long\n";
00307 assert(0);
00308 }
00309 istrstream linestr(line);
00310 string opcode;
00311 char subcode;
00312 linestr >> opcode >> subcode;
00313
00314 if( opcode == "end" )
00315 break;
00316
00317 else if( opcode == "add" )
00318 {
00319 switch (subcode)
00320 {
00321 case 'p':
00322 {
00323 string classname;
00324 linestr >> classname;
00325
00326 if (classname == "Collision" )
00327 continue;
00328
00329 string name;
00330 int type;
00331 float mass, width, cut, charge, spin, lifetime;
00332
00333 linestr >> name >> type;
00334 linestr >> mass >> width >> cut >> charge;
00335 linestr >> spin >> lifetime;
00336
00337 charge /= 3.0;
00338 if (classname != "Meson")
00339 spin /= 2.0;
00340
00341
00342 if (lifetime > 0.0 && width < 1e-10)
00343 width = HBARC / (lifetime/10.0);
00344
00345 addParticle(name.c_str(), lundId(type), spin, charge, mass, width, cut);
00346 break;
00347 }
00348
00349 case 'c':
00350 {
00351 int ptype, nchild;
00352 float bf;
00353 string decayer;
00354
00355 linestr >> ptype >> bf >> decayer >> nchild;
00356 PdtEntry *parent = lookup(lundId(ptype));
00357 if (parent == 0) continue;
00358
00359 vector<PdtEntry*> *kids = new vector<PdtEntry*>;
00360
00361 int i;
00362 for(i=0; i<nchild; i++ )
00363 {
00364 int ctype;
00365 linestr >> ctype;
00366 PdtEntry* secondary = lookup(lundId(ctype));
00367 if( secondary ==0 ) break;
00368 kids->push_back(secondary);
00369 }
00370
00371 parent->addDecay(bf, kids );
00372 break;
00373 }
00374
00375 case 'd':
00376 break;
00377
00378 default:
00379 cerr << "Pdt.readMCppTable : unknown subcode '" << subcode
00380 << "' for operation add\n";
00381 break;
00382 }
00383 }
00384 }
00385
00386
00387
00388
00389 fclose(ifl);
00390
00391 }
00392
00393 void Pdt::deleteAll(){
00394
00395 PdtEntry *p;
00396
00397 map<int, PdtEntry*>::iterator iter=_entries->begin();
00398
00399 while ( iter != _entries->end()) {
00400 p = iter->second;
00401 assert(0!=p);
00402 delete p;
00403 ++iter;
00404 }
00405
00406 _entries->clear();
00407 _entriesn->clear();
00408 _entriesg->clear();
00409 _entriesl->clear();
00410 }
00411
00412
00413 PdtLund::LundType Pdt::lundId(const PdtGeant::GeantType gId){
00414
00415 if (gId == PdtGeant::deuteron ||
00416 gId == PdtGeant::tritium ||
00417 gId == PdtGeant::alpha ||
00418 gId == PdtGeant::geantino ||
00419 gId == PdtGeant::He3 ||
00420 gId == PdtGeant::Cerenkov)
00421 { return PdtLund::null; }
00422 else
00423 { return PdtGeant::_lundId[gId - PdtGeant::_firstGeantId]; }
00424
00425 }
00426
00427 PdtGeant::GeantType Pdt::geantId(const PdtLund::LundType lId){
00428
00429
00430
00431
00432
00433 if (lId == PdtLund::nu_e ||
00434 lId == PdtLund::anti_nu_e ||
00435 lId == PdtLund::nu_mu ||
00436 lId == PdtLund::anti_nu_mu ||
00437 lId == PdtLund::nu_tau ||
00438 lId == PdtLund::anti_nu_tau )
00439 { return PdtGeant::nu_e; }
00440 else
00441 {
00442 int i;
00443 for(i = 0; i< PdtGeant::_nGeantId; i++)
00444 { if (PdtGeant::_lundId[i] == lId) return geantId(i); }
00445 return PdtGeant::null;
00446 }
00447 }
00448
00449 PdtPid::PidType Pdt::pidId(const PdtLund::LundType lId){
00450
00451 switch (lId)
00452 {
00453 case PdtLund::e_minus:
00454 case PdtLund::e_plus:
00455 return PdtPid::electron;
00456 case PdtLund::mu_minus:
00457 case PdtLund::mu_plus:
00458 return PdtPid::muon;
00459 case PdtLund::pi_minus:
00460 case PdtLund::pi_plus:
00461 return PdtPid::pion;
00462 case PdtLund::K_minus:
00463 case PdtLund::K_plus:
00464 return PdtPid::kaon;
00465 case PdtLund::p_plus:
00466 case PdtLund::anti_p_minus:
00467 return PdtPid::proton;
00468 default:
00469 return PdtPid::null;
00470 }
00471 }
00472
00473 PdtPid::PidNeutralType Pdt::pidNeutId(const PdtLund::LundType lId){
00474
00475 switch (lId)
00476 {
00477 case PdtLund::gamma:
00478 return PdtPid::gamma;
00479 case PdtLund::pi0:
00480 return PdtPid::pi0;
00481 case PdtLund::K_L0:
00482 return PdtPid::K0L;
00483 case PdtLund::anti_n0:
00484 return PdtPid::anti_neutron;
00485 case PdtLund::n0:
00486 return PdtPid::neutron;
00487 default:
00488 return PdtPid::none;
00489 }
00490 }
00491
00492 PdtPid::PidType Pdt::pidId(const PdtGeant::GeantType gId){
00493 return pidId(lundId(gId));
00494 }
00495
00496 PdtPid::PidNeutralType Pdt::pidNeutId(const PdtGeant::GeantType gId){
00497 return pidNeutId(lundId(gId));
00498 }
00499
00500 PdtGeant::GeantType Pdt::geantId(const PdtPid::PidType pId, int charge){
00501 return geantId(lundId(pId, charge));
00502 }
00503
00504 PdtGeant::GeantType Pdt::geantId(const PdtPid::PidNeutralType pId, int charge){
00505 return geantId(lundId(pId, charge));
00506 }
00507
00508 PdtLund::LundType Pdt::lundId (const PdtPid::PidType pId, int charge){
00509
00510 if(charge == -1)
00511 {
00512 switch(pId)
00513 {
00514 case PdtPid::electron:
00515 return PdtLund::e_minus;
00516 case PdtPid::muon:
00517 return PdtLund::pi_minus;
00518 case PdtPid::pion:
00519 return PdtLund::pi_minus;
00520 case PdtPid::kaon:
00521 return PdtLund::K_minus;
00522 case PdtPid::proton:
00523 return PdtLund::anti_p_minus;
00524 default:
00525 return PdtLund::null;
00526 }
00527 }
00528 else if (charge == 1)
00529 {
00530 switch(pId)
00531 {
00532 case PdtPid::electron:
00533 return PdtLund::e_plus;
00534 case PdtPid::muon:
00535 return PdtLund::pi_plus;
00536 case PdtPid::pion:
00537 return PdtLund::pi_plus;
00538 case PdtPid::kaon:
00539 return PdtLund::K_plus;
00540 case PdtPid::proton:
00541 return PdtLund::p_plus;
00542 default:
00543 return PdtLund::null;
00544 }
00545 }
00546 else
00547 {
00548 cerr << "error: wrong charge; must be +1 or -1" << endl;
00549 return PdtLund::null;
00550 }
00551 }
00552
00553 PdtLund::LundType Pdt::lundId (const PdtPid::PidNeutralType pId, int charge){
00554
00555 if (charge == 0)
00556 {
00557 switch(pId)
00558 {
00559 case PdtPid::gamma:
00560 return PdtLund::gamma;
00561 case PdtPid::pi0:
00562 return PdtLund::pi0;
00563 case PdtPid::K0L:
00564 return PdtLund::K_L0;
00565 case PdtPid::neutron:
00566 return PdtLund::n0;
00567 case PdtPid::anti_neutron:
00568 return PdtLund::anti_n0;
00569 default:
00570 return PdtLund::null;
00571 }
00572 }
00573 else
00574 {
00575 cerr << "error: wrong charge; must be +1 or -1" << endl;
00576 return PdtLund::null;
00577 }
00578 }
00579
00580 PdtLund::LundType Pdt::lundId( const PdtPdg::PdgType pdg ) {
00581
00582 int ret;
00583 switch ( (int) pdg ) {
00584 case 551:
00585 ret = 10551 ;
00586 case 10443:
00587 ret = 20443 ;
00588 case 20443:
00589 ret = 30443 ;
00590 case 30443:
00591 ret = 40443 ;
00592 case 20553:
00593 ret = 30553 ;
00594 case 4322:
00595 ret = 4232 ;
00596 case -4322:
00597 ret = -4232 ;
00598 case 4232:
00599 ret = 4322 ;
00600 case -4232:
00601 ret = -4322 ;
00602 default:
00603 ret = (int) pdg;
00604 }
00605 return lundId(ret);
00606 }
00607
00608 PdtPdg::PdgType Pdt::pdgId( const PdtLund::LundType lund ) {
00609
00610 int ret;
00611 switch ( (int) lund ) {
00612 case 551:
00613 ret = 0 ;
00614 case 10551:
00615 ret = 551 ;
00616 case 10443:
00617 ret = 0 ;
00618 case 20443:
00619 ret = 10443 ;
00620 case 30443:
00621 ret = 20443 ;
00622 case 40443:
00623 ret = 30443 ;
00624 case 30553:
00625 ret = 20553 ;
00626 case 4232:
00627 ret = 4322 ;
00628 case -4232:
00629 ret = -4322 ;
00630 case 4322:
00631 ret = 4232 ;
00632 case -4322:
00633 ret = -4232 ;
00634 default:
00635 ret = (int) lund ;
00636 }
00637 return pdgId(ret);
00638 }
00639
00640
00641
00642 unsigned Pdt::PdtNumHash(const int& num){
00643
00644
00645
00646 if (num>0) return num%PdtHashSize;
00647 return (1-num)%PdtHashSize;
00648 }
00649
00650 unsigned Pdt::PdtNameHash(const string& c){
00651
00652 int i;
00653 int hash=0;
00654
00655 for(i=0;i<(int)c.length();i++){
00656 hash=(153*hash+(int)c[i])%PdtHashSize;
00657 }
00658
00659 return hash;
00660
00661 }
00662
00663 bool Pdt::sameOrConj( PdtLund::LundType id1, PdtLund::LundType id2) {
00664 return Pdt::sameOrConj(Pdt::lookup(id1),Pdt::lookup(id2));
00665 }
00666
00667 bool Pdt::sameOrConj( PdtPdg::PdgType id1, PdtPdg::PdgType id2) {
00668 return Pdt::sameOrConj(Pdt::lookup(id1),Pdt::lookup(id2));
00669 }
00670
00671 bool Pdt::sameOrConj( PdtEntry* pdt1, PdtEntry* pdt2) {
00672 if ( (pdt1==0) || (pdt2==0) ) return false;
00673 if ( pdt1->pdgId() == pdt2->pdgId() ) return true;
00674 if ( pdt1->conjugate()->pdgId() == pdt2->pdgId() ) return true;
00675 return false;
00676 }
00677
00678
00679
00680 const PdtEntry*
00681 Pdt:: conjugate( const PdtEntry* pdt )
00682 {
00683 if( pdt==0 ) return 0;
00684 PdtEntry* conjPdt = (PdtEntry*) pdt;
00685
00686
00687 int key = (int) pdt->pdgId();
00688
00689
00690 key *= -1;
00691
00692
00693
00694 map<int, PdtEntry*>::iterator iter=_entries->find(key);
00695
00696 if ( iter!=_entries->end() ) {
00697
00698 conjPdt=iter->second;
00699 }
00700
00701
00702
00703 return conjPdt;
00704 }
00705
00706
00707
00708