Commit 44afc942 authored by ext-calvo_p's avatar ext-calvo_p

Merge branch 'OPAL-2.0' into 'OPAL-2.0'

Improving beam stripping

See merge request !177
parents c8bc3087 48127b02
......@@ -225,7 +225,9 @@ ParallelCyclotronTracker::~ParallelCyclotronTracker() {
double ParallelCyclotronTracker::initializeTimeStep() {
double E_mean = 0.0;
// Time step changes in each step accordng to energy
double E_mean = 0.0;
double dT = 0.0;
double dti = itsBunch_m->getdTinit();
......@@ -235,24 +237,17 @@ double ParallelCyclotronTracker::initializeTimeStep() {
double Ef = itsBunch_m->getEfinal();
if(dti > 0) {
// *gmsg << "* Time step changes in each step accordng to energy " << endl;
switch (mode_m)
{
case MODE::SEO:
{
case MODE::SEO: {
E_mean = (sqrt(1.0 + dot(itsBunch_m->P[0], itsBunch_m->P[0])) - 1) * itsBunch_m->getM();
break;
}
case MODE::SINGLE:
{
break;}
case MODE::SINGLE: {
E_mean = (sqrt(1.0 + dot(itsBunch_m->P[0], itsBunch_m->P[0])) - 1) * itsBunch_m->getM();
break;
}
case MODE::BUNCH:
{
break;}
case MODE::BUNCH: {
E_mean = itsBunch_m->get_meanKineticEnergy() * 1.0e6;
break;
}
break; }
case MODE::UNDEFINED:
default:
throw OpalException("ParallelCyclotronTracker::initializeTimeStep",
......@@ -269,9 +264,8 @@ double ParallelCyclotronTracker::initializeTimeStep() {
// *gmsg << "* dT = " << dT << endl;
itsBunch_m->setdT(dT);
double harm = getHarmonicNumber();
dT *= 1.0e9 * harm;
double harm = getHarmonicNumber();
dT *= 1.0e9 * harm;
return dT;
}
......@@ -343,7 +337,7 @@ void ParallelCyclotronTracker::bgf_main_collision_test() {
if(res >= 0) {
lossDs_m->addParticle(itsBunch_m->R[i]*1000, itsBunch_m->P[i], itsBunch_m->ID[i], itsBunch_m->getT()*1e9, turnnumber_m);
itsBunch_m->Bin[i] = -1;
INFOMSG(level2 << "* Particle " << itsBunch_m->ID[i] << " lost on boundary geometry" << endl;);
*gmsg << level4 << "* Particle " << itsBunch_m->ID[i] << " lost on boundary geometry" << endl;
}
}
}
......@@ -672,22 +666,27 @@ void ParallelCyclotronTracker::visitBeamStripping(const BeamStripping &bstp) {
BeamStripping* elptr = dynamic_cast<BeamStripping *>(bstp.clone());
myElements.push_back(elptr);
double pressure = elptr->getPressure();
*gmsg << "* Pressure = " << pressure << " [mbar]" << endl;
double BcParameter[8];
for(int i = 0; i < 6; i++)
BcParameter[i] = 0.0;
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);
double BcParameter[8];
for(int i = 0; i < 4; i++)
BcParameter[i] = 0.0;
BcParameter[0] = pressure;
elptr->initialise(itsBunch_m, elptr->getPScale());
BcParameter[1] = temperature;
BcParameter[2] = stop;
......@@ -1170,7 +1169,9 @@ void ParallelCyclotronTracker::visitStripper(const Stripper &stripper) {
Stripper *elptr = dynamic_cast<Stripper *>(stripper.clone());
myElements.push_back(elptr);
*gmsg << "* Name = " << elptr->getName() << endl;
double xstart = elptr->getXstart();
*gmsg << "* Xstart = " << xstart << " [mm]" << endl;
......@@ -1633,7 +1634,7 @@ bool ParallelCyclotronTracker::checkGapCross(Vector_t Rold, Vector_t Rnew,
double distOld = (Rold[0] * sinx - Rold[1] * cosx) - PerpenDistance;
if(distOld > 0.0 && distNew <= 0.0) flag = true;
// This parameter is used correct cavity phase
Dold = 1.0e3 * distOld; // m --> mm
Dold = 1.0e3 * distOld;
return flag;
}
......@@ -2317,7 +2318,7 @@ void ParallelCyclotronTracker::applyPluginElements(const double dt) {
-> checkStripper(itsBunch_m, turnnumber_m, itsBunch_m->getT() * 1e9 /*[ns]*/, dt);
if(flag_stripper) {
itsBunch_m->updateNumTotal();
//INFOMSG("* Total number of particles collide with stripper = " << itsBunch_m->getTotalNum() << endl;);
//INFOMSG(level2 << "* Total number of particles collide with stripper = " << itsBunch_m->getTotalNum() << endl;);
}
}
......
This diff is collapsed.
......@@ -3,15 +3,11 @@
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: BeamStripping
// Defines the abstract interface for a beam BeamStripping.
// *** MISSING *** BeamStripping interface is still incomplete.
// Defines the abstract interface for a beam BeamStripping
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
......@@ -32,6 +28,35 @@
class BeamlineVisitor;
class LossDataSink;
struct PFieldData {
std::string filename;
// known from file: field and three theta derivatives
std::vector<double> pfld; //Pz
// Grid-Size
//need to be read from inputfile.
int nrad, ntet;
// one more grid line is stored in azimuthal direction:
int ntetS;
// total grid points number.
int ntot;
};
struct PPositions {
// these 4 parameters are need to be read from field file.
double rmin, delr;
double tetmin, dtet;
// Radii and step width of initial Grid
std::vector<double> rarr;
double Pfact; // MULTIPLICATION FACTOR FOR PRESSURE MAP
};
// Class BeamStripping
// ------------------------------------------------------------------------
......@@ -51,7 +76,7 @@ public:
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B);
virtual bool checkBeamStripping(PartBunchBase<double, 3> *bunch, Cyclotron* cycl, const int turnnumber, const double t, const double tstep);
......@@ -59,7 +84,7 @@ public:
virtual void initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField);
virtual void initialise(PartBunchBase<double, 3> *bunch);
virtual void initialise(PartBunchBase<double, 3> *bunch, const double &scaleFactor);
virtual void finalise();
......@@ -75,15 +100,15 @@ public:
void print();
string getBeamStrippingShape();
void setOutputFN(string fn);
string getOutputFN();
std::string getBeamStrippingShape();
void setOutputFN(std::string fn);
std::string getOutputFN();
unsigned int getLosses() const;
int checkPoint(const double &x, const double &y, const double &z);
// --------Cyclotron beam stripping
double checkPressure(const double &x, const double &y);
void setPressure(double pressure) ;
double getPressure() const;
......@@ -91,32 +116,58 @@ public:
void setTemperature(double temperature) ;
double getTemperature() const;
void setPressureMapFN(std::string pmapfn);
virtual std::string getPressureMapFN() const;
void setPScale(double ps);
virtual double getPScale() const;
void setResidualGas(std::string gas);
virtual std::string getResidualGas() const;
void setStop(bool stopflag);
virtual bool getStop() const;
protected:
private:
void initR(double rmin, double dr, int nrad);
void getPressureFromFile(const double &scaleFactor);
inline int idx(int irad, int ktet) {return (ktet + PField.ntetS * irad);}
// Not implemented.
void operator=(const BeamStripping &);
private:
string filename_m; /**< The name of the outputfile*/
std::string filename_m; /**< The name of the outputfile*/
bool informed_m;
///parameters for BeamStripping
double pressure_m; ///< mbar
///@{ parameters for BeamStripping
std::string gas_m;
double pressure_m;
std::string pmapfn_m; /// stores the filename of the pressure map
double pscale_m; /// a scale factor for the P-field
double temperature_m; ///< K
double stop_m;
///@}
///@{ size limits took from cyclotron
double minr_m;
double maxr_m;
double minz_m;
double maxz_m;
///@}
unsigned int losses_m;
std::unique_ptr<LossDataSink> lossDs_m;
ParticleMatterInteractionHandler *parmatint_m;
protected:
// object of Matrices including pressure field map and its derivates
PFieldData PField;
// object of parameters about the map grid
PPositions PP;
};
inline
......
......@@ -297,21 +297,25 @@ bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t &
if (zpos > maxz_m || zpos < minz_m || rpos > maxr_m || rpos < minr_m){
flagNeedUpdate = true;
Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
gmsgALL << getName() << ": particle "<< id <<" out of the global aperture of cyclotron!"<< endl;
gmsgALL << getName() << ": Coords: "<< RefPartBunch_m->R[id] << endl;
gmsgALL << level2 << getName()
<< ": particle " << id
<< " out of the global aperture of cyclotron!" << endl;
gmsgALL << level2 << getName()
<< ": Coords: "<< RefPartBunch_m->R[id] << endl;
} else{
flagNeedUpdate = apply(RefPartBunch_m->R[id], RefPartBunch_m->P[id], t, E, B);
if(flagNeedUpdate){
Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
gmsgALL << getName() << ": particle "<< id <<" out of the field map boundary!"<< endl;
gmsgALL << getName() << ": Coords: "<< RefPartBunch_m->R[id] << endl;
gmsgALL << level2 << getName()
<< ": particle "<< id
<< " out of the field map boundary!" << endl;
gmsgALL << level2 << getName()
<< ": Coords: "<< RefPartBunch_m->R[id] << endl;
}
}
if (flagNeedUpdate) {
lossDs_m->addParticle(RefPartBunch_m->R[id], RefPartBunch_m->P[id], id, t);
lossDs_m->addParticle(RefPartBunch_m->R[id]*1000, RefPartBunch_m->P[id], id, t);
RefPartBunch_m->Bin[id] = -1;
}
......
......@@ -291,12 +291,16 @@ bool Stripper::checkStripper(PartBunchBase<double, 3> *bunch, const int turnnum
strippoint(1) = (A_m*A_m*bunch->R[i](1) - A_m*B_m*bunch->R[i](0)-B_m*C_m)/(R_m*R_m);
strippoint(2) = bunch->R[i](2);
lossDs_m->addParticle(strippoint, bunch->P[i], bunch->ID[i], t+dt, turnnumber);
*gmsg << level4 << "* Particle " << bunch->ID[i] << " is deleted by stripper " << getName() << endl;
*gmsg << level4 << "* strippoint " << strippoint << endl;
if (stop_m) {
bunch->Bin[i] = -1;
flagNeedUpdate = true;
//*gmsg << level4 << "* Particle " << bunch->ID[i] << " is deleted by stripper " << getName() << endl;
}else{
//*gmsg << level4 << "* Total number of particles after stripper " << getName() << " = " << bunch->getTotalNum() << endl;
flagNeedUpdate = true;
// change charge and mass of PartData when the reference particle hits the stripper.
if(bunch->ID[i] == 0)
......
......@@ -84,8 +84,6 @@ void FM3DH5Block::readMap() {
#if defined (NDEBUG)
(void)h5err;
#endif
INFOMSG (level3 << typeset_msg("reading fieldmap '" + Filename_m + "'", "info") << "\n"
<< endl);
auto props = H5CreateFileProp ();
auto comm = Ippl::getComm();
......@@ -129,8 +127,7 @@ void FM3DH5Block::readMap() {
h5err = H5CloseFile(file);
assert (h5err != H5_ERR);
INFOMSG(level3 << typeset_msg("fieldmap '" + Filename_m + "' read", "info") << "\n"
<< endl);
INFOMSG(level3 << typeset_msg("fieldmap '" + Filename_m + "' read", "info") << endl);
}
......@@ -143,8 +140,7 @@ void FM3DH5Block::freeMap() {
FieldstrengthHy_m.clear();
FieldstrengthHz_m.clear();
INFOMSG(level3 << typeset_msg("freed fieldmap '" + Filename_m + "'", "info") << "\n"
<< endl);
INFOMSG(level3 << typeset_msg("freed fieldmap '" + Filename_m + "'", "info") << endl);
}
}
......
#ifndef BEAMSTRIPPINGPHYSICS_HH
#define BEAMSTRIPPINGPHYSICS_HH
//Class:BeamStrippingPhysics
// Defines the beam stripping physics models
// ------------------------------------------------------------------------
//
// Class: BeamStrippingPhysics
// Defines the beam stripping physics models
//
// ------------------------------------------------------------------------
// Class category:
// Class category: ParticleMatterInteraction
// ------------------------------------------------------------------------
// $Date: 2018/11 $
// $Author: PedroCalvo$
......@@ -41,13 +45,14 @@ public:
BeamStrippingPhysics(const std::string &name, ElementBase *element);
~BeamStrippingPhysics();
void setCyclotron(Cyclotron* cycl) { cycl_m = cycl; };
void apply(PartBunchBase<double, 3> *bunch,
const std::pair<Vector_t, double> &boundingSphere,
size_t numParticlesInSimulation = 0);
virtual const std::string getType() const;
void print(Inform& msg);
bool stillActive();
bool stillAlive(PartBunchBase<double, 3> *bunch);
......@@ -58,25 +63,27 @@ public:
unsigned getRediffused() {return rediffusedStat_m;}
unsigned int getNumEntered() {return bunchToMatStat_m;}
inline void doPhysics(PartBunchBase<double, 3> *bunch);
void setCyclotron(Cyclotron* cycl) { cycl_m = cycl; };
private:
void Material();
void CrossSection(const Vector_t &R, double &Eng);
void MolecularDensity(const double &pressure, const double &temperature, int &iComp);
void CrossSection(double &Eng);
double CSAnalyticFunction(double Eng, double Eth,
double a1, double a2, double a3, double a4, double a5, double a6);
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 CSChebyshevFitting(double Eng, double Emin, double Emax);
bool GasStripping(double &deltas);
bool LorentzStripping(double &gamma, double &E);
void SecondaryParticles(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(double &Eng,
const double deltat,
......@@ -84,42 +91,58 @@ private:
Cyclotron *cycl_m;
BeamStripping *bstp_m;
gsl_rng * r_m;
std::string material_m;
ElementBase::ElementType bstpshape_m;
double T_m;
gsl_rng * r_m;
double T_m;
double dT_m;
double mass_m;
double charge_m;
double m_h;
int NbComponents;
// double totalmolecularDensity_m;
double totalmolecularDensity_m;
double molecularDensity[3];
std::string gas_m;
double pressure_m;
///@{ size limits took from cyclotron
double rhoinit_m;
double rhofinal_m;
std::unique_ptr<LossDataSink> lossDs_m;
double NCS_single_all;
double NCS_double_all;
double NCS_total_all;
double NCS_a;
double NCS_b;
double NCS_c;
double NCS_total;
unsigned bunchToMatStat_m;
unsigned stoppedPartStat_m;
unsigned rediffusedStat_m;
size_t locPartsInMat_m;
static const double fMolarFraction[3];
static const double CSCoefSingle_Hminus[3][7];
static const double CSCoefDouble_Hminus[3][7];
static const double CSCoefSingle_Hminus[3][9];
static const double CSCoefDouble_Hminus[3][9];
static const double CSCoefSingle_Hplus[3][9];
static const double CSCoefDouble_Hplus[3][9];
static const double CSCoefSingleLoss_H[3][7];
static const double CSCoefSingleLoss_H[3][9];
static const double CSCoefSingleCapt_H[3][9];
static const double CSCoefHminusProduction_H_Tabata[13];
static const double CSCoefProtonProduction_H_Tabata[9];
static const double CSCoefProtonProduction_H2plus_Tabata[11];
static const double CSCoefH3plusProduction_H2plus_Tabata[7];
static const double CSCoefSingle_Hminus_Chebyshev[11];
static const double CSCoefDouble_Hminus_Chebyshev[11];
static const double CSCoefSingle_Hplus_Chebyshev[11];
static const double CSCoefDouble_Hplus_Chebyshev[11];
static const double CSCoefHydrogenProduction_H2plus_Chebyshev[11];
static double a_m[9];
static double b_m[3][9];
};
#endif //BEAMSTRIPPINGPHYSICS_HH
// ------------------------------------------------------------------------
// $RCSfile: OpalBeamStripping.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
......@@ -29,18 +27,24 @@ OpalBeamStripping::OpalBeamStripping():
OpalElement(SIZE, "BEAMSTRIPPING",
"The \"BEAMSTRIPPING\" element defines a beam stripping interaction"),
parmatint_m(NULL) {
itsAttr[PRESSURE] = Attributes::makeReal
("PRESSURE", " Pressure os the accelerator, [mbar]");
itsAttr[TEMPERATURE] = Attributes::makeReal
itsAttr[PRESSURE] = Attributes::makeReal
("PRESSURE", " Pressure in the accelerator, [mbar]");
itsAttr[TEMPERATURE] = Attributes::makeReal
("TEMPERATURE", " Temperature of the accelerator, [K]");
itsAttr[STOP] = Attributes::makeBool
itsAttr[PMAPFN] = Attributes::makeString
("PMAPFN", "Filename for the Pressure fieldmap");
itsAttr[PSCALE] = Attributes::makeReal
("PSCALE", "Scale factor for the P-field", 1.0);
itsAttr[GAS] = Attributes::makeString
("GAS", "The composition of residual gas");
itsAttr[STOP] = Attributes::makeBool
("STOP", "Option Whether stop tracking after beam stripping. Default: true", true);
itsAttr[OUTFN] = Attributes::makeString
("OUTFN", "BeamStripping output filename");
registerRealAttribute("PRESSURE");
registerRealAttribute("TEMPERATURE");
registerStringAttribute("OUTFN");
registerStringAttribute("PMAPFN");
registerRealAttribute("PSCALE");
registerStringAttribute("GAS");
registerOwnership();
......@@ -78,13 +82,17 @@ void OpalBeamStripping::update() {
double pressure = Attributes::getReal(itsAttr[PRESSURE]);
double temperature = Attributes::getReal(itsAttr[TEMPERATURE]);
std::string pmap = Attributes::getString(itsAttr[PMAPFN]);
double pscale = Attributes::getReal(itsAttr[PSCALE]);
bool stop = Attributes::getBool(itsAttr[STOP]);
std::string gas = Attributes::getString(itsAttr[GAS]);
bstp->setPressure(pressure);
bstp->setTemperature(temperature);
bstp->setPressureMapFN(pmap);
bstp->setPScale(pscale);
bstp->setStop(stop);
bstp->setOutputFN(Attributes::getString(itsAttr[OUTFN]));
bstp->setResidualGas(gas);
if(itsAttr[PARTICLEMATTERINTERACTION] && parmatint_m == NULL) {
parmatint_m = (ParticleMatterInteraction::find(Attributes::getString(itsAttr[PARTICLEMATTERINTERACTION])))->clone(getOpalName() + std::string("_parmatint"));
......
#ifndef OPAL_OpalBeamStripping_HH
#define OPAL_OpalBeamStripping_HH
// ------------------------------------------------------------------------
// $RCSfile: OpalSlit.h,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: OpalBeamStripping
// ------------------------------------------------------------------------
// $Date: 2018/11 $
......@@ -28,10 +23,12 @@ public:
/// The attributes of class OpalBeamStripping.
enum {
PRESSURE = COMMON,
TEMPERATURE,
STOP,
OUTFN,
PRESSURE = COMMON, // Pressure in mbar
TEMPERATURE, // Temperature of residual gas
PMAPFN, // The filename of the mid-plane pressure map
PSCALE, // A scalar to scale the P-field
GAS, // Type of gas for residual vaccum
STOP, // whether the secondary particles are tracked
SIZE
};
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment