diff --git a/src/Algorithms/MultiBunchHandler.cpp b/src/Algorithms/MultiBunchHandler.cpp
index 0d12afb95376806a15da55652221d5660b5d8176..516107b65e0d9db0c9b978ef8304e7974c22628b 100644
--- a/src/Algorithms/MultiBunchHandler.cpp
+++ b/src/Algorithms/MultiBunchHandler.cpp
@@ -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;
}
diff --git a/src/Algorithms/ParallelCyclotronTracker.cpp b/src/Algorithms/ParallelCyclotronTracker.cpp
index 983ed6cbcd6041b837f7ef88f467fcf9528d6362..229d92bcce672e8dede9ea63a0ef44273fdf98a8 100644
--- a/src/Algorithms/ParallelCyclotronTracker.cpp
+++ b/src/Algorithms/ParallelCyclotronTracker.cpp
@@ -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);
}
}
diff --git a/src/Algorithms/ParallelTTracker.cpp b/src/Algorithms/ParallelTTracker.cpp
index 3565273eaf2cfc02711feb8bb8115e8e288d8596..62affcfd093b94bc6b5f59041afe3db3935a99fc 100644
--- a/src/Algorithms/ParallelTTracker.cpp
+++ b/src/Algorithms/ParallelTTracker.cpp
@@ -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
+}
diff --git a/src/Classic/AbsBeamline/BeamStripping.cpp b/src/Classic/AbsBeamline/BeamStripping.cpp
index 90f384fe91d3ff3ed50011f9024521f74f37dd6a..b2617619ed8ab67260155be79f79238004c4a3d7 100644
--- a/src/Classic/AbsBeamline/BeamStripping.cpp
+++ b/src/Classic/AbsBeamline/BeamStripping.cpp
@@ -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
diff --git a/src/Classic/AbsBeamline/BeamStripping.h b/src/Classic/AbsBeamline/BeamStripping.h
index c96ba3ddd84b53356e35b4332214bb130d882a39..3d4325386af22eeb07bb1a21d5adf5a892718043 100644
--- a/src/Classic/AbsBeamline/BeamStripping.h
+++ b/src/Classic/AbsBeamline/BeamStripping.h
@@ -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
diff --git a/src/Classic/AbsBeamline/CCollimator.cpp b/src/Classic/AbsBeamline/CCollimator.cpp
index a1007272edf6349a809b229bb19fe447026c0031..3ce0cca889b98ad1a50d54d3e1df4a721c0192b9 100644
--- a/src/Classic/AbsBeamline/CCollimator.cpp
+++ b/src/Classic/AbsBeamline/CCollimator.cpp
@@ -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
+}
diff --git a/src/Classic/AbsBeamline/Degrader.cpp b/src/Classic/AbsBeamline/Degrader.cpp
index 13930d0bcb082581531d3ec566fd941b9d922cb8..7128202ff00b0fb90b5c17ad2c6d3f8d063640c3 100644
--- a/src/Classic/AbsBeamline/Degrader.cpp
+++ b/src/Classic/AbsBeamline/Degrader.cpp
@@ -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
+}
diff --git a/src/Classic/AbsBeamline/Stripper.cpp b/src/Classic/AbsBeamline/Stripper.cpp
index 53d326e61151ae1d0e742c715be036abf4e166e0..f0440ca4b1cf3baa839e803692b164caf6da0fd8 100644
--- a/src/Classic/AbsBeamline/Stripper.cpp
+++ b/src/Classic/AbsBeamline/Stripper.cpp
@@ -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
+}
diff --git a/src/Classic/Algorithms/PartBunchBase.h b/src/Classic/Algorithms/PartBunchBase.h
index be9af282f1a973c900efb1b834b76eb33386bd30..3dfd20502bd5fb3221c13b44b1cb5abdce043ad4 100644
--- a/src/Classic/Algorithms/PartBunchBase.h
+++ b/src/Classic/Algorithms/PartBunchBase.h
@@ -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);
double getdT() const;
@@ -240,7 +249,6 @@ public:
void set_sPos(double s);
double get_gamma() const;
-
double get_meanKineticEnergy() const;
Vector_t get_origin() const;
Vector_t get_maxExtent() const;
@@ -258,7 +266,6 @@ public:
double get_Dx() const;
double get_Dy() const;
-
double get_DDx() const;
double get_DDy() const;
@@ -270,7 +277,6 @@ public:
void get_PBounds(Vector_t &min, Vector_t &max) const;
void calcBeamParameters();
-
void calcBeamParametersInitial(); // Calculate initial beam parameters before emission.
double getCouplingConstant() const;
@@ -344,21 +350,24 @@ public:
double getM() const;
double getP() const;
double getE() const;
- ParticleType::type getPType() const;
+ ParticleOrigin getPOrigin() const;
+ ParticleType getPType() const;
+ std::string getPTypeString();
double getInitialBeta() const;
double getInitialGamma() const;
///@}
///@{ Set reference data
void resetQ(double q);
void resetM(double m);
- void setPType(ParticleType::type);
+ void setPOrigin(ParticleOrigin);
+ void setPType(std::string type);
///@}
double getdE() const;
virtual double getGamma(int i);
virtual double getBeta(int i);
virtual void actT();
- const PartData *getReference() const;
+ const PartData* getReference() const;
double getEmissionDeltaT();
@@ -373,7 +382,7 @@ public:
void calcEMean();
- Inform &print(Inform &os);
+ Inform& print(Inform& os);
/*
* (Pure) virtual member functions
@@ -385,14 +394,14 @@ public:
virtual void resetInterpolationCache(bool clearCache = false);
- /** \brief calculates back the max/min of the efield on the grid */
+ //brief calculates back the max/min of the efield on the grid
virtual VectorPair_t getEExtrema() = 0;
virtual double getRho(int x, int y, int z) = 0;
virtual void computeSelfFields() = 0;
- /** /brief used for self fields with binned distribution */
+ //brief used for self fields with binned distribution
virtual void computeSelfFields(int bin) = 0;
virtual void computeSelfFields_cycl(double gamma) = 0;
@@ -414,7 +423,7 @@ public:
// virtual Mesh_t &getMesh() = 0;
// virtual void setFieldLayout(FieldLayout_t* fLayout) = 0;
- virtual FieldLayout_t &getFieldLayout() = 0;
+ virtual FieldLayout_t& getFieldLayout() = 0;
virtual void resizeMesh() { };
@@ -433,7 +442,7 @@ public:
unsigned int getMinimumNumberOfParticlesPerCore() const;
void setMinimumNumberOfParticlesPerCore(unsigned int n);
- ParticleLayout<T, Dim> & getLayout();
+ ParticleLayout<T, Dim>& getLayout();
const ParticleLayout<T, Dim>& getLayout() const;
bool getUpdateFlag(UpdateFlags_t f) const;
@@ -485,38 +494,36 @@ public:
/*
* Bunch attributes
*/
-
-
ParticlePos_t& R;
ParticleIndex_t& ID;
-
// Particle container attributes
- ParticleAttrib< Vector_t > P; // particle momentum // ParticleSpatialLayout<double, 3>::ParticlePos_t P;
- ParticleAttrib< double > Q; // charge per simulation particle, unit: C.
- ParticleAttrib< double > M; // mass per simulation particle, for multi-species particle tracking, unit:GeV/c^2.
- ParticleAttrib< double > Phi; // the electric potential
- ParticleAttrib< Vector_t > Ef; // e field vector
- ParticleAttrib< Vector_t > Eftmp; // e field vector for gun simulations
-
- ParticleAttrib< Vector_t > Bf; // b field vector
- ParticleAttrib< int > Bin; // holds the bin in which the particle is in, if zero particle is marked for deletion
- ParticleAttrib< double > dt; // holds the dt timestep for particle
-
- ParticleAttrib< short > PType; // we can distinguish dark current particles from primary particle
- ParticleAttrib< int > TriID; // holds the ID of triangle that the particle hit. Only for BoundaryGeometry case.
- ParticleAttrib< short > cavityGapCrossed; ///< particle just crossed cavity gap (for ParallelCyclotronTracker)
-
- ParticleAttrib< short > bunchNum; // bunch number to which particle belongs (multi-bunch mode)
-
+ ParticleAttrib< Vector_t > P; // particle momentum // ParticleSpatialLayout<double, 3>::ParticlePos_t P;
+ ParticleAttrib< double > Q; // charge per simulation particle, unit: C.
+ ParticleAttrib< double > M; // mass per simulation particle, for multi-species particle tracking, unit:GeV/c^2.
+ ParticleAttrib< double > Phi; // the electric potential
+ ParticleAttrib< Vector_t > Ef; // e field vector
+ ParticleAttrib< Vector_t > Eftmp; // e field vector for gun simulations
+
+ ParticleAttrib< Vector_t > Bf; // b field vector
+ ParticleAttrib< int > Bin; // holds the bin in which the particle is in, if zero particle is marked for deletion
+ ParticleAttrib< double > dt; // holds the dt timestep for particle
+ ParticleAttrib< ParticleType > PType; // particle names
+ ParticleAttrib< ParticleOrigin > POrigin; // we can distinguish dark current particles from primary particle
+ ParticleAttrib< int > TriID; // holds the ID of triangle that the particle hit. Only for BoundaryGeometry case.
+ ParticleAttrib< short > cavityGapCrossed; // particle just crossed cavity gap (for ParallelCyclotronTracker)
+ ParticleAttrib< short > bunchNum; // bunch number to which particle belongs (multi-bunch mode)
Vector_t RefPartR_m;
Vector_t RefPartP_m;
- ParticleType::type refPType_m;
+
CoordinateSystemTrafo toLabTrafo_m;
+ ParticleOrigin refPOrigin_m;
+ ParticleType refPType_m;
+
// The structure for particle binning
- PartBins *pbin_m;
+ PartBins* pbin_m;
/// timer for IC, can not be in Distribution.h
IpplTimings::TimerRef distrReload_m;
@@ -548,13 +555,11 @@ protected:
/// timer for selfField calculation
IpplTimings::TimerRef selfFieldTimer_m;
+ const PartData* reference;
- const PartData *reference;
-
-
-// /*
-// Member variables starts here
-// */
+ /*
+ Member variables starts here
+ */
// unit state of PartBunch
UnitState_t unit_state_;
@@ -601,7 +606,6 @@ protected:
/// rms emittance (not normalized)
Vector_t eps_m;
-
/// rms normalized emittance
Vector_t eps_norm_m;
@@ -613,7 +617,6 @@ protected:
/// dispersion x & y
double Dx_m;
double Dy_m;
-
/// derivative of the dispersion
double DDx_m;
double DDy_m;
@@ -624,7 +627,7 @@ protected:
Vektor<int, 3> nr_m;
/// stores the used field solver
- FieldSolver *fs_m;
+ FieldSolver* fs_m;
double couplingConstant_m;
@@ -674,8 +677,7 @@ protected:
std::unique_ptr<size_t[]> globalPartPerNode_m;
-
- Distribution *dist_m;
+ Distribution* dist_m;
// flag to tell if we are a DC-beam
bool dcBeam_m;
@@ -685,4 +687,4 @@ protected:
#include "PartBunchBase.hpp"
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Algorithms/PartBunchBase.hpp b/src/Classic/Algorithms/PartBunchBase.hpp
index 4559a3d35181d8189f8ec3cc0d22c0710ccd3a5b..44a493be1b8228c378ca5f2dfe529b2e1aa3261d 100644
--- a/src/Classic/Algorithms/PartBunchBase.hpp
+++ b/src/Classic/Algorithms/PartBunchBase.hpp
@@ -2,7 +2,7 @@
// Class PartBunchBase
// Base class for representing particle bunches.
//
-// Copyright (c) 2008 - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
+// Copyright (c) 2008 - 2021, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
@@ -23,7 +23,7 @@
#include "Distribution/Distribution.h"
-#include "AbstractObjects/OpalData.h" // OPAL file
+#include "AbstractObjects/OpalData.h"
#include "Algorithms/PartBins.h"
#include "Algorithms/PartBinsCyc.h"
#include "Algorithms/PartData.h"
@@ -35,10 +35,10 @@
#include "Utilities/SwitcherError.h"
#include "Utilities/Util.h"
-extern Inform *gmsg;
+extern Inform* gmsg;
template <class T, unsigned Dim>
-PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData *ref)
+PartBunchBase<T, Dim>::PartBunchBase(AbstractParticle<T, Dim>* pb, const PartData* ref)
: R(*(pb->R_p)),
ID(*(pb->ID_p)),
pbin_m(nullptr),
@@ -159,7 +159,7 @@ template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setEnergyBins(int numberOfEnergyBins) {
bingamma_m = std::unique_ptr<double[]>(new double[numberOfEnergyBins]);
binemitted_m = std::unique_ptr<size_t[]>(new size_t[numberOfEnergyBins]);
- for(int i = 0; i < numberOfEnergyBins; i++)
+ for (int i = 0; i < numberOfEnergyBins; i++)
binemitted_m[i] = 0;
}
@@ -177,16 +177,16 @@ bool PartBunchBase<T, Dim>::weHaveEnergyBins() {
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::switchToUnitlessPositions(bool use_dt_per_particle) {
- if(unit_state_ == unitless)
+ if (unit_state_ == unitless)
throw SwitcherError("PartBunch::switchToUnitlessPositions",
"Cannot make a unitless PartBunch unitless");
bool hasToReset = false;
- if(!R.isDirty()) hasToReset = true;
+ if (!R.isDirty()) hasToReset = true;
- for(size_t i = 0; i < getLocalNum(); i++) {
+ for (size_t i = 0; i < getLocalNum(); i++) {
double dt = getdT();
- if(use_dt_per_particle)
+ if (use_dt_per_particle)
dt = this->dt[i];
R[i] /= Vector_t(Physics::c * dt);
@@ -194,24 +194,23 @@ void PartBunchBase<T, Dim>::switchToUnitlessPositions(bool use_dt_per_particle)
unit_state_ = unitless;
- if(hasToReset) R.resetDirtyFlag();
+ if (hasToReset) R.resetDirtyFlag();
}
-
//FIXME: unify methods, use convention that all particles have own dt
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::switchOffUnitlessPositions(bool use_dt_per_particle) {
- if(unit_state_ == units)
+ if (unit_state_ == units)
throw SwitcherError("PartBunch::switchOffUnitlessPositions",
"Cannot apply units twice to PartBunch");
bool hasToReset = false;
- if(!R.isDirty()) hasToReset = true;
+ if (!R.isDirty()) hasToReset = true;
- for(size_t i = 0; i < getLocalNum(); i++) {
+ for (size_t i = 0; i < getLocalNum(); i++) {
double dt = getdT();
- if(use_dt_per_particle)
+ if (use_dt_per_particle)
dt = this->dt[i];
R[i] *= Vector_t(Physics::c * dt);
@@ -219,7 +218,7 @@ void PartBunchBase<T, Dim>::switchOffUnitlessPositions(bool use_dt_per_particle)
unit_state_ = units;
- if(hasToReset) R.resetDirtyFlag();
+ if (hasToReset) R.resetDirtyFlag();
}
@@ -230,13 +229,11 @@ void PartBunchBase<T, Dim>::do_binaryRepart() {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::setDistribution(Distribution *d,
- std::vector<Distribution *> addedDistributions,
- size_t &np)
-{
+void PartBunchBase<T, Dim>::setDistribution(Distribution* d,
+ std::vector<Distribution*> addedDistributions,
+ size_t& np) {
Inform m("setDistribution " );
dist_m = d;
-
dist_m->createOpalT(this, addedDistributions, np);
// if (Options::cZero)
@@ -281,8 +278,7 @@ bool PartBunchBase<T, Dim>::doEmission() {
template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::weHaveBins() const {
-
- if(pbin_m != NULL)
+ if (pbin_m != NULL)
return pbin_m->weHaveBins();
else
return false;
@@ -290,7 +286,7 @@ bool PartBunchBase<T, Dim>::weHaveBins() const {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::setPBins(PartBins *pbin) {
+void PartBunchBase<T, Dim>::setPBins(PartBins* pbin) {
pbin_m = pbin;
*gmsg << *pbin_m << endl;
setEnergyBins(pbin_m->getNBins());
@@ -298,7 +294,7 @@ void PartBunchBase<T, Dim>::setPBins(PartBins *pbin) {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::setPBins(PartBinsCyc *pbin) {
+void PartBunchBase<T, Dim>::setPBins(PartBinsCyc* pbin) {
pbin_m = pbin;
setEnergyBins(pbin_m->getNBins());
}
@@ -329,7 +325,7 @@ void PartBunchBase<T, Dim>::rebin() {
template <class T, unsigned Dim>
int PartBunchBase<T, Dim>::getLastemittedBin() {
- if(pbin_m != NULL)
+ if (pbin_m != NULL)
return pbin_m->getLastemittedBin();
else
return 0;
@@ -338,7 +334,7 @@ int PartBunchBase<T, Dim>::getLastemittedBin() {
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setLocalBinCount(size_t num, int bin) {
- if(pbin_m != NULL) {
+ if (pbin_m != NULL) {
pbin_m->setPartNum(bin, num);
}
}
@@ -351,18 +347,18 @@ void PartBunchBase<T, Dim>::calcGammas() {
size_t s = 0;
- for(int i = 0; i < emittedBins; i++)
+ for (int i = 0; i < emittedBins; i++)
bingamma_m[i] = 0.0;
- for(unsigned int n = 0; n < getLocalNum(); n++)
+ for (unsigned int n = 0; n < getLocalNum(); n++)
bingamma_m[this->Bin[n]] += std::sqrt(1.0 + dot(this->P[n], this->P[n]));
std::unique_ptr<size_t[]> particlesInBin(new size_t[emittedBins]);
reduce(bingamma_m.get(), bingamma_m.get() + emittedBins, bingamma_m.get(), OpAddAssign());
reduce(binemitted_m.get(), binemitted_m.get() + emittedBins, particlesInBin.get(), OpAddAssign());
- for(int i = 0; i < emittedBins; i++) {
+ for (int i = 0; i < emittedBins; i++) {
size_t &pInBin = particlesInBin[i];
- if(pInBin != 0) {
+ if (pInBin != 0) {
bingamma_m[i] /= pInBin;
INFOMSG(level2 << "Bin " << std::setw(3) << i
<< " gamma = " << std::setw(8) << std::scientific
@@ -378,12 +374,12 @@ void PartBunchBase<T, Dim>::calcGammas() {
particlesInBin.reset();
- if(s != getTotalNum() && !OpalData::getInstance()->hasGlobalGeometry())
+ if (s != getTotalNum() && !OpalData::getInstance()->hasGlobalGeometry())
ERRORMSG("sum(Bins)= " << s << " != sum(R)= " << getTotalNum() << endl;);
- if(emittedBins >= 2) {
- for(int i = 1; i < emittedBins; i++) {
- if(binemitted_m[i - 1] != 0 && binemitted_m[i] != 0)
+ if (emittedBins >= 2) {
+ for (int i = 1; i < emittedBins; i++) {
+ if (binemitted_m[i - 1] != 0 && binemitted_m[i] != 0)
INFOMSG(level2 << "d(gamma)= " << 100.0 * std::abs(bingamma_m[i - 1] - bingamma_m[i]) / bingamma_m[i] << " [%] "
<< "between bin " << i - 1 << " and " << i << endl);
}
@@ -396,9 +392,9 @@ void PartBunchBase<T, Dim>::calcGammas_cycl() {
const int emittedBins = pbin_m->getLastemittedBin();
- for(int i = 0; i < emittedBins; i++)
+ for (int i = 0; i < emittedBins; i++)
bingamma_m[i] = 0.0;
- for(unsigned int n = 0; n < getLocalNum(); n++) {
+ for (unsigned int n = 0; n < getLocalNum(); n++) {
if ( this->Bin[n] > -1 ) {
bingamma_m[this->Bin[n]] += std::sqrt(1.0 + dot(this->P[n], this->P[n]));
}
@@ -406,11 +402,12 @@ void PartBunchBase<T, Dim>::calcGammas_cycl() {
allreduce(*bingamma_m.get(), emittedBins, std::plus<double>());
- for(int i = 0; i < emittedBins; i++) {
- if(pbin_m->getTotalNumPerBin(i) > 0)
+ for (int i = 0; i < emittedBins; i++) {
+ if (pbin_m->getTotalNumPerBin(i) > 0) {
bingamma_m[i] /= pbin_m->getTotalNumPerBin(i);
- else
+ } else {
bingamma_m[i] = 0.0;
+ }
INFOMSG("Bin " << i << " : particle number = " << pbin_m->getTotalNumPerBin(i)
<< " gamma = " << bingamma_m[i] << endl);
}
@@ -442,7 +439,7 @@ size_t PartBunchBase<T, Dim>::calcNumPartsOutside(Vector_t x) {
std::size_t localnum = 0;
const Vector_t meanR = get_rmean();
- for(unsigned long k = 0; k < getLocalNum(); ++ k)
+ for (unsigned long k = 0; k < getLocalNum(); ++ k)
if (std::abs(R[k](0) - meanR(0)) > x(0) ||
std::abs(R[k](1) - meanR(1)) > x(1) ||
std::abs(R[k](2) - meanR(2)) > x(2)) {
@@ -471,9 +468,8 @@ size_t PartBunchBase<T, Dim>::calcNumPartsOutside(Vector_t x) {
*/
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::calcLineDensity(unsigned int nBins,
- std::vector<double> &lineDensity,
- std::pair<double, double> &meshInfo)
-{
+ std::vector<double>& lineDensity,
+ std::pair<double, double>& meshInfo) {
Vector_t rmin, rmax;
get_bounds(rmin, rmax);
@@ -509,7 +505,6 @@ void PartBunchBase<T, Dim>::calcLineDensity(unsigned int nBins,
}
-
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::boundp() {
/*
@@ -517,12 +512,12 @@ void PartBunchBase<T, Dim>::boundp() {
*/
IpplTimings::startTimer(boundpTimer_m);
- //if(!R.isDirty() && stateOfLastBoundP_ == unit_state_) return;
+ //if (!R.isDirty() && stateOfLastBoundP_ == unit_state_) return;
if ( !(R.isDirty() || ID.isDirty() ) && stateOfLastBoundP_ == unit_state_) return; //-DW
stateOfLastBoundP_ = unit_state_;
- if(!isGridFixed()) {
+ if (!isGridFixed()) {
const int dimIdx = (dcBeam_m? 2: 3);
/**
@@ -538,7 +533,7 @@ void PartBunchBase<T, Dim>::boundp() {
Vector_t len = rmax_m - rmin_m;
double volume = 1.0;
- for(int i = 0; i < dimIdx; i++) {
+ for (int i = 0; i < dimIdx; i++) {
double length = std::abs(rmax_m[i] - rmin_m[i]);
if (length < 1e-10) {
rmax_m[i] += 1e-10;
@@ -579,7 +574,7 @@ void PartBunchBase<T, Dim>::boundp() {
}
//INFOMSG("It is a full boundp hz= " << hr_m << " rmax= " << rmax_m << " rmin= " << rmin_m << endl);
- if(hr_m[0] * hr_m[1] * hr_m[2] <= 0) {
+ if (hr_m[0] * hr_m[1] * hr_m[2] <= 0) {
throw GeneralClassicException("boundp() ", "h<0, can not build a mesh");
}
@@ -605,15 +600,14 @@ void PartBunchBase<T, Dim>::boundp_destroy() {
IpplTimings::startTimer(boundpTimer_m);
std::unique_ptr<size_t[]> countLost;
- if(weHaveBins()) {
+ if (weHaveBins()) {
const int tempN = pbin_m->getLastemittedBin();
countLost = std::unique_ptr<size_t[]>(new size_t[tempN]);
- for(int ii = 0; ii < tempN; ii++) countLost[ii] = 0;
+ for (int ii = 0; ii < tempN; ii++) countLost[ii] = 0;
}
this->updateDomainLength(nr_m);
-
IpplTimings::startTimer(boundpBoundsTimer_m);
get_bounds(rmin_m, rmax_m);
IpplTimings::stopTimer(boundpBoundsTimer_m);
@@ -626,20 +620,20 @@ void PartBunchBase<T, Dim>::boundp_destroy() {
if (checkfactor != 0) {
//INFOMSG("checkfactor= " << checkfactor << endl);
// check the bunch if its full size is larger than checkfactor times of its rms size
- if(checkfactor < 0) {
+ if (checkfactor < 0) {
checkfactor *= -1;
if (len[0] > checkfactor * rrms_m[0] ||
len[1] > checkfactor * rrms_m[1] ||
- len[2] > checkfactor * rrms_m[2])
- {
- for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
+ len[2] > checkfactor * rrms_m[2] ) {
+
+ for (unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
/* delete the particle if the distance to the beam center
* is larger than 8 times of beam's rms size
*/
if (std::abs(R[ii](0) - rmean_m(0)) > checkfactor * rrms_m[0] ||
std::abs(R[ii](1) - rmean_m(1)) > checkfactor * rrms_m[1] ||
- std::abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2])
- {
+ std::abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2] ) {
+
// put particle onto deletion list
destroy(1, ii);
//update bin parameter
@@ -651,18 +645,17 @@ void PartBunchBase<T, Dim>::boundp_destroy() {
}
}
}
- }
- else {
+ } else {
if (len[0] > checkfactor * rrms_m[0] ||
- len[2] > checkfactor * rrms_m[2])
- {
- for(unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
+ len[2] > checkfactor * rrms_m[2] ) {
+
+ for (unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
/* delete the particle if the distance to the beam center
* is larger than 8 times of beam's rms size
*/
if (std::abs(R[ii](0) - rmean_m(0)) > checkfactor * rrms_m[0] ||
- std::abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2])
- {
+ std::abs(R[ii](2) - rmean_m(2)) > checkfactor * rrms_m[2] ) {
+
// put particle onto deletion list
destroy(1, ii);
//update bin parameter
@@ -677,7 +670,7 @@ void PartBunchBase<T, Dim>::boundp_destroy() {
}
}
- for(int i = 0; i < dimIdx; i++) {
+ for (int i = 0; i < dimIdx; i++) {
double length = std::abs(rmax_m[i] - rmin_m[i]);
rmax_m[i] += dh_m * length;
rmin_m[i] -= dh_m * length;
@@ -687,7 +680,7 @@ void PartBunchBase<T, Dim>::boundp_destroy() {
// rescale mesh
this->updateFields(hr_m, rmin_m);
- if(weHaveBins()) {
+ if (weHaveBins()) {
pbin_m->updatePartInBin_cyc(countLost.get());
}
@@ -718,21 +711,22 @@ size_t PartBunchBase<T, Dim>::boundp_destroyT() {
size_t ne = 0;
const size_t localNum = getLocalNum();
- if(weHaveEnergyBins()) {
+ if (weHaveEnergyBins()) {
tmpbinemitted = std::unique_ptr<size_t[]>(new size_t[getNumberOfEnergyBins()]);
- for(int i = 0; i < getNumberOfEnergyBins(); i++) {
+ for (int i = 0; i < getNumberOfEnergyBins(); i++) {
tmpbinemitted[i] = 0;
}
- for(unsigned int i = 0; i < localNum; i++) {
+ for (unsigned int i = 0; i < localNum; i++) {
if (Bin[i] < 0) {
ne++;
destroy(1, i);
- } else
+ } else {
tmpbinemitted[Bin[i]]++;
+ }
}
} else {
- for(unsigned int i = 0; i < localNum; i++) {
- if((Bin[i] < 0) && ((localNum - ne) > minNumParticlesPerCore)) { // need in minimum x particles per node
+ for (unsigned int i = 0; i < localNum; i++) {
+ if ((Bin[i] < 0) && ((localNum - ne) > minNumParticlesPerCore)) { // need in minimum x particles per node
ne++;
destroy(1, i);
}
@@ -751,9 +745,9 @@ size_t PartBunchBase<T, Dim>::boundp_destroyT() {
calcBeamParameters();
gatherLoadBalanceStatistics();
- if(weHaveEnergyBins()) {
+ if (weHaveEnergyBins()) {
const int lastBin = dist_m->getLastEmittedEnergyBin();
- for(int i = 0; i < lastBin; i++) {
+ for (int i = 0; i < lastBin; i++) {
binemitted_m[i] = tmpbinemitted[i];
}
}
@@ -772,12 +766,12 @@ size_t PartBunchBase<T, Dim>::destroyT() {
const size_t totalNum = getTotalNum();
size_t ne = 0;
- if(weHaveEnergyBins()) {
+ if (weHaveEnergyBins()) {
tmpbinemitted = std::unique_ptr<size_t[]>(new size_t[getNumberOfEnergyBins()]);
- for(int i = 0; i < getNumberOfEnergyBins(); i++) {
+ for (int i = 0; i < getNumberOfEnergyBins(); i++) {
tmpbinemitted[i] = 0;
}
- for(unsigned int i = 0; i < localNum; i++) {
+ for (unsigned int i = 0; i < localNum; i++) {
if (Bin[i] < 0) {
destroy(1, i);
++ ne;
@@ -786,8 +780,8 @@ size_t PartBunchBase<T, Dim>::destroyT() {
}
} else {
Inform dmsg("destroy: ", INFORM_ALL_NODES);
- for(size_t i = 0; i < localNum; i++) {
- if((Bin[i] < 0)) {
+ for (size_t i = 0; i < localNum; i++) {
+ if ((Bin[i] < 0)) {
if ((localNum - ne) > minNumParticlesPerCore) { // need in minimum x particles per node
ne++;
destroy(1, i);
@@ -807,7 +801,7 @@ size_t PartBunchBase<T, Dim>::destroyT() {
if (weHaveEnergyBins()) {
const int lastBin = dist_m->getLastEmittedEnergyBin();
- for(int i = 0; i < lastBin; i++) {
+ for (int i = 0; i < lastBin; i++) {
binemitted_m[i] = tmpbinemitted[i];
}
}
@@ -819,64 +813,56 @@ size_t PartBunchBase<T, Dim>::destroyT() {
return totalNum - newTotalNum;
}
+
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getPx(int /*i*/) {
return 0.0;
}
-
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getPy(int) {
return 0.0;
}
-
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getPz(int) {
return 0.0;
}
-
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getPx0(int) {
return 0.0;
}
-
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getPy0(int) {
return 0;
}
-
//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getX(int i) {
return this->R[i](0);
}
-
//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getY(int i) {
return this->R[i](1);
}
-
//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getZ(int i) {
return this->R[i](2);
}
-
//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getX0(int /*i*/) {
return 0.0;
}
-
//ff
template <class T, unsigned Dim>
double PartBunchBase<T, Dim>::getY0(int /*i*/) {
@@ -890,7 +876,7 @@ void PartBunchBase<T, Dim>::setZ(int /*i*/, double /*zcoo*/) {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::get_bounds(Vector_t &rmin, Vector_t &rmax) {
+void PartBunchBase<T, Dim>::get_bounds(Vector_t& rmin, Vector_t& rmax) {
this->getLocalBounds(rmin, rmax);
@@ -911,7 +897,7 @@ void PartBunchBase<T, Dim>::get_bounds(Vector_t &rmin, Vector_t &rmax) {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::getLocalBounds(Vector_t &rmin, Vector_t &rmax) {
+void PartBunchBase<T, Dim>::getLocalBounds(Vector_t& rmin, Vector_t& rmax) {
const size_t localNum = getLocalNum();
if (localNum == 0) {
double maxValue = 1e8;
@@ -1013,12 +999,12 @@ OpalParticle PartBunchBase<T, Dim>::get_part(int ii) {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::maximumAmplitudes(const FMatrix<double, 6, 6> &D,
- double &axmax, double &aymax) {
+void PartBunchBase<T, Dim>::maximumAmplitudes(const FMatrix<double, 6, 6>& D,
+ double& axmax, double& aymax) {
axmax = aymax = 0.0;
OpalParticle part;
- for(unsigned int ii = 0; ii < getLocalNum(); ii++) {
+ for (unsigned int ii = 0; ii < getLocalNum(); ii++) {
part = get_part(ii);
@@ -1210,7 +1196,7 @@ void PartBunchBase<T, Dim>::set_meshEnlargement(double dh) {
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::gatherLoadBalanceStatistics() {
- for(int i = 0; i < Ippl::getNodes(); i++)
+ for (int i = 0; i < Ippl::getNodes(); i++)
globalPartPerNode_m[i] = 0;
std::size_t localnum = getLocalNum();
@@ -1243,8 +1229,8 @@ void PartBunchBase<T, Dim>::calcBeamParameters() {
get_bounds(rmin_m, rmax_m);
- if(totalNum == 0) {
- for(unsigned int i = 0 ; i < Dim; i++) {
+ if (totalNum == 0) {
+ for (unsigned int i = 0 ; i < Dim; i++) {
rmean_m(i) = 0.0;
pmean_m(i) = 0.0;
rrms_m(i) = 0.0;
@@ -1261,19 +1247,19 @@ void PartBunchBase<T, Dim>::calcBeamParameters() {
const size_t intN = calcMoments();
const double N = static_cast<double>(intN);
- for(unsigned int i = 0 ; i < Dim; i++) {
+ for (unsigned int i = 0 ; i < Dim; i++) {
rmean_m(i) = centroid_m[2 * i] / N;
pmean_m(i) = centroid_m[(2 * i) + 1] / N;
rsqsum(i) = moments_m(2 * i, 2 * i) - N * rmean_m(i) * rmean_m(i);
psqsum(i) = moments_m((2 * i) + 1, (2 * i) + 1) - N * pmean_m(i) * pmean_m(i);
- if(psqsum(i) < 0)
+ if (psqsum(i) < 0)
psqsum(i) = 0;
rpsum(i) = moments_m((2 * i), (2 * i) + 1) - N * rmean_m(i) * pmean_m(i);
}
eps2 = (rsqsum * psqsum - rpsum * rpsum) / (N * N);
rpsum /= N;
- for(unsigned int i = 0 ; i < Dim; i++) {
+ for (unsigned int i = 0 ; i < Dim; i++) {
rrms_m(i) = std::sqrt(rsqsum(i) / N);
prms_m(i) = std::sqrt(psqsum(i) / N);
eps_norm_m(i) = std::sqrt(std::max(eps2(i), zero));
@@ -1291,7 +1277,7 @@ void PartBunchBase<T, Dim>::calcBeamParameters() {
// Find unnormalized emittance.
double gamma = 0.0;
- for(size_t i = 0; i < locNp; i++)
+ for (size_t i = 0; i < locNp; i++)
gamma += std::sqrt(1.0 + dot(P[i], P[i]));
allreduce(&gamma, 1, std::plus<double>());
@@ -1322,14 +1308,14 @@ void PartBunchBase<T, Dim>::calcBeamParametersInitial() {
const double N = static_cast<double>(getTotalNum());
- if(N == 0) {
+ if (N == 0) {
rmean_m = Vector_t(0.0);
pmean_m = Vector_t(0.0);
rrms_m = Vector_t(0.0);
prms_m = Vector_t(0.0);
eps_m = Vector_t(0.0);
} else {
- if(Ippl::myNode() == 0) {
+ if (Ippl::myNode() == 0) {
// fixme: the following code is crap!
// Only use one node as this function will get called only once before
// particles have been emitted (at least in principle).
@@ -1339,19 +1325,19 @@ void PartBunchBase<T, Dim>::calcBeamParametersInitial() {
const double N = static_cast<double>(pbin_m->getNp());
calcMomentsInitial();
- for(unsigned int i = 0 ; i < Dim; i++) {
+ for (unsigned int i = 0 ; i < Dim; i++) {
rmean_m(i) = centroid_m[2 * i] / N;
pmean_m(i) = centroid_m[(2 * i) + 1] / N;
rsqsum(i) = moments_m(2 * i, 2 * i) - N * rmean_m(i) * rmean_m(i);
psqsum(i) = moments_m((2 * i) + 1, (2 * i) + 1) - N * pmean_m(i) * pmean_m(i);
- if(psqsum(i) < 0)
+ if (psqsum(i) < 0)
psqsum(i) = 0;
rpsum(i) = moments_m((2 * i), (2 * i) + 1) - N * rmean_m(i) * pmean_m(i);
}
eps2 = (rsqsum * psqsum - rpsum * rpsum) / (N * N);
rpsum /= N;
- for(unsigned int i = 0 ; i < Dim; i++) {
+ for (unsigned int i = 0 ; i < Dim; i++) {
rrms_m(i) = std::sqrt(rsqsum(i) / N);
prms_m(i) = std::sqrt(psqsum(i) / N);
@@ -1380,7 +1366,7 @@ void PartBunchBase<T, Dim>::setCouplingConstant(double c) {
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setCharge(double q) {
qi_m = q;
- if(getTotalNum() != 0)
+ if (getTotalNum() != 0)
Q = qi_m;
else
WARNMSG("Could not set total charge in PartBunch::setCharge based on getTotalNum" << endl);
@@ -1420,7 +1406,7 @@ void PartBunchBase<T, Dim>::setSolver(FieldSolver *fs) {
CAN not re-inizialize ParticleLayout
this is an IPPL issue
*/
- if(!OpalData::getInstance()->hasBunchAllocated()) {
+ if (!OpalData::getInstance()->hasBunchAllocated()) {
this->initialize(fs_m->getFieldLayout());
// this->setMesh(fs_m->getMesh());
// this->setFieldLayout(fs_m->getFieldLayout());
@@ -1430,7 +1416,7 @@ void PartBunchBase<T, Dim>::setSolver(FieldSolver *fs) {
template <class T, unsigned Dim>
bool PartBunchBase<T, Dim>::hasFieldSolver() {
- if(fs_m)
+ if (fs_m)
return fs_m->hasValidSolver();
else
return false;
@@ -1440,7 +1426,7 @@ bool PartBunchBase<T, Dim>::hasFieldSolver() {
/// \brief Return the fieldsolver type if we have a fieldsolver
template <class T, unsigned Dim>
std::string PartBunchBase<T, Dim>::getFieldSolverType() const {
- if(fs_m)
+ if (fs_m)
return fs_m->getFieldSolverType();
else
return "";
@@ -1603,21 +1589,20 @@ double PartBunchBase<T, Dim>::calcMeanPhi() {
double py[emittedBins];
double meanPhi = 0.0;
- for(int ii = 0; ii < emittedBins; ii++) {
+ for (int ii = 0; ii < emittedBins; ii++) {
phi[ii] = 0.0;
px[ii] = 0.0;
py[ii] = 0.0;
}
- for(unsigned int ii = 0; ii < getLocalNum(); ii++) {
-
+ for (unsigned int ii = 0; ii < getLocalNum(); ii++) {
px[Bin[ii]] += P[ii](0);
py[Bin[ii]] += P[ii](1);
}
reduce(px, px + emittedBins, px, OpAddAssign());
reduce(py, py + emittedBins, py, OpAddAssign());
- for(int ii = 0; ii < emittedBins; ii++) {
+ for (int ii = 0; ii < emittedBins; ii++) {
phi[ii] = calculateAngle(px[ii], py[ii]);
meanPhi += phi[ii];
INFOMSG("Bin " << ii << " mean phi = " << phi[ii] * 180.0 / Physics::pi - 90.0 << "[degree]" << endl);
@@ -1642,31 +1627,31 @@ bool PartBunchBase<T, Dim>::resetPartBinID2(const double eta) {
calcGammas_cycl();
int maxbin = pbin_m->getNBins();
size_t partInBin[maxbin];
- for(int ii = 0; ii < maxbin; ii++) partInBin[ii] = 0;
+ for (int ii = 0; ii < maxbin; ii++) partInBin[ii] = 0;
double pMin0 = 1.0e9;
double pMin = 0.0;
double maxbinIndex = 0;
- for(unsigned long int n = 0; n < getLocalNum(); n++) {
+ for (unsigned long int n = 0; n < getLocalNum(); n++) {
double temp_betagamma = std::sqrt(std::pow(P[n](0), 2) + std::pow(P[n](1), 2));
- if(pMin0 > temp_betagamma)
+ if (pMin0 > temp_betagamma)
pMin0 = temp_betagamma;
}
reduce(pMin0, pMin, OpMinAssign());
INFOMSG("minimal beta*gamma = " << pMin << endl);
double asinh0 = std::asinh(pMin);
- for(unsigned long int n = 0; n < getLocalNum(); n++) {
+ for (unsigned long int n = 0; n < getLocalNum(); n++) {
double temp_betagamma = std::sqrt(std::pow(P[n](0), 2) + std::pow(P[n](1), 2));
int itsBinID = std::floor((std::asinh(temp_betagamma) - asinh0) / eta + 1.0E-6);
Bin[n] = itsBinID;
- if(maxbinIndex < itsBinID) {
+ if (maxbinIndex < itsBinID) {
maxbinIndex = itsBinID;
}
- if(itsBinID >= maxbin) {
+ if (itsBinID >= maxbin) {
ERRORMSG("The bin number limit is " << maxbin << ", please increase the energy interval and try again" << endl);
return false;
} else
@@ -1682,7 +1667,6 @@ bool PartBunchBase<T, Dim>::resetPartBinID2(const double eta) {
calcGammas_cycl();
return true;
-
}
@@ -1733,15 +1717,97 @@ double PartBunchBase<T, Dim>::getE() const {
template <class T, unsigned Dim>
-ParticleType::type PartBunchBase<T, Dim>::getPType() const {
+ParticleOrigin PartBunchBase<T, Dim>::getPOrigin() const {
+ return refPOrigin_m;
+}
+
+template <class T, unsigned Dim>
+void PartBunchBase<T, Dim>::setPOrigin(ParticleOrigin origin) {
+ refPOrigin_m = origin;
+}
+
+
+template <class T, unsigned Dim>
+ParticleType PartBunchBase<T, Dim>::getPType() const {
return refPType_m;
}
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::setPType(ParticleType::type type) {
- refPType_m = type;
+std::string PartBunchBase<T, Dim>::getPTypeString() {
+
+ switch (refPType_m) {
+ case ParticleType::ELECTRON:
+ return "ELECTRON";
+ case ParticleType::PROTON:
+ return "PROTON";
+ case ParticleType::POSITRON:
+ return "POSITRON";
+ case ParticleType::ANTIPROTON:
+ return "ANTIPROTON";
+ case ParticleType::CARBON:
+ return "CARBON";
+ case ParticleType::HMINUS:
+ return "HMINUS";
+ case ParticleType::URANIUM:
+ return "URANIUM";
+ case ParticleType::MUON:
+ return "MUON";
+ case ParticleType::DEUTERON:
+ return "DEUTERON";
+ case ParticleType::XENON:
+ return "XENON";
+ case ParticleType::H2P:
+ return "H2P";
+ case ParticleType::ALPHA:
+ return "ALPHA";
+ case ParticleType::HYDROGEN:
+ return "HYDROGEN";
+ case ParticleType::H3P:
+ return "H3P";
+ default:
+ INFOMSG("Unknown type for OPAL particles" << endl);
+ return "UNNAMED";
+ }
}
+template <class T, unsigned Dim>
+void PartBunchBase<T, Dim>::setPType(std::string type) {
+
+ if (type == "ELECTRON") {
+ refPType_m = ParticleType::ELECTRON;
+ } else if (type == "PROTON") {
+ refPType_m = ParticleType::PROTON;
+ } else if (type == "POSITRON") {
+ refPType_m = ParticleType::POSITRON;
+ } else if (type == "ANTIPROTON") {
+ refPType_m = ParticleType::ANTIPROTON;
+ } else if (type == "CARBON") {
+ refPType_m = ParticleType::CARBON;
+ } else if (type == "HMINUS") {
+ refPType_m = ParticleType::HMINUS;
+ } else if (type == "URANIUM") {
+ refPType_m = ParticleType::URANIUM;
+ } else if (type == "MUON") {
+ refPType_m = ParticleType::MUON;
+ } else if (type == "DEUTERON") {
+ refPType_m = ParticleType::DEUTERON;
+ } else if (type == "XENON") {
+ refPType_m = ParticleType::XENON;
+ } else if (type == "H2P") {
+ refPType_m = ParticleType::H2P;
+ } else if (type == "ALPHA") {
+ refPType_m = ParticleType::ALPHA;
+ } else if (type == "HYDROGEN") {
+ refPType_m = ParticleType::HYDROGEN;
+ } else if (type == "H3P") {
+ refPType_m = ParticleType::H3P;
+ } else {
+ INFOMSG("Unknown type for OPAL particles. Default ParticleType::UNNAMED is set" << endl);
+ refPType_m = ParticleType::UNNAMED;
+ }
+}
+
+
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::resetQ(double q) {
const_cast<PartData *>(reference)->setQ(q);
@@ -1785,14 +1851,13 @@ double PartBunchBase<T, Dim>::getBeta(int /*i*/) {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::actT()
-{
+void PartBunchBase<T, Dim>::actT() {
// do nothing here
};
template <class T, unsigned Dim>
-const PartData *PartBunchBase<T, Dim>::getReference() const {
+const PartData* PartBunchBase<T, Dim>::getReference() const {
return reference;
}
@@ -1817,7 +1882,7 @@ void PartBunchBase<T, Dim>::calcEMean() {
eKin_m = 0.0;
- for(unsigned int k = 0; k < locNp; k++) {
+ for (unsigned int k = 0; k < locNp; k++) {
eKin_m += std::sqrt(dot(P[k], P[k]) + 1.0);
}
@@ -1830,9 +1895,9 @@ void PartBunchBase<T, Dim>::calcEMean() {
}
template <class T, unsigned Dim>
-Inform &PartBunchBase<T, Dim>::print(Inform &os) {
+Inform& PartBunchBase<T, Dim>::print(Inform& os) {
- if(getTotalNum() != 0) { // to suppress Nans
+ if (getTotalNum() != 0) { // to suppress Nans
Inform::FmtFlags_t ff = os.flags();
double pathLength = get_sPos();
@@ -1887,12 +1952,10 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
std::vector<double> loc_moments(4 * Dim + Dim * ( 2 * Dim + 1 ));
long int totalNum = this->getTotalNum();
- if (!Options::amr && OpalData::getInstance()->isInOPALCyclMode())
- {
- /* FIXME After issue 287 is resolved this shouldn't be necessary
- * anymore
- */
- for(unsigned long k = 0; k < localNum; ++ k) {
+ if (!Options::amr && OpalData::getInstance()->isInOPALCyclMode()) {
+ //FIXME After issue 287 is resolved this shouldn't be necessary anymore
+
+ for (unsigned long k = 0; k < localNum; ++ k) {
if (ID[k] == 0) {
part[1] = P[k](0);
part[3] = P[k](1);
@@ -1904,7 +1967,7 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
unsigned int l = 2 * Dim;
for (unsigned int i = 0; i < 2 * Dim; ++i) {
loc_moments[i] -= part[i];
- for(unsigned int j = 0; j <= i; j++) {
+ for (unsigned int j = 0; j <= i; j++) {
loc_moments[l++] -= part[i] * part[j];
}
}
@@ -1922,7 +1985,7 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
allreduce(totalNum, 1, std::less<long int>());
}
- for(unsigned long k = 0; k < localNum; ++ k) {
+ for (unsigned long k = 0; k < localNum; ++ k) {
part[1] = P[k](0);
part[3] = P[k](1);
part[5] = P[k](2);
@@ -1933,7 +1996,7 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
unsigned int l = 2 * Dim;
for (unsigned int i = 0; i < 2 * Dim; ++i) {
loc_moments[i] += part[i];
- for(unsigned int j = 0; j <= i; j++) {
+ for (unsigned int j = 0; j <= i; j++) {
loc_moments[l++] += part[i] * part[j];
}
}
@@ -1953,7 +2016,7 @@ size_t PartBunchBase<T, Dim>::calcMoments() {
unsigned int l = 2 * Dim;
for (unsigned int i = 0; i < 2 * Dim; ++i) {
- for(unsigned int j = 0; j <= i; j++) {
+ for (unsigned int j = 0; j <= i; j++) {
moments_m(i, j) = loc_moments[l++];
moments_m(j, i) = moments_m(i, j);
}
@@ -1987,19 +2050,19 @@ void PartBunchBase<T, Dim>::calcMomentsInitial() {
double part[2 * Dim];
- for(unsigned int i = 0; i < 2 * Dim; ++i) {
+ for (unsigned int i = 0; i < 2 * Dim; ++i) {
centroid_m[i] = 0.0;
- for(unsigned int j = 0; j <= i; ++j) {
+ for (unsigned int j = 0; j <= i; ++j) {
moments_m(i, j) = 0.0;
moments_m(j, i) = moments_m(i, j);
}
}
- for(size_t k = 0; k < pbin_m->getNp(); k++) {
- for(int binNumber = 0; binNumber < pbin_m->getNBins(); binNumber++) {
+ for (size_t k = 0; k < pbin_m->getNp(); k++) {
+ for (int binNumber = 0; binNumber < pbin_m->getNBins(); binNumber++) {
std::vector<double> p;
- if(pbin_m->getPart(k, binNumber, p)) {
+ if (pbin_m->getPart(k, binNumber, p)) {
part[0] = p.at(0);
part[1] = p.at(3);
part[2] = p.at(1);
@@ -2007,9 +2070,9 @@ void PartBunchBase<T, Dim>::calcMomentsInitial() {
part[4] = p.at(2);
part[5] = p.at(5);
- for(unsigned int i = 0; i < 2 * Dim; ++i) {
+ for (unsigned int i = 0; i < 2 * Dim; ++i) {
centroid_m[i] += part[i];
- for(unsigned int j = 0; j <= i; ++j) {
+ for (unsigned int j = 0; j <= i; ++j) {
moments_m(i, j) += part[i] * part[j];
}
}
@@ -2017,8 +2080,8 @@ void PartBunchBase<T, Dim>::calcMomentsInitial() {
}
}
- for(unsigned int i = 0; i < 2 * Dim; ++i) {
- for(unsigned int j = 0; j < i; ++j) {
+ for (unsigned int i = 0; i < 2 * Dim; ++i) {
+ for (unsigned int j = 0; j < i; ++j) {
moments_m(j, i) = moments_m(i, j);
}
}
@@ -2035,7 +2098,7 @@ double PartBunchBase<T, Dim>::calculateAngle(double x, double y) {
template <class T, unsigned Dim>
-Inform &operator<<(Inform &os, PartBunchBase<T, Dim> &p) {
+Inform& operator<<(Inform &os, PartBunchBase<T, Dim>& p) {
return p.print(os);
}
@@ -2070,6 +2133,7 @@ void PartBunchBase<T, Dim>::swap(unsigned int i, unsigned int j) {
std::swap(Bin[i], Bin[j]);
std::swap(dt[i], dt[j]);
std::swap(PType[i], PType[j]);
+ std::swap(POrigin[i], POrigin[j]);
std::swap(TriID[i], TriID[j]);
std::swap(cavityGapCrossed[i], cavityGapCrossed[j]);
std::swap(bunchNum[i], bunchNum[j]);
@@ -2095,26 +2159,26 @@ void PartBunchBase<T, Dim>::setBCForDCBeam() {
template <class T, unsigned Dim>
-void PartBunchBase<T, Dim>::updateFields(const Vector_t &/*hr*/, const Vector_t &/*origin*/) {
+void PartBunchBase<T, Dim>::updateFields(const Vector_t& /*hr*/, const Vector_t& /*origin*/) {
}
template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::setup(AbstractParticle<T, Dim>* pb) {
+
pb->addAttribute(P);
pb->addAttribute(Q);
pb->addAttribute(M);
pb->addAttribute(Phi);
pb->addAttribute(Ef);
pb->addAttribute(Eftmp);
-
pb->addAttribute(Bf);
pb->addAttribute(Bin);
pb->addAttribute(dt);
pb->addAttribute(PType);
+ pb->addAttribute(POrigin);
pb->addAttribute(TriID);
pb->addAttribute(cavityGapCrossed);
-
pb->addAttribute(bunchNum);
boundpTimer_m = IpplTimings::getTimer("Boundingbox");
@@ -2128,7 +2192,6 @@ void PartBunchBase<T, Dim>::setup(AbstractParticle<T, Dim>* pb) {
distrCreate_m = IpplTimings::getTimer("Create Distr");
distrReload_m = IpplTimings::getTimer("Load Distr");
-
globalPartPerNode_m = std::unique_ptr<size_t[]>(new size_t[Ippl::getNodes()]);
pmsg_m.release();
@@ -2213,7 +2276,7 @@ template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::update() {
try {
pbase->update();
- } catch (const IpplException &ex) {
+ } catch (const IpplException& ex) {
throw OpalException(ex.where(), ex.what());
}
}
@@ -2222,7 +2285,7 @@ template <class T, unsigned Dim>
void PartBunchBase<T, Dim>::update(const ParticleAttrib<char>& canSwap) {
try {
pbase->update(canSwap);
- } catch (const IpplException &ex) {
+ } catch (const IpplException& ex) {
throw OpalException(ex.where(), ex.what());
}
}
@@ -2282,4 +2345,4 @@ FMatrix<double, 2 * Dim, 2 * Dim> PartBunchBase<T, Dim>::getSigmaMatrix() {
return sigmaMatrix;
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Air.h b/src/Classic/Physics/Air.h
index 139b5c17d9f6e070880956c38fd31a3fe89ee0a6..ee45271c8227589287ff3b75857f86c02801bcf6 100644
--- a/src/Classic/Physics/Air.h
+++ b/src/Classic/Physics/Air.h
@@ -1,3 +1,24 @@
+//
+// Class Air
+// Defines Air as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef AIR_H
#define AIR_H
@@ -9,11 +30,12 @@ namespace Physics {
Air():
Material(7,
14,
- 0.0012,
- 37.99,
+ 1.205e-3,
+ 36.62,
85.7,
- {{3.350, 1.683e3, 1.900e3, 2.513e-2}})
+ {{2.954, 3.350, 1.683e3, 1.900e3, 2.513e-2,
+ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/AluminaAL2O3.h b/src/Classic/Physics/AluminaAL2O3.h
index ee2609d7398d4ef019de7e975ff0928882a8a14d..2288756f2ff65da4a9df743559b41142155ca4e7 100644
--- a/src/Classic/Physics/AluminaAL2O3.h
+++ b/src/Classic/Physics/AluminaAL2O3.h
@@ -1,3 +1,24 @@
+//
+// Class AluminaAL2O3
+// Defines Alumina as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef ALUMINAAL2O3_H
#define ALUMINAAL2O3_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
AluminaAL2O3():
Material(50,
- 101.96,
+ 101.9600768,
3.97,
27.94,
145.2,
- {{7.227, 1.121e4, 3.864e2, 4.474e-3}})
+ {{1.187e1, 1.343e1, 1.069e4, 7.723e2, 2.153e-2,
+ 5.4, 0.53, 103.1, 3.931, 7.767}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Aluminum.h b/src/Classic/Physics/Aluminum.h
index 5e31fb46eeceeff7da8feaa8fe34e408d7bbe608..fbdf7b36b26ef6e7a12e7abff170a39c9e56406e 100644
--- a/src/Classic/Physics/Aluminum.h
+++ b/src/Classic/Physics/Aluminum.h
@@ -1,3 +1,24 @@
+//
+// Class Aluminum
+// Defines Aluminum as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef ALUMINUM_H
#define ALUMINUM_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
Aluminum():
Material(13,
- 26.98,
- 2.7,
+ 26.9815384,
+ 2.699,
24.01,
166.0,
- {{4.739, 2.766e3, 1.645e2, 2.023e-2}})
+ {{4.154, 4.739, 2.766e3, 1.645e2, 2.023e-2,
+ 2.5, 0.625, 45.7, 0.1, 4.359}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Beryllium.h b/src/Classic/Physics/Beryllium.h
index d7175db2ade1009bd54317830d588388a8def01d..cf4b7df36f1fc788a008b92cf0f172f9d66b8fb8 100644
--- a/src/Classic/Physics/Beryllium.h
+++ b/src/Classic/Physics/Beryllium.h
@@ -1,3 +1,24 @@
+//
+// Class Beryllium
+// Defines Beryllium as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef BERYLLIUM_H
#define BERYLLIUM_H
@@ -8,12 +29,13 @@ namespace Physics {
public:
Beryllium():
Material(4,
- 9.012,
+ 9.0121831,
1.848,
65.19,
63.7,
- {{2.590, 9.660e2, 1.538e2, 3.475e-2}})
+ {{2.248, 2.590, 9.660e2, 1.538e2, 3.475e-2,
+ 2.1895, 0.47183, 7.2362, 134.30, 197.96}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/BoronCarbide.h b/src/Classic/Physics/BoronCarbide.h
index 6c07c03e7a7284d6e11456a8403fde41180045df..27ec86d9f94fe4a090bfffcec5063ce5c643f42e 100644
--- a/src/Classic/Physics/BoronCarbide.h
+++ b/src/Classic/Physics/BoronCarbide.h
@@ -1,3 +1,24 @@
+//
+// Class BoronCarbide
+// Defines BoronCarbide as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef BORONCARBIDE_H
#define BORONCARBIDE_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
BoronCarbide():
Material(26,
- 55.25,
- 2.48,
- 50.14,
+ 55.251,
+ 2.52,
+ 50.13,
84.7,
- {{3.963, 6065.0, 1243.0, 7.782e-3}})
+ {{3.519, 3.963, 6065.0, 1243.0, 7.782e-3,
+ 5.013, 0.4707, 85.8, 16.55, 3.211}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Copper.h b/src/Classic/Physics/Copper.h
index 9e3c70697284fffc94e0d50e3677c1e93ff09d22..583127a48ae1dbf5ffcc73d920921a199169882d 100644
--- a/src/Classic/Physics/Copper.h
+++ b/src/Classic/Physics/Copper.h
@@ -1,3 +1,24 @@
+//
+// Class Copper
+// Defines Copper as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef COPPER_H
#define COPPER_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
Copper():
Material(29,
- 63.54,
+ 63.546,
8.96,
12.86,
322.0,
- {{4.194, 4.649e3, 8.113e1, 2.242e-2}})
+ {{3.969, 4.194, 4.649e3, 8.113e1, 2.242e-2,
+ 3.114, 0.5236, 76.67, 7.62, 6.385}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Gold.h b/src/Classic/Physics/Gold.h
index 506e54cd11a0df332a966e35025da60ef4ebeccb..b742fddb129dc31786a637d03317c9e949a699a1 100644
--- a/src/Classic/Physics/Gold.h
+++ b/src/Classic/Physics/Gold.h
@@ -1,3 +1,24 @@
+//
+// Class Gold
+// Defines Gold as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef GOLD_H
#define GOLD_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
Gold():
Material(79,
- 197,
- 19.3,
+ 196.966570,
+ 19.32,
6.46,
790.0,
- {{5.458, 7.852e3, 9.758e2, 2.077e-2}})
+ {{4.844, 5.458, 7.852e3, 9.758e2, 2.077e-2,
+ 3.223, 0.5883, 232.7, 2.954, 1.05}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Graphite.h b/src/Classic/Physics/Graphite.h
index e1cbf0cb21bba5a018dc83d692d7f53e061fd447..2113ed11554a20add7cb6ea0b2ea026cb177b218 100644
--- a/src/Classic/Physics/Graphite.h
+++ b/src/Classic/Physics/Graphite.h
@@ -1,3 +1,24 @@
+//
+// Class Graphite
+// Defines Graphite as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef GRAPHITE_H
#define GRAPHITE_H
@@ -12,8 +33,9 @@ namespace Physics {
2.210,
42.70,
78.0,
- {{2.601, 1.701e3, 1.279e3, 1.638e-2}})
+ {{0.0, 2.601, 1.701e3, 1.279e3, 1.638e-2,
+ 3.80133, 0.41590, 12.9966, 117.83, 242.28}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/GraphiteR6710.h b/src/Classic/Physics/GraphiteR6710.h
index 4cddd6fe18166bcc5a8844dc31670d2151bcde3b..7a2cad7ab2657b442e348fefcacf8dee9c888fef 100644
--- a/src/Classic/Physics/GraphiteR6710.h
+++ b/src/Classic/Physics/GraphiteR6710.h
@@ -1,3 +1,24 @@
+//
+// Class GraphiteR6710
+// Defines GraphiteR6710 as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef GRAPHITER6710_H
#define GRAPHITER6710_H
@@ -12,8 +33,9 @@ namespace Physics {
1.88,
42.70,
78.0,
- {{2.601, 1.701e3, 1.279e3, 1.638e-2}})
+ {{0.0, 2.601, 1.701e3, 1.279e3, 1.638e-2,
+ 3.80133, 0.41590, 12.9966, 117.83, 242.28}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Kapton.h b/src/Classic/Physics/Kapton.h
index 9331da245537a838c98c8def68b0a9cd3d2e3811..2c93455e9ff3b07a2fc9e5aca92ea156c750a1be 100644
--- a/src/Classic/Physics/Kapton.h
+++ b/src/Classic/Physics/Kapton.h
@@ -1,3 +1,24 @@
+//
+// Class Kapton
+// Defines Kapton as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef KAPTON_H
#define KAPTON_H
@@ -9,11 +30,12 @@ namespace Physics {
Kapton():
Material(6,
12,
- 1.4,
- 39.95,
- 79.6, // source: pstar from NIST
- {{2.601, 1.701e3, 1.279e3, 1.638e-2}})
+ 1.420,
+ 40.58,
+ 79.6,
+ {{0.0, 2.601, 1.701e3, 1.279e3, 1.638e-2,
+ 3.83523, 0.42993, 12.6125, 227.41, 188.97}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Material.cpp b/src/Classic/Physics/Material.cpp
index 74235251e0c82fd4ab2f22c9e09dbc68dbb08cc2..e3d2a6767305535458c878fa9790c18e46e5143c 100644
--- a/src/Classic/Physics/Material.cpp
+++ b/src/Classic/Physics/Material.cpp
@@ -1,3 +1,20 @@
+//
+// Class Material
+// Base class for representing materials
+//
+// Copyright (c) 2019 - 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 "Physics/Material.h"
#include "Physics/Air.h"
#include "Physics/AluminaAL2O3.h"
@@ -22,7 +39,7 @@ using namespace Physics;
std::map<std::string, std::shared_ptr<Material> > Material::protoTable_sm;
-std::shared_ptr<Material> Material::addMaterial(const std::string &name,
+std::shared_ptr<Material> Material::addMaterial(const std::string& name,
std::shared_ptr<Material> mat_ptr) {
std::string nameUp = Util::toUpper(name);
if (protoTable_sm.find(nameUp) != protoTable_sm.end())
@@ -33,7 +50,7 @@ std::shared_ptr<Material> Material::addMaterial(const std::string &name,
return mat_ptr;
}
-std::shared_ptr<Material> Material::getMaterial(const std::string &name) {
+std::shared_ptr<Material> Material::getMaterial(const std::string& name) {
std::string nameUp = Util::toUpper(name);
if (protoTable_sm.find(nameUp) != protoTable_sm.end()) return protoTable_sm[nameUp];
@@ -73,4 +90,4 @@ namespace {
titanium);
auto water = Material::addMaterial("Water",
std::shared_ptr<Material>(new Water()));
-}
\ No newline at end of file
+}
diff --git a/src/Classic/Physics/Material.h b/src/Classic/Physics/Material.h
index c9c7a1d4219c2101e5e6a5bc5e5e018a395d3ca9..41f20e78fd7846ccc5bf0dcf038cc569ebed3f9a 100644
--- a/src/Classic/Physics/Material.h
+++ b/src/Classic/Physics/Material.h
@@ -1,3 +1,20 @@
+//
+// Class Material
+// Base class for representing materials
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef MATERIAL_H
#define MATERIAL_H
@@ -11,10 +28,16 @@ namespace Physics {
class Material {
public:
enum FitCoeffs {
- A2 = 0,
+ A1 = 0,
+ A2,
A3,
A4,
- A5
+ A5,
+ B1,
+ B2,
+ B3,
+ B4,
+ B5
};
Material(double atomicNumber,
@@ -22,7 +45,7 @@ namespace Physics {
double massDensity,
double radiationLength,
double meanExcitationEnergy,
- std::array<double, 4> fitCoefficients):
+ std::array<double, 10> fitCoefficients):
atomicNumber_m(atomicNumber),
atomicMass_m(atomicMass),
massDensity_m(massDensity),
@@ -38,8 +61,8 @@ namespace Physics {
double getMeanExcitationEnergy() const; // [eV]
double getStoppingPowerFitCoefficients(FitCoeffs n) const;
- static std::shared_ptr<Material> getMaterial(const std::string &name);
- static std::shared_ptr<Material> addMaterial(const std::string &name,
+ static std::shared_ptr<Material> getMaterial(const std::string& name);
+ static std::shared_ptr<Material> addMaterial(const std::string& name,
std::shared_ptr<Material> mat_ptr);
private:
static
@@ -50,7 +73,7 @@ namespace Physics {
const double massDensity_m;
const double radiationLength_m;
const double meanExcitationEnergy_m;
- const std::array<double,4> stoppingPowerFitCoefficients_m;
+ const std::array<double,10> stoppingPowerFitCoefficients_m;
};
inline
diff --git a/src/Classic/Physics/Molybdenum.h b/src/Classic/Physics/Molybdenum.h
index 8dda9e5e9a38dbec6a6c76fd4ea5d10da85922c3..346bc43127e7549abfc8acdb1e7521540418eef3 100644
--- a/src/Classic/Physics/Molybdenum.h
+++ b/src/Classic/Physics/Molybdenum.h
@@ -1,3 +1,24 @@
+//
+// Class Molybdenum
+// Defines Molybdenum as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef MOLYBDENUM_H
#define MOLYBDENUM_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
Molybdenum():
Material(42,
- 95.94,
+ 95.95,
10.22,
- 9.8,
+ 9.80,
424.0,
- {{7.248, 9.545e3, 4.802e2, 5.376e-3}})
+ {{6.424, 7.248, 9.545e3, 4.802e2, 5.376e-3,
+ 9.276, 0.418, 157.1, 8.038, 1.29}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Mylar.h b/src/Classic/Physics/Mylar.h
index 6fa7aa629178af61063adbab2d46dd3006629657..8b98133411f13895f7a94a085227b5274abcc944 100644
--- a/src/Classic/Physics/Mylar.h
+++ b/src/Classic/Physics/Mylar.h
@@ -1,3 +1,24 @@
+//
+// Class Mylar
+// Defines Mylar as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef MYLAR_H
#define MYLAR_H
@@ -9,11 +30,12 @@ namespace Physics {
Mylar():
Material(6.702,
12.88,
- 1.4,
+ 1.400,
39.95,
78.7,
- {{3.350, 1683, 1900, 2.513e-02}})
+ {{2.954, 3.350, 1683, 1900, 2.513e-02,
+ 1.9259, 0.5550, 27.15125, 26.0665, 6.2768}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Physics.h b/src/Classic/Physics/Physics.h
index 922bea72559bd4506e26f232d620dfedd9bb59d7..00f85a6af77d1d2c349e27feb3927bc6a1c44541 100644
--- a/src/Classic/Physics/Physics.h
+++ b/src/Classic/Physics/Physics.h
@@ -2,7 +2,7 @@
// Class Physics
// A namespace defining various mathematical and physical constants.
//
-// Copyright (c) 2015-2019, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin
+// Copyright (c) 2015-2021, Christof Metzger-Kraus, Helmholtz-Zentrum Berlin
// Pedro Calvo, CIEMAT, Spain
// All rights reserved
//
@@ -24,9 +24,6 @@
#define CLASSIC_Physics_HH
-// Class Physics
-// ------------------------------------------------------------------------
-
namespace Physics {
/// The value of \f[ \pi \f]
@@ -128,6 +125,9 @@ namespace Physics {
/// The xenon rest mass in GeV
constexpr double m_xe = 124 * amu;
+ /// The alpha particle rest mass in GeV
+ constexpr double m_alpha = 4.001506179127 * amu;
+
/// The hydrogen atom rest mass in GeV
constexpr double m_h = 1.00782503224 * amu;
@@ -137,7 +137,7 @@ namespace Physics {
/// The H3+ rest mass in GeV
constexpr double m_h3p = 3.02293 * amu;
- constexpr double PMASS = 1.67262192369e-27; // kg
+ constexpr double PMASS = 1.67262192369e-27; // kg
constexpr double EMASS = 9.1093837015e-31; // kg
@@ -149,7 +149,6 @@ namespace Physics {
constexpr double e0m = 1.75882001076e+11;
// e/mc
constexpr double e0mc = e0m / c;
-
};
-#endif // CLASSIC_Physics_HH
\ No newline at end of file
+#endif // CLASSIC_Physics_HH
diff --git a/src/Classic/Physics/Titanium.h b/src/Classic/Physics/Titanium.h
index 0b82127264a7200406c08e5335717bfaee6a24ba..a28b8f501f73210362fdac2b04868f7c7f65fe66 100644
--- a/src/Classic/Physics/Titanium.h
+++ b/src/Classic/Physics/Titanium.h
@@ -1,3 +1,24 @@
+//
+// Class Titanium
+// Defines Titanium as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef TITANIUM_H
#define TITANIUM_H
@@ -10,12 +31,13 @@ namespace Physics {
public:
Titanium():
Material(22,
- 47.8,
- 4.54,
+ 47.867,
+ 4.540,
16.16,
233.0,
- {{5.489, 5.260e3, 6.511e2, 8.930e-3}})
+ {{4.858, 5.489, 5.260e3, 6.511e2, 8.930e-3,
+ 4.71, 0.5087, 65.28, 8.806, 5.948}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Physics/Water.h b/src/Classic/Physics/Water.h
index 5022243b9061da68ac28ad1ba8074201a58bb1d9..7fcf1dcd3a9581f1e9469a3edcbc2247bc58a17e 100644
--- a/src/Classic/Physics/Water.h
+++ b/src/Classic/Physics/Water.h
@@ -1,3 +1,24 @@
+//
+// Class Water
+// Defines Water as a material for particle-matter interactions
+//
+// Data taken from Standard Atomic Weight 2019 (https://www.qmul.ac.uk/sbcs/iupac/AtWt/),
+// atomic properties from PDG database (https://pdg.lbl.gov/2020/AtomicNuclearProperties)
+// and fit coefficients from ICRU Report 49.
+//
+// Copyright (c) 2019 - 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/>.
+//
#ifndef WATER_H
#define WATER_H
@@ -8,12 +29,13 @@ namespace Physics {
public:
Water():
Material(10,
- 18,
+ 18.0152,
1,
36.08,
- 75.0, // source: pstar from NIST
- {{2.199, 2.393e3, 2.699e3, 1.568e-2}})
+ 75.0,
+ {{4.015, 4.542, 3.955e3, 4.847e2, 7.904e-3,
+ 2.9590, 0.53255, 34.247, 60.655, 15.153}})
{ }
};
}
-#endif
\ No newline at end of file
+#endif
diff --git a/src/Classic/Solvers/BeamStrippingPhysics.cpp b/src/Classic/Solvers/BeamStrippingPhysics.cpp
index bb4cafcb6f11d84cd35af762d3bb7ad8d92e0d76..c03bc202eda16a176004191bb1340b71632e38ec 100644
--- a/src/Classic/Solvers/BeamStrippingPhysics.cpp
+++ b/src/Classic/Solvers/BeamStrippingPhysics.cpp
@@ -1,9 +1,9 @@
//
// Class BeamStrippingPhysics
-// This class provides beam stripping physical processes as
-// particle matter interaction type.
+// Defines the physical processes of residual gas
+// interactions and Lorentz stripping
//
-// 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
@@ -20,7 +20,6 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-
#include "Solvers/BeamStrippingPhysics.hh"
#include "AbsBeamline/BeamStripping.h"
@@ -29,6 +28,7 @@
#include "Algorithms/PartData.h"
#include "Physics/Physics.h"
#include "Structure/LossDataSink.h"
+#include "Utilities/GeneralClassicException.h"
#include "Utilities/Options.h"
#include "Utility/Inform.h"
@@ -40,37 +40,31 @@
#include <sys/time.h>
#include <boost/math/special_functions/chebyshev.hpp>
-using Physics::q_e;
-using Physics::m_e;
-using Physics::m_hm;
-using Physics::m_p;
-using Physics::m_h;
-using Physics::m_h2p;
namespace {
struct InsideTester {
virtual ~InsideTester() {}
- virtual bool checkHit(const Vector_t &R) = 0;
- virtual double getPressure(const Vector_t &R) = 0;
+ virtual bool checkHit(const Vector_t& R) = 0;
+ virtual double getPressure(const Vector_t& R) = 0;
};
struct BeamStrippingInsideTester: public InsideTester {
- explicit BeamStrippingInsideTester(ElementBase * el) {
+ explicit BeamStrippingInsideTester(ElementBase* el) {
bstp_m = static_cast<BeamStripping*>(el);
}
- virtual bool checkHit(const Vector_t &R) {
+ virtual bool checkHit(const Vector_t& R) {
return bstp_m->checkPoint(R(0), R(1), R(2));
}
- double getPressure(const Vector_t &R) {
+ double getPressure(const Vector_t& R) {
return bstp_m->checkPressure(R(0), R(1));
}
private:
- BeamStripping *bstp_m;
+ BeamStripping* bstp_m;
};
}
-BeamStrippingPhysics::BeamStrippingPhysics(const std::string &name, ElementBase *element):
+BeamStrippingPhysics::BeamStrippingPhysics(const std::string& name, ElementBase* element):
ParticleMatterInteractionHandler(name, element),
T_m(0.0),
dT_m(0.0),
@@ -86,11 +80,11 @@ BeamStrippingPhysics::BeamStrippingPhysics(const std::string &name, ElementBase
rediffusedStat_m(0),
locPartsInMat_m(0)
{
- bstp_m = dynamic_cast<BeamStripping *>(getElement());
+ bstp_m = dynamic_cast<BeamStripping* >(getElement());
bstpshape_m = element_ref_m->getType();
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(element_ref_m->getName(), !Options::asciidump));
- const gsl_rng_type * T;
+ const gsl_rng_type* T;
gsl_rng_env_setup();
struct timeval tv; // Seed generation based on time
gettimeofday(&tv,0);
@@ -98,8 +92,6 @@ BeamStrippingPhysics::BeamStrippingPhysics(const std::string &name, ElementBase
T = gsl_rng_default; // Generator setup
r_m = gsl_rng_alloc (T);
gsl_rng_set(r_m, mySeed);
-
- gas_m = bstp_m->getResidualGas();
}
@@ -112,30 +104,30 @@ const std::string BeamStrippingPhysics::getType() const {
return "BeamStrippingPhysics";
}
-void BeamStrippingPhysics::apply(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &/*boundingSphere*/) {
+void BeamStrippingPhysics::apply(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& /*boundingSphere*/) {
dT_m = bunch->getdT();
- double mass = bunch->getM()*1E-9;
+ ParticleType pType = bunch->getPType();
+ if (pType == ParticleType::PROTON ||
+ pType == ParticleType::HMINUS ||
+ pType == ParticleType::H2P ||
+ pType == ParticleType::HYDROGEN ) {
- if (bunch->get_sPos() != 0) {
- if (std::abs(mass-m_hm) < 1E-6 ||
- std::abs(mass-m_p) < 1E-6 ||
- std::abs(mass-m_h2p) < 1E-6 ||
- std::abs(mass-m_h) < 1E-6) {
+ if (bunch->get_sPos() != 0) {
doPhysics(bunch);
}
- else {
- Inform gmsgALL("OPAL", INFORM_ALL_NODES);
- gmsgALL << getName() << ": Unsupported type of particle for residual gas interactions!"<< endl;
- gmsgALL << getName() << "-> Beam Stripping Physics not apply"<< endl;
- }
+ } else {
+ throw GeneralClassicException(
+ "BeamStrippingPhysics::apply",
+ "Particle " + bunch->getPTypeString() +
+ " is not supported for residual stripping interactions!");
}
}
-void BeamStrippingPhysics::doPhysics(PartBunchBase<double, 3> *bunch) {
+void BeamStrippingPhysics::doPhysics(PartBunchBase<double, 3>* bunch) {
/*
Do physics if
-- particle in material
@@ -148,7 +140,7 @@ void BeamStrippingPhysics::doPhysics(PartBunchBase<double, 3> *bunch) {
bool stop = bstp_m->getStop();
- InsideTester *tester;
+ InsideTester* tester;
tester = new BeamStrippingInsideTester(element_ref_m);
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
@@ -167,10 +159,10 @@ void BeamStrippingPhysics::doPhysics(PartBunchBase<double, 3> *bunch) {
double beta = std::sqrt(1.0 - 1.0 / (gamma * gamma));
double deltas = dT_m * beta * Physics::c;
- crossSection(Eng);
+ crossSection(bunch, i, Eng);
pdead_GS = gasStripping(deltas);
- if (std::abs(mass_m-m_hm) < 1E-6 && charge_m == -q_e) {
+ if (bunch->PType[i] == ParticleType::HMINUS) {
cycl_m->apply(bunch->R[i], bunch->P[i], T_m, extE, extB);
double B = 0.1 * std::sqrt(extB[0]*extB[0] + extB[1]*extB[1] + extB[2]*extB[2]); //T
double E = gamma * beta * Physics::c * B;
@@ -185,10 +177,9 @@ void BeamStrippingPhysics::doPhysics(PartBunchBase<double, 3> *bunch) {
stoppedPartStat_m++;
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i]
<< " is deleted by beam stripping" << endl;
- }
- else {
+ } else {
secondaryParticles(bunch, i, pdead_LS);
- bunch->updateNumTotal();
+ //bunch->updateNumTotal();
gmsgALL << level4 << getName()
<< ": Total number of particles after beam stripping = "
<< bunch->getTotalNum() << endl;
@@ -200,12 +191,13 @@ void BeamStrippingPhysics::doPhysics(PartBunchBase<double, 3> *bunch) {
}
-void BeamStrippingPhysics::crossSection(double Eng){
+void BeamStrippingPhysics::crossSection(PartBunchBase<double, 3>* bunch, size_t& i, double Eng){
- const double temperature = bstp_m->getTemperature(); // K
+ const double temperature = bstp_m->getTemperature(); // K
+ const ParticleType& pType = bunch->PType[i];
+ const ResidualGas& gas = bstp_m->getResidualGas();
Eng *=1E6; //keV
-
double Eth = 0.0;
double a1 = 0.0;
double a2 = 0.0;
@@ -222,9 +214,9 @@ void BeamStrippingPhysics::crossSection(double Eng){
double molecularDensity[3] ={};
- if (gas_m == "H2") {
+ if (gas == ResidualGas::H2) {
- molecularDensity[0] = 100 * pressure_m / (Physics::kB * q_e * temperature);
+ molecularDensity[0] = 100 * pressure_m / (Physics::kB * Physics::q_e * temperature);
double Emin = 0.0;
double Emax = 0.0;
@@ -233,37 +225,37 @@ void BeamStrippingPhysics::crossSection(double Eng){
double CS_c = 0.0;
double CS_total = 0.0;
- if (std::abs(mass_m-m_hm) < 1E-6 && charge_m == -q_e) {
+ if (pType == ParticleType::HMINUS) {
// Single-electron detachment - Hydrogen Production
Emin = csCoefSingle_Hminus_Chebyshev[0];
Emax = csCoefSingle_Hminus_Chebyshev[1];
for (int i = 0; i < 9; ++i)
a_m[i] = {csCoefSingle_Hminus_Chebyshev[i+2]};
CS_a = csChebyshevFitting(Eng*1E3, Emin, Emax);
-
+
// Double-electron detachment - Proton Production
Emin = csCoefDouble_Hminus_Chebyshev[0];
Emax = csCoefDouble_Hminus_Chebyshev[1];
for (int i = 0; i < 9; ++i)
a_m[i] = {csCoefDouble_Hminus_Chebyshev[i+2]};
CS_b = csChebyshevFitting(Eng*1E3, Emin, Emax);
- }
- else if (std::abs(mass_m-m_p) < 1E-6 && charge_m == q_e) {
+
+ } else if (pType == ParticleType::PROTON) {
// Single-electron capture - Hydrogen Production
Emin = csCoefSingle_Hplus_Chebyshev[0];
Emax = csCoefSingle_Hplus_Chebyshev[1];
for (int i = 0; i < 9; ++i)
a_m[i] = {csCoefSingle_Hplus_Chebyshev[i+2]};
CS_a = csChebyshevFitting(Eng*1E3, Emin, Emax);
-
+
// Double-electron capture - Hminus Production
Emin = csCoefDouble_Hplus_Chebyshev[0];
Emax = csCoefDouble_Hplus_Chebyshev[1];
for (int i = 0; i < 9; ++i)
a_m[i] = {csCoefDouble_Hplus_Chebyshev[i+2]};
CS_b = csChebyshevFitting(Eng*1E3, Emin, Emax);
- }
- else if (std::abs(mass_m-m_h) < 1E-6 && charge_m == 0.0) {
+
+ } else if (pType == ParticleType::HYDROGEN) {
// Single-electron detachment - Proton Production
Eth = csCoefProtonProduction_H_Tabata[0];
a1 = csCoefProtonProduction_H_Tabata[1];
@@ -274,7 +266,7 @@ void BeamStrippingPhysics::crossSection(double Eng){
a6 = csCoefProtonProduction_H_Tabata[6];
CS_a = csAnalyticFunctionTabata(Eng, Eth, a1, a2, 0.0, 0.0, a3, a4) +
a5 * csAnalyticFunctionTabata(Eng/a6, Eth/a6, a1, a2, 0.0, 0.0, a3, a4);
-
+
// Single-electron capture - Hminus Production
Eth = csCoefHminusProduction_H_Tabata[0];
a1 = csCoefHminusProduction_H_Tabata[1];
@@ -291,8 +283,8 @@ void BeamStrippingPhysics::crossSection(double Eng){
a12 = csCoefHminusProduction_H_Tabata[12];
CS_b = csAnalyticFunctionTabata(Eng, Eth, a1, a2, a3, a4, a5, a6) +
csAnalyticFunctionTabata(Eng, Eth, a7, a8, a9, a10, a11, a12);
- }
- else if (std::abs(mass_m-m_h2p) < 1E-6 && charge_m == q_e) {
+
+ } else if (pType == ParticleType::H2P) {
// Proton production
Eth = csCoefProtonProduction_H2plus_Tabata[0];
a1 = csCoefProtonProduction_H2plus_Tabata[1];
@@ -308,14 +300,14 @@ void BeamStrippingPhysics::crossSection(double Eng){
CS_a = csAnalyticFunctionTabata(Eng, Eth, a1, a2, 0.0, 0.0, a3, a4) +
csAnalyticFunctionTabata(Eng, Eth, a5, a6, 0.0, 0.0, a7, a8) +
a9 * csAnalyticFunctionTabata(Eng/a10, Eth/a10, a5, a6, 0.0, 0.0, a7, a8);
-
+
// Hydrogen production
Emin = csCoefHydrogenProduction_H2plus_Chebyshev[0];
Emax = csCoefHydrogenProduction_H2plus_Chebyshev[1];
for (int i = 0; i < 9; ++i)
a_m[i] = {csCoefHydrogenProduction_H2plus_Chebyshev[i+2]};
CS_b = csChebyshevFitting(Eng*1E3, Emin, Emax);
-
+
// H3+, H
Eth = csCoefH3plusProduction_H2plus_Tabata[0];
a1 = csCoefH3plusProduction_H2plus_Tabata[1];
@@ -327,14 +319,13 @@ void BeamStrippingPhysics::crossSection(double Eng){
CS_c = csAnalyticFunctionTabata(Eng, Eth, a1, a2, a3, a4, a5, a6);
}
CS_total = CS_a + CS_b + CS_c;
-
+
NCS_a = CS_a * 1E-4 * molecularDensity[0];
NCS_b = CS_b * 1E-4 * molecularDensity[0];
NCS_c = CS_c * 1E-4 * molecularDensity[0];
NCS_total = CS_total * 1E-4 * molecularDensity[0];
- }
- else if (gas_m == "AIR") {
+ } else if (gas == ResidualGas::AIR) {
static const double fMolarFraction[3] = {
// Nitrogen
@@ -342,8 +333,8 @@ void BeamStrippingPhysics::crossSection(double Eng){
//Oxygen
20.947/100,
//Argon
- 0.934/100
- };
+ 0.934/100 };
+
double CS_single[3];
double CS_double[3];
double CS_total[3];
@@ -356,9 +347,9 @@ void BeamStrippingPhysics::crossSection(double Eng){
for (int i = 0; i < 3; ++i) {
- molecularDensity[i] = 100 * pressure_m * fMolarFraction[i] / (Physics::kB * q_e * temperature);
+ molecularDensity[i] = 100 * pressure_m * fMolarFraction[i] / (Physics::kB * Physics::q_e * temperature);
- if (std::abs(mass_m-m_hm) < 1E-6 && charge_m == -q_e) {
+ if (pType == ParticleType::HMINUS) {
// Single-electron detachment - Hydrogen Production
Eth = csCoefSingle_Hminus[i][0];
for (int j = 0; j < 8; ++j)
@@ -370,8 +361,8 @@ void BeamStrippingPhysics::crossSection(double Eng){
for (int j = 0; j < 8; ++j)
b_m[i][j] = csCoefDouble_Hminus[i][j+1];
CS_double[i] = csAnalyticFunctionNakai(Eng, Eth, i);
- }
- else if (std::abs(mass_m-m_p) < 1E-6 && charge_m == q_e) {
+
+ } else if (pType == ParticleType::PROTON) {
// Single-electron capture - Hydrogen Production
Eth = csCoefSingle_Hplus[i][0];
for (int j = 0; j < 8; ++j)
@@ -386,11 +377,11 @@ void BeamStrippingPhysics::crossSection(double Eng){
if (b_m[i][7] != 0) {
CS_double[i] = csAnalyticFunctionNakai(Eng, Eth, i) +
b_m[i][6] * csAnalyticFunctionNakai(Eng/b_m[i][7], Eth/b_m[i][7], i);
- }
- else
+ } else {
CS_double[i] = csAnalyticFunctionNakai(Eng, Eth, i);
- }
- else if (std::abs(mass_m-m_h) < 1E-6 && charge_m == 0.0) {
+ }
+
+ } else if (pType == ParticleType::HYDROGEN) {
// Single-electron detachment - Proton Production
Eth = csCoefSingleLoss_H[i][0];
for (int j = 0; j < 8; ++j)
@@ -404,11 +395,10 @@ void BeamStrippingPhysics::crossSection(double Eng){
if (b_m[i][7] != 0) {
CS_double[i] = csAnalyticFunctionNakai(Eng, Eth, i) +
b_m[i][6] * csAnalyticFunctionNakai(Eng/b_m[i][7], Eth/b_m[i][7], i);
- }
- else
+ } else {
CS_double[i] = csAnalyticFunctionNakai(Eng, Eth, i);
- }
- else {
+ }
+ } else {
CS_single[i] = {0.0};
CS_double[i] = {0.0};
}
@@ -426,13 +416,12 @@ void BeamStrippingPhysics::crossSection(double Eng){
NCS_b = NCS_double_all;
NCS_total = NCS_total_all;
}
-
}
double BeamStrippingPhysics::csAnalyticFunctionNakai(double Eng, double Eth, int &i) {
- const double E_R = Physics::E_ryd*1e6 * m_h / m_e; //keV
+ const double E_R = Physics::E_ryd*1e6 * Physics::m_h / Physics::m_e; //keV
const double sigma_0 = 1E-16;
double sigma = 0.0; //cm2
if (Eng > Eth) {
@@ -447,8 +436,8 @@ double BeamStrippingPhysics::csAnalyticFunctionNakai(double Eng, double Eth, int
}
double BeamStrippingPhysics::csAnalyticFunctionTabata(double Eng, double Eth,
- double a1, double a2, double a3, double a4, double a5, double a6) {
-
+ double a1, double a2, double a3,
+ double a4, double a5, double a6) {
const double sigma_0 = 1E-16;
double sigma = 0.0; //cm2
if (Eng > Eth) {
@@ -479,7 +468,7 @@ double BeamStrippingPhysics::csChebyshevFitting(double Eng, double Emin, double
}
-bool BeamStrippingPhysics::gasStripping(double &deltas) {
+bool BeamStrippingPhysics::gasStripping(double& deltas) {
double xi = gsl_rng_uniform(r_m);
double fg = 1-std::exp(-NCS_total*deltas);
@@ -487,7 +476,7 @@ bool BeamStrippingPhysics::gasStripping(double &deltas) {
return (fg >= xi);
}
-bool BeamStrippingPhysics::lorentzStripping(double &gamma, double &E) {
+bool BeamStrippingPhysics::lorentzStripping(double& gamma, double& E) {
double xi = gsl_rng_uniform(r_m);
@@ -497,93 +486,92 @@ bool BeamStrippingPhysics::lorentzStripping(double &gamma, double &E) {
// tau = (A1/E) * std::exp(A2/E);
//Theoretical
- const double eps0 = 0.75419 * q_e;
- const double hbar = Physics::h_bar*1E9*q_e;
- const double me = m_e*1E9*q_e/(Physics::c*Physics::c);
+ const double eps0 = 0.75419 * Physics::q_e;
+ const double hbar = Physics::h_bar * 1E9 * Physics::q_e;
+ const double me = Physics::m_e * 1E9 * Physics::q_e / (Physics::c*Physics::c);
const double p = 0.0126;
const double S0 = 0.783;
const double a = 2.01407/Physics::a0;
const double k0 = std::sqrt(2 * me * eps0)/hbar;
const double N = (std::sqrt(2 * k0 * (k0+a) * (2*k0+a)))/a;
- double zT = eps0 / (q_e * E);
+ double zT = eps0 / (Physics::q_e * E);
double tau = (4 * me * zT)/(S0 * N * N * hbar * (1+p)*(1+p) * (1-1/(2*k0*zT))) * std::exp(4*k0*zT/3);
double fL = 1 - std::exp( - dT_m / (gamma * tau));
return (fL >= xi);
}
-void BeamStrippingPhysics::secondaryParticles(PartBunchBase<double, 3> *bunch, size_t &i, bool pdead_LS) {
+void BeamStrippingPhysics::secondaryParticles(PartBunchBase<double, 3>* bunch, size_t& i, bool pdead_LS) {
double r = gsl_rng_uniform(r_m);
+ const ParticleType& pType = bunch->PType[i];
+
// change the mass_m and charge_m
- if (std::abs(mass_m-m_hm) < 1E-6 && charge_m == -q_e) {
- if (pdead_LS == true)
+ if (pType == ParticleType::HMINUS) {
+ if (pdead_LS == true) {
transformToHydrogen(bunch, i);
- else {
+ } else {
if (r > NCS_b/NCS_total)
transformToHydrogen(bunch, i);
else
transformToProton(bunch, i);
}
- }
- else if (std::abs(mass_m-m_p) < 1E-6 && charge_m == q_e) {
+ } else if (pType == ParticleType::PROTON) {
if (r > NCS_b/NCS_total)
transformToHydrogen(bunch, i);
else
transformToHminus(bunch, i);
- }
- else if (std::abs(mass_m-m_h) < 1E-6 && charge_m == 0.0) {
+ } else if (pType == ParticleType::HYDROGEN) {
if (r > NCS_b/NCS_total)
transformToProton(bunch, i);
else
transformToHminus(bunch, i);
- }
- else if (std::abs(mass_m-m_h2p) < 1E-6 && charge_m == q_e) {
- if (NCS_c>NCS_b && NCS_b>NCS_a){
+ } else if (pType == ParticleType::H2P) {
+ if (NCS_c>NCS_b && NCS_b>NCS_a) {
if (r > (NCS_a+NCS_b)/NCS_total)
transformToH3plus(bunch, i);
else if (r > NCS_a/NCS_total)
transformToHydrogen(bunch, i);
else
transformToProton(bunch, i);
- }
- else if (NCS_a>NCS_b && NCS_b>NCS_c) {
+
+ } else if (NCS_a>NCS_b && NCS_b>NCS_c) {
if (r > (NCS_c+NCS_b)/NCS_total)
transformToProton(bunch, i);
else if (r > NCS_c/NCS_total)
transformToHydrogen(bunch, i);
else
transformToH3plus(bunch, i);
- }
- else if (NCS_a>NCS_b && NCS_c>NCS_a) {
+
+ } else if (NCS_a>NCS_b && NCS_c>NCS_a) {
if (r > (NCS_a+NCS_b)/NCS_total)
transformToH3plus(bunch, i);
else if (r > NCS_b/NCS_total)
transformToProton(bunch, i);
else
transformToHydrogen(bunch, i);
- }
- else if (NCS_a>NCS_c && NCS_c>NCS_b) {
+
+ } else if (NCS_a>NCS_c && NCS_c>NCS_b) {
if (r > (NCS_c+NCS_b)/NCS_total)
transformToProton(bunch, i);
else if (r > NCS_b/NCS_total)
transformToH3plus(bunch, i);
else
transformToHydrogen(bunch, i);
- }
- else if (NCS_b>NCS_c && NCS_b>NCS_a && NCS_a>NCS_c) {
+
+ } else if (NCS_b>NCS_c && NCS_b>NCS_a && NCS_a>NCS_c) {
if (r > (NCS_c+NCS_a)/NCS_total)
transformToHydrogen(bunch, i);
else if (r > NCS_c/NCS_total)
transformToProton(bunch, i);
else
transformToH3plus(bunch, i);
- }
- else {
+
+ } else {
if (r > (NCS_c+NCS_a)/NCS_total)
transformToHydrogen(bunch, i);
else if (r > NCS_a/NCS_total)
@@ -593,39 +581,43 @@ void BeamStrippingPhysics::secondaryParticles(PartBunchBase<double, 3> *bunch, s
}
}
- bunch->PType[i] = ParticleType::SECONDARY;
+ bunch->POrigin[i] = ParticleOrigin::SECONDARY;
if (bunch->weHaveBins())
bunch->Bin[bunch->getLocalNum()-1] = bunch->Bin[i];
}
-void BeamStrippingPhysics::transformToProton(PartBunchBase<double, 3> *bunch, size_t &i) {
+void BeamStrippingPhysics::transformToProton(PartBunchBase<double, 3>* bunch, size_t& i) {
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i]
<< " is transformed to proton" << endl;
- bunch->M[i] = m_p;
- bunch->Q[i] = q_e;
+ bunch->M[i] = Physics::m_p;
+ bunch->Q[i] = Physics::q_e;
+ bunch->PType[i] = ParticleType::PROTON;
}
-void BeamStrippingPhysics::transformToHydrogen(PartBunchBase<double, 3> *bunch, size_t &i) {
+void BeamStrippingPhysics::transformToHydrogen(PartBunchBase<double, 3>* bunch, size_t& i) {
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i]
<< " is transformed to neutral hydrogen" << endl;
- bunch->M[i] = m_h;
+ bunch->M[i] = Physics::m_h;
bunch->Q[i] = 0.0;
+ bunch->PType[i] = ParticleType::HYDROGEN;
}
-void BeamStrippingPhysics::transformToHminus(PartBunchBase<double, 3> *bunch, size_t &i) {
+void BeamStrippingPhysics::transformToHminus(PartBunchBase<double, 3>* bunch, size_t& i) {
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i]
<< " is transformed to negative hydrogen ion" << endl;
- bunch->M[i] = m_hm;
- bunch->Q[i] = -q_e;
+ bunch->M[i] = Physics::m_hm;
+ bunch->Q[i] = -Physics::q_e;
+ bunch->PType[i] = ParticleType::HMINUS;
}
-void BeamStrippingPhysics::transformToH3plus(PartBunchBase<double, 3> *bunch, size_t &i) {
+void BeamStrippingPhysics::transformToH3plus(PartBunchBase<double, 3>* bunch, size_t& i) {
Inform gmsgALL("OPAL", INFORM_ALL_NODES);
gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i]
<< " is transformed to H3+" << endl;
bunch->M[i] = Physics::m_h3p;
- bunch->Q[i] = q_e;
+ bunch->Q[i] = Physics::q_e;
+ bunch->PType[i] = ParticleType::H3P;
}
void BeamStrippingPhysics::print(Inform& /*msg*/) {
@@ -701,4 +693,4 @@ const double BeamStrippingPhysics::csCoefHydrogenProduction_H2plus_Chebyshev[11]
2.00E+03, 1.00E+05, -70.670173645, -0.632612288, -0.6065212488, -0.0915143117, -0.0121710282, 0.0168179292, 0.0104796877, 0, 0
};
double BeamStrippingPhysics::a_m[9] = {};
-double BeamStrippingPhysics::b_m[3][9] = {};
\ No newline at end of file
+double BeamStrippingPhysics::b_m[3][9] = {};
diff --git a/src/Classic/Solvers/BeamStrippingPhysics.hh b/src/Classic/Solvers/BeamStrippingPhysics.hh
index 8aa88ae3b72cc416dd63bb961fbaa830c633822a..feed19f14f7167f70c862366456572faa2a050cc 100644
--- a/src/Classic/Solvers/BeamStrippingPhysics.hh
+++ b/src/Classic/Solvers/BeamStrippingPhysics.hh
@@ -1,9 +1,9 @@
//
// Class BeamStrippingPhysics
-// This class provides beam stripping physical processes as
-// particle matter interaction type.
+// Defines the physical processes of residual gas
+// interactions and Lorentz stripping
//
-// 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
@@ -43,13 +43,13 @@ class BeamStrippingPhysics: public ParticleMatterInteractionHandler {
public:
- BeamStrippingPhysics(const std::string &name, ElementBase *element);
+ 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);
+ virtual void apply(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& boundingSphere);
virtual const std::string getType() const;
virtual void print(Inform& msg);
@@ -60,46 +60,47 @@ public:
virtual size_t getParticlesInMat() {return locPartsInMat_m;}
virtual unsigned getRediffused() {return rediffusedStat_m;}
virtual unsigned int getNumEntered() {return bunchToMatStat_m;}
- inline void doPhysics(PartBunchBase<double, 3> *bunch);
+
+ inline void doPhysics(PartBunchBase<double, 3>* bunch);
private:
- void crossSection(double Eng);
+ void crossSection(PartBunchBase<double, 3>* bunch, size_t& i, double Eng);
- double csAnalyticFunctionNakai(double Eng, double Eth, int &i);
+ double csAnalyticFunctionNakai(double Eng, double Eth, int& i);
double csAnalyticFunctionTabata(double Eng, double Eth,
- double a1, double a2, double a3, double a4, double a5, double a6);
+ double a1, double a2, double a3,
+ double a4, double a5, double a6);
double csChebyshevFitting(double Eng, double Emin, double Emax);
- bool gasStripping(double &deltas);
-
- bool lorentzStripping(double &gamma, double &E);
+ bool gasStripping(double& deltas);
+ bool lorentzStripping(double& gamma, double& E);
- void secondaryParticles(PartBunchBase<double, 3> *bunch, size_t &i, bool pdead_LS);
- void transformToProton(PartBunchBase<double, 3> *bunch, size_t &i);
- void transformToHydrogen(PartBunchBase<double, 3> *bunch, size_t &i);
- void transformToHminus(PartBunchBase<double, 3> *bunch, size_t &i);
- void transformToH3plus(PartBunchBase<double, 3> *bunch, size_t &i);
+ void secondaryParticles(PartBunchBase<double, 3>* bunch, size_t& i, bool pdead_LS);
+ void transformToProton(PartBunchBase<double, 3>* bunch, size_t& i);
+ void transformToHydrogen(PartBunchBase<double, 3>* bunch, size_t& i);
+ void transformToHminus(PartBunchBase<double, 3>* bunch, size_t& i);
+ void transformToH3plus(PartBunchBase<double, 3>* bunch, size_t& i);
- bool computeEnergyLoss(Vector_t &/*P*/, const double /*deltat*/, bool /*includeFluctuations*/) const {
+ bool computeEnergyLoss(PartBunchBase<double, 3>* /*bunch*/,
+ Vector_t& /*P*/,
+ const double /*deltat*/,
+ bool /*includeFluctuations*/) const {
return false;
}
- Cyclotron *cycl_m;
- BeamStripping *bstp_m;
+ Cyclotron* cycl_m;
+ BeamStripping* bstp_m;
ElementBase::ElementType bstpshape_m;
- gsl_rng * r_m;
-
- double T_m; /// s
- double dT_m; /// s
+ gsl_rng* r_m;
+ double T_m; // s
+ double dT_m; // s
double mass_m;
double charge_m;
-
- std::string gas_m;
double pressure_m;
std::unique_ptr<LossDataSink> lossDs_m;
@@ -136,4 +137,4 @@ private:
static double b_m[3][9];
};
-#endif //BEAMSTRIPPINGPHYSICS_HH
\ No newline at end of file
+#endif //BEAMSTRIPPINGPHYSICS_HH
diff --git a/src/Classic/Solvers/CollimatorPhysics.cpp b/src/Classic/Solvers/CollimatorPhysics.cpp
index 009b98410387117f9241454f529b1379bca91074..59a43f190b7c2beb912ee410d4eaa582016dc159 100644
--- a/src/Classic/Solvers/CollimatorPhysics.cpp
+++ b/src/Classic/Solvers/CollimatorPhysics.cpp
@@ -1,9 +1,9 @@
//
// Class CollimatorPhysics
+// Defines the physical processes of beam scattering
+// and energy loss by heavy charged particles
//
-// Defines the collimator physics models
-//
-// Copyright (c) 2009 - 2020, Bi, Yang, Stachel, Adelmann
+// Copyright (c) 2009 - 2021, Bi, Yang, Stachel, Adelmann
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
@@ -17,7 +17,6 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-
#include "Solvers/CollimatorPhysics.hh"
#include "Physics/Physics.h"
#include "Physics/Material.h"
@@ -39,63 +38,65 @@
#include <gsl/gsl_randist.h>
+#include <algorithm>
#include <cmath>
-#include <iostream>
#include <fstream>
-#include <algorithm>
+#include <iostream>
#include <sys/time.h>
namespace {
struct DegraderInsideTester: public InsideTester {
- explicit DegraderInsideTester(ElementBase * el) {
+ explicit DegraderInsideTester(ElementBase* el) {
deg_m = static_cast<Degrader*>(el);
}
virtual
- bool checkHit(const Vector_t &R) override {
+ bool checkHit(const Vector_t& R) override {
return deg_m->isInMaterial(R(2));
}
private:
- Degrader *deg_m;
+ Degrader* deg_m;
};
struct CollimatorInsideTester: public InsideTester {
- explicit CollimatorInsideTester(ElementBase * el) {
+ explicit CollimatorInsideTester(ElementBase* el) {
col_m = static_cast<CCollimator*>(el);
}
virtual
- bool checkHit(const Vector_t &R) override {
+ bool checkHit(const Vector_t& R) override {
return col_m->checkPoint(R(0), R(1));
}
private:
- CCollimator *col_m;
+ CCollimator* col_m;
};
struct FlexCollimatorInsideTester: public InsideTester {
- explicit FlexCollimatorInsideTester(ElementBase * el) {
+ explicit FlexCollimatorInsideTester(ElementBase* el) {
col_m = static_cast<FlexibleCollimator*>(el);
}
virtual
- bool checkHit(const Vector_t &R) override {
+ bool checkHit(const Vector_t& R) override {
return col_m->isStopped(R);
}
private:
- FlexibleCollimator *col_m;
+ FlexibleCollimator* col_m;
};
}
-CollimatorPhysics::CollimatorPhysics(const std::string &name,
- ElementBase *element,
- std::string &material,
+CollimatorPhysics::CollimatorPhysics(const std::string& name,
+ ElementBase* element,
+ std::string& material,
bool enableRutherford,
double lowEnergyThr):
ParticleMatterInteractionHandler(name, element),
T_m(0.0),
dT_m(0.0),
+ mass_m(0.0),
+ charge_m(0.0),
material_m(material),
hitTester_m(nullptr),
Z_m(0),
@@ -103,6 +104,7 @@ CollimatorPhysics::CollimatorPhysics(const std::string &name,
rho_m(0.0),
X0_m(0.0),
I_m(0.0),
+ A1_c(0.0),
A2_c(0.0),
A3_c(0.0),
A4_c(0.0),
@@ -151,8 +153,8 @@ CollimatorPhysics::CollimatorPhysics(const std::string &name,
lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
- DegraderApplyTimer_m = IpplTimings::getTimer("DegraderApply");
- DegraderLoopTimer_m = IpplTimings::getTimer("DegraderLoop");
+ DegraderApplyTimer_m = IpplTimings::getTimer("DegraderApply");
+ DegraderLoopTimer_m = IpplTimings::getTimer("DegraderLoop");
DegraderDestroyTimer_m = IpplTimings::getTimer("DegraderDestroy");
}
@@ -177,14 +179,20 @@ void CollimatorPhysics::configureMaterialParameters() {
rho_m = material->getMassDensity();
X0_m = material->getRadiationLength();
I_m = material->getMeanExcitationEnergy();
+ A1_c = material->getStoppingPowerFitCoefficients(Physics::Material::A1);
A2_c = material->getStoppingPowerFitCoefficients(Physics::Material::A2);
A3_c = material->getStoppingPowerFitCoefficients(Physics::Material::A3);
A4_c = material->getStoppingPowerFitCoefficients(Physics::Material::A4);
A5_c = material->getStoppingPowerFitCoefficients(Physics::Material::A5);
+ B1_c = material->getStoppingPowerFitCoefficients(Physics::Material::B1);
+ B2_c = material->getStoppingPowerFitCoefficients(Physics::Material::B2);
+ B3_c = material->getStoppingPowerFitCoefficients(Physics::Material::B3);
+ B4_c = material->getStoppingPowerFitCoefficients(Physics::Material::B4);
+ B5_c = material->getStoppingPowerFitCoefficients(Physics::Material::B5);
}
-void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere) {
+void CollimatorPhysics::apply(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& boundingSphere) {
IpplTimings::startTimer(DegraderApplyTimer_m);
/*
@@ -208,6 +216,9 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
dT_m = bunch->getdT();
T_m = bunch->getT();
+ mass_m = bunch->getM();
+ charge_m = bunch->getQ();
+
bool onlyOneLoopOverParticles = ! (allParticleInMat_m);
do {
@@ -215,15 +226,13 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
push();
setTimeStepForLeavingParticles();
- /*
- if we are not looping copy newly arrived particles
- */
- if (onlyOneLoopOverParticles)
+ // if we are not looping copy newly arrived particles
+ if (onlyOneLoopOverParticles) {
copyFromBunch(bunch, boundingSphere);
-
+ }
addBackToBunch(bunch);
- computeInteraction();
+ computeInteraction(bunch);
push();
resetTimeStep();
@@ -244,100 +253,139 @@ void CollimatorPhysics::apply(PartBunchBase<double, 3> *bunch,
IpplTimings::stopTimer(DegraderApplyTimer_m);
}
-void CollimatorPhysics::computeInteraction() {
- /***
+void CollimatorPhysics::computeInteraction(PartBunchBase<double, 3>* bunch) {
+ /*
Do physics if
+ -- correct type of particle
-- particle not stopped (locParts_m[i].label != -1.0)
Absorbed particle i: locParts_m[i].label = -1.0;
*/
-
- for (size_t i = 0; i < locParts_m.size(); ++i) {
- if (locParts_m[i].label != -1) {
- Vector_t &R = locParts_m[i].Rincol;
- Vector_t &P = locParts_m[i].Pincol;
- double &dt = locParts_m[i].DTincol;
-
- if (hitTester_m->checkHit(R)) {
- bool pdead = computeEnergyLoss(P, dt);
- if (!pdead) {
- /*
- Now scatter and transport particle in material.
- The checkHit call just above will detect if the
- particle is rediffused from the material into vacuum.
- */
-
- computeCoulombScattering(R, P, dt);
- } else {
- // The particle is stopped in the material, set label to -1
- locParts_m[i].label = -1.0;
- ++ stoppedPartStat_m;
- lossDs_m->addParticle(R, P, -locParts_m[i].IDincol);
+ ParticleType pType = bunch->getPType();
+ if (pType == ParticleType::PROTON ||
+ pType == ParticleType::DEUTERON ||
+ pType == ParticleType::HMINUS ||
+ pType == ParticleType::MUON ||
+ pType == ParticleType::H2P ||
+ pType == ParticleType::ALPHA ) {
+
+ for (size_t i = 0; i < locParts_m.size(); ++i) {
+ if (locParts_m[i].label != -1) {
+ Vector_t& R = locParts_m[i].Rincol;
+ Vector_t& P = locParts_m[i].Pincol;
+ double& dt = locParts_m[i].DTincol;
+
+ if (hitTester_m->checkHit(R)) {
+ bool pdead = computeEnergyLoss(bunch, P, dt);
+ if (!pdead) {
+ /*
+ Now scatter and transport particle in material.
+ The checkHit call just above will detect if the
+ particle is rediffused from the material into vacuum.
+ */
+
+ computeCoulombScattering(R, P, dt);
+ } else {
+ // The particle is stopped in the material, set label to -1
+ locParts_m[i].label = -1.0;
+ ++ stoppedPartStat_m;
+ lossDs_m->addParticle(R, P, -locParts_m[i].IDincol);
+ }
}
}
}
+ } else {
+ throw GeneralClassicException(
+ "CollimatorPhysics::computeInteraction",
+ "Particle " + bunch->getPTypeString() +
+ " is not supported for scattering interactions!");
}
-
- /*
- delete absorbed particles
- */
+ // delete absorbed particles
deleteParticleFromLocalVector();
}
-/// Energy Loss: using the Bethe-Bloch equation.
-/// Energy straggling: For relatively thick absorbers such that the number of collisions is large,
-/// the energy loss distribution is shown to be Gaussian in form.
+/// Energy Loss: using the Bethe-Bloch equation.
+/// In low-energy region use Andersen-Ziegler fitting (only for protons)
+/// Energy straggling: For relatively thick absorbers such that the number of collisions
+/// is large, the energy loss distribution is shown to be Gaussian in form.
/// See Particle Physics Booklet, chapter 'Passage of particles through matter' or
/// Review of Particle Physics, DOI: 10.1103/PhysRevD.86.010001, page 329 ff
// -------------------------------------------------------------------------
-bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
+bool CollimatorPhysics::computeEnergyLoss(PartBunchBase<double, 3>* bunch,
+ Vector_t& P,
const double deltat,
bool includeFluctuations) const {
+ ParticleType pType = bunch->getPType();
constexpr double m2cm = 1e2;
+ constexpr double eV2keV = 1e-3;
constexpr double GeV2keV = 1e6;
+
constexpr double massElectron_keV = Physics::m_e * GeV2keV;
- constexpr double massProton_keV = Physics::m_p * GeV2keV;
- constexpr double massProton_amu = Physics::m_p / Physics::amu;
- constexpr double chargeProton = Physics::z_p;
+
+ const double mass_keV = mass_m * eV2keV;
+
constexpr double K = 4.0 * Physics::pi * Physics::Avo * Physics::r_e * m2cm * Physics::r_e * m2cm * massElectron_keV;
double gamma = Util::getGamma(P);
const double gammaSqr = std::pow(gamma, 2);
const double betaSqr = 1.0 - 1.0 / gammaSqr;
double beta = std::sqrt(betaSqr);
- double Ekin = (gamma - 1) * massProton_keV;
+ double Ekin = (gamma - 1) * mass_keV;
double dEdx = 0.0;
+ double epsilon = 0.0;
const double deltas = deltat * beta * Physics::c;
const double deltasrho = deltas * m2cm * rho_m;
- if (Ekin >= 600) {
- constexpr double eV2keV = 1e-3;
- constexpr double massRatio = Physics::m_e / Physics::m_p;
- double Tmax = (2.0 * massElectron_keV * betaSqr * gammaSqr /
- (std::pow(gamma + massRatio, 2) - (gammaSqr - 1.0)));
-
- dEdx = (-K * std::pow(chargeProton, 2) * Z_m / (A_m * betaSqr) *
- (0.5 * std::log(2 * massElectron_keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
-
- Ekin += deltasrho * dEdx;
- } else if (Ekin > 10) {
- const double Ts = Ekin / massProton_amu;
- const double epsilon_low = A2_c * pow(Ts, 0.45);
- const double epsilon_high = (A3_c / Ts) * log(1 + (A4_c / Ts) + (A5_c * Ts));
- const double epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
- dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
-
- Ekin += deltasrho * dEdx;
+ const double massRatio = massElectron_keV / mass_keV;
+ double Tmax = (2.0 * massElectron_keV * betaSqr * gammaSqr /
+ (std::pow(gamma + massRatio, 2) - (gammaSqr - 1.0)));
+
+ if (pType != ParticleType::ALPHA) {
+
+ if (Ekin >= 600 && Ekin < 1e7) {
+ dEdx = (-K * std::pow(charge_m, 2) * Z_m / (A_m * betaSqr) *
+ (0.5 * std::log(2 * massElectron_keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
+ } else if (pType == ParticleType::PROTON && Ekin < 600) {
+ constexpr double massProton_amu = Physics::m_p / Physics::amu;
+ const double Ts = Ekin / massProton_amu;
+ if (Ekin > 10) {
+ const double epsilon_low = A2_c * std::pow(Ts, 0.45);
+ const double epsilon_high = (A3_c / Ts) * std::log(1 + (A4_c / Ts) + (A5_c * Ts));
+ epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
+ } else if (Ekin > 1) {
+ epsilon = A1_c * std::pow(Ts, 0.5);
+ }
+ dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
+ } else {
+ INFOMSG(level4 << "Particle energy out of the valid range "
+ "for energy loss calculation." << endl);
+ }
+ } else {
+ if (Ekin > 10000 && Ekin < 1e6) {
+ dEdx = (-K * std::pow(charge_m, 2) * Z_m / (A_m * betaSqr) *
+ (0.5 * std::log(2 * massElectron_keV * betaSqr * gammaSqr * Tmax / std::pow(I_m * eV2keV, 2)) - betaSqr));
+ } else if (Ekin > 1 && Ekin <= 10000) {
+ const double T = Ekin * 1e-3; //MeV
+ const double epsilon_low = B1_c * std::pow(1000*T, B2_c);
+ const double epsilon_high = (B3_c / T) * std::log(1 + (B4_c / T) + (B5_c * T));
+ epsilon = (epsilon_low * epsilon_high) / (epsilon_low + epsilon_high);
+ dEdx = -epsilon / (1e18 * (A_m / Physics::Avo));
+ } else {
+ INFOMSG(level4 << "Particle energy out of the valid range "
+ "for energy loss calculation." << endl);
+ }
}
+ Ekin += deltasrho * dEdx;
+
if (includeFluctuations) {
double sigma_E = std::sqrt(K * massElectron_keV * rho_m * (Z_m / A_m) * deltas * m2cm);
Ekin += gsl_ran_gaussian(rGen_m, sigma_E);
}
- gamma = Ekin / massProton_keV + 1.0;
+ gamma = Ekin / mass_keV + 1.0;
beta = std::sqrt(1.0 - 1.0 / std::pow(gamma, 2));
P = gamma * beta * P / euclidean_norm(P);
@@ -347,38 +395,38 @@ bool CollimatorPhysics::computeEnergyLoss(Vector_t &P,
// Implement the rotation in 2 dimensions here
-// For details see:
-// J. Beringer et al. (Particle Data Group), "Passage of particles through matter"
-// Phys. Rev. D 86, 010001 (2012),
-void CollimatorPhysics::applyRotation(Vector_t &P,
- Vector_t &R,
+// For details see: J. Beringer et al. (Particle Data Group),
+// "Passage of particles through matter", Phys. Rev. D 86, 010001 (2012)
+// -------------------------------------------------------------------------
+void CollimatorPhysics::applyRotation(Vector_t& P,
+ Vector_t& R,
double shift,
double thetacou) {
// Calculate the angle between the transverse and longitudinal component of the momentum
double Psixz = std::fmod(std::atan2(P(0), P(2)) + Physics::two_pi, Physics::two_pi);
- R(0) = R(0) + shift * cos(Psixz);
- R(2) = R(2) - shift * sin(Psixz);
+ R(0) = R(0) + shift * std::cos(Psixz);
+ R(2) = R(2) - shift * std::sin(Psixz);
// Apply the rotation about the random angle thetacou
double Px = P(0);
- P(0) = Px * cos(thetacou) + P(2) * sin(thetacou);
- P(2) = -Px * sin(thetacou) + P(2) * cos(thetacou);
+ P(0) = Px * std::cos(thetacou) + P(2) * std::sin(thetacou);
+ P(2) = -Px * std::sin(thetacou) + P(2) * std::cos(thetacou);
}
-void CollimatorPhysics::applyRandomRotation(Vector_t &P, double theta0) {
+void CollimatorPhysics::applyRandomRotation(Vector_t& P, double theta0) {
double thetaru = 2.5 / std::sqrt(gsl_rng_uniform(rGen_m)) * 2.0 * theta0;
double phiru = Physics::two_pi * gsl_rng_uniform(rGen_m);
double normPtrans = std::sqrt(P(0) * P(0) + P(1) * P(1));
double Theta = std::atan(normPtrans / std::abs(P(2)));
- double CosT = cos(Theta);
- double SinT = sin(Theta);
+ double CosT = std::cos(Theta);
+ double SinT = std::sin(Theta);
- Vector_t X(cos(phiru)*sin(thetaru),
- sin(phiru)*sin(thetaru),
- cos(thetaru));
+ Vector_t X(std::cos(phiru)*std::sin(thetaru),
+ std::sin(phiru)*std::sin(thetaru),
+ std::cos(thetaru));
X *= euclidean_norm(P);
Vector_t W(-P(1), P(0), 0.0);
@@ -391,21 +439,18 @@ void CollimatorPhysics::applyRandomRotation(Vector_t &P, double theta0) {
/// Coulomb Scattering: Including Multiple Coulomb Scattering and large angle Rutherford Scattering.
/// Using the distribution given in Classical Electrodynamics, by J. D. Jackson.
//--------------------------------------------------------------------------
-void CollimatorPhysics::computeCoulombScattering(Vector_t &R,
- Vector_t &P,
+void CollimatorPhysics::computeCoulombScattering(Vector_t& R,
+ Vector_t& P,
double dt) {
- constexpr double GeV2eV = 1e9;
- constexpr double massProton_eV = Physics::m_p * GeV2eV;
- constexpr double chargeProton = Physics::z_p;
constexpr double sqrtThreeInv = 0.57735026918962576451; // sqrt(1.0 / 3.0)
const double normP = euclidean_norm(P);
const double gammaSqr = std::pow(normP, 2) + 1.0;
const double beta = std::sqrt(1.0 - 1.0 / gammaSqr);
const double deltas = dt * beta * Physics::c;
- const double theta0 = (13.6e6 / (beta * normP * massProton_eV) *
- chargeProton * std::sqrt(deltas / X0_m) *
- (1.0 + 0.038 * log(deltas / X0_m)));
+ const double theta0 = (13.6e6 / (beta * normP * mass_m) *
+ charge_m * std::sqrt(deltas / X0_m) *
+ (1.0 + 0.038 * std::log(deltas / X0_m)));
double phi = Physics::two_pi * gsl_rng_uniform(rGen_m);
for (unsigned int i = 0; i < 2; ++ i) {
@@ -438,7 +483,7 @@ void CollimatorPhysics::computeCoulombScattering(Vector_t &R,
}
}
-void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
+void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3>* bunch) {
const size_t nL = locParts_m.size();
if (nL == 0) return;
@@ -447,7 +492,7 @@ void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
const double elementLength = element_ref_m->getElementLength();
for (size_t i = 0; i < nL; ++ i) {
- Vector_t &R = locParts_m[i].Rincol;
+ Vector_t& R = locParts_m[i].Rincol;
if (R[2] >= elementLength) {
@@ -457,7 +502,6 @@ void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
Binincol is still <0, but now the particle is rediffused
from the material and hence this is not a "lost" particle anymore
*/
-
bunch->Bin[numLocalParticles] = 1;
bunch->R[numLocalParticles] = R;
bunch->P[numLocalParticles] = locParts_m[i].Pincol;
@@ -478,14 +522,12 @@ void CollimatorPhysics::addBackToBunch(PartBunchBase<double, 3> *bunch) {
}
}
- /*
- delete particles that went to the bunch
- */
+ // delete particles that went to the bunch
deleteParticleFromLocalVector();
}
-void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere)
+void CollimatorPhysics::copyFromBunch(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& boundingSphere)
{
const size_t nL = bunch->getLocalNum();
if (nL == 0) return;
@@ -555,8 +597,8 @@ void CollimatorPhysics::print(Inform &msg) {
Inform::FmtFlags_t ff = msg.flags();
if (totalPartsInMat_m > 0 ||
- bunchToMatStat_m > 0 ||
- rediffusedStat_m > 0 ||
+ bunchToMatStat_m > 0 ||
+ rediffusedStat_m > 0 ||
stoppedPartStat_m > 0) {
OPALTimer::Timer time;
@@ -586,8 +628,8 @@ namespace {
void CollimatorPhysics::deleteParticleFromLocalVector() {
/*
- the particle to be deleted (label < 0) are all at the end of
- the vector.
+ the particle to be deleted (label < 0) are all
+ at the end of the vector.
*/
sort(locParts_m.begin(), locParts_m.end(), myCompF);
@@ -611,8 +653,8 @@ void CollimatorPhysics::deleteParticleFromLocalVector() {
void CollimatorPhysics::push() {
for (size_t i = 0; i < locParts_m.size(); ++i) {
- Vector_t &R = locParts_m[i].Rincol;
- Vector_t &P = locParts_m[i].Pincol;
+ Vector_t& R = locParts_m[i].Rincol;
+ Vector_t& P = locParts_m[i].Pincol;
double gamma = Util::getGamma(P);
R += 0.5 * dT_m * Physics::c * P / gamma;
@@ -629,9 +671,9 @@ void CollimatorPhysics::setTimeStepForLeavingParticles() {
const double elementLength = element_ref_m->getElementLength();
for (size_t i = 0; i < locParts_m.size(); ++i) {
- Vector_t &R = locParts_m[i].Rincol;
- Vector_t &P = locParts_m[i].Pincol;
- double &dt = locParts_m[i].DTincol;
+ Vector_t& R = locParts_m[i].Rincol;
+ Vector_t& P = locParts_m[i].Pincol;
+ double& dt = locParts_m[i].DTincol;
double gamma = Util::getGamma(P);
Vector_t stepLength = dT_m * Physics::c * P / gamma;
@@ -651,7 +693,7 @@ void CollimatorPhysics::setTimeStepForLeavingParticles() {
void CollimatorPhysics::resetTimeStep() {
for (size_t i = 0; i < locParts_m.size(); ++i) {
- double &dt = locParts_m[i].DTincol;
+ double& dt = locParts_m[i].DTincol;
dt = dT_m;
}
diff --git a/src/Classic/Solvers/CollimatorPhysics.hh b/src/Classic/Solvers/CollimatorPhysics.hh
index c4a8fa73dbd761f3c5f69d1d00061d3f896a3b82..db5da27a595a1febae2d5e908e1ef5749ab36ae3 100644
--- a/src/Classic/Solvers/CollimatorPhysics.hh
+++ b/src/Classic/Solvers/CollimatorPhysics.hh
@@ -1,9 +1,9 @@
//
// Class CollimatorPhysics
+// Defines the physical processes of beam scattering
+// and energy loss by heavy charged particles
//
-// Defines the collimator physics models
-//
-// Copyright (c) 2009 - 2020, Bi, Yang, Stachel, Adelmann
+// Copyright (c) 2009 - 2021, Bi, Yang, Stachel, Adelmann
// Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved.
//
@@ -17,7 +17,6 @@
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
-
#ifndef COLLIMATORPHYSICS_HH
#define COLLIMATORPHYSICS_HH
@@ -25,13 +24,14 @@
#include "AbsBeamline/ElementBase.h"
#include "Algorithms/Vektor.h"
+
#include <gsl/gsl_rng.h>
#include "Utility/IpplTimings.h"
#include <memory>
-#include <utility>
#include <string>
+#include <utility>
#include <vector>
template <class T, unsigned Dim>
@@ -64,15 +64,15 @@ struct InsideTester {
class CollimatorPhysics: public ParticleMatterInteractionHandler {
public:
- CollimatorPhysics(const std::string &name,
- ElementBase *element,
- std::string &mat,
+ CollimatorPhysics(const std::string& name,
+ ElementBase* element,
+ std::string& mat,
bool enableRutherford,
double lowEnergyThr);
~CollimatorPhysics();
- virtual void apply(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere);
+ virtual void apply(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& boundingSphere);
virtual const std::string getType() const;
@@ -84,28 +84,29 @@ public:
virtual size_t getParticlesInMat();
virtual unsigned getRediffused();
virtual unsigned int getNumEntered();
- void computeInteraction();
+ void computeInteraction(PartBunchBase<double, 3>* bunch);
- virtual bool computeEnergyLoss(Vector_t &P,
+ virtual bool computeEnergyLoss(PartBunchBase<double, 3>* bunch,
+ Vector_t& P,
const double deltat,
bool includeFluctuations = true) const;
private:
void configureMaterialParameters();
- void computeCoulombScattering(Vector_t &R,
- Vector_t &P,
+ void computeCoulombScattering(Vector_t& R,
+ Vector_t& P,
double dt);
- void applyRotation(Vector_t &P,
- Vector_t &R,
+ void applyRotation(Vector_t& P,
+ Vector_t& R,
double xplane,
double thetacou);
- void applyRandomRotation(Vector_t &P, double theta0);
+ void applyRandomRotation(Vector_t& P, double theta0);
- void copyFromBunch(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere);
- void addBackToBunch(PartBunchBase<double, 3> *bunch);
+ void copyFromBunch(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& boundingSphere);
+ void addBackToBunch(PartBunchBase<double, 3>* bunch);
void deleteParticleFromLocalVector();
@@ -119,6 +120,9 @@ private:
double T_m; // own time, maybe larger than in the bunch object
double dT_m; // dt from bunch
+ double mass_m; // mass from bunch (eV)
+ double charge_m; // charge from bunch (elementary charges)
+
gsl_rng *rGen_m; // random number generator
std::string material_m; // type of material e.g. aluminum
@@ -133,10 +137,21 @@ private:
double X0_m; // the radiation length in [m]
double I_m; // the mean excitation energy [eV]
- double A2_c; // coefficients to fit model to measurement data
- double A3_c; // see e.g. page 16-20 in H. H. Andersen, J. F. Ziegler,
- double A4_c; // "Hydrogen Stopping Powers and Ranges in All Elements",
- double A5_c; // Pergamon Press, 1977
+ /*
+ coefficients to fit model to measurement data according to Andersen-Ziegler formulae.
+ see ICRU-49, "Stopping Powers and Ranges for Protons and Alpha Particles",
+ chapter 'Electronic (Collision) Stopping Powers in the Low-Energy Region'
+ */
+ double A1_c;
+ double A2_c;
+ double A3_c;
+ double A4_c;
+ double A5_c;
+ double B1_c;
+ double B2_c;
+ double B3_c;
+ double B4_c;
+ double B5_c;
// number of particles that enter the material in current step (count for single step)
unsigned int bunchToMatStat_m;
diff --git a/src/Classic/Solvers/ParticleMatterInteractionHandler.hh b/src/Classic/Solvers/ParticleMatterInteractionHandler.hh
index af4fc75ad42382060a326fa9746be44156b08d85..ed81c3673fed42642389ffbce508bbf25930c53e 100644
--- a/src/Classic/Solvers/ParticleMatterInteractionHandler.hh
+++ b/src/Classic/Solvers/ParticleMatterInteractionHandler.hh
@@ -1,3 +1,20 @@
+//
+// Class ParticleMatterInteractionHandler
+// Defines the handler for particle-matter interactions in the elements
+//
+// 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/>.
+//
#ifndef PARTICLEMATTERINTERACTIONHANDLER_HH
#define PARTICLEMATTERINTERACTIONHANDLER_HH
@@ -11,11 +28,12 @@ class PartBunchBase;
class Inform;
class ParticleMatterInteractionHandler {
+
public:
- ParticleMatterInteractionHandler(std::string name, ElementBase *elref);
+ ParticleMatterInteractionHandler(std::string name, ElementBase* elref);
virtual ~ParticleMatterInteractionHandler() { };
- virtual void apply(PartBunchBase<double, 3> *bunch,
- const std::pair<Vector_t, double> &boundingSphere) = 0;
+ virtual void apply(PartBunchBase<double, 3>* bunch,
+ const std::pair<Vector_t, double>& boundingSphere) = 0;
virtual const std::string getType() const = 0;
virtual void print(Inform& os) = 0;
virtual bool stillActive() = 0;
@@ -26,29 +44,30 @@ public:
virtual unsigned int getNumEntered() = 0;
void setFlagAllParticlesIn(bool p);
bool getFlagAllParticlesIn() const;
- void updateElement(ElementBase *newref);
+ void updateElement(ElementBase* newref);
ElementBase* getElement();
- virtual bool computeEnergyLoss(Vector_t &P,//double &Eng,
+ virtual bool computeEnergyLoss(PartBunchBase<double, 3>* bunch,
+ Vector_t& P,
const double deltat,
bool includeFluctuations = true) const = 0;
protected:
- ElementBase *element_ref_m;
+ ElementBase* element_ref_m;
bool allParticleInMat_m; ///< if all particles are in matter stay inside the particle matter interaction
+
private:
const std::string name_m;
-
};
inline
-ParticleMatterInteractionHandler::ParticleMatterInteractionHandler(std::string name, ElementBase *elref):
+ParticleMatterInteractionHandler::ParticleMatterInteractionHandler(std::string name, ElementBase* elref):
element_ref_m(elref),
allParticleInMat_m(false),
name_m(name)
{}
inline
-void ParticleMatterInteractionHandler::updateElement(ElementBase *newref) {
+void ParticleMatterInteractionHandler::updateElement(ElementBase* newref) {
element_ref_m = newref;
}
@@ -66,4 +85,4 @@ inline
bool ParticleMatterInteractionHandler::getFlagAllParticlesIn() const {
return allParticleInMat_m;
}
-#endif // PARTICLEMATTERINTERACTION_HH
\ No newline at end of file
+#endif // PARTICLEMATTERINTERACTION_HH
diff --git a/src/Distribution/Distribution.cpp b/src/Distribution/Distribution.cpp
index 6a720620fd4c5ca9afadea33210e4fee05945a17..d98c1abaf9b580a38bf69dc9bfdd3db883a1606a 100644
--- a/src/Distribution/Distribution.cpp
+++ b/src/Distribution/Distribution.cpp
@@ -1635,7 +1635,8 @@ size_t Distribution::emitParticles(PartBunchBase<double, 3> *beam, double eZ) {
beam->Q[numberOfEmittedParticles] = beam->getChargePerParticle();
beam->Ef[numberOfEmittedParticles] = Vector_t(0.0);
beam->Bf[numberOfEmittedParticles] = Vector_t(0.0);
- beam->PType[numberOfEmittedParticles] = ParticleType::REGULAR;
+ beam->PType[numberOfEmittedParticles] = beam->getPType();
+ beam->POrigin[numberOfEmittedParticles] = ParticleOrigin::REGULAR;
beam->TriID[numberOfEmittedParticles] = 0;
numberOfEmittedParticles++;
@@ -2751,7 +2752,8 @@ void Distribution::injectBeam(PartBunchBase<double, 3> *beam) {
beam->Q[partIndex] = beam->getChargePerParticle();
beam->Ef[partIndex] = Vector_t(0.0);
beam->Bf[partIndex] = Vector_t(0.0);
- beam->PType[partIndex] = ParticleType::REGULAR;
+ beam->PType[partIndex] = beam->getPType();
+ beam->POrigin[partIndex] = ParticleOrigin::REGULAR;
beam->TriID[partIndex] = 0;
if (numberOfEnergyBins_m > 0) {
@@ -4451,4 +4453,4 @@ void Distribution::adjustPhaseSpace(double massIneV) {
tOrZDist_m[i] -= avrg[4];
pzDist_m[i] += eps;
}
-}
\ No newline at end of file
+}
diff --git a/src/Structure/Beam.cpp b/src/Structure/Beam.cpp
index deb647b362b6ba361f154e03ee04f50fff5bebf3..8f4b70171acdadf6e7bd3538605c67700b879aca 100644
--- a/src/Structure/Beam.cpp
+++ b/src/Structure/Beam.cpp
@@ -1,22 +1,24 @@
-// ------------------------------------------------------------------------
-// $RCSfile: Beam.cpp,v $
-// ------------------------------------------------------------------------
-// $Revision: 1.3.4.1 $
-// ------------------------------------------------------------------------
-// Copyright: see Copyright.readme
-// ------------------------------------------------------------------------
//
-// Class: Beam
+// 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
//
-// $Date: 2003/08/11 22:09:00 $
-// $Author: dbruhwil $
+// 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"
@@ -30,72 +32,62 @@
using namespace Expressions;
-// Class Beam
-// ------------------------------------------------------------------------
// The attributes of class Beam.
namespace {
enum {
- // DESCRIPTION OF SINGLE PARTICLE:
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
-
- // BEAM CURRENT AND EMITTANCES:
BCURRENT, // Beam current in A
-
- // BEAM FREQUENCY
- BFREQ, // Beam frequency in MHz
-
- // DESCRIPTION OF BUNCHES:
+ BFREQ, // Beam frequency in MHz
NPART, // Number of particles per bunch
SIZE
};
}
-
const double Beam::energy_scale = 1.0e9;
-
Beam::Beam():
Definition(SIZE, "BEAM",
"The \"BEAM\" statement defines data for the particles "
"in a beam."),
reference(1.0, Physics::m_p *energy_scale, 1.0 * energy_scale) {
- // DESCRIPTION OF SINGLE PARTICLE:
itsAttr[PARTICLE] = Attributes::makeString
("PARTICLE", "Name of particle to be used");
- itsAttr[MASS] = Attributes::makeReal
- ("MASS", "Particle rest mass in GeV");
- itsAttr[CHARGE] = Attributes::makeReal
- ("CHARGE", "Particle charge in proton charges");
- itsAttr[ENERGY] = Attributes::makeReal
- ("ENERGY", "Particle energy in GeV");
- itsAttr[PC] = Attributes::makeReal
- ("PC", "Particle momentum in GeV/c");
+
+ 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]");
PtrToScalar<double> expr = new SRefExpr<double>("P0", "");
itsAttr[PC].set(new SAutomatic<double>(expr));
- itsAttr[GAMMA] = Attributes::makeReal
- ("GAMMA", "ENERGY / MASS");
- // BEAM CURRENT AND EMITTANCES:
+ itsAttr[GAMMA] = Attributes::makeReal
+ ("GAMMA", "ENERGY / MASS");
+
itsAttr[BCURRENT] = Attributes::makeReal
- ("BCURRENT", "Beam current in A (all bunches)");
+ ("BCURRENT", "Beam current [A] (all bunches)");
- // BEAM FREQUENCY
- itsAttr[BFREQ] = Attributes::makeReal
- ("BFREQ", "Beam frequency in MHz (all bunches)");
+ itsAttr[BFREQ] = Attributes::makeReal
+ ("BFREQ", "Beam frequency [MHz] (all bunches)");
- // DESCRIPTION OF BUNCHES:
- itsAttr[NPART] = Attributes::makeReal
- ("NPART", "Number of particles in bunch");
+ itsAttr[NPART] = Attributes::makeReal
+ ("NPART", "Number of particles in bunch");
// Set up default beam.
- Beam *defBeam = clone("UNNAMED_BEAM");
+ Beam* defBeam = clone("UNNAMED_BEAM");
defBeam->builtin = true;
try {
@@ -109,7 +101,7 @@ Beam::Beam():
}
-Beam::Beam(const std::string &name, Beam *parent):
+Beam::Beam(const std::string& name, Beam* parent):
Definition(name, parent),
reference(parent->reference)
{}
@@ -119,13 +111,13 @@ Beam::~Beam()
{}
-bool Beam::canReplaceBy(Object *object) {
+bool Beam::canReplaceBy(Object* object) {
// Can replace only by another BEAM.
- return dynamic_cast<Beam *>(object) != 0;
+ return dynamic_cast<Beam*>(object) != 0;
}
-Beam *Beam::clone(const std::string &name) {
+Beam* Beam::clone(const std::string& name) {
return new Beam(name, this);
}
@@ -135,10 +127,10 @@ void Beam::execute() {
}
-Beam *Beam::find(const std::string &name) {
- Beam *beam = dynamic_cast<Beam *>(OpalData::getInstance()->find(name));
+Beam* Beam::find(const std::string& name) {
+ Beam* beam = dynamic_cast<Beam*>(OpalData::getInstance()->find(name));
- if(beam == 0) {
+ if (beam == 0) {
throw OpalException("Beam::find()", "Beam \"" + name + "\" not found.");
}
@@ -149,9 +141,9 @@ size_t Beam::getNumberOfParticles() const {
return (size_t)Attributes::getReal(itsAttr[NPART]);
}
-const PartData &Beam::getReference() const {
+const PartData& Beam::getReference() const {
// Cast away const, to allow logically constant Beam to update.
- const_cast<Beam *>(this)->update();
+ const_cast<Beam*>(this)->update();
return reference;
}
@@ -187,9 +179,20 @@ double Beam::getMassPerParticle() const {
void Beam::update() {
// Find the particle name.
- if(itsAttr[PARTICLE]) {
+ if (itsAttr[PARTICLE]) {
static const char *names[] = {
- "ELECTRON", "PROTON", "POSITRON", "ANTIPROTON", "CARBON", "HMINUS", "URANIUM", "MUON", "DEUTERON", "XENON", "H2P"
+ "ELECTRON",
+ "PROTON",
+ "POSITRON",
+ "ANTIPROTON",
+ "CARBON",
+ "HMINUS",
+ "URANIUM",
+ "MUON",
+ "DEUTERON",
+ "XENON",
+ "H2P",
+ "ALPHA"
};
static const double masses[] = {
@@ -203,17 +206,18 @@ void Beam::update() {
Physics::m_mu,
Physics::m_d,
Physics::m_xe,
- Physics::m_h2p
+ Physics::m_h2p,
+ Physics::m_alpha
};
static const double charges[] = {
- -1.0, 1.0, 1.0, -1.0, 12.0, -1.0, 35.0, -1.0, 1.0, 20.0, 1.0
+ -1.0, 1.0, 1.0, -1.0, 12.0, -1.0, 35.0, -1.0, 1.0, 20.0, 1.0, 2.0
};
const unsigned int numParticleNames = std::end(names) - std::begin(names);
std::string pName = Attributes::getString(itsAttr[PARTICLE]);
- for(unsigned int i = 0; i < numParticleNames; ++ i) {
- if(pName == names[i]) {
+ for (unsigned int i = 0; i < numParticleNames; ++ i) {
+ if (pName == names[i]) {
Attributes::setReal(itsAttr[MASS], masses[i]);
Attributes::setReal(itsAttr[CHARGE], charges[i]);
break;
@@ -228,25 +232,25 @@ void Beam::update() {
reference = PartData(charge, mass, 1.0);
// Checks
- if(itsAttr[GAMMA]) {
+ if (itsAttr[GAMMA]) {
double gamma = Attributes::getReal(itsAttr[GAMMA]);
- if(gamma > 1.0) {
+ if (gamma > 1.0) {
reference.setGamma(gamma);
} else {
throw OpalException("Beam::update()",
"\"GAMMA\" should be greater than 1.");
}
- } else if(itsAttr[ENERGY]) {
+ } else if (itsAttr[ENERGY]) {
double energy = Attributes::getReal(itsAttr[ENERGY]) * energy_scale;
- if(energy > reference.getM()) {
+ if (energy > reference.getM()) {
reference.setE(energy);
} else {
throw OpalException("Beam::update()",
"\"ENERGY\" should be greater than \"MASS\".");
}
- } else if(itsAttr[PC]) {
+ } else if (itsAttr[PC]) {
double pc = Attributes::getReal(itsAttr[PC]) * energy_scale;
- if(pc > 0.0) {
+ if (pc > 0.0) {
reference.setP(pc);
} else {
throw OpalException("Beam::update()",
@@ -255,7 +259,7 @@ void Beam::update() {
}
// Set default name.
- if(getOpalName().empty()) setOpalName("UNNAMED_BEAM");
+ if (getOpalName().empty()) setOpalName("UNNAMED_BEAM");
}
@@ -270,7 +274,7 @@ double Beam::getPC() const { //obtain value for PC
}
-void Beam::print(std::ostream &os) const {
+void Beam::print(std::ostream& os) const {
double charge = Attributes::getReal(itsAttr[CHARGE]);
os << "* ************* B E A M ************************************************************ " << std::endl;
os << "* BEAM " << getOpalName() << '\n'
diff --git a/src/Structure/Beam.h b/src/Structure/Beam.h
index 5f221a50c236208451f6fe61b85e80ffe4a1825c..a5b26b2c2dd35e34e7f44091383184522114866a 100644
--- a/src/Structure/Beam.h
+++ b/src/Structure/Beam.h
@@ -1,22 +1,24 @@
-#ifndef OPAL_Beam_HH
-#define OPAL_Beam_HH
-
-// ------------------------------------------------------------------------
-// $RCSfile: Beam.h,v $
-// ------------------------------------------------------------------------
-// $Revision: 1.1.1.1 $
-// ------------------------------------------------------------------------
-// Copyright: see Copyright.readme
-// ------------------------------------------------------------------------
//
-// Class: Beam
+// 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.
//
-// $Date: 2000/03/27 09:33:44 $
-// $Author: Andreas Adelmann $
+// 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 OPAL_Beam_HH
+#define OPAL_Beam_HH
#include "AbstractObjects/Definition.h"
#include "Algorithms/PartData.h"
@@ -26,12 +28,6 @@
class Inform;
-// Class Beam
-// ------------------------------------------------------------------------
-/// The BEAM definition.
-// A BEAM definition is used by most physics commands to define the
-// particle charge and the reference momentum, together with some other
-// data.
class Beam: public Definition {
@@ -44,16 +40,16 @@ public:
/// Test if replacement is allowed.
// Can replace only by another BEAM.
- virtual bool canReplaceBy(Object *object);
+ virtual bool canReplaceBy(Object* object);
/// Make clone.
- virtual Beam *clone(const std::string &name);
+ virtual Beam* clone(const std::string& name);
/// Check the BEAM data.
virtual void execute();
/// Find named BEAM.
- static Beam *find(const std::string &name);
+ static Beam* find(const std::string& name);
//ff => get gamma value
double getGamma() const;
@@ -65,7 +61,7 @@ public:
size_t getNumberOfParticles() const;
/// Return the embedded CLASSIC PartData.
- const PartData &getReference() const;
+ const PartData& getReference() const;
/// Return the beam current in A
double getCurrent() const;
@@ -91,16 +87,16 @@ public:
/// Update the BEAM data.
virtual void update();
- void print(std::ostream &os) const;
+ void print(std::ostream& os) const;
private:
// Not implemented.
- Beam(const Beam &);
- void operator=(const Beam &);
+ Beam(const Beam&);
+ void operator=(const Beam&);
// Clone constructor.
- Beam(const std::string &name, Beam *parent);
+ Beam(const std::string& name, Beam* parent);
// The particle reference data.
PartData reference;
@@ -109,10 +105,9 @@ private:
static const double energy_scale;
};
-inline std::ostream &operator<<(std::ostream &os, const Beam &b) {
+inline std::ostream &operator<<(std::ostream& os, const Beam& b) {
b.print(os);
return os;
}
-
#endif // OPAL_Beam_HH
diff --git a/src/Structure/H5PartWrapperForPT.cpp b/src/Structure/H5PartWrapperForPT.cpp
index 7215ce182c086856d14e1619615a92cd84f9191a..55756fbb41ef1f17225b6dac595773a7dd5c49b4 100644
--- a/src/Structure/H5PartWrapperForPT.cpp
+++ b/src/Structure/H5PartWrapperForPT.cpp
@@ -1,7 +1,20 @@
//
-// Copyright & License: See Copyright.readme in src directory
+// Class H5PartWrapperForPT
+// A class that manages all calls to H5Part for the Parallel-T tracker.
+//
+// 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/H5PartWrapperForPT.h"
#include "OPALconfig.h"
@@ -14,14 +27,15 @@
#include "h5core/h5_types.h"
+#include <cmath>
#include <sstream>
#include <set>
-H5PartWrapperForPT::H5PartWrapperForPT(const std::string &fileName, h5_int32_t flags):
+H5PartWrapperForPT::H5PartWrapperForPT(const std::string& fileName, h5_int32_t flags):
H5PartWrapper(fileName, flags)
{ }
-H5PartWrapperForPT::H5PartWrapperForPT(const std::string &fileName, int restartStep, std::string sourceFile, h5_int32_t flags):
+H5PartWrapperForPT::H5PartWrapperForPT(const std::string& fileName, int restartStep, std::string sourceFile, h5_int32_t flags):
H5PartWrapper(fileName, restartStep, sourceFile, flags)
{
if (restartStep == -1) {
@@ -73,7 +87,7 @@ void H5PartWrapperForPT::readHeader() {
H5ReadFileAttribInt64(file_m, "nAutoPhaseCavities", &numAutoPhaseCavities) != H5_SUCCESS) {
numAutoPhaseCavities = 0;
} else {
- for(long i = 0; i < numAutoPhaseCavities; ++ i) {
+ for (long i = 0; i < numAutoPhaseCavities; ++ i) {
std::string elementName = "Cav-" + std::to_string(i + 1) + "-name";
std::string elementPhase = "Cav-" + std::to_string(i + 1) + "-value";
@@ -123,9 +137,9 @@ void H5PartWrapperForPT::readStepHeader(PartBunchBase<double, 3>* bunch) {
Vector_t TaitBryant;
READSTEPATTRIB(Float64, file_m, "TaitBryantAngles", (h5_float64_t *)&TaitBryant);
- Quaternion rotTheta(cos(0.5 * TaitBryant[0]), 0, sin(0.5 * TaitBryant[0]), 0);
- Quaternion rotPhi(cos(0.5 * TaitBryant[1]), sin(0.5 * TaitBryant[1]), 0, 0);
- Quaternion rotPsi(cos(0.5 * TaitBryant[2]), 0, 0, sin(0.5 * TaitBryant[2]));
+ Quaternion rotTheta(std::cos(0.5 * TaitBryant[0]), 0, std::sin(0.5 * TaitBryant[0]), 0);
+ Quaternion rotPhi(std::cos(0.5 * TaitBryant[1]), std::sin(0.5 * TaitBryant[1]), 0, 0);
+ Quaternion rotPsi(std::cos(0.5 * TaitBryant[2]), 0, 0, std::sin(0.5 * TaitBryant[2]));
Quaternion rotation = rotTheta * (rotPhi * rotPsi);
bunch->toLabTrafo_m = CoordinateSystemTrafo(-rotation.conjugate().rotate(RefPartR), rotation);
}
@@ -143,47 +157,47 @@ void H5PartWrapperForPT::readStepData(PartBunchBase<double, 3>* bunch, h5_ssize_
std::vector<char> buffer(numParticles * sizeof(h5_float64_t));
char* buffer_ptr = Util::c_data(buffer);
- h5_float64_t *f64buffer = reinterpret_cast<h5_float64_t*>(buffer_ptr);
- h5_int32_t *i32buffer = reinterpret_cast<h5_int32_t*>(buffer_ptr);
+ h5_float64_t* f64buffer = reinterpret_cast<h5_float64_t*>(buffer_ptr);
+ h5_int32_t* i32buffer = reinterpret_cast<h5_int32_t*>(buffer_ptr);
READDATA(Float64, file_m, "x", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->R[n](0) = f64buffer[n];
bunch->Bin[n] = 0;
}
READDATA(Float64, file_m, "y", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->R[n](1) = f64buffer[n];
}
READDATA(Float64, file_m, "z", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->R[n](2) = f64buffer[n];
}
READDATA(Float64, file_m, "px", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->P[n](0) = f64buffer[n];
}
READDATA(Float64, file_m, "py", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->P[n](1) = f64buffer[n];
}
READDATA(Float64, file_m, "pz", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->P[n](2) = f64buffer[n];
}
READDATA(Float64, file_m, "q", f64buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->Q[n] = f64buffer[n];
}
READDATA(Int32, file_m, "id", i32buffer);
- for(long int n = 0; n < numParticles; ++ n) {
+ for (long int n = 0; n < numParticles; ++ n) {
bunch->ID[n] = i32buffer[n];
}
@@ -206,6 +220,7 @@ void H5PartWrapperForPT::writeHeader() {
WRITESTRINGFILEATTRIB(file_m, "pzUnit", "#beta#gamma");
WRITESTRINGFILEATTRIB(file_m, "ptypeUnit", "1");
+ WRITESTRINGFILEATTRIB(file_m, "poriginUnit", "1");
WRITESTRINGFILEATTRIB(file_m, "qUnit", "C");
if (Options::ebDump) {
@@ -273,7 +288,7 @@ void H5PartWrapperForPT::writeHeader() {
WRITEFILEATTRIB(Float64, file_m, "dPhiGlobal", &dphi, 1);
}
-void H5PartWrapperForPT::writeStep(PartBunchBase<double, 3>* bunch, const std::map<std::string, double> &additionalStepAttributes) {
+void H5PartWrapperForPT::writeStep(PartBunchBase<double, 3>* bunch, const std::map<std::string, double>& additionalStepAttributes) {
if (bunch->getTotalNum() == 0) return;
open(H5_O_APPENDONLY);
@@ -282,7 +297,7 @@ void H5PartWrapperForPT::writeStep(PartBunchBase<double, 3>* bunch, const std::m
close();
}
-void H5PartWrapperForPT::writeStepHeader(PartBunchBase<double, 3>* bunch, const std::map<std::string, double> &additionalStepAttributes) {
+void H5PartWrapperForPT::writeStepHeader(PartBunchBase<double, 3>* bunch, const std::map<std::string, double>& additionalStepAttributes) {
bunch->calcBeamParameters();
double actPos = bunch->get_sPos();
@@ -324,26 +339,26 @@ void H5PartWrapperForPT::writeStepHeader(PartBunchBase<double, 3>* bunch, const
REPORTONERROR(H5SetStep(file_m, numSteps_m));
- char const *OPALFlavour = "opal-t";
+ char const* OPALFlavour = "opal-t";
WRITESTRINGSTEPATTRIB(file_m, "OPAL_flavour", OPALFlavour);
WRITESTEPATTRIB(Float64, file_m, "SPOS", &actPos, 1);
- WRITESTEPATTRIB(Float64, file_m, "RefPartR", (h5_float64_t *)&RefPartR, 3);
- WRITESTEPATTRIB(Float64, file_m, "centroid", (h5_float64_t *)¢roid, 3);
- WRITESTEPATTRIB(Float64, file_m, "RMSX", (h5_float64_t *)&xsigma, 3);
+ WRITESTEPATTRIB(Float64, file_m, "RefPartR", (h5_float64_t*)&RefPartR, 3);
+ WRITESTEPATTRIB(Float64, file_m, "centroid", (h5_float64_t*)¢roid, 3);
+ WRITESTEPATTRIB(Float64, file_m, "RMSX", (h5_float64_t*)&xsigma, 3);
- WRITESTEPATTRIB(Float64, file_m, "RefPartP", (h5_float64_t *)&RefPartP, 3);
- WRITESTEPATTRIB(Float64, file_m, "MEANP", (h5_float64_t *)&pmean, 3);
- WRITESTEPATTRIB(Float64, file_m, "RMSP", (h5_float64_t *)&psigma, 3);
- WRITESTEPATTRIB(Float64, file_m, "TaitBryantAngles", (h5_float64_t *)&TaitBryant, 3);
+ WRITESTEPATTRIB(Float64, file_m, "RefPartP", (h5_float64_t*)&RefPartP, 3);
+ WRITESTEPATTRIB(Float64, file_m, "MEANP", (h5_float64_t*)&pmean, 3);
+ WRITESTEPATTRIB(Float64, file_m, "RMSP", (h5_float64_t*)&psigma, 3);
+ WRITESTEPATTRIB(Float64, file_m, "TaitBryantAngles", (h5_float64_t*)&TaitBryant, 3);
- WRITESTEPATTRIB(Float64, file_m, "#varepsilon", (h5_float64_t *)&vareps, 3);
- WRITESTEPATTRIB(Float64, file_m, "#varepsilon-geom", (h5_float64_t *)&geomvareps, 3);
+ WRITESTEPATTRIB(Float64, file_m, "#varepsilon", (h5_float64_t*)&vareps, 3);
+ WRITESTEPATTRIB(Float64, file_m, "#varepsilon-geom", (h5_float64_t*)&geomvareps, 3);
- WRITESTEPATTRIB(Float64, file_m, "minX", (h5_float64_t *)&rmin, 3);
- WRITESTEPATTRIB(Float64, file_m, "maxX", (h5_float64_t *)&rmax, 3);
+ WRITESTEPATTRIB(Float64, file_m, "minX", (h5_float64_t*)&rmin, 3);
+ WRITESTEPATTRIB(Float64, file_m, "maxX", (h5_float64_t*)&rmax, 3);
- WRITESTEPATTRIB(Float64, file_m, "minP", (h5_float64_t *)&minP, 3);
- WRITESTEPATTRIB(Float64, file_m, "maxP", (h5_float64_t *)&maxP, 3);
+ WRITESTEPATTRIB(Float64, file_m, "minP", (h5_float64_t*)&minP, 3);
+ WRITESTEPATTRIB(Float64, file_m, "maxP", (h5_float64_t*)&maxP, 3);
WRITESTEPATTRIB(Int64, file_m, "Step", &numSteps_m, 1);
WRITESTEPATTRIB(Int64, file_m, "LocalTrackStep", &localTrackStep, 1);
@@ -373,8 +388,8 @@ void H5PartWrapperForPT::writeStepHeader(PartBunchBase<double, 3>* bunch, const
additionalStepAttributes.at("E-ref_z"),
additionalStepAttributes.at("E-ref_y"));
- WRITESTEPATTRIB(Float64, file_m, "B-ref", (h5_float64_t *)&referenceB, 3);
- WRITESTEPATTRIB(Float64, file_m, "E-ref", (h5_float64_t *)&referenceE, 3);
+ WRITESTEPATTRIB(Float64, file_m, "B-ref", (h5_float64_t*)&referenceB, 3);
+ WRITESTEPATTRIB(Float64, file_m, "E-ref", (h5_float64_t*)&referenceE, 3);
} catch (std::out_of_range & m) {
ERRORMSG(m.what() << endl);
@@ -392,75 +407,79 @@ void H5PartWrapperForPT::writeStepData(PartBunchBase<double, 3>* bunch) {
std::vector<char> buffer(numLocalParticles * sizeof(h5_float64_t));
char* buffer_ptr = Util::c_data(buffer);
- h5_float64_t *f64buffer = reinterpret_cast<h5_float64_t*>(buffer_ptr);
- h5_int64_t *i64buffer = reinterpret_cast<h5_int64_t*>(buffer_ptr);
- h5_int32_t *i32buffer = reinterpret_cast<h5_int32_t*>(buffer_ptr);
+ h5_float64_t* f64buffer = reinterpret_cast<h5_float64_t*>(buffer_ptr);
+ h5_int64_t* i64buffer = reinterpret_cast<h5_int64_t*>(buffer_ptr);
+ h5_int32_t* i32buffer = reinterpret_cast<h5_int32_t*>(buffer_ptr);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->R[i](0);
WRITEDATA(Float64, file_m, "x", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->R[i](1);
WRITEDATA(Float64, file_m, "y", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->R[i](2);
WRITEDATA(Float64, file_m, "z", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->P[i](0);
WRITEDATA(Float64, file_m, "px", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->P[i](1);
WRITEDATA(Float64, file_m, "py", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->P[i](2);
WRITEDATA(Float64, file_m, "pz", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Q[i];
WRITEDATA(Float64, file_m, "q", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
i64buffer[i] = bunch->ID[i];
WRITEDATA(Int64, file_m, "id", i64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
i32buffer[i] = (h5_int32_t) bunch->PType[i];
WRITEDATA(Int32, file_m, "ptype", i32buffer);
- if(Options::ebDump) {
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
+ i32buffer[i] = (h5_int32_t) bunch->POrigin[i];
+ WRITEDATA(Int32, file_m, "porigin", i32buffer);
+
+ if (Options::ebDump) {
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Ef[i](0);
WRITEDATA(Float64, file_m, "Ex", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Ef[i](1);
WRITEDATA(Float64, file_m, "Ey", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Ef[i](2);
WRITEDATA(Float64, file_m, "Ez", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Bf[i](0);
WRITEDATA(Float64, file_m, "Bx", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Bf[i](1);
WRITEDATA(Float64, file_m, "By", f64buffer);
- for(size_t i = 0; i < numLocalParticles; ++ i)
+ for (size_t i = 0; i < numLocalParticles; ++ i)
f64buffer[i] = bunch->Bf[i](2);
WRITEDATA(Float64, file_m, "Bz", f64buffer);
}
/// Write space charge field map if asked for.
- if(Options::rhoDump) {
+ if (Options::rhoDump) {
NDIndex<3> idx = bunch->getFieldLayout().getLocalNDIndex();
NDIndex<3> elem;
h5_err_t herr = H5Block3dSetView(
@@ -476,9 +495,9 @@ void H5PartWrapperForPT::writeStepData(PartBunchBase<double, 3>* bunch) {
// h5block uses the fortran convention of storing data:
// INTEGER, DIMENSION(2,3) :: a
// => {a(1,1), a(2,1), a(1,2), a(2,2), a(1,3), a(2,3)}
- for(int i = idx[2].min(); i <= idx[2].max(); ++ i) {
- for(int j = idx[1].min(); j <= idx[1].max(); ++ j) {
- for(int k = idx[0].min(); k <= idx[0].max(); ++ k) {
+ for (int i = idx[2].min(); i <= idx[2].max(); ++ i) {
+ for (int j = idx[1].min(); j <= idx[1].max(); ++ j) {
+ for (int k = idx[0].min(); k <= idx[0].max(); ++ k) {
data[ii] = bunch->getRho(k, j, i);
++ ii;
}
@@ -501,4 +520,4 @@ void H5PartWrapperForPT::writeStepData(PartBunchBase<double, 3>* bunch) {
reportOnError(herr, __FILE__, __LINE__);
}
-}
\ No newline at end of file
+}
diff --git a/src/Structure/H5PartWrapperForPT.h b/src/Structure/H5PartWrapperForPT.h
index 80022896068a7e9569cbd41f25ba52b2f61ecf04..08b9935d5d4266362826cdae3d7ca1802c0b6c93 100644
--- a/src/Structure/H5PartWrapperForPT.h
+++ b/src/Structure/H5PartWrapperForPT.h
@@ -1,13 +1,22 @@
-#ifndef OPAL_H5PARTWRAPPERFORPT_H
-#define OPAL_H5PARTWRAPPERFORPT_H
-
//
-// Copyright & License: See Copyright.readme in src directory
+// Class H5PartWrapperForPT
+// A class that manages all calls to H5Part for the Parallel-T tracker.
//
-
-/*!
- H5PartWrapperForPT: a class that manages all calls to H5Part for the Parallel-T tracker
-*/
+// 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/>.
+//
+#ifndef OPAL_H5PARTWRAPPERFORPT_H
+#define OPAL_H5PARTWRAPPERFORPT_H
#include "Structure/H5PartWrapper.h"
@@ -15,15 +24,15 @@
class H5PartWrapperForPT: public H5PartWrapper {
public:
- H5PartWrapperForPT(const std::string &fileName, h5_int32_t flags = H5_O_WRONLY);
- H5PartWrapperForPT(const std::string &fileName, int restartStep, std::string sourceFile, h5_int32_t flags = H5_O_RDWR);
+ H5PartWrapperForPT(const std::string& fileName, h5_int32_t flags = H5_O_WRONLY);
+ H5PartWrapperForPT(const std::string& fileName, int restartStep, std::string sourceFile, h5_int32_t flags = H5_O_RDWR);
virtual ~H5PartWrapperForPT();
virtual void readHeader();
virtual void readStep(PartBunchBase<double, 3>*, h5_ssize_t firstParticle, h5_ssize_t lastParticle);
virtual void writeHeader();
- virtual void writeStep(PartBunchBase<double, 3>*, const std::map<std::string, double> &additionalStepAttributes);
+ virtual void writeStep(PartBunchBase<double, 3>*, const std::map<std::string, double>& additionalStepAttributes);
virtual bool predecessorIsSameFlavour() const;
@@ -31,7 +40,7 @@ private:
void readStepHeader(PartBunchBase<double, 3>*);
void readStepData(PartBunchBase<double, 3>*, h5_ssize_t, h5_ssize_t);
- void writeStepHeader(PartBunchBase<double, 3>*, const std::map<std::string, double> &);
+ void writeStepHeader(PartBunchBase<double, 3>*, const std::map<std::string, double>&);
void writeStepData(PartBunchBase<double, 3>*);
};
@@ -40,4 +49,4 @@ bool H5PartWrapperForPT::predecessorIsSameFlavour() const {
return (predecessorOPALFlavour_m == "opal-t");
}
-#endif //OPAL_H5PARTWRAPPERFORPT_H
\ No newline at end of file
+#endif //OPAL_H5PARTWRAPPERFORPT_H
diff --git a/src/Track/TrackRun.cpp b/src/Track/TrackRun.cpp
index fc046a22974beb346340e05494acf6a561010e17..ef3b14ed560a0d5cb8947642247789df97985d97 100644
--- a/src/Track/TrackRun.cpp
+++ b/src/Track/TrackRun.cpp
@@ -323,6 +323,8 @@ void TrackRun::setupTTracker(){
Beam *beam = Beam::find(Attributes::getString(itsAttr[BEAM]));
Track::block->bunch->setBeamFrequency(beam->getFrequency() * 1e6);
+ Track::block->bunch->setPType(beam->getParticleName());
+
if (Attributes::getString(itsAttr[BOUNDARYGEOMETRY]) != "NONE") {
// Ask the dictionary if BoundaryGeometry is allocated.
// If it is allocated use the allocated BoundaryGeometry
@@ -436,7 +438,8 @@ void TrackRun::setupCyclotronTracker(){
setupFieldsolver();
- Track::block->bunch->PType = ParticleType::REGULAR;
+ Track::block->bunch->setPType(beam->getParticleName());
+ Track::block->bunch->POrigin = ParticleOrigin::REGULAR;
std::vector<std::string> distr_str = Attributes::getStringArray(itsAttr[DISTRIBUTION]);
if (distr_str.size() == 0) {