TEvtGen/EvtGenBase/EvtPDL.cxx

   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: EvtPDL.cc
  12 //
  13 // Description: routines to store particle properties in EvtPDL structure.
  14 //
  15 // Modification history:
  16 //
  17 //    DJL/RYD     September 25, 1996         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <iostream>
  23 #include <fstream>
  24 #include <stdlib.h>
  25 #include <ctype.h>
  26 #include <string.h>
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenBase/EvtPartProp.hh"
  29 #include "EvtGenBase/EvtId.hh"
  30 #include "EvtGenBase/EvtParticle.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32 using std::endl;
  33 using std::fstream;
  34 using std::ifstream;
  35
  36 static int first=1;
  37
  38 unsigned int EvtPDL::_firstAlias;
  39 int EvtPDL::_nentries;
  40
  41 std::map<std::string,int> EvtPDL::_particleNameLookup;
  42
  43 EvtPDL::EvtPDL() {
  44
  45   if (first!=0) {
  46     first=0;
  47     _nentries=0;
  48     _firstAlias=999999;
  49   }
  50
  51 }
  52
  53
  54 EvtPDL::~EvtPDL(){
  55
  56 }
  57
  58 void EvtPDL::read(const char* fname)
  59 {
  60   readPDT(fname);
  61 }
  62
  63 void EvtPDL::readPDT(const std::string fname){
  64
  65
  66   ifstream indec;
  67
  68   indec.open(fname.c_str());
  69
  70   char cmnd[100];
  71   char xxxx[100];
  72
  73   char pname[100];
  74   int  stdhepid;
  75   double mass;
  76   double pwidth;
  77   double pmaxwidth;
  78   int    chg3;
  79   int    spin2;
  80   double ctau;
  81   int    lundkc;
  82   EvtId i;
  83
  84   if (!indec) {
  85     report(ERROR,"EvtGen") << "Could not open:"<<fname.c_str()<<"EvtPDL"<<endl;
  86     return;
  87   }
  88
  89   do{
  90
  91     char ch,ch1;
  92
  93     do{
  94
  95       indec.get(ch);
  96       if (ch=='\n') indec.get(ch);
  97       if (ch!='*') {
  98         indec.putback(ch);
  99       }
 100       else{
 101         while (indec.get(ch1),ch1!='\n');
 102       }
 103     } while(ch=='*');
 104
 105     indec >> cmnd;
 106
 107     if (strcmp(cmnd,"end")){
 108
 109       if (!strcmp(cmnd,"add")){
 110
 111         indec >> xxxx;
 112         indec >> xxxx;
 113         indec >> pname;
 114         indec >> stdhepid;
 115         indec >> mass;
 116         indec >> pwidth;
 117         indec >> pmaxwidth;
 118         indec >> chg3;
 119         indec >> spin2;
 120         indec >> ctau;
 121         indec >> lundkc;
 122
 123
 124         i=EvtId(_nentries,_nentries);
 125
 126         EvtPartProp tmp;
 127
 128         tmp.setSpinType(EvtSpinType::SCALAR);
 129
 130
 131         if (spin2==0) tmp.setSpinType(EvtSpinType::SCALAR);
 132         if (spin2==1) tmp.setSpinType(EvtSpinType::DIRAC);
 133         if (spin2==2) tmp.setSpinType(EvtSpinType::VECTOR);
 134         if (spin2==3) tmp.setSpinType(EvtSpinType::RARITASCHWINGER);
 135         if (spin2==4) tmp.setSpinType(EvtSpinType::TENSOR);
 136         if (spin2==5) tmp.setSpinType(EvtSpinType::SPIN5HALF);
 137         if (spin2==6) tmp.setSpinType(EvtSpinType::SPIN3);
 138         if (spin2==7) tmp.setSpinType(EvtSpinType::SPIN7HALF);
 139         if (spin2==8) tmp.setSpinType(EvtSpinType::SPIN4);
 140         if (spin2==2 && mass < 0.0001 ) tmp.setSpinType(EvtSpinType::PHOTON);
 141         if (spin2==1 && mass < 0.0001 ) tmp.setSpinType(EvtSpinType::NEUTRINO);
 142
 143
 144         if (!strcmp(pname,"string")){
 145           tmp.setSpinType(EvtSpinType::STRING);
 146         }
 147
 148         if (!strcmp(pname,"vpho")){
 149           tmp.setSpinType(EvtSpinType::VECTOR);
 150         }
 151
 152
 153         tmp.setId(i);
 154         tmp.setIdChgConj(EvtId(-1,-1));
 155         tmp.setStdHep(stdhepid);
 156         tmp.setLundKC(lundkc);
 157         tmp.setName(pname);
 158         if (_particleNameLookup.find(std::string(pname))!=
 159             _particleNameLookup.end()) {
 160             report(ERROR,"EvtGen")<<"The particle name:"<<pname<<" is already defined."<<endl;
 161             report(ERROR,"EvtGen") << "Will terminate execution.";
 162             ::abort();
 163         }
 164         _particleNameLookup[std::string(pname)]=_nentries;
 165         tmp.setctau(ctau);
 166         tmp.setChg3(chg3);
 167
 168         tmp.initLineShape(mass,pwidth,pmaxwidth);
 169
 170
 171         partlist().push_back(tmp);
 172         _nentries++;
 173
 174       }
 175
 176       // if find a set read information and discard it
 177
 178       if (!strcmp(cmnd,"set")){
 179
 180         indec >> xxxx;
 181         indec >> xxxx;
 182         indec >> xxxx;
 183         indec >> xxxx;
 184       }
 185
 186     }
 187
 188   }while(strcmp(cmnd,"end"));
 189
 190   setUpConstsPdt();
 191
 192 }
 193
 194
 195 void EvtPDL::aliasChgConj(EvtId a,EvtId abar){
 196
 197   if (EvtPDL::chargeConj(EvtId(a.getId(),a.getId()))!=
 198                  EvtId(abar.getId(),abar.getId())) {
 199
 200     report(ERROR,"EvtGen")<<"Can't charge conjugate the two aliases:"
 201                           <<EvtPDL::name(a).c_str()<<" and "<<EvtPDL::name(abar).c_str()<<endl;
 202
 203     ::abort();
 204
 205   }
 206
 207   partlist()[a.getAlias()].setIdChgConj(abar);
 208   partlist()[abar.getAlias()].setIdChgConj(a);
 209
 210 }
 211
 212 EvtId EvtPDL::chargeConj(EvtId id){
 213   EvtId idchg=partlist()[id.getAlias()].getIdChgConj();
 214
 215   if (idchg!=EvtId(-1,-1)) return idchg;
 216
 217   if (id.getId()!=id.getAlias()){
 218     if (chargeConj(EvtId(id.getId(),id.getId()))==EvtId(id.getId(),id.getId())){
 219
 220       partlist()[id.getAlias()].setIdChgConj(id);
 221       return id;
 222     }
 223   }
 224
 225   if (id.getAlias()!=id.getId()) {
 226
 227     report(ERROR,"EvtGen")<<"Trying to charge conjugate alias particle:"
 228                           <<name(id).c_str()<<" without defining the alias!"<<endl;
 229
 230     ::abort();
 231
 232   }
 233
 234   for (size_t i=0;i<partlist().size();i++){
 235     if (partlist()[i].getStdHep()==-partlist()[id.getId()].getStdHep()){
 236       partlist()[id.getId()].setIdChgConj(partlist()[i].getId());
 237       return partlist()[i].getId();
 238     }
 239   }
 240
 241   partlist()[id.getId()].setIdChgConj(id);
 242   return id;
 243
 244 }
 245
 246 EvtId EvtPDL::evtIdFromStdHep(int stdhep){
 247
 248   for (size_t i=0;i<partlist().size();i++){
 249     if (partlist()[i].getStdHep()==stdhep)
 250       return partlist()[i].getId();
 251   }
 252
 253   return EvtId(-1,-1);
 254
 255 }
 256
 257
 258
 259 void EvtPDL::alias(EvtId num,const std::string& newname){
 260
 261   if ( _firstAlias < partlist().size() ) {
 262     for(size_t i=_firstAlias;i<partlist().size();i--){
 263       if (newname==partlist()[i].getName()){
 264         report(WARNING,"EvtGen")<<"Redefining alias:"<<newname.c_str()<<" will be ignored!"<<endl;
 265         return;
 266       }
 267     }
 268   }
 269   else{
 270     _firstAlias=partlist().size();
 271   }
 272
 273   partlist().push_back(partlist()[num.getId()]);
 274   int entry=partlist().size()-1;
 275   partlist()[entry].setName(newname);
 276   if (_particleNameLookup.find(std::string(newname))!=
 277       _particleNameLookup.end()){
 278             report(ERROR,"EvtGen")<<"The particle name:"<<newname<<" is already defined."<<endl;
 279             report(ERROR,"EvtGen") << "Will terminate execution.";
 280             ::abort();
 281   }
 282   _particleNameLookup[std::string(newname)]=entry;
 283   partlist()[entry].setId(EvtId(num.getId(),entry));
 284   //Lange - Dec7, 2003. Unset the charge conjugate.
 285   partlist()[entry].setIdChgConj(EvtId(-1,-1));
 286
 287 }
 288
 289 EvtId EvtPDL::getId(const std::string& name ){
 290
 291   std::map<std::string,int>::iterator it=_particleNameLookup.find(std::string(name));
 292   if (it==_particleNameLookup.end()) return EvtId(-1,-1);
 293
 294   return partlist()[it->second].getId();
 295
 296 }
 297
 298 void EvtPDL::setUpConstsPdt(){
 299
 300 }
 301
 302
 303 // Function to get EvtId from LundKC ( == Pythia Hep Code , KF )
 304 EvtId EvtPDL::evtIdFromLundKC(int pythiaId){
 305
 306   unsigned int i;
 307
 308   for (i=0;i<partlist().size();i++){
 309     if (partlist()[i].getLundKC()==pythiaId)
 310       return partlist()[i].getId();
 311   }
 312
 313   return EvtId(-1,-1);
 314
 315 }
 316
 317 double EvtPDL::getMeanMass(EvtId i ) {
 318   return partlist()[i.getId()].getMass();
 319 }
 320
 321 double EvtPDL::getMass(EvtId i ) {
 322   return partlist()[i.getId()].rollMass();
 323 }
 324
 325 double EvtPDL::getRandMass(EvtId i, EvtId *parId, int nDaug, EvtId *dauId,
 326                            EvtId *othDaugId,double maxMass,
 327                            double *dauMasses ) {
 328   return partlist()[i.getId()].getRandMass(parId,nDaug,dauId,
 329                                           othDaugId,maxMass,dauMasses);
 330 }
 331
 332 double EvtPDL::getMassProb(EvtId i, double mass, double massPar, int nDaug,
 333                            double *massDau) {
 334   return partlist()[i.getId()].getMassProb(mass,massPar,nDaug,massDau);
 335 }
 336
 337 double EvtPDL::getMaxMass(EvtId i ) {
 338   return partlist()[i.getId()].getMassMax();
 339 }
 340
 341 double EvtPDL::getMinMass(EvtId i ) {
 342   return partlist()[i.getId()].getMassMin();
 343 }
 344
 345 double EvtPDL::getMaxRange(EvtId i ) {
 346   return partlist()[i.getId()].getMaxRange();
 347 }
 348
 349 double EvtPDL::getWidth(EvtId i ) {
 350   return partlist()[i.getId()].getWidth();
 351 }
 352
 353 double EvtPDL::getctau(EvtId i ) {
 354   return partlist()[i.getId()].getctau();
 355 }
 356
 357 int EvtPDL::getStdHep(EvtId id ) {
 358   return partlist()[id.getId()].getStdHep();
 359 }
 360
 361 int EvtPDL::getLundKC(EvtId id ) {
 362   return partlist()[id.getId()].getLundKC();
 363 }
 364
 365 int EvtPDL::chg3(EvtId i ) {
 366   return partlist()[i.getId()].getChg3();
 367 }
 368
 369 EvtSpinType::spintype EvtPDL::getSpinType(EvtId i ) {
 370   return partlist()[i.getId()].getSpinType();
 371 }
 372
 373 std::string EvtPDL::name(EvtId i) {
 374   return partlist()[i.getAlias()].getName();
 375 }
 376
 377 size_t EvtPDL::entries() {
 378   return partlist().size();
 379 }
 380
 381 EvtId EvtPDL::getEntry(int i) {
 382   return partlist()[i].getId();
 383 }
 384
 385 void EvtPDL::reSetMass(EvtId i, double mass) {
 386   partlist()[i.getId()].reSetMass(mass);
 387 }
 388
 389 void EvtPDL::reSetWidth(EvtId i, double width) {
 390   partlist()[i.getId()].reSetWidth(width);
 391 }
 392
 393 void EvtPDL::reSetMassMin(EvtId i, double mass) {
 394   partlist()[i.getId()].reSetMassMin(mass);
 395 }
 396
 397 void EvtPDL::reSetMassMax(EvtId i,double mass) {
 398   partlist()[i.getId()].reSetMassMax(mass);
 399 }
 400
 401 void EvtPDL::reSetBlatt(EvtId i,double blatt) {
 402   partlist()[i.getId()].reSetBlatt(blatt);
 403 }
 404
 405 void EvtPDL::includeBirthFactor(EvtId i,bool yesno) {
 406   partlist()[i.getId()].includeBirthFactor(yesno);
 407 }
 408
 409 void EvtPDL::includeDecayFactor(EvtId i,bool yesno) {
 410   partlist()[i.getId()].includeDecayFactor(yesno);
 411 }
 412
 413 void EvtPDL::changeLS(EvtId i, std::string &newLS ) {
 414   partlist()[i.getId()].newLineShape(newLS);
 415 }
 416
 417 void EvtPDL::setPWForDecay(EvtId i, int spin, EvtId d1, EvtId d2) {
 418   partlist()[i.getId()].setPWForDecay(spin,d1,d2);
 419 }
 420
 421 void EvtPDL::setPWForBirthL(EvtId i, int spin, EvtId par, EvtId othD) {
 422   partlist()[i.getId()].setPWForBirthL(spin,par,othD);
 423 }