G4INCLParticleSpecies.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
// INCL++ intra-nuclear cascade model
// Alain Boudard, CEA-Saclay, France
// Joseph Cugnon, University of Liege, Belgium
// Jean-Christophe David, CEA-Saclay, France
// Pekka Kaitaniemi, CEA-Saclay, France, and Helsinki Institute of Physics, Finland
// Sylvie Leray, CEA-Saclay, France
// Davide Mancusi, CEA-Saclay, France
//
#define INCLXX_IN_GEANT4_MODE 1
 
#include "globals.hh"
 
/*
 * G4INCLParticleSpecies.cc
 *
 *  \date Nov 25, 2011
 * \author Davide Mancusi
 */
 
#include "G4INCLParticleSpecies.hh"
#include "G4INCLParticleTable.hh"
#include <algorithm>
#include <cctype>
#include <sstream>
#include <algorithm>
 
namespace G4INCL {
 
  ParticleSpecies::ParticleSpecies(std::string const &pS) {
    // Normalise the string to lower case
    if(pS=="p" || pS=="proton") {
      theA = 1;
      theZ = 1;
      theS = 0;
      theType = G4INCL::Proton;
    } else if(pS=="n" || pS=="neutron") {
      theA = 1;
      theZ = 0;
      theS = 0;
      theType = G4INCL::Neutron;
    } else if(pS=="delta++" || pS=="deltaplusplus") {
      theA = 1;
      theZ = 2;
      theS = 0;
      theType = G4INCL::DeltaPlusPlus;
    } else if(pS=="delta+" || pS=="deltaplus") {
      theA = 1;
      theZ = 1;
      theS = 0;
      theType = G4INCL::DeltaPlus;
    } else if(pS=="delta0" || pS=="deltazero") {
      theA = 1;
      theZ = 0;
      theS = 0;
      theType = G4INCL::DeltaZero;
    } else if(pS=="delta-" || pS=="deltaminus") {
      theA = 1;
      theZ = -1;
      theS = 0;
      theType = G4INCL::DeltaMinus;
    } else if(pS=="pi+" || pS=="pion+" || pS=="piplus" || pS=="pionplus") {
      theA = 0;
      theZ = 1;
      theS = 0;
      theType = G4INCL::PiPlus;
    } else if(pS=="pi0" || pS=="pion0" || pS=="pizero" || pS=="pionzero") {
      theA = 0;
      theZ = 0;
      theS = 0;
      theType = G4INCL::PiZero;
    } else if(pS=="pi-" || pS=="pion-" || pS=="piminus" || pS=="pionminus") {
      theA = 0;
      theZ = -1;
      theS = 0;
      theType = G4INCL::PiMinus;
    } else if(pS=="lambda" || pS=="l" || pS=="l0")  {
      theA = 1;
      theZ = 0;
      theS = -1;
      theType = G4INCL::Lambda;
    } else if(pS=="s+" || pS=="sigma+" || pS=="sigmaplus")  {
      theA = 1;
      theZ = 1;
      theS = -1;
      theType = G4INCL::SigmaPlus;
    } else if(pS=="s0" || pS=="sigma0" || pS=="sigmazero")  {
      theA = 1;
      theZ = 0;
      theS = -1;
      theType = G4INCL::SigmaZero;
    } else if(pS=="s-" || pS=="sigma-" || pS=="sigmaminus")  { //Sm = Samarium
      theA = 1;
      theZ = -1;
      theS = -1;
      theType = G4INCL::SigmaMinus; 
    } else if(pS=="xi-" || pS=="x-") {
      theA = 1;
      theZ = -1;
      theS = -2;
      theType = G4INCL::XiMinus; 
    } else if(pS=="xi0" || pS=="x0") {
      theA = 1;
      theZ = 0;
      theS = -2;
      theType = G4INCL::XiZero;  
    } else if(pS=="pb" || pS=="pbar" || pS=="antiproton") {
      theA = -1;
      theZ = -1;
      theS = 0;
      theType = G4INCL::antiProton;   
    } else if(pS=="nb" || pS=="nbar" || pS=="antineutron") {
      theA = -1;
      theZ = 0;
      theS = 0;
      theType = G4INCL::antiNeutron;
    } else if(pS=="s+b" || pS=="antisigma+" || pS=="antisigmaplus")  {
      theA = -1;
      theZ = -1;
      theS = 1;
      theType = G4INCL::antiSigmaPlus;
    } else if(pS=="s0b" || pS=="antisigma0" || pS=="antisigmazero")  {
      theA = -1;
      theZ = 0;
      theS = 1;
      theType = G4INCL::antiSigmaZero;
    } else if(pS=="s-b" || pS=="antisigma-" || pS=="antisigmaminus")  { //Sm = Samarium; Whats wrong with the sign?
      theA = -1;
      theZ = 1;
      theS = 1;
      theType = G4INCL::antiSigmaMinus;
    } else if(pS=="antilambda" || pS=="lb" || pS=="l0b")  {
      theA = -1;
      theZ = 0;
      theS = 1;
      theType = G4INCL::antiLambda;
    } else if(pS=="antixi-" || pS=="x-b") {
      theA = -1;
      theZ = 1;
      theS = 2;
      theType = G4INCL::antiXiMinus; 
    } else if(pS=="antixi0" || pS=="x0b") {
      theA = -1;
      theZ = 0;
      theS = 2;
      theType = G4INCL::antiXiZero;  
    } else if(pS=="k+" || pS=="kaon+" || pS=="kplus" || pS=="kaonplus") {
      theA = 0;
      theZ = 1;
      theS = 1;
      theType = G4INCL::KPlus;
    } else if(pS=="k0" || pS=="kaon0" || pS=="kzero" || pS=="kaonzero") {
      theA = 0;
      theZ = 0;
      theS = 1;
      theType = G4INCL::KZero;
    } else if(pS=="k0b" || pS=="kzb" || pS=="kaon0bar" || pS=="kzerobar" || pS=="kaonzerobar") {
      theA = 0;
      theZ = 0;
      theS = -1;
      theType = G4INCL::KZeroBar;
    } else if(pS=="k-" || pS=="kaon-" || pS=="kminus" || pS=="kaonminus") {
      theA = 0;
      theZ = -1;
      theS = -1;
      theType = G4INCL::KMinus;
    } else if(pS=="k0s" || pS=="kshort" || pS=="ks" || pS=="kaonshort") {
      theA = 0;
      theZ = 0;
//    theS not defined
      theType = G4INCL::KShort;
    } else if(pS=="k0l" || pS=="klong" || pS=="kl" || pS=="kaonlong") {
      theA = 0;
      theZ = 0;
//    theS not defined
      theType = G4INCL::KLong;
    } else if(pS=="d" || pS=="deuteron") {
      theA = 2;
      theZ = 1;
      theS = 0;
      theType = G4INCL::Composite;
    } else if(pS=="t" || pS=="triton") {
      theA = 3;
      theZ = 1;
      theS = 0;
      theType = G4INCL::Composite;
    } else if(pS=="a" || pS=="alpha") {
      theA = 4;
      theZ = 2;
      theS = 0;
      theType = G4INCL::Composite;
    } else if(pS=="eta") {
      theA = 0;
      theZ = 0;
      theS = 0;
      theType = G4INCL::Eta;
    } else if(pS=="omega") {
      theA = 0;
      theZ = 0;
      theS = 0;
      theType = G4INCL::Omega;
    } else if(pS=="etaprime" || pS=="etap") {
      theA = 0;
      theZ = 0;
      theS = 0;
      theType = G4INCL::EtaPrime;
    } else if(pS=="photon") {
      theA = 0;
      theZ = 0;
      theS = 0;
      theType = G4INCL::Photon;    
    } else if (pS=="db" || pS=="dbar" || pS=="antideuteron"){
      theA = -2;
      theZ = -1;
      theS = 0;
      theType = G4INCL::antiComposite;     
    } else
      parseNuclide(pS);
  }
 
  ParticleSpecies::ParticleSpecies(ParticleType const t) :
    theType(t),
    theA(ParticleTable::getMassNumber(theType)),
    theZ(ParticleTable::getChargeNumber(theType)),
    theS(ParticleTable::getStrangenessNumber(theType))
  {}
 
  ParticleSpecies::ParticleSpecies(const int A, const int Z){
    if (A>=0){
      theType = Composite;
      theA = A;
      theZ = Z;
      theS = 0;
    } else {
      theType = antiComposite;
      theA = A;
      theZ = Z;
      theS = 0;
    }
  }
 
  ParticleSpecies::ParticleSpecies(const int A, const int Z, const int S){
    if (A>=0){
      theType = Composite;
      theA = A;
      theZ = Z;
      theS = S;
    } else {
      theType = antiComposite;
      theA = A;
      theZ = Z;
      theS = S;
    }
  }
 
  void ParticleSpecies::parseNuclide(std::string const &pS) {
    theType = Composite;
    theS = 0; // no hypernuclei projectile or target for now
 
    // Allowed characters
    const std::string separators("-_");
    std::string allowed("0123456789abcdefghijklmnopqrstuvwxyzABCDEFGHIJKLMNOPQRSTUVWXYZ");
    allowed += separators;
 
    // There must be at least one character
    if(pS.find_first_not_of(allowed)!=std::string::npos) {
      // Malformed input string
      // Setting unknown particle species
      (*this) = ParticleSpecies(UnknownParticle);
      return;
    }
    if(pS.size()<1) {
      // Malformed input string
      // Setting unknown particle species
      (*this) = ParticleSpecies(UnknownParticle);
      return;
    }
 
    std::size_t firstSeparator = pS.find_first_of(separators);
    std::size_t lastSeparator = pS.find_last_of(separators);
    if(firstSeparator!=std::string::npos && firstSeparator!=lastSeparator) {
      // Several separators in malformed input string
      // Setting unknown particle species
      (*this) = ParticleSpecies(UnknownParticle);
      return;
    }
 
    // Identify the type of the first character
    G4int (*predicate)(G4int);
    G4bool startsWithAlpha = std::isalpha(pS.at(0));
    if(startsWithAlpha) {
      predicate=std::isdigit;
    } else if(std::isdigit(pS.at(0))) {
      predicate=std::isalpha;
    } else {
      // Non-alphanumeric character in string
      // Setting unknown particle species
      (*this) = ParticleSpecies(UnknownParticle);
      return;
    }
 
    G4bool hasIsotope = true;
    size_t endFirstSection, beginSecondSection;
    if(firstSeparator==std::string::npos) {
      // No separator, Fe56 or 56Fe style
      // Identify the end of the first section
 
      // Find the first character that is not of the same type as the first one
      beginSecondSection = std::find_if(pS.begin()+1, pS.end(), predicate) - pS.begin();
 
      if(beginSecondSection>=pS.size()) {
        if(startsWithAlpha) {
          // Only alphabetic characters are present -- must be an element name
          hasIsotope = false;
        } else {
          // Only numeric characters in the string
          // Setting unknown particle species
          (*this) = ParticleSpecies(UnknownParticle);
          return;
        }
      }
 
      endFirstSection = beginSecondSection;
 
    } else {
      // One separator, Fe-56 or 56-Fe style or hypercluster style: Fe56-1 (iron 56 including 1 lambda)
      endFirstSection = firstSeparator;
      beginSecondSection = firstSeparator+1;
    }
 
    std::string firstSection(pS.substr(0,endFirstSection));
    std::string secondSection(pS.substr(beginSecondSection,std::string::npos));
    std::stringstream parsingStream;
 
    if(std::isalpha(firstSection.at(0)) && std::isdigit(firstSection.at(endFirstSection-1))) { // Hypernucleus, must be Fe56-1 style
      std::stringstream parseStrangeness;
      parseStrangeness.str(secondSection);
      parseStrangeness >> theS;
      if(parsingStream.fail()) {
        // Couldn't parse the strange charge section
        // Setting unknown particle species
        (*this) = ParticleSpecies(UnknownParticle);
        return;
      }
      theS *= (-1);
      beginSecondSection = std::find_if(pS.begin()+1, pS.end(), predicate) - pS.begin(); // predicate == std::isdigit(G4int) in this case
      firstSection = pS.substr(0, beginSecondSection);
      secondSection = pS.substr(beginSecondSection, endFirstSection);
    }
 
    // Parse the sections
    G4bool success;
    if(startsWithAlpha) {
      parsingStream.str(secondSection);
      success = parseElement(firstSection);
    } else {
      parsingStream.str(firstSection);
      success = parseElement(secondSection);
    }
    if(!success) {
      // Couldn't parse the element section
      // Setting unknown particle species
      (*this) = ParticleSpecies(UnknownParticle);
      return;
    }
 
    if(hasIsotope) {
      parsingStream >> theA;
      if(parsingStream.fail()) {
        // Couldn't parse the mass section
        // Setting unknown particle species
        (*this) = ParticleSpecies(UnknownParticle);
        return;
      }
    } else
      theA = 0;
 
    // Check that Z<=A
    if(theZ>theA && hasIsotope) {
      // Setting unknown particle species
      (*this) = ParticleSpecies(UnknownParticle);
      return;
    }
 
    // Special particle type for protons
    if(theZ==1 && theA==1)
      theType = Proton;
  }
 
  G4bool ParticleSpecies::parseElement(std::string const &s) {
    theZ = ParticleTable::parseElement(s);
 
    if(theZ<0)
      theZ = ParticleTable::parseIUPACElement(s);
 
    if(theZ<0)
      return false;
    else
      return true;
  }
 
  G4bool ParticleSpecies::parseIUPACElement(std::string const &s) {
    theZ = ParticleTable::parseIUPACElement(s);
    if(theZ==0)
      return false;
    else
      return true;
  }
    
  G4int ParticleSpecies::getPDGCode() const {
    switch (theType) {
        case Proton:
            return 2212;
            break;
        case Neutron:
            return 2112;
            break;
        case DeltaPlusPlus:
            return 2224;
            break;
        case DeltaPlus:
            return 2214;
            break;
        case DeltaZero:
            return 2114;
            break;
        case DeltaMinus:
            return 1114;
            break;
        case PiPlus:
            return 211;
            break;
        case PiZero:
            return 111;
            break;
        case PiMinus:
            return -211;
            break;
        case Eta:
            return 221;
            break;
        case Omega:
            return 223;
            break;
        case EtaPrime:
            return 331;
            break;
        case Photon:
            return 22;
            break;
        case Lambda:
            return 3122;
            break;
        case SigmaPlus:
            return 3222;
            break;
        case SigmaZero:
            return 3212;
            break;
        case SigmaMinus:
            return 3112;
            break;          
      case antiProton:
        return -2212;
      break;
        case XiMinus:
            return 3312;
            break;
        case XiZero:
            return 3322;
            break;          
        case antiNeutron:
            return -2112;
            break;
        case antiLambda:
            return -3122;
            break;
        case antiSigmaPlus:
            return -3222;
            break;
        case antiSigmaZero:
            return -3212;
            break;
        case antiSigmaMinus:
            return -3112;
            break;
        case antiXiMinus:
            return -3312;
            break;
        case antiXiZero:
            return -3322;
            break;              
        case KPlus:
            return 321;
            break;
        case KZero:
            return 311;
            break;
        case KZeroBar:
            return -311;
            break;
        case KShort:
            return 310;
            break;
        case KLong:
            return 130;
            break;
        case KMinus:
            return -321;
            break;      
        case Composite:
            if(theA == 1 && theZ == 1 && theS == 0) return 2212;
            else if(theA == 1 && theZ == 0 && theS == 0) return 2112;
            else if(theA == 1 && theZ == 0 && theS == -1) return 3122;
            else return theA+theZ*1000-theS*1e6; // Here -theS because hyper-nucleus -> theS < 0
            break;
    case antiComposite:
        if(theA == 1 && theZ == 1 && theS == 0) return -2212;
        else if(theA == 1 && theZ == 0 && theS == 0) return -2112;
        else if(theA == 1 && theZ == 0 && theS == -1) return -3122;
        else return -(theA + theZ*1000 - theS*1e6);
      break;
        default:
            INCL_ERROR("ParticleSpecies::getPDGCode: Unknown particle type." << '\n');
            return 0;
            break;
    }   
  }
}