Commit eceb7728 authored by Jianjun Yang's avatar Jianjun Yang
Browse files

- implement the global aperture of cyclotron by 4 arguments in cyclotron element:

		MINR, MAXR, MINZ, MAXZ. 
- if a particle is out of cyclotron aperture or out of the field map, 
   it will be deleted after this step and its parameters are recorded in the *.loss file
- change data type of particle id from int to size_t
- many minor cosmetics of display massage
parent 1f2a3f2c
......@@ -48,11 +48,11 @@ void BeamBeam::accept(BeamlineVisitor &visitor) const {
visitor.visitBeamBeam(*this);
}
bool BeamBeam::apply(const int &i, const double &t, double E[], double B[]) {
bool BeamBeam::apply(const size_t &i, const double &t, double E[], double B[]) {
return false;
}
bool BeamBeam::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool BeamBeam::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
......
......@@ -63,9 +63,9 @@ public:
// Units are metres.
virtual const Vector3D &getBunchDisplacement() const = 0;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -156,12 +156,12 @@ void Collimator::accept(BeamlineVisitor &visitor) const {
visitor.visitCollimator(*this);
}
bool Collimator::apply(const int &i, const double &t, double E[], double B[]) {
bool Collimator::apply(const size_t &i, const double &t, double E[], double B[]) {
Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
return apply(i, t, Ev, Bv);
}
bool Collimator::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Collimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
const Vector_t &R = RefPartBunch_m->R[i] - Vector_t(dx_m, dy_m, ds_m); // including the missaligment
const Vector_t &P = RefPartBunch_m->P[i];
const double recpgamma = Physics::c * RefPartBunch_m->getdT() / sqrt(1.0 + dot(P, P));
......
......@@ -69,9 +69,9 @@ public:
/// Return the vertical half-aperture.
virtual double getYsize() const {return b_m;}
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -102,9 +102,9 @@ public:
virtual void addKT(int i, double t, Vector_t &K) {};
virtual bool apply(const int &i, const double &t, double E[], double B[]) = 0;
virtual bool apply(const size_t &i, const double &t, double E[], double B[]) = 0;
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B) = 0;
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) = 0;
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) = 0;
......
......@@ -48,11 +48,11 @@ void Corrector::accept(BeamlineVisitor &visitor) const {
visitor.visitCorrector(*this);
}
bool Corrector::apply(const int &i, const double &t, double E[], double B[]) {
bool Corrector::apply(const size_t &i, const double &t, double E[], double B[]) {
return false;
}
bool Corrector::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Corrector::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
......
......@@ -77,9 +77,9 @@ public:
/// Return the plane on which the corrector acts.
virtual Plane getPlane() const = 0;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -22,6 +22,7 @@
#include "Algorithms/PartBunch.h"
#include "AbsBeamline/BeamlineVisitor.h"
#include "Physics/Physics.h"
#include "Structure/LossDataSink.h"
#include "Fields/Fieldmap.hh"
#include "Utilities/OpalException.h"
#include <fstream>
......@@ -61,6 +62,10 @@ Cyclotron::Cyclotron(const Cyclotron &right):
tcr2_m(right.tcr2_m),
mbtc_m(right.mbtc_m),
slptc_m(right.slptc_m),
minr_m(right.minr_m),
maxr_m(right.maxr_m),
minz_m(right.minz_m),
maxz_m(right.maxz_m),
RFfilename_m(right.RFfilename_m) {
}
......@@ -241,9 +246,10 @@ double Cyclotron::getMaxZ() const {
return maxz_m;
}
bool Cyclotron::apply(const int &i, const double &t, double E[], double B[]) {
bool Cyclotron::apply(const size_t &id, const double &t, double E[], double B[]) {
Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
if(apply(RefPartBunch_m->R[i], Vector_t(0.0), t, Ev, Bv)) return true;
if(apply(id, t, Ev, Bv)) return true;
E[0] = Ev(0);
E[1] = Ev(1);
......@@ -255,8 +261,33 @@ bool Cyclotron::apply(const int &i, const double &t, double E[], double B[]) {
return false;
}
bool Cyclotron::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
return apply(RefPartBunch_m->R[i], Vector_t(0.0), t, E, B);
bool Cyclotron::apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B) {
bool flagNeedUpdate = false;
const double rpos = sqrt(RefPartBunch_m->R[id](0) * RefPartBunch_m->R[id](0)
+ RefPartBunch_m->R[id](1) * RefPartBunch_m->R[id](1));
const double zpos = RefPartBunch_m->R[id](2);
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;
} else{
flagNeedUpdate = apply(RefPartBunch_m->R[id], Vector_t(0.0), t, E, B);
if(flagNeedUpdate){
Inform gmsgALL("OPAL ", INFORM_ALL_NODES);
gmsgALL<<getName() << ": particle "<< id <<" out of the field map boundary!"<< endl;
}
}
if (flagNeedUpdate) {
lossDs_m->addParticle(RefPartBunch_m->R[id], RefPartBunch_m->P[id],id);
lossDs_m->save(getName());
RefPartBunch_m->Bin[id] = -1;
}
return flagNeedUpdate;
}
bool Cyclotron::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {
......@@ -415,12 +446,9 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &centroid, const double
B[0] = br * cos(tet_rad) - bt * sin(tet_rad);
B[1] = br * sin(tet_rad) + bt * cos(tet_rad);
B[2] = bz;
} else {
ERRORMSG("Error!" << getName() << ".getFieldstrength(): out of boudaries (z,r) = (" << R[0] << "," << R[1] << "," << R[2] << ")" << endl
<< " rad=" << rad << " [mm], theta=" << tet << " [deg], it=" << it << ", ir=" << ir << endl);
return true;
return true;
}
if(myBFieldType_m == BANDRF) {
//The RF field is suppose to be sampled on a cartesian grid
......@@ -431,40 +459,34 @@ bool Cyclotron::apply(const Vector_t &R, const Vector_t &centroid, const double
double xBegin(0), xEnd(0), yBegin(0), yEnd(0), zBegin(0), zEnd(0);
int fcount = 0;
for(; fi != RFfields_m.end(); ++fi, ++rffi, ++rfphii, ++escali, ++fcount) {
(*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
if (R(0) >= xBegin && R(0) <= xEnd && R(1) >= yBegin && R(1) <= yEnd && R(2) >= zBegin && R(2) <= zEnd) {
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if(!(*fi)->getFieldstrength(R, tmpE, tmpB)) {
double phase = 2.0 * pi * 1E-3 * (*rffi) * t + *rfphii;
double ebscale = *escali;
E += ebscale * cos(phase) * tmpE;
B -= ebscale * sin(phase) * tmpB;
//INFOMSG("Field " << fcount << " BANDRF E= " << tmpE << " R= " << R << " phase " << phase << endl);
}
if (!superpose_m) break;
}
(*fi)->getFieldDimensions(xBegin, xEnd, yBegin, yEnd, zBegin, zEnd);
if (R(0) >= xBegin && R(0) <= xEnd && R(1) >= yBegin && R(1) <= yEnd && R(2) >= zBegin && R(2) <= zEnd) {
Vector_t tmpE(0.0, 0.0, 0.0), tmpB(0.0, 0.0, 0.0);
if(!(*fi)->getFieldstrength(R, tmpE, tmpB)) {
double phase = 2.0 * pi * 1E-3 * (*rffi) * t + *rfphii;
double ebscale = *escali;
E += ebscale * cos(phase) * tmpE;
B -= ebscale * sin(phase) * tmpB;
//INFOMSG("Field " << fcount << " BANDRF E= " << tmpE << " R= " << R << " phase " << phase << endl);
}
if (!superpose_m) break;
}
}
}
return false;
}
void Cyclotron::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
RefPartBunch_m = bunch;
online_m = true;
}
void Cyclotron::finalise() {
online_m = false;
*gmsg << "Finalize cyclotron" << endl;
if(lossDs_m)
delete lossDs_m;
}
bool Cyclotron::bends() const {
return true;
}
// calculate derivatives with 5-point lagrange's formula.
double Cyclotron::gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr)
......@@ -837,9 +859,17 @@ void Cyclotron::initR(double rmin, double dr, int nrad) {
BP.delr = dr;
}
void Cyclotron::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
RefPartBunch_m = bunch;
online_m = true;
}
void Cyclotron::initialise(PartBunch *bunch, const int &fieldflag, const double &scaleFactor) {
RefPartBunch_m = bunch;
bool doH5 = false;
lossDs_m = new LossDataSink(bunch->getTotalNum(), doH5);
// PSIBF, AVFEQBF, ANSYSBF, FFAGBF
// for your own format field, you should add your own getFieldFromFile() function by yourself.
......
......@@ -26,6 +26,7 @@
#include "Fields/BMultipoleField.h"
class Fieldmap;
class LossDataSink;
enum BFieldType {PSIBF, CARBONBF,ANSYSBF,AVFEQBF, FFAGBF,BANDRF};
......@@ -186,9 +187,9 @@ public:
virtual double getMaxZ() const;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &id, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &id, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......@@ -250,6 +251,9 @@ private:
std::vector<Fieldmap *> RFfields_m;
std::vector<std::string> RFfilename_m;
// handling for store the particle out of region
LossDataSink *lossDs_m;
void getdiffs();
double gutdf5d(double *f, double dx, const int kor, const int krl, const int lpr);
......
......@@ -102,7 +102,7 @@ bool CyclotronValley::getFast() const {
return fast_m;
}
bool CyclotronValley::apply(const int &i, const double &t, double E[], double B[]) {
bool CyclotronValley::apply(const size_t &i, const double &t, double E[], double B[]) {
Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
Vector_t Rt(RefPartBunch_m->getX(i), RefPartBunch_m->getY(i), RefPartBunch_m->getZ(i));
if(apply(Rt, Vector_t(0.0), t, Ev, Bv)) return true;
......@@ -117,7 +117,7 @@ bool CyclotronValley::apply(const int &i, const double &t, double E[], double B[
return false;
}
bool CyclotronValley::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool CyclotronValley::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
const Vector_t tmpR(RefPartBunch_m->getX(i) - dx_m, RefPartBunch_m->getY(i) - dy_m , RefPartBunch_m->getZ(i) - startField_m - ds_m);
......
......@@ -58,9 +58,9 @@ public:
bool getFast() const;
const string &getType() const;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -48,11 +48,11 @@ void Diagnostic::accept(BeamlineVisitor &visitor) const {
visitor.visitDiagnostic(*this);
}
bool Diagnostic::apply(const int &i, const double &t, double E[], double B[]) {
bool Diagnostic::apply(const size_t &i, const double &t, double E[], double B[]) {
return false;
}
bool Diagnostic::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Diagnostic::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
......
......@@ -43,9 +43,9 @@ public:
/// Apply visitor to Diagnostic.
virtual void accept(BeamlineVisitor &) const;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -50,11 +50,11 @@ void Drift::accept(BeamlineVisitor &visitor) const {
visitor.visitDrift(*this);
}
bool Drift::apply(const int &i, const double &t, double E[], double B[]) {
bool Drift::apply(const size_t &i, const double &t, double E[], double B[]) {
return false;
}
bool Drift::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Drift::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
......
......@@ -44,9 +44,9 @@ public:
/// Apply visitor to Drift.
virtual void accept(BeamlineVisitor &) const;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -48,11 +48,11 @@ void Lambertson::accept(BeamlineVisitor &visitor) const {
visitor.visitLambertson(*this);
}
bool Lambertson::apply(const int &i, const double &t, double E[], double B[]) {
bool Lambertson::apply(const size_t &i, const double &t, double E[], double B[]) {
return false;
}
bool Lambertson::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Lambertson::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
......
......@@ -44,9 +44,9 @@ public:
/// Apply visitor to Lambertson.
virtual void accept(BeamlineVisitor &) const;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -49,11 +49,11 @@ void Marker::accept(BeamlineVisitor &visitor) const {
visitor.visitMarker(*this);
}
bool Marker::apply(const int &i, const double &t, double E[], double B[]) {
bool Marker::apply(const size_t &i, const double &t, double E[], double B[]) {
return false;
}
bool Marker::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Marker::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
return false;
}
......
......@@ -43,9 +43,9 @@ public:
/// Apply visitor to Marker.
virtual void accept(BeamlineVisitor &) const;
virtual bool apply(const int &i, const double &t, double E[], double B[]);
virtual bool apply(const size_t &i, const double &t, double E[], double B[]);
virtual bool apply(const int &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B);
virtual bool apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B);
......
......@@ -90,12 +90,12 @@ void Monitor::accept(BeamlineVisitor &visitor) const {
visitor.visitMonitor(*this);
}
bool Monitor::apply(const int &i, const double &t, double E[], double B[]) {
bool Monitor::apply(const size_t &i, const double &t, double E[], double B[]) {
Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
return apply(i, t, Ev, Bv);
}
bool Monitor::apply(const int &i, const double &t, Vector_t &E, Vector_t &B) {
bool Monitor::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
const Vector_t &R = RefPartBunch_m->R[i];
const Vector_t &P = RefPartBunch_m->P[i];
const double recpgamma = Physics::c * RefPartBunch_m->getdT() / sqrt(1.0 + dot(P, P));
......
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