Beam.cpp

Go to the documentation of this file.

//
// Class Beam
//   The class for the OPAL BEAM command.
//   A BEAM definition is used by most physics commands to define the
//   particle charge and the reference momentum, together with some other data.
//
// Copyright (c) 200x - 2021, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Structure/Beam.h"
 
#include "AbstractObjects/Expressions.h"
#include "AbstractObjects/OpalData.h"
#include "Attributes/Attributes.h"
#include "Expressions/SAutomatic.h"
#include "Expressions/SRefExpr.h"
#include "Physics/ParticleProperties.h"
#include "Physics/Physics.h"
#include "Physics/Units.h"
#include "Utilities/OpalException.h"
#include "Utilities/Options.h"
 
#include <cmath>
#include <iterator>
 
using namespace Expressions;
 
// The attributes of class Beam.
namespace {
    constexpr const char* photonParticleName = "PHOTON";
 
    enum {
        PARTICLE,         // The particle name
        MASS,             // The particle rest mass in GeV
        CHARGE,           // The particle charge in proton charges
        ENERGY,           // The particle energy in GeV
        PC,               // The particle momentum in GeV/c
        GAMMA,            // ENERGY / MASS
        BCURRENT,         // Legacy, unused in OPALX (holdover from OPALCycl)
        BFREQ,            // Legacy, unused in OPALX (holdover from OPALCycl)
        BCHARGE,          // Bunch charge in C
        NALLOC,           // Allocation size (macroparticles) for this beam
        SOURCES,          // Name of EMISSIONSOURCELIST
        GLOBALPROCESSES,  // Global physics processes active for this beam
        DAUGHTERBEAM,     // Name of the beam that receives decay daughter particles
        POLARIZATION,     // Initial polarization vector P (rest-frame, lab-frame axes)
        SIZE
    };
}  // namespace
 
Beam::Beam()
    : Definition(
              SIZE, "BEAM",
              "The \"BEAM\" statement defines data for the particles "
              "in a beam."),
      reference(1.0, Physics::m_p * Units::GeV2eV, 1.0 * Units::GeV2eV) {
    itsAttr[PARTICLE] = Attributes::makePredefinedString(
            "PARTICLE", "Name of particle to be used",
            {"PHOTON", "ELECTRON", "POSITRON", "MUON", "PION", "PROTON", "ANTIPROTON", "DEUTERON",
             "HMINUS", "H2P", "ALPHA", "CARBON", "XENON", "URANIUM"});
 
    itsAttr[MASS] = Attributes::makeReal("MASS", "Particle rest mass [GeV]");
 
    itsAttr[CHARGE] = Attributes::makeReal("CHARGE", "Particle charge in proton charges");
 
    itsAttr[ENERGY] = Attributes::makeReal("ENERGY", "Particle energy [GeV]");
 
    itsAttr[PC] = Attributes::makeReal("PC", "Particle momentum [GeV/c]");
 
    itsAttr[GAMMA] = Attributes::makeReal("GAMMA", "ENERGY / MASS");
 
    itsAttr[BCURRENT] =
            Attributes::makeReal("BCURRENT", "Legacy, unused in OPALX. Use BCHARGE instead.");
 
    itsAttr[BFREQ] = Attributes::makeReal("BFREQ", "Legacy, unused in OPALX. Use BCHARGE instead.");
 
    itsAttr[BCHARGE] = Attributes::makeReal("BCHARGE", "Bunch charge [C]");
 
    itsAttr[NALLOC] =
            Attributes::makeReal("NALLOC", "Allocation size (macroparticles) for this beam");
 
    itsAttr[SOURCES] = Attributes::makeString(
            "SOURCES", "Name of the emission sources list (EMISSIONSOURCELIST).");
 
    itsAttr[GLOBALPROCESSES] = Attributes::makeUpperCaseStringArray(
            "GLOBALPROCESSES", "Global physics processes active for this beam.");
 
    itsAttr[DAUGHTERBEAM] = Attributes::makeString(
            "DAUGHTERBEAM", "Name of the BEAM that receives decay daughter particles.");
 
    itsAttr[POLARIZATION] = Attributes::makeRealArray(
            "POLARIZATION",
            "Initial polarization vector P = {Px, Py, Pz} for this beam "
            "(rest-frame Pauli expectation values along lab-frame axes). "
            "Must be length-3 with |P| in [0, 1]. Only valid for PARTICLE=MUON; "
            "setting POLARIZATION on any other species is rejected. Setting "
            "POLARIZATION (even {0,0,0}) enables per-particle spin tracking for this "
            "beam (Thomas-BMT integration); leaving it unset disables spin tracking. "
            "For muons produced by an upstream decay (e.g. PionDecay) the per-particle "
            "value is overwritten by the decay, but POLARIZATION must still be set to "
            "enable the spin storage that receives it.");
 
    // Set up default beam.
    Beam* defBeam    = clone("UNNAMED_BEAM");
    defBeam->builtin = true;
 
    try {
        defBeam->update();
        OpalData::getInstance()->define(defBeam);
    } catch (...) {
        delete defBeam;
    }
 
    registerOwnership(AttributeHandler::STATEMENT);
}
 
Beam::Beam(const std::string& name, Beam* parent)
    : Definition(name, parent), reference(parent->reference) {}
 
Beam::~Beam() {}
 
bool Beam::canReplaceBy(Object* object) {
    // Can replace only by another BEAM.
    return dynamic_cast<Beam*>(object) != 0;
}
 
Beam* Beam::clone(const std::string& name) { return new Beam(name, this); }
 
void Beam::execute() {
    const bool photon = itsAttr[PARTICLE] && getParticleName() == photonParticleName;
 
    if (photon) {
        if (!itsAttr[ENERGY]) {
            throw OpalException("Beam::execute()", "\"ENERGY\" must be set for PARTICLE=PHOTON.");
        }
        if (itsAttr[MASS]) {
            throw OpalException(
                    "Beam::execute()",
                    "\"MASS\" is not allowed for PARTICLE=PHOTON. Use \"ENERGY\".");
        }
        if (itsAttr[CHARGE]) {
            throw OpalException(
                    "Beam::execute()", "\"CHARGE\" is not allowed for PARTICLE=PHOTON.");
        }
        if (itsAttr[PC]) {
            throw OpalException(
                    "Beam::execute()",
                    "\"PC\" is not allowed for PARTICLE=PHOTON. Use \"ENERGY\".");
        }
        if (itsAttr[GAMMA]) {
            throw OpalException(
                    "Beam::execute()",
                    "\"GAMMA\" is not allowed for PARTICLE=PHOTON. Use \"ENERGY\".");
        }
        if (itsAttr[SOURCES]) {
            throw OpalException(
                    "Beam::execute()", "\"SOURCES\" is not allowed for PARTICLE=PHOTON.");
        }
    }
 
    if (itsAttr[BCURRENT] || itsAttr[BFREQ]) {
        throw OpalException(
                "Beam::execute()",
                "\"BCURRENT\" and \"BFREQ\" are no longer used in OPALX. "
                "Use \"BCHARGE\" [C] to specify the bunch charge directly.");
    }
 
    update();
 
    if (!(itsAttr[PARTICLE]) && (!itsAttr[MASS] || !(itsAttr[CHARGE]))) {
        throw OpalException(
                "Beam::execute()",
                "The beam particle hasn't been set. "
                "Set either \"PARTICLE\" or \"MASS\" and \"CHARGE\".");
    }
 
    if (!(itsAttr[NALLOC])) {
        throw OpalException("Beam::execute()", "\"NALLOC\" must be set.");
    }
 
    if (photon) {
        const double energy = Attributes::getReal(itsAttr[ENERGY]);
        if (energy <= 0.0) {
            throw OpalException(
                    "Beam::execute()", "\"ENERGY\" should be greater than 0 for PARTICLE=PHOTON.");
        }
        return;
    }
 
    // Beam-only validation: each non-photon beam must specify its EMISSIONSOURCELIST.
    (void)getEmissionSourceListName();
 
    // Currently supported global process names (extend as new processes are implemented).
    for (const std::string& name : getGlobalProcessNames()) {
        if (name != "DECAY") {
            throw OpalException(
                    "Beam::execute()", "Unsupported entry in \"GLOBALPROCESSES\": \"" + name
                                               + "\". Supported values: DECAY.");
        }
    }
 
    validatePolarization();
}
 
void Beam::validatePolarization() const {
    if (!itsAttr[POLARIZATION]) {
        return;
    }
    const std::vector<double> pol = Attributes::getRealArray(itsAttr[POLARIZATION]);
    if (pol.empty()) {
        // Attribute exists in the registry but the user did not specify a value;
        // treat as "not set" and skip all polarization-specific validation.
        return;
    }
    if (pol.size() != 3) {
        throw OpalException(
                "Beam::execute()", "\"POLARIZATION\" must be a length-3 vector {Px, Py, Pz}.");
    }
    // Only muons are currently spin-tracked. Reject POLARIZATION on any other
    // species regardless of magnitude — setting it (even to zero) is meaningless
    // for spin-0 particles and for species whose g-factor isn't wired into T-BMT.
    if (!itsAttr[PARTICLE]) {
        throw OpalException(
                "Beam::execute()",
                "\"POLARIZATION\" can only be set when \"PARTICLE\" is also set.");
    }
    const ParticleType pType = ParticleProperties::getParticleType(getParticleName());
    if (pType != ParticleType::MUON) {
        throw OpalException(
                "Beam::execute()",
                "\"POLARIZATION\" is only supported for PARTICLE=MUON; got \""
                        + getParticleName()
                        + "\". For muons produced by PionDecay, leave POLARIZATION unset — "
                          "the daughter polarization is determined by the pion decay.");
    }
    const double pMag2 = pol[0] * pol[0] + pol[1] * pol[1] + pol[2] * pol[2];
    if (pMag2 > 1.0 + 1.0e-12) {
        throw OpalException(
                "Beam::execute()", "\"POLARIZATION\" magnitude must be in [0, 1]; got |P| = "
                                           + std::to_string(std::sqrt(pMag2)) + ".");
    }
}
 
std::string Beam::getEmissionSourceListName() const {
    if (!itsAttr[SOURCES]) {
        throw OpalException(
                "Beam::getEmissionSourceListName()",
                "\"SOURCES\" must be set for a beam (name of EMISSIONSOURCELIST).");
    }
 
    const std::string name = Attributes::getString(itsAttr[SOURCES]);
    if (name.empty()) {
        throw OpalException(
                "Beam::getEmissionSourceListName()",
                "\"SOURCES\" must not be empty for a beam (name of EMISSIONSOURCELIST).");
    }
    return name;
}
 
std::vector<std::string> Beam::getGlobalProcessNames() const {
    return Attributes::getStringArray(itsAttr[GLOBALPROCESSES]);
}
 
std::string Beam::getDaughterBeamName() const {
    return Attributes::getString(itsAttr[DAUGHTERBEAM]);
}
 
std::vector<double> Beam::getPolarization() const {
    std::vector<double> pol = Attributes::getRealArray(itsAttr[POLARIZATION]);
    if (pol.empty()) {
        return {0.0, 0.0, 0.0};
    }
    return pol;
}
 
bool Beam::hasPolarization() const {
    return !Attributes::getRealArray(itsAttr[POLARIZATION]).empty();
}
 
Beam* Beam::find(const std::string& name) {
    Beam* beam = dynamic_cast<Beam*>(OpalData::getInstance()->find(name));
 
    if (beam == 0) {
        throw OpalException("Beam::find()", "Beam \"" + name + "\" not found.");
    }
 
    return beam;
}
 
size_t Beam::getNumAlloc() const {
    if (Attributes::getReal(itsAttr[NALLOC]) > 0) {
        return (size_t)Attributes::getReal(itsAttr[NALLOC]);
    } else {
        throw OpalException(
                "Beam::getNumAlloc()",
                "Wrong allocation size for beam! \"NALLOC\" must be positive");
    }
}
 
const PartData& Beam::getReference() const {
    // Cast away const, to allow logically constant Beam to update.
    const_cast<Beam*>(this)->update();
    return reference;
}
 
double Beam::getCurrent() const { return Attributes::getReal(itsAttr[BCURRENT]); }
 
double Beam::getBunchCharge() const { return Attributes::getReal(itsAttr[BCHARGE]); }
 
double Beam::getCharge() const { return Attributes::getReal(itsAttr[CHARGE]); }
 
double Beam::getMass() const { return Attributes::getReal(itsAttr[MASS]); }
 
double Beam::getMomentum() const { return reference.getP() / 1.e9; }
 
std::string Beam::getParticleName() const { return Attributes::getString(itsAttr[PARTICLE]); }
 
bool Beam::isPhoton() const { return itsAttr[PARTICLE] && getParticleName() == photonParticleName; }
 
double Beam::getFrequency() const { return Attributes::getReal(itsAttr[BFREQ]); }
 
bool Beam::hasExplicitEnergy() const { return itsAttr[GAMMA] || itsAttr[ENERGY] || itsAttr[PC]; }
 
double Beam::getChargePerParticle() const {
    return std::copysign(1.0, getCharge()) * getBunchCharge() / getNumAlloc();
}
 
double Beam::getMassPerParticle() const {
    return getMass() * getChargePerParticle() / (getCharge() * Physics::q_e);
}
 
void Beam::update() {
    if (itsAttr[PARTICLE]) {
        std::string pName  = getParticleName();
        ParticleType pType = ParticleProperties::getParticleType(pName);
        if (!itsAttr[MASS]) {
            Attributes::setReal(itsAttr[MASS], ParticleProperties::getParticleMass(pType));
        }
        if (!itsAttr[CHARGE]) {
            Attributes::setReal(itsAttr[CHARGE], ParticleProperties::getParticleCharge(pType));
        }
    }
 
    if (isPhoton()) {
        reference = PartData();
        reference.setQ(0.0);
        reference.setM(0.0);
        return;
    }
 
    // Set up particle reference; convert all to eV for OPALX.
    double mass   = (itsAttr[MASS] ? getMass() : Physics::m_p) * Units::GeV2eV;
    double charge = itsAttr[CHARGE] ? getCharge() : 1.0;
 
    reference = PartData(charge, mass, 1.0);
 
    // Magnetic moment anomaly G = (g-2)/2 for the species (used by T-BMT).
    if (itsAttr[PARTICLE]) {
        const ParticleType pType = ParticleProperties::getParticleType(getParticleName());
        reference.setAnomaly(ParticleProperties::getParticleAnomaly(pType));
    }
 
    if (itsAttr[GAMMA]) {
        double gamma = Attributes::getReal(itsAttr[GAMMA]);
        if (gamma > 1.0) {
            reference.setGamma(gamma);
        } else {
            throw OpalException("Beam::update()", "\"GAMMA\" should be greater than 1.");
        }
    } else if (itsAttr[ENERGY]) {
        double energy = Attributes::getReal(itsAttr[ENERGY]) * Units::GeV2eV;
        if (energy > reference.getM()) {
            reference.setE(energy);
        } else {
            throw OpalException("Beam::update()", "\"ENERGY\" should be greater than \"MASS\".");
        }
    } else if (itsAttr[PC]) {
        double pc = Attributes::getReal(itsAttr[PC]) * Units::GeV2eV;
        if (pc > 0.0) {
            reference.setP(pc);
        } else {
            throw OpalException("Beam::update()", "\"PC\" should be greater than 0.");
        }
    }
 
    // Set default name.
    if (getOpalName().empty()) setOpalName("UNNAMED_BEAM");
}
 
void Beam::print(std::ostream& os) const {
    double charge = Attributes::getReal(itsAttr[CHARGE]);
    os << "* ************* B E A M ************************************************************ "
       << std::endl;
    os << "* BEAM        " << getOpalName() << '\n'
       << "* PARTICLE    " << Attributes::getString(itsAttr[PARTICLE]) << '\n'
       << "* REST MASS   " << Attributes::getReal(itsAttr[MASS]) << " [GeV]\n"
       << "* CHARGE      " << (charge > 0 ? '+' : '-') << "e * " << std::abs(charge) << " \n"
       << "* MOMENTUM    " << reference.getP() << " [eV/c]\n"
       << "* MOMENTUM    " << Attributes::getReal(itsAttr[PC]) << " [GeV/c]\n"
       << "* BCHARGE     " << Attributes::getReal(itsAttr[BCHARGE]) << " [C]\n"
       << "* NALLOC      " << Attributes::getReal(itsAttr[NALLOC]) << '\n';
    os << "* ********************************************************************************** "
       << std::endl;
}