BeamStrippingPhysics.h

Go to the documentation of this file.

//
// Class BeamStrippingPhysics
//   Defines the physical processes of residual gas 
//   interactions and Lorentz stripping
//
// Copyright (c) 2018 - 2021, Pedro Calvo, CIEMAT, Spain
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Optimizing the radioisotope production of the novel AMIT
// superconducting weak focusing cyclotron"
//
// 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/>.
//
#ifndef BEAMSTRIPPINGPHYSICS_H
#define BEAMSTRIPPINGPHYSICS_H

#include "AbsBeamline/Component.h"
#include "Algorithms/PartBunchBase.h"
#include "Algorithms/Vektor.h"
#include "Solvers/ParticleMatterInteractionHandler.h"

#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>

#include <memory>
#include <string>
#include <utility>

template <class T, unsigned Dim>
class PartBunchBase;

class LossDataSink;
class Inform;
class Cyclotron;
class Vacuum;

class BeamStrippingPhysics: public ParticleMatterInteractionHandler {

public:
    BeamStrippingPhysics(const std::string& name, ElementBase* element);
    ~BeamStrippingPhysics();

    void setCyclotron(Cyclotron* cycl) { cycl_m = cycl; };

    virtual void apply(PartBunchBase<double, 3>* bunch,
                       const std::pair<Vector_t, double>& boundingSphere) override;

    virtual const std::string getType() const override;
    virtual void print(Inform& msg) override;
    virtual bool stillActive() override;

    virtual double getTime() override;
    virtual std::string getName() override;
    virtual size_t getParticlesInMat() override;
    virtual unsigned getRediffused() override;
    virtual unsigned int getNumEntered() override;

    void doPhysics(PartBunchBase<double, 3>* bunch);

private:

    void computeCrossSection(double energy);

    double computeCrossSectionNakai(double energy, double energyThreshold, int& i);
    
    double computeCrossSectionTabata(double energy, double energyThreshold,
                                     double a1, double a2, double a3,
                                     double a4, double a5, double a6);
                              
    double computeCrossSectionChebyshev(double energy, double energyMin, double energyMax);

    double computeCrossSectionBohr(double energy, int zTarget, double massInAmu);

    bool evalGasStripping(double& deltas);
    bool evalLorentzStripping(double& gamma, double& eField);

    void getSecondaryParticles(PartBunchBase<double, 3>* bunch, size_t& i, bool pdead_LS);

    void transformToSecondary(PartBunchBase<double, 3>* bunch, size_t& i, ParticleType type);

    bool computeEnergyLoss(PartBunchBase<double, 3>* /*bunch*/,
                           Vector_t& /*P*/,
                           const double /*deltat*/,
                           bool /*includeFluctuations*/) const override {
        return false;
    }

    void gatherStatistics();

    ParticleType pType_m;

    Vacuum* vac_m;
    Cyclotron* cycl_m;

    gsl_rng* r_m;

    double T_m;  // s
    double dT_m; // s
    double mass_m; // GeV/c2
    double pressure_m; // mbar
    double temperature_m; // K

    std::unique_ptr<LossDataSink> lossDs_m;

    double nCSA_m;
    double nCSB_m;
    double nCSC_m;
    double nCSTotal_m;

    unsigned bunchToMatStat_m;
    unsigned stoppedPartStat_m;
    unsigned rediffusedStat_m;
    unsigned totalPartsInMat_m;

    static const double csCoefSingle_Hminus[3][9];
    static const double csCoefDouble_Hminus[3][9];
    static const double csCoefSingle_Hplus[3][9];
    static const double csCoefDouble_Hplus[3][9];
    static const double csCoefSingleLoss_H[3][9];
    static const double csCoefSingleCapt_H[3][9];

    static const double csCoefSingle_Hplus_Tabata[10];
    static const double csCoefHminusProduction_H_Tabata[13];
    static const double csCoefProtonProduction_H_Tabata[9];
    static const double csCoefProtonProduction_H2plus_Tabata[11];
    static const double csCoefH3plusProduction_H2plus_Tabata[7];

    static const double csCoefSingle_Hminus_Chebyshev[11];
    static const double csCoefDouble_Hminus_Chebyshev[11];
    static const double csCoefSingle_Hplus_Chebyshev[11];
    static const double csCoefDouble_Hplus_Chebyshev[11];
    static const double csCoefHydrogenProduction_H2plus_Chebyshev[11];
    static const double csCoefProtonProduction_H2plus_Chebyshev[11];
    static const double energyRangeH2plusinH2[2];
    static double a_m[9];
    static double b_m[3][9];
};


inline
double BeamStrippingPhysics::getTime() {
    return T_m;
}

inline
std::string BeamStrippingPhysics::getName() {
    return (element_ref_m->getName() + "_" + name_m);
}

inline
size_t BeamStrippingPhysics::getParticlesInMat() {
    return totalPartsInMat_m;
}

inline
unsigned int BeamStrippingPhysics::getRediffused() {
    return rediffusedStat_m;
}

inline
unsigned int BeamStrippingPhysics::getNumEntered() {
    return bunchToMatStat_m;
}

inline
bool BeamStrippingPhysics::stillActive() {
    return totalPartsInMat_m != 0;
}

inline
const std::string BeamStrippingPhysics::getType() const {
    return "BeamStrippingPhysics";
}

#endif //BEAMSTRIPPINGPHYSICS_H