Commit 88cb9258 authored by ext-calvo_p's avatar ext-calvo_p
Browse files

Implementation beam stripping physics into master

parent 8ac6b29e
......@@ -108,6 +108,9 @@ void LieMapper::visitBeamBeam(const BeamBeam &map) {
// *** MISSING *** Map for beam-beam.
}
void LieMapper::visitBeamStripping(const BeamStripping &map) {
// *** MISSING *** Map for beam stripping.
}
void LieMapper::visitCCollimator(const CCollimator &coll) {
applyDrift(flip_s * coll.getElementLength());
......@@ -526,4 +529,4 @@ void LieMapper::applyTransform(const Euclid3D &euclid, double refLength) {
theMap.assign(F);
itsMap = itsMap.catenate(theMap);
}
}
}
\ No newline at end of file
......@@ -93,6 +93,9 @@ public:
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a BeamStripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a Collimator.
virtual void visitCCollimator(const CCollimator &);
......@@ -188,4 +191,4 @@ private:
int itsOrder;
};
#endif // OPAL_LieMapper_HH
#endif // OPAL_LieMapper_HH
\ No newline at end of file
......@@ -18,6 +18,7 @@ template <class T, unsigned Dim>
class PartBunchBase;
class AlignWrapper;
class BeamBeam;
class BeamStripping;
class CCollimator;
class Corrector;
class CyclotronValley;
......@@ -55,6 +56,7 @@ public:
NIL_VISITELEMENT(AlignWrapper)
NIL_VISITELEMENT(Beamline)
NIL_VISITELEMENT(BeamBeam)
NIL_VISITELEMENT(BeamStripping)
NIL_VISITELEMENT(CCollimator)
NIL_VISITELEMENT(Corrector)
NIL_VISITELEMENT(CyclotronValley)
......
......@@ -27,6 +27,7 @@
#include "AbstractObjects/Element.h"
#include "AbstractObjects/OpalData.h"
#include "AbsBeamline/BeamStripping.h"
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Cyclotron.h"
......@@ -384,12 +385,12 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << "* -------------------------- Adding Cyclotron ---------------------------- *" << endl;
Cyclotron *elptr = dynamic_cast<Cyclotron *>(cycl.clone());
myElements.push_back(elptr);
cycl_m = dynamic_cast<Cyclotron *>(cycl.clone());
myElements.push_back(cycl_m);
// Is this a Spiral Inflector Simulation? If yes, we'll give the user some
// useful information
spiral_flag = elptr->getSpiralFlag();
spiral_flag = cycl_m->getSpiralFlag();
if(spiral_flag) {
......@@ -415,11 +416,11 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
// Get reference values from cyclotron element
// For now, these are still stored in mm. should be the only ones. -DW
referenceR = elptr->getRinit();
referenceTheta = elptr->getPHIinit();
referenceZ = elptr->getZinit();
referencePr = elptr->getPRinit();
referencePz = elptr->getPZinit();
referenceR = cycl_m->getRinit();
referenceTheta = cycl_m->getPHIinit();
referenceZ = cycl_m->getZinit();
referencePr = cycl_m->getPRinit();
referencePz = cycl_m->getPZinit();
if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
throw OpalException("Error in ParallelCyclotronTracker::visitCyclotron",
......@@ -500,29 +501,29 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
*gmsg << "* Reference axial momentum (Pz) = " << referencePz * 1000.0 << " [MCU]" << endl;
*gmsg << endl;
double sym = elptr->getSymmetry();
double sym = cycl_m->getSymmetry();
*gmsg << "* " << sym << "-fold field symmetry " << endl;
// ckr: this just returned the default value as defined in Component.h
// double rff = elptr->getRfFrequ();
// double rff = cycl_m->getRfFrequ();
// *gmsg << "* Rf frequency= " << rff << " [MHz]" << endl;
std::string fmfn = elptr->getFieldMapFN();
std::string fmfn = cycl_m->getFieldMapFN();
*gmsg << "* Field map file name = " << fmfn << " " << endl;
std::string type = elptr->getCyclotronType();
std::string type = cycl_m->getCyclotronType();
*gmsg << "* Type of cyclotron = " << type << " " << endl;
double rmin = elptr->getMinR();
double rmax = elptr->getMaxR();
double rmin = cycl_m->getMinR();
double rmax = cycl_m->getMaxR();
*gmsg << "* Radial aperture = " << rmin << " ... " << rmax<<" [m] "<< endl;
double zmin = elptr->getMinZ();
double zmax = elptr->getMaxZ();
double zmin = cycl_m->getMinZ();
double zmax = cycl_m->getMaxZ();
*gmsg << "* Vertical aperture = " << zmin << " ... " << zmax<<" [m]"<< endl;
double h = elptr->getCyclHarm();
*gmsg << "* Number of trimcoils = " << elptr->getNumberOfTrimcoils() << endl;
double h = cycl_m->getCyclHarm();
*gmsg << "* Number of trimcoils = " << cycl_m->getNumberOfTrimcoils() << endl;
*gmsg << "* Harmonic number h = " << h << " " << endl;
/**
......@@ -536,18 +537,18 @@ void ParallelCyclotronTracker::visitCyclotron(const Cyclotron &cycl) {
* fieldflag = 6, readin both median plane B field map and 3D E field map of RF cavity for compact cyclotron
* fieldflag = 7, read in fields for Daniel's synchrocyclotron simulations
*/
int fieldflag = elptr->getFieldFlag(type);
int fieldflag = cycl_m->getFieldFlag(type);
// Read in cyclotron field maps (midplane + 3D fields if desired).
elptr->initialise(itsBunch_m, fieldflag, elptr->getBScale());
cycl_m->initialise(itsBunch_m, fieldflag, cycl_m->getBScale());
double BcParameter[8] = {};
BcParameter[0] = 0.001 * elptr->getRmin();
BcParameter[1] = 0.001 * elptr->getRmax();
BcParameter[0] = 0.001 * cycl_m->getRmin();
BcParameter[1] = 0.001 * cycl_m->getRmax();
// Store inner radius and outer radius of cyclotron field map in the list
buildupFieldList(BcParameter, ElementBase::CYCLOTRON, elptr);
buildupFieldList(BcParameter, ElementBase::CYCLOTRON, cycl_m);
}
/**
......@@ -558,6 +559,39 @@ void ParallelCyclotronTracker::visitBeamBeam(const BeamBeam &) {
*gmsg << "In BeamBeam tracker is missing " << endl;
}
void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) {
*gmsg << "* ------------------------------ Beam Stripping ------------------------------" << endl;
BeamStripping* elptr = dynamic_cast<BeamStripping *>(bstp.clone());
myElements.push_back(elptr);
double BcParameter[8] = {};
if(elptr->getPressureMapFN() == "") {
double pressure = elptr->getPressure();
*gmsg << "* Pressure = " << pressure << " [mbar]" << endl;
BcParameter[0] = pressure;
}
double temperature = elptr->getTemperature();
*gmsg << "* Temperature = " << temperature << " [K]" << endl;
std::string gas = elptr->getResidualGas();
*gmsg << "* Residual gas = " << gas << endl;
bool stop = elptr->getStop();
*gmsg << std::boolalpha << "* Particles stripped will be deleted after interaction -> " << stop << endl;
elptr->initialise(itsBunch_m, elptr->getPScale());
BcParameter[1] = temperature;
BcParameter[2] = stop;
buildupFieldList(BcParameter, ElementBase::BEAMSTRIPPING, elptr);
}
/**
*
*
......@@ -1183,7 +1217,7 @@ void ParallelCyclotronTracker::execute() {
opalRing_m->lockRing();
// Display the selected elements
*gmsg << "* -------------------------------------" << endl;
*gmsg << "* ---------------------------------------------------" << endl;
*gmsg << "* The selected Beam line elements are :" << endl;
for(auto fd : FieldDimensions) {
......@@ -1202,7 +1236,7 @@ void ParallelCyclotronTracker::execute() {
}
}
*gmsg << "* -------------------------------------" << endl;
*gmsg << "* ---------------------------------------------------" << endl;
// Get BoundaryGeometry that is already initialized
bgf_m = OpalData::getInstance()->getGlobalGeometry();
......@@ -2074,6 +2108,14 @@ bool ParallelCyclotronTracker::applyPluginElements(const double dt) {
// Plugin Elements are all defined in mm, change beam to mm before applying
itsBunch_m->R *= Vector_t(1000.0);
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);
}
}
bool flag = false;
for (PluginElement* element : pluginElements_m) {
bool tmp = element->check(itsBunch_m,
......@@ -2084,7 +2126,8 @@ bool ParallelCyclotronTracker::applyPluginElements(const double dt) {
if ( tmp ) {
itsBunch_m->updateNumTotal();
*gmsg << "* Total number of particles = " << itsBunch_m->getTotalNum() << endl;
*gmsg << "* Total number of particles after PluginElement= "
<< itsBunch_m->getTotalNum() << endl;
}
}
......@@ -2173,11 +2216,13 @@ bool ParallelCyclotronTracker::deleteParticle(bool flagNeedUpdate){
reduce(locNumLost, globLostParticleNum, 1, std::plus<size_t>());
*gmsg << "At step " << step_m
<< ", lost " << globLostParticleNum << " particles "
<< "on stripper, collimator, septum, geometry, "
<< "out of cyclotron aperture or beam stripping"
<< endl;
if ( globLostParticleNum > 0 ) {
*gmsg << "At step " << step_m
<< ", lost " << globLostParticleNum << " particles "
<< "on stripper, collimator, septum, geometry, "
<< "out of cyclotron aperture or beam stripping"
<< endl;
}
if (totalnum[bunchCount] == 0) {
IpplTimings::stopTimer(DelParticleTimer_m);
......@@ -2509,7 +2554,6 @@ void ParallelCyclotronTracker::singleParticleDump() {
} else {
for(size_t i = 0; i < itsBunch_m->getLocalNum(); i++) {
if(itsBunch_m->ID[i] == 0 || itsBunch_m->ID[i] == 1) {
outfTrackOrbit_m << "ID" << itsBunch_m->ID[i] << " ";
......@@ -2675,6 +2719,7 @@ void ParallelCyclotronTracker::bunchDumpPhaseSpaceData() {
}
double const betagamma_temp = sqrt(dot(meanP, meanP));
double const E = itsBunch_m->get_meanKineticEnergy();
// Bunch (global) angle w.r.t. x-axis (cylinder coordinates)
......@@ -3529,4 +3574,4 @@ void ParallelCyclotronTracker::initPathLength() {
// we need to reset the path length of each bunch
itsDataSink->setMultiBunchInitialPathLengh(mbHandler_m.get());
}
}
}
\ No newline at end of file
......@@ -54,7 +54,6 @@ public:
SEO = 1,
BUNCH = 2
};
typedef std::vector<double> dvector_t;
typedef std::vector<int> ivector_t;
typedef std::pair<double[8], Component *> element_pair;
......@@ -91,6 +90,9 @@ public:
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a Beam Stripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a collimator.
virtual void visitCCollimator(const CCollimator &);
......@@ -222,6 +224,8 @@ private:
BoundaryGeometry *bgf_m;
Cyclotron *cycl_m;
/// The maximal number of steps the system is integrated
int maxSteps_m;
......@@ -576,4 +580,4 @@ bool ParallelCyclotronTracker::hasMultiBunch() const {
return ( isMultiBunch() && mbHandler_m->getNumBunch() > 1 );
}
#endif // OPAL_ParallelCyclotronTracker_HH
#endif // OPAL_ParallelCyclotronTracker_HH
\ No newline at end of file
......@@ -32,6 +32,7 @@
#include "Algorithms/IndexMap.h"
#include "AbsBeamline/AlignWrapper.h"
#include "AbsBeamline/BeamBeam.h"
#include "AbsBeamline/BeamStripping.h"
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/Diagnostic.h"
......@@ -105,6 +106,9 @@ public:
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a BeamStripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a collimator.
virtual void visitCCollimator(const CCollimator &);
......@@ -325,6 +329,9 @@ inline void ParallelTTracker::visitBeamBeam(const BeamBeam &bb) {
itsOpalBeamline_m.visit(bb, *this, itsBunch_m);
}
inline void ParallelTTracker::visitBeamStripping(const BeamStripping &bstp) {
itsOpalBeamline_m.visit(bstp, *this, itsBunch_m);
}
inline void ParallelTTracker::visitCCollimator(const CCollimator &coll) {
itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
......
......@@ -92,6 +92,9 @@ void ThickMapper::visitBeamBeam(const BeamBeam &) {
// *** MISSING *** Map for beam-beam.
}
void ThickMapper::visitBeamStripping(const BeamStripping &) {
// *** MISSING *** Map for beam Stripping.
}
void ThickMapper::visitCCollimator(const CCollimator &coll) {
applyDrift(flip_s * coll.getElementLength());
......@@ -763,4 +766,4 @@ namespace {
//std::cerr << "==> Leaving implicitIntStep(zin,f,gradf,ds,nx)" << std::endl;
return zf;
}
};
};
\ No newline at end of file
......@@ -98,6 +98,9 @@ public:
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a BeamStripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a collimator.
virtual void visitCCollimator(const CCollimator &);
......@@ -176,4 +179,4 @@ private:
};
#endif // OPAL_ThickMapper_HH
#endif // OPAL_ThickMapper_HH
\ No newline at end of file
......@@ -22,7 +22,7 @@
#include "Algorithms/IndexMap.h"
#include "AbsBeamline/BeamBeam.h"
#include "AbsBeamline/BeamStripping.h"
#include "AbsBeamline/CCollimator.h"
#include "AbsBeamline/Corrector.h"
#include "AbsBeamline/CyclotronValley.h"
......@@ -149,6 +149,9 @@ public:
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a BeamStripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a collimator.
virtual void visitCCollimator(const CCollimator &);
......@@ -391,6 +394,9 @@ inline void ThickTracker::visitBeamBeam(const BeamBeam &bb) {
this->throwElementError_m("BeamBeam");
}
inline void ThickTracker::visitBeamStripping(const BeamStripping &bstp) {
itsOpalBeamline_m.visit(bstp, *this, itsBunch_m);
}
inline void ThickTracker::visitCCollimator(const CCollimator &coll) {
// itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
......@@ -582,4 +588,4 @@ inline void ThickTracker::visitCyclotronValley(const CyclotronValley &cv) {
this->throwElementError_m("CyclotronValley");
}
#endif // OPAL_ThickTracker_HH
#endif // OPAL_ThickTracker_HH
\ No newline at end of file
......@@ -120,6 +120,9 @@ void TransportMapper::visitBeamBeam(const BeamBeam &) {
// *** MISSING *** Map for beam-beam.
}
void TransportMapper::visitBeamStripping(const BeamStripping &) {
// *** MISSING *** Map for beam Stripping.
}
void TransportMapper::visitCCollimator(const CCollimator &coll) {
applyDrift(flip_s * coll.getElementLength());
......@@ -757,4 +760,4 @@ buildSBendVectorPotential(const BMultipoleField &field, double h) {
}
return As;
}
}
\ No newline at end of file
......@@ -102,6 +102,9 @@ public:
/// Apply the algorithm to a BeamBeam.
virtual void visitBeamBeam(const BeamBeam &);
/// Apply the algorithm to a BeamStripping.
virtual void visitBeamStripping(const BeamStripping &);
/// Apply the algorithm to a collimator.
virtual void visitCCollimator(const CCollimator &);
......@@ -225,4 +228,4 @@ private:
TransportMap <double, 6> itsMap;
};
#endif // OPAL_TransportMapper_HH
#endif // OPAL_TransportMapper_HH
\ No newline at end of file
// ------------------------------------------------------------------------
// $RCSfile: BeamStripping.cpp,v $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: BeamStripping
// Defines the abstract interface for a beam stripping for cyclotrons.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
// $Date: 2018/11 $
// $Author: PedroCalvo$
// ------------------------------------------------------------------------
#include "AbsBeamline/BeamlineVisitor.h"
#include "AbsBeamline/BeamStripping.h"
#include "AbsBeamline/Cyclotron.h"
#include "Algorithms/PartBunchBase.h"
#include "Fields/Fieldmap.h"
#include "Physics/Physics.h"
#include "Solvers/BeamStrippingPhysics.hh"
#include "Solvers/ParticleMatterInteractionHandler.hh"
#include "Utilities/LogicalError.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "Utilities/GeneralClassicException.h"
#include <fstream>
#include <memory>
#include <fstream>
#include <iostream>
#include <fstream>
#include <algorithm>
#include "Ippl.h"
#include "gsl/gsl_spline.h"
#include "gsl/gsl_interp.h"
using Physics::q_e;
using Physics::pi;
#define CHECK_BSTP_FSCANF_EOF(arg) if(arg == EOF)\
throw GeneralClassicException("BeamStripping::getPressureFromFile",\
"fscanf returned EOF at " #arg);
extern Inform *gmsg;
using namespace std;
// Class BeamStripping
// ------------------------------------------------------------------------
BeamStripping::BeamStripping():
BeamStripping("")
{}
BeamStripping::BeamStripping(const BeamStripping &right):
Component(right),
gas_m(right.gas_m),
pressure_m(right.pressure_m),
pmapfn_m(right.pmapfn_m),
pscale_m(right.pscale_m),
temperature_m(right.temperature_m),
stop_m(right.stop_m),
minr_m(right.minr_m),
maxr_m(right.maxr_m),
minz_m(right.minz_m),
maxz_m(right.maxz_m),
parmatint_m(NULL) {
}
BeamStripping::BeamStripping(const std::string &name):
Component(name),
gas_m(""),
pressure_m(0.0),
pmapfn_m(""),
pscale_m(0.0),
temperature_m(0.0),
stop_m(true),
minr_m(0.0),
maxr_m(0.0),
minz_m(0.0),
maxz_m(0.0),
parmatint_m(NULL) {
}
BeamStripping::~BeamStripping() {
if (online_m)
goOffline();
}
void BeamStripping::accept(BeamlineVisitor &visitor) const {
visitor.visitBeamStripping(*this);
}
void BeamStripping::setPressure(double pressure) {
pressure_m = pressure;
}
double BeamStripping::getPressure() const {
if(pressure_m > 0.0)
return pressure_m;
else {
throw LogicalError("BeamStripping::getPressure",
"Pressure must not be zero");
}
}
void BeamStripping::setTemperature(double temperature) {
temperature_m = temperature;
}
double BeamStripping::getTemperature() const {
if(temperature_m > 0.0)
return temperature_m;
else {
throw LogicalError("BeamStripping::getTemperature",
"Temperature must not be zero");
}
}
void BeamStripping::setPressureMapFN(std::string p) {
pmapfn_m = p;
}
string BeamStripping::getPressureMapFN() const {
return pmapfn_m;
}
void BeamStripping::setPScale(double ps) {
pscale_m = ps;
}
double BeamStripping::getPScale() const {
return pscale_m;
}
void BeamStripping::setResidualGas(std::string gas) {
gas_m = gas;
}
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");
}
}
void BeamStripping::setStop(bool stopflag) {
stop_m = stopflag;
}
bool BeamStripping::getStop() const {
return stop_m;
}
bool BeamStripping::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
bool BeamStripping::apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
bool BeamStripping::checkBeamStripping(Vector_t r, Vector_t rmin, Vector_t rmax) {
int pflag = checkPoint(r(0), r(1), r(2));
bool isDead = (pflag != 0);
return isDead;
}
bool BeamStripping::checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotron* cycl,
const int turnnumber, const double t, const double tstep) {