Skip to content

GitLab

Commit 06d6e6d5 authored by ext-calvo_p's avatar ext-calvo_p
Browse files

Resolve "Reviewing physics behind particle matter interaction models"

parent 49fc5b28
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 637 additions and 345 deletions
src/Algorithms/MultiBunchHandler.cpp
View file @ 06d6e6d5
......@@ -112,7 +112,7 @@ void MultiBunchHandler::saveBunch(PartBunchBase<double, 3> *beam)
Ppos_t coord, momentum;
ParticleAttrib<double> mass, charge;
ParticleAttrib<short> ptype;
ParticleAttrib<ParticleOrigin> porigin;
std::size_t localNum = beam->getLocalNum();
......@@ -128,8 +128,8 @@ void MultiBunchHandler::saveBunch(PartBunchBase<double, 3> *beam)
charge.create(localNum);
charge = beam->Q;
ptype.create(localNum);
ptype = beam->PType;
porigin.create(localNum);
porigin = beam->POrigin;
std::map<std::string, double> additionalAttributes = {
std::make_pair("REFPR", 0.0),
......@@ -234,7 +234,7 @@ bool MultiBunchHandler::readBunch(PartBunchBase<double, 3> *beam,
beam->P[localNum] = tmpBunch->P[ii];
beam->M[localNum] = tmpBunch->M[ii];
beam->Q[localNum] = tmpBunch->Q[ii];
beam->PType[localNum] = ParticleType::REGULAR;
beam->POrigin[localNum] = ParticleOrigin::REGULAR;
beam->Bin[localNum] = bunchNum;
beam->bunchNum[localNum] = bunchNum;
}
......
src/Algorithms/ParallelCyclotronTracker.cpp
View file @ 06d6e6d5
......@@ -565,7 +565,7 @@ void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) {
double temperature = elptr->getTemperature();
*gmsg << "* Temperature = " << temperature << " [K]" << endl;
std::string gas = elptr->getResidualGas();
std::string gas = elptr->getResidualGasName();
*gmsg << "* Residual gas = " << gas << endl;
bool stop = elptr->getStop();
......@@ -2047,8 +2047,7 @@ bool ParallelCyclotronTracker::applyPluginElements(const double dt) {
for(beamline_list::iterator sindex = ++(FieldDimensions.begin()); sindex != FieldDimensions.end(); ++sindex) {
if(((*sindex)->first) == ElementBase::BEAMSTRIPPING) {
BeamStripping *bstp = static_cast<BeamStripping *>(((*sindex)->second).second);
bstp->checkBeamStripping(itsBunch_m, cycl_m, turnnumber_m,
itsBunch_m->getT()*1e9 /*[ns]*/, dt);
bstp->checkBeamStripping(itsBunch_m, cycl_m);
}
}
......
src/Algorithms/ParallelTTracker.cpp
View file @ 06d6e6d5
......@@ -1341,8 +1341,8 @@ void ParallelTTracker::evenlyDistributeParticles() {
send_msgbuf.insert(send_msgbuf.end(), buffer, buffer + sizeof(double));
buffer = reinterpret_cast<const char*>(&(itsBunch_m->dt[idx]));
send_msgbuf.insert(send_msgbuf.end(), buffer, buffer + sizeof(double));
buffer = reinterpret_cast<const char*>(&(itsBunch_m->PType[idx]));
send_msgbuf.insert(send_msgbuf.end(), buffer, buffer + sizeof(short));
buffer = reinterpret_cast<const char*>(&(itsBunch_m->POrigin[idx]));
send_msgbuf.insert(send_msgbuf.end(), buffer, buffer + sizeof(ParticleOrigin));
buffer = reinterpret_cast<const char*>(&(itsBunch_m->TriID[idx]));
send_msgbuf.insert(send_msgbuf.end(), buffer, buffer + sizeof(int));
buffer = reinterpret_cast<const char*>(&(itsBunch_m->ID[idx]));
......@@ -1384,10 +1384,10 @@ void ParallelTTracker::evenlyDistributeParticles() {
j += 9 * sizeof(double);
{
const short *buffer = reinterpret_cast<const short*>(recvbuf + j);
itsBunch_m->PType[idx] = buffer[0];
const ParticleOrigin *buffer = reinterpret_cast<const ParticleOrigin*>(recvbuf + j);
itsBunch_m->POrigin[idx] = buffer[0];
}
j += sizeof(short);
j += sizeof(ParticleOrigin);
{
const int *buffer = reinterpret_cast<const int*>(recvbuf + j);
......@@ -1408,4 +1408,4 @@ void ParallelTTracker::evenlyDistributeParticles() {
if (requests.size() > 0) {
MPI_Waitall(requests.size(), &(requests[0]), MPI_STATUSES_IGNORE);
}
}
\ No newline at end of file
}
src/Classic/AbsBeamline/BeamStripping.cpp
View file @ 06d6e6d5
......@@ -3,7 +3,7 @@
// Defines the abstract interface environment for
// beam stripping physics.
//
// Copyright (c) 2018-2019, Pedro Calvo, CIEMAT, Spain
// Copyright (c) 2018-2021, Pedro Calvo, CIEMAT, Spain
// All rights reserved
//
// Implemented as part of the PhD thesis
......@@ -44,17 +44,14 @@
throw GeneralClassicException("BeamStripping::getPressureFromFile",\
"fscanf returned EOF at " #arg);
extern Inform *gmsg;
extern Inform* gmsg;
// Class BeamStripping
// ------------------------------------------------------------------------
BeamStripping::BeamStripping():
BeamStripping("")
{}
BeamStripping::BeamStripping(const BeamStripping &right):
BeamStripping::BeamStripping(const BeamStripping& right):
Component(right),
gas_m(right.gas_m),
pressure_m(right.pressure_m),
......@@ -69,9 +66,9 @@ BeamStripping::BeamStripping(const BeamStripping &right):
parmatint_m(NULL) {
}
BeamStripping::BeamStripping(const std::string &name):
BeamStripping::BeamStripping(const std::string& name):
Component(name),
gas_m(""),
gas_m(ResidualGas::NOGAS),
pressure_m(0.0),
pmapfn_m(""),
pscale_m(0.0),
......@@ -84,13 +81,13 @@ BeamStripping::BeamStripping(const std::string &name):
parmatint_m(NULL) {
}
BeamStripping::~BeamStripping() {
if (online_m)
goOffline();
}
void BeamStripping::accept(BeamlineVisitor &visitor) const {
void BeamStripping::accept(BeamlineVisitor& visitor) const {
visitor.visitBeamStripping(*this);
}
......@@ -142,16 +139,27 @@ double BeamStripping::getPScale() const {
}
void BeamStripping::setResidualGas(std::string gas) {
gas_m = gas;
if (gas == "AIR") {
gas_m = ResidualGas::AIR;
} else if (gas == "H2") {
gas_m = ResidualGas::H2;
}
}
std::string BeamStripping::getResidualGas() const {
if (gas_m == "H2" || gas_m == "AIR")
return gas_m;
else {
throw GeneralClassicException("BeamStripping::getResidualGas",
"Residual gas " + gas_m + " not found");
}
ResidualGas BeamStripping::getResidualGas() const {
return gas_m;
}
std::string BeamStripping::getResidualGasName() {
switch (gas_m) {
case ResidualGas::AIR:
return "AIR";
case ResidualGas::H2:
return "H2";
default:
throw GeneralClassicException("BeamStripping::getResidualGasName",
"Residual gas not found");
}
}
void BeamStripping::setStop(bool stopflag) {
......@@ -163,8 +171,7 @@ bool BeamStripping::getStop() const {
}
bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotron* cycl,
const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {
bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3>* bunch, Cyclotron* cycl) {
bool flagNeedUpdate = false;
......@@ -194,12 +201,12 @@ bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotro
}
void BeamStripping::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
void BeamStripping::initialise(PartBunchBase<double, 3>* bunch, double& startField, double& endField) {
endField = startField + getElementLength();
initialise(bunch, pscale_m);
}
void BeamStripping::initialise(PartBunchBase<double, 3> *bunch, const double &scaleFactor) {
void BeamStripping::initialise(PartBunchBase<double, 3>* bunch, const double& scaleFactor) {
RefPartBunch_m = bunch;
......@@ -214,7 +221,6 @@ void BeamStripping::initialise(PartBunchBase<double, 3> *bunch, const double &sc
goOnline(-1e6);
// change the mass and charge to simulate real particles
//*gmsg << "* Mass and charge have been reseted for beam stripping " << endl;
for (size_t i = 0; i < bunch->getLocalNum(); ++i) {
bunch->M[i] = bunch->getM()*1E-9;
bunch->Q[i] = bunch->getQ() * Physics::q_e;
......@@ -240,7 +246,7 @@ bool BeamStripping::bends() const {
return false;
}
void BeamStripping::getDimensions(double &/*zBegin*/, double &/*zEnd*/) const {}
void BeamStripping::getDimensions(double& /*zBegin*/, double& /*zEnd*/) const {}
ElementBase::ElementType BeamStripping::getType() const {
return BEAMSTRIPPING;
......@@ -250,17 +256,18 @@ std::string BeamStripping::getBeamStrippingShape() {
return "BeamStripping";
}
int BeamStripping::checkPoint(const double &x, const double &y, const double &z) {
int BeamStripping::checkPoint(const double& x, const double& y, const double& z) {
int cn;
double rpos = std::sqrt(x * x + y * y);
if (z >= maxz_m || z <= minz_m || rpos >= maxr_m || rpos <= minr_m)
if (z >= maxz_m || z <= minz_m || rpos >= maxr_m || rpos <= minr_m) {
cn = 0;
else
} else {
cn = 1;
}
return (cn); // 0 if out, and 1 if in
}
double BeamStripping::checkPressure(const double &x, const double &y) {
double BeamStripping::checkPressure(const double& x, const double& y) {
double pressure = 0.0;
......@@ -296,7 +303,7 @@ double BeamStripping::checkPressure(const double &x, const double &y) {
// include zero degree point
it++;
double epsilon = 0.06;
if (tet > 360 - epsilon && tet < 360 + epsilon) it = 0;
if (tet > 360 - epsilon && tet < 360 + epsilon) it = 0;
int r1t1, r2t1, r1t2, r2t2;
// r1t1 : the index of the "min angle, min radius" point in the 2D field array.
......@@ -317,18 +324,15 @@ double BeamStripping::checkPressure(const double &x, const double &y) {
*gmsg << level4 << getName() << ": Take constant value through BeamStripping::getPressure" << endl;
pressure = getPressure();
}
}
else if (ir >= PField.nrad) {
} else if (ir >= PField.nrad) {
*gmsg << level4 << getName() << ": Particle out of maximum radial position of pressure field map." << endl;
*gmsg << level4 << getName() << ": Take constant value through BeamStripping::getPressure" << endl;
pressure = getPressure();
}
else {
} else {
throw GeneralClassicException("BeamStripping::checkPressure",
"Pressure data not found");
}
}
else {
} else {
pressure = getPressure();
}
......@@ -345,9 +349,9 @@ void BeamStripping::initR(double rmin, double dr, int nrad) {
}
// Read pressure map from external file.
void BeamStripping::getPressureFromFile(const double &scaleFactor) {
void BeamStripping::getPressureFromFile(const double& scaleFactor) {
FILE *f = NULL;
FILE* f = NULL;
*gmsg << "* " << endl;
*gmsg << "* Reading pressure field map " << pmapfn_m << endl;
......@@ -356,7 +360,8 @@ void BeamStripping::getPressureFromFile(const double &scaleFactor) {
if ((f = std::fopen(pmapfn_m.c_str(), "r")) == NULL) {
throw GeneralClassicException("BeamStripping::getPressureFromFile",
"failed to open file '" + pmapfn_m + "', please check if it exists");
"failed to open file '" + pmapfn_m +
"', please check if it exists");
}
CHECK_BSTP_FSCANF_EOF(std::fscanf(f, "%lf", &PP.rmin));
......@@ -403,4 +408,4 @@ void BeamStripping::getPressureFromFile(const double &scaleFactor) {
std::fclose(f);
}
#undef CHECK_BSTP_FSCANF_EOF
\ No newline at end of file
#undef CHECK_BSTP_FSCANF_EOF
src/Classic/AbsBeamline/BeamStripping.h
View file @ 06d6e6d5
......@@ -58,53 +58,55 @@ struct PPositions {
double Pfact; // MULTIPLICATION FACTOR FOR PRESSURE MAP
};
enum class ResidualGas:short {AIR,
H2,
NOGAS};
// Class BeamStripping
// ------------------------------------------------------------------------
class BeamStripping: public Component {
public:
/// Constructor with given name.
explicit BeamStripping(const std::string &name);
explicit BeamStripping(const std::string& name);
BeamStripping();
BeamStripping(const BeamStripping &rhs);
BeamStripping(const BeamStripping& rhs);
virtual ~BeamStripping();
/// Apply visitor to BeamStripping.
virtual void accept(BeamlineVisitor &) const;
virtual void accept(BeamlineVisitor&) const;
virtual bool checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotron* cycl, const int turnnumber, const double t, const double tstep);
virtual bool checkBeamStripping(PartBunchBase<double, 3>* bunch,
Cyclotron* cycl);
virtual void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField);
virtual void initialise(PartBunchBase<double, 3>* bunch,
double& startField, double& endField);
virtual void initialise(PartBunchBase<double, 3> *bunch, const double &scaleFactor);
virtual void initialise(PartBunchBase<double, 3>* bunch,
const double& scaleFactor);
virtual void finalise();
virtual bool bends() const;
virtual void goOnline(const double &kineticEnergy);
virtual void goOnline(const double& kineticEnergy);
virtual void goOffline();
virtual ElementBase::ElementType getType() const;
virtual void getDimensions(double &zBegin, double &zEnd) const;
virtual void getDimensions(double& zBegin, double& zEnd) const;
std::string getBeamStrippingShape();
int checkPoint(const double &x, const double &y, const double &z);
int checkPoint(const double& x, const double& y, const double& z);
double checkPressure(const double &x, const double &y);
double checkPressure(const double& x, const double& y);
void setPressure(double pressure) ;
void setPressure(double pressure);
double getPressure() const;
void setTemperature(double temperature) ;
void setTemperature(double temperature);
double getTemperature() const;
void setPressureMapFN(std::string pmapfn);
......@@ -114,7 +116,8 @@ public:
virtual double getPScale() const;
void setResidualGas(std::string gas);
virtual std::string getResidualGas() const;
ResidualGas getResidualGas() const;
std::string getResidualGasName();
void setStop(bool stopflag);
virtual bool getStop() const;
......@@ -130,7 +133,7 @@ protected:
private:
///@{ parameters for BeamStripping
std::string gas_m;
ResidualGas gas_m;
double pressure_m; /// mbar
std::string pmapfn_m; /// stores the filename of the pressure map
double pscale_m; /// a scale factor for the P-field
......@@ -145,7 +148,7 @@ private:
double maxz_m; /// mm
///@}
ParticleMatterInteractionHandler *parmatint_m;
ParticleMatterInteractionHandler* parmatint_m;
protected:
// object of Matrices including pressure field map and its derivates
......@@ -155,4 +158,4 @@ protected:
PPositions PP;
};
#endif // CLASSIC_BeamStripping_HH
\ No newline at end of file
#endif // CLASSIC_BeamStripping_HH
src/Classic/AbsBeamline/CCollimator.cpp
View file @ 06d6e6d5
......@@ -82,7 +82,7 @@ bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumbe
int pflag = 0;
// now check each particle in bunch
for (unsigned int i = 0; i < tempnum; ++i) {
if (bunch->PType[i] == ParticleType::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
if (bunch->POrigin[i] == ParticleOrigin::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
// only now careful check in r
pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
/// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
......@@ -172,4 +172,4 @@ void CCollimator::doSetGeom() {
// *gmsg << "point " << i << " ( " << geom_m[i].x << ", " << geom_m[i].y << ")" << endl;
// }
// *gmsg << "rmin " << rmin_m << " rmax " << rmax_m << endl;
}
\ No newline at end of file
}
src/Classic/AbsBeamline/Degrader.cpp
View file @ 06d6e6d5
......@@ -111,7 +111,7 @@ bool Degrader::applyToReferenceParticle(const Vector_t &R,
if (!isInMaterial(R(2))) return false;
Vector_t updatedP = P;
bool isDead = getParticleMatterInteraction()->computeEnergyLoss(updatedP, RefPartBunch_m->getdT(), false);
bool isDead = getParticleMatterInteraction()->computeEnergyLoss(RefPartBunch_m, updatedP, RefPartBunch_m->getdT(), false);
double deltaP = euclidean_norm(updatedP) - euclidean_norm(P);
E(2) += deltaP * RefPartBunch_m->getM() / (RefPartBunch_m->getdT() * RefPartBunch_m->getQ() * Physics::c);
......@@ -188,4 +188,4 @@ ElementBase::ElementType Degrader::getType() const {
std::string Degrader::getDegraderShape() {
return "DEGRADER";
}
\ No newline at end of file
}
src/Classic/AbsBeamline/Stripper.cpp
View file @ 06d6e6d5
......@@ -117,7 +117,7 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
for (unsigned int i = 0; i < tempnum; ++i) {
if (bunch->PType[i] != ParticleType::REGULAR) continue;
if (bunch->POrigin[i] != ParticleOrigin::REGULAR) continue;
double tangle = calculateIncidentAngle(bunch->P[i](0), bunch->P[i](1));
changeWidth(bunch, i, tstep, tangle);
......@@ -142,12 +142,12 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i] << " collide in stripper " << getName() << endl;
// change charge and mass of PartData when the reference particle hits the stripper.
if (bunch->ID[i] == 0)
bunch->setPType(ParticleType::STRIPPED);
bunch->setPOrigin(ParticleOrigin::STRIPPED);
// change the mass and charge
bunch->M[i] = opmass_m;
bunch->Q[i] = opcharge_m * Physics::q_e;
bunch->PType[i] = ParticleType::STRIPPED;
bunch->POrigin[i] = ParticleOrigin::STRIPPED;
int j = 1;
//create new particles
......@@ -158,8 +158,8 @@ bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, co
bunch->P[index] = bunch->P[i];
bunch->Q[index] = bunch->Q[i];
bunch->M[index] = bunch->M[i];
// once the particle is stripped, change PType from 0 to 1 as a flag so as to avoid repetitive stripping.
bunch->PType[index] = ParticleType::STRIPPED;
// once the particle is stripped, change POrigin from 0 to 1 as a flag so as to avoid repetitive stripping.
bunch->POrigin[index] = ParticleOrigin::STRIPPED;
if (bunch->weHaveBins())
bunch->Bin[index] = bunch->Bin[i];
......@@ -176,7 +176,7 @@ bool Stripper::doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpd
if (!stop_m){
// change charge and mass of PartData when the reference particle hits the stripper.
if (bunch->getPType() == ParticleType::STRIPPED) {
if (bunch->getPOrigin() == ParticleOrigin::STRIPPED) {
bunch->resetM(opmass_m * 1.0e9); // GeV -> eV
bunch->resetQ(opcharge_m); // elementary charge
}
......@@ -187,4 +187,4 @@ bool Stripper::doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpd
ElementBase::ElementType Stripper::getType() const {
return STRIPPER;
}
\ No newline at end of file
}
src/Classic/Algorithms/PartBunchBase.h
View file @ 06d6e6d5
......@@ -41,13 +41,25 @@ class PartBins;
class PartBinsCyc;
class PartData;
namespace ParticleType {
enum type { REGULAR,
FIELDEMISSION,
SECONDARY,
NEWSECONDARY,
STRIPPED};
}
enum class ParticleOrigin:short {REGULAR,
SECONDARY,
STRIPPED};
enum class ParticleType:short {ELECTRON,
PROTON,
POSITRON,
ANTIPROTON,
CARBON,
HMINUS,
URANIUM,
MUON,
DEUTERON,
XENON,
H2P,
ALPHA,
HYDROGEN,
H3P,
UNNAMED};
template <class T, unsigned Dim>
class PartBunchBase
......@@ -65,19 +77,18 @@ public:
enum UnitState_t { units = 0, unitless = 1 };
public:
virtual ~PartBunchBase() { }
PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData *ref);
PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData* ref);
PartBunchBase(const PartBunchBase &rhs) = delete; // implement if needed
PartBunchBase(const PartBunchBase& rhs) = delete; // implement if needed
/*
* Bunch common member functions
*/
// This is required since we initialize the Layout and the RegionLayout with default constructor
virtual void initialize(FieldLayout_t *fLayout) = 0;
virtual void initialize(FieldLayout_t* fLayout) = 0;
bool getIfBeamEmitting();
......@@ -99,9 +110,9 @@ public:
//FIXME: unify methods, use convention that all particles have own dt
void switchOffUnitlessPositions(bool use_dt_per_particle = false);
void setDistribution(Distribution *d,
std::vector<Distribution *> addedDistributions,
size_t &np);
void setDistribution(Distribution* d,
std::vector<Distribution*> addedDistributions,
size_t& np);
bool isGridFixed() const;
......@@ -120,9 +131,9 @@ public:
bool weHaveBins() const;
void setPBins(PartBins *pbin);
void setPBins(PartBins* pbin);
void setPBins(PartBinsCyc *pbin);
void setPBins(PartBinsCyc* pbin);
/** \brief Emit particles in the given bin
i.e. copy the particles from the bin structure into the
......@@ -155,8 +166,8 @@ public:
/** \brief returns the number of particles outside of a box defined by x */
size_t calcNumPartsOutside(Vector_t x);
void calcLineDensity(unsigned int nBins, std::vector<double> &lineDensity,
std::pair<double, double> &meshInfo);
void calcLineDensity(unsigned int nBins, std::vector<double>& lineDensity,
std::pair<double, double>& meshInfo);
void setBeamFrequency(double v);
......@@ -177,7 +188,6 @@ public:
/*
Read out coordinates
*/
virtual double getPx(int i);
virtual double getPy(int i);
virtual double getPz(int i);
......@@ -194,9 +204,9 @@ public:
virtual void setZ(int i, double zcoo);
void get_bounds(Vector_t &rmin, Vector_t &rmax);
void get_bounds(Vector_t& rmin, Vector_t& rmax);
void getLocalBounds(Vector_t &rmin, Vector_t &rmax);
void getLocalBounds(Vector_t& rmin, Vector_t& rmax);
std::pair<Vector_t, double> getBoundingSphere();
......@@ -206,7 +216,6 @@ public:
/*
Compatibility function push_back
*/
void push_back(OpalParticle p);
void set_part(FVector<double, 6> z, int ii);
......@@ -219,8 +228,8 @@ public:
// The matrix [b]D[/b] is used to normalise the first two modes.
// The maximum normalised amplitudes for these modes are stored
// in [b]axmax[/b] and [b]aymax[/b].
void maximumAmplitudes(const FMatrix<double, 6, 6> &D,
double &axmax, double &aymax);
void maximumAmplitudes(const FMatrix<double, 6, 6>& D,
double& axmax, double& aymax);
void setdT(double dt);