---

HepPDT.cxx

Go to the documentation of this file.

// -*- C++ -*-  $Id: HepPDT.cxx,v 1.2 2000/01/23 03:37:53 pfkeb Exp $
//
// This file is part of Gismo 2
//
#include "HepPDT.h"
#include "gismo/Interactor.h"
#include "facilities/error.h"

#include <fstream>
#include <strstream>

#include <cstring>
#include <cstdio>
#include <cstdlib>

HepPDT::HepPDT(Interactor* emint, Interactor* hadint)
{
   if(emint)
          addInteractor(emint);
        if(hadint)
          addInteractor(hadint);

}
static char * fullName(char* name, int charge, int anti=0)
{
        static char* cname[]={"--","-","0","+","++"};
        static char buf1[20],buf2[20];    // required to keep these both
        static int which=0;
        which = ! which;
        char * buf = which ? buf1 : buf2;
        strcpy(buf,name);
        if( anti ) strcat(buf,"~");
        if(abs(charge)<3) strcat(buf,cname[charge+2]);
        return buf;
}

void
HepPDT::createTheParticles(char* name,  long* id,int *charge,float mass, float width)
{

        for( int i=0; i< 4 && id[i]>0; i++)
        {
                unsigned digit[5];
                for(long j=0, ii=id[0]; j< 5; j++, ii/=10)
                        digit[j] = unsigned(ii%10);

                enum particleType  {gauge, lepton, meson, baryon, diquark} type;

                if(      digit[3]    )  type = id[0]&1 ? diquark: baryon;
                else if( digit[2]    )  type = meson;
                else if( digit[1]==1 )  type = lepton;
                else                    type = gauge;


                int spin = digit[0];    // last digit sometimes a spin index
                switch ( type )
                {
                 case gauge:
                        if( id[i]==22 ) setInteraction("electromagnetic");
                        else            setInteraction("weak");
                        addParticle(name, id[i],  spin, charge[i], mass, width);
                        break;

                 case lepton:
                        setInteraction("weak");
                        if( charge[i]==0 )          // neutrino its own anti particle here
                         addParticle(name,id[i],1,0,mass,width);
                        else
                        {  if( id[i]==11 ) setInteraction("electromagnetic");   // electron
                                addParticle(fullName(name,-1),fullName(name,1),id[i],1,charge[i],mass,width);
                        }
                        break;

                 case meson:
                  setInteraction("hadronic");
                  if( charge[i]==0 )
                  {
                         if( digit[1]==digit[2])
                         {
                                if( digit[1]==1 )  //pi0, rho ...
                                 addParticle(fullName(name,0) ,id[i], spin,0, mass,width);
                          else   //  eta, phi, psi, ...
                                 addParticle(name ,id[i], spin,0, mass,width);
                         }
                         else
                         {  // K0 , K~0
                                addParticle(fullName(name,0),fullName(name,0,1),
                                        id[i],  spin, charge[i],  mass, width);
                         }
                  }
                  else // charged
                          addParticle(fullName(name,charge[i]),fullName(name,-charge[i]),
                                        id[i],  spin, charge[i],  mass, width);

                  break;


                case baryon:
                 setInteraction("hadronic");
                 spin--;
                 if( id[i]==2212 || id[i]==2112 )   // no charge for p,n
                 {
                         addParticle(name, fullName(name,99,1),
                                id[i], spin, charge[i],  mass, width );
                 }
                 else
                 {
                         addParticle(fullName(name,charge[i]),fullName(name,-charge[i],1),
                                        id[i], spin, charge[i], mass, width);
                 }
                 break;
                default:
                 FATAL(" unexpected particle type in PDT table ");
          }
        }

}
HepPDT&
HepPDT::readPDGMasses(const char* pdgMassFile)
{
  std::ifstream in(pdgMassFile);
  if( in.bad() )
  {
         std::ostrstream error;
         error << "could not open pdg mass file: " << pdgMassFile << '\0';
         FATAL(error.str());
  }
  HepPDT& ret = readPDGMasses(in);
  in.close();
  return ret;
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
//      member function that takes an existing PDT, and reads info from a
//      table of masses and widths
//
//      Format of table is:
//* Particle ID(s)           Value (GeV)    Errors (GeV)     Name          Charges
//   113   213              7.681E-01      +5.0E-04-5.0E-04 rho(770)          0,+
//          1         2         3         4         5         6         7
//01234567890123456789012345678901234567890123456789012345678901234567890123456789                                                              1          1             1


HepPDT&
HepPDT::readPDGMasses(std::istream& in)
{
        for(;;)  // read until eof
  {
         char  buf[82];
         long  id[4];         // array of ID's
         int   charge[4];     // array of charges
         float mass = -1.0;
         char  name[16];      // particle name
         char  lastLine = buf[0];  // keep track of previous line

         in.getline(buf,sizeof(buf));
         if( in.bad() )
                FATAL("file error, reading from pdg mass file");
         if( lastLine=='M')
         {
                 float width = (!in.eof() && buf[0]=='W') ? atof(buf+26) : 0;
                 createTheParticles(name,id,charge,mass,width);
         }
         if (in.eof() ) break;

         if( buf[0]=='M' )
         {  // first line: decode everything
                 id[0] = atol(buf+2);
                 id[1] =  buf[12]!=' '?  atol(buf+8) : 0;
                 id[2] =  buf[18]!=' '?  atol(buf+14): 0;
                 id[3] =  buf[24]!=' '?  atol(buf+20): 0;
                 mass = atof(buf+26);

                 // float error[2]={atof(buf+42),atof(buf+50)}; // the errors
                 char chargeList[9];
                 sscanf(buf+59,"%15s %8s",name,chargeList);  // name and list of charges

                 // decode the list of charges

                 int* ch = charge;
                 for(char* token = strtok(chargeList, ",");
                        token;
                        token = strtok(0, ","))
                 {
                        if( strcmp(token,"0")==0 )
                                *ch++ = 0;
                        else if (strcmp(token,"+")==0)
                                *ch++ = 1;
                        else if (strcmp(token,"-")==0)
                                *ch++ = -1;
                        else if (strcmp(token,"++")==0)
                                *ch++ = 2;
                        else
                                WARNING("bad charge found");
                 }
          }
  }
  return *this;
}

---