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

Change the arguments of CCOLLIMATOR element to XSTART, XEND, YSTART, YEND,...

Change the arguments of CCOLLIMATOR element to XSTART, XEND, YSTART, YEND, WIDTH so as to model non-radial oriented collimator
parent e82c1657
......@@ -53,11 +53,11 @@ Collimator::Collimator():
b_m(0.0),
x0_m(0.0),
y0_m(0.0),
as_m(0.0),
rs_m(0.0),
ae_m(0.0),
re_m(0.0),
w_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
isAPepperPot_m(false),
isASlit_m(false),
isARColl_m(false),
......@@ -89,11 +89,11 @@ Collimator::Collimator(const Collimator &right):
b_m(right.b_m),
x0_m(right.x0_m),
y0_m(right.y0_m),
as_m(right.as_m),
rs_m(right.rs_m),
ae_m(right.ae_m),
re_m(right.re_m),
w_m(right.w_m),
xstart_m(right.xstart_m),
xend_m(right.xend_m),
ystart_m(right.ystart_m),
yend_m(right.yend_m),
width_m(right.width_m),
isAPepperPot_m(right.isAPepperPot_m),
isASlit_m(right.isASlit_m),
isARColl_m(right.isARColl_m),
......@@ -127,11 +127,11 @@ Collimator::Collimator(const string &name):
b_m(0.0),
x0_m(0.0),
y0_m(0.0),
as_m(0.0),
rs_m(0.0),
ae_m(0.0),
re_m(0.0),
w_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
isAPepperPot_m(false),
isASlit_m(false),
isARColl_m(false),
......@@ -244,27 +244,29 @@ bool Collimator::checkCollimator(PartBunch &bunch, const int turnnumber, const d
bool flagNeedUpdate = false;
Vector_t rmin, rmax;
bunch.get_bounds(rmin, rmax);
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
if(r1 > rs_m - 50.0 && r1 < re_m + 50.0 ){
size_t tempnum = bunch.getLocalNum();
int pflag = 0;
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch.PType[i] == 0) {
pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
if(pflag != 0) {
lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
bunch.Bin[i] = -1;
flagNeedUpdate = true;
}
}
}
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
size_t tempnum = bunch.getLocalNum();
int pflag = 0;
for(unsigned int i = 0; i < tempnum; ++i) {
if(bunch.PType[i] == 0) {
pflag = checkPoint(bunch.R[i](0), bunch.R[i](1));
if(pflag != 0) {
lossDs_m->addParticle(bunch.R[i], bunch.P[i], bunch.ID[i]);
bunch.Bin[i] = -1;
flagNeedUpdate = true;
}
}
}
}
reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
if(flagNeedUpdate) lossDs_m->save(getName());
return flagNeedUpdate;
reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
if(flagNeedUpdate) lossDs_m->save(getName());
return flagNeedUpdate;
}
......@@ -318,7 +320,7 @@ void Collimator::goOnline() {
else if(isARColl_m)
msg << "RCollimator a= " << getXsize() << " b= " << b_m << " start= " << position_m << " fn= " << filename_m << " ny= " << nHolesY_m << " pitch= " << pitch_m << endl;
else if(isACColl_m)
msg << "CCollimator angstart= " << as_m << " angend " << ae_m << " rstart " << rs_m << " rend " << re_m << endl;
msg << "CCollimator angstart= " << xstart_m << " angend " << ystart_m << " rstart " << xend_m << " rend " << yend_m << endl;
else if(isAWire_m)
msg << "Wire x= " << x0_m << " y= " << y0_m << endl;
else
......@@ -539,44 +541,44 @@ double Collimator::getYpos() {
// --------Cyclotron collimator
}
void Collimator::setAngStart(double as) {
as_m = as;
void Collimator::setXStart(double xstart) {
xstart_m = xstart;
}
void Collimator::setRStart(double rs) {
rs_m = rs;
void Collimator::setXEnd(double xend) {
xend_m = xend;
}
void Collimator::setAngEnd(double ae) {
ae_m = ae;
void Collimator::setYStart(double ystart) {
ystart_m = ystart;
}
void Collimator::setREnd(double re) {
re_m = re;
void Collimator::setYEnd(double yend) {
yend_m = yend;
}
void Collimator::setWidth(double w) {
w_m = w;
void Collimator::setWidth(double width) {
width_m = width;
}
double Collimator::getAngStart() {
return as_m;
double Collimator::getXStart() {
return xstart_m;
}
double Collimator::getRStart() {
return rs_m;
double Collimator::getXEnd() {
return xend_m;
}
double Collimator::getAngEnd() {
return ae_m;
double Collimator::getYStart() {
return ystart_m;
}
double Collimator::getREnd() {
return re_m;
double Collimator::getYEnd() {
return yend_m;
}
double Collimator::getWidth() {
return w_m;
return width_m;
}
//-------------------------------
......@@ -644,22 +646,30 @@ string Collimator::getCollimatorShape() {
void Collimator::setGeom() {
geom_m[0].x = rs_m*cos(as_m);
geom_m[0].y = rs_m*sin(as_m);
double slope;
if (xend_m == xstart_m)
slope = 1.0e12;
else
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
double coeff2 = sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
double halfdist = width_m / 2.0;
geom_m[0].x = xstart_m - halfdist * coeff1;
geom_m[0].y = ystart_m + halfdist / coeff2;
geom_m[1].x = re_m*cos(as_m);
geom_m[1].y = re_m*sin(as_m);
geom_m[1].x = xstart_m + halfdist * coeff1;
geom_m[1].y = ystart_m - halfdist / coeff2;
geom_m[2].x = re_m*cos(ae_m);
geom_m[2].y = re_m*sin(ae_m);
geom_m[2].x = xend_m + halfdist * coeff1;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[3].x = rs_m*cos(ae_m);
geom_m[3].y = rs_m*sin(ae_m);
geom_m[3].x = xend_m - halfdist * coeff1;
geom_m[3].y = yend_m + halfdist / coeff2;
geom_m[4].x = geom_m[0].x;
geom_m[4].y = geom_m[0].y;
}
......
......@@ -115,23 +115,23 @@ public:
// --------Cyclotron collimator
void setAngStart(double as) ;
void setXStart(double xstart) ;
void setRStart(double rs) ;
void setYStart(double ystart) ;
void setAngEnd(double ae) ;
void setXEnd(double xend) ;
void setREnd(double re) ;
void setYEnd(double yend) ;
void setWidth(double w) ;
void setWidth(double width) ;
double getAngStart() ;
double getXStart() ;
double getRStart() ;
double getYStart() ;
double getAngEnd() ;
double getXEnd() ;
double getREnd() ;
double getYEnd() ;
double getWidth() ;
......@@ -170,11 +170,13 @@ private:
double b_m;
double x0_m;
double y0_m;
double as_m;
double rs_m;
double ae_m;
double re_m;
double w_m;
//parameters for CCollimator
double xstart_m;
double xend_m;
double ystart_m;
double yend_m;
double width_m;
/** This defines a pepperpot */
......
......@@ -187,17 +187,20 @@ void Probe::setGeom(const double dist) {
else
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
geom_m[0].x = xstart_m - dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[0].y = ystart_m + dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
double coeff2 = sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
double halfdist = dist / 2.0;
geom_m[0].x = xstart_m - halfdist * coeff1;
geom_m[0].y = ystart_m + halfdist / coeff2;
geom_m[1].x = xstart_m + dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[1].y = ystart_m - dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
geom_m[1].x = xstart_m + halfdist * coeff1;
geom_m[1].y = ystart_m - halfdist / coeff2;
geom_m[2].x = xend_m + dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[2].y = yend_m - dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
geom_m[2].x = xend_m + halfdist * coeff1;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[3].x = xend_m - dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[3].y = yend_m + dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
geom_m[3].x = xend_m - halfdist * coeff1;
geom_m[3].y = yend_m + halfdist / coeff2;
geom_m[4].x = geom_m[0].x;
geom_m[4].y = geom_m[0].y;
......@@ -211,7 +214,7 @@ bool Probe::checkProbe(PartBunch &bunch, const int turnnumber, const double t,
double r_start = sqrt(xstart_m * xstart_m + ystart_m * ystart_m);
double r_end = sqrt(xend_m * xend_m + yend_m * yend_m);
double r1 = sqrt(rmax(0) * rmax(0) + rmax(1) * rmax(1));
if( r1 > r_start - 10.0 && r1 < r_end + 10.0 ){
size_t tempnum = bunch.getLocalNum();
......
......@@ -99,26 +99,30 @@ Stripper::Stripper(const string &name):
void Stripper::setGeom(const double dist) {
double slope;
double slope;
if (xend_m == xstart_m)
slope = 1.0e12;
else
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
geom_m[0].x = xstart_m - dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[0].y = ystart_m + dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
double coeff2 = sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
double halfdist = dist / 2.0;
geom_m[0].x = xstart_m - halfdist * coeff1;
geom_m[0].y = ystart_m + halfdist / coeff2;
geom_m[1].x = xstart_m + dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[1].y = ystart_m - dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
geom_m[1].x = xstart_m + halfdist * coeff1;
geom_m[1].y = ystart_m - halfdist / coeff2;
geom_m[2].x = xend_m + dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[2].y = yend_m - dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
geom_m[2].x = xend_m + halfdist * coeff1;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[3].x = xend_m - dist / 2.0 * slope / sqrt(1 + slope * slope);
geom_m[3].y = yend_m + dist / 2.0 * 1.0 / sqrt(1 + slope * slope);
geom_m[3].x = xend_m - halfdist * coeff1;
geom_m[3].y = yend_m + halfdist / coeff2;
geom_m[4].x = geom_m[0].x;
geom_m[4].y = geom_m[0].y;
}
......
......@@ -4,7 +4,7 @@
// Class category:
// ------------------------------------------------------------------------
// $Date: 2009/07/20 09:32:31 $
// $Author: bi $
// $Author: Bi, Yang $
//-------------------------------------------------------------------------
#include "Solvers/CollimatorPhysics.hh"
#include "Physics/Physics.h"
......@@ -36,10 +36,10 @@ CollimatorPhysics::CollimatorPhysics(const string &name, ElementBase *element, c
b_m(minor),
xp_m(0.0),
yp_m(0.0),
angstart_m(0.0),
angend_m(0.0),
rstart_m(0.0),
rend_m(0.0),
xstart_m(0.0),
xend_m(0.0),
ystart_m(0.0),
yend_m(0.0),
width_m(0.0),
Begin_m(0.0),
End_m(0.0),
......@@ -64,11 +64,12 @@ CollimatorPhysics::CollimatorPhysics(const string &name, ElementBase *element, c
collshape_m = coll->getCollimatorShape();
xp_m = coll->getXpos();
yp_m = coll->getYpos();
angstart_m = coll->getAngStart();
rstart_m = coll->getRStart();
angend_m = coll->getAngEnd();
rend_m = coll->getREnd();
xstart_m = coll->getXStart();
ystart_m = coll->getYStart();
xend_m = coll->getXEnd();
yend_m = coll->getYEnd();
width_m = coll->getWidth();
setCColimatorGeom();
} else if(dynamic_cast<Drift *>(element_ref_m)) {
Drift *drf = dynamic_cast<Drift *>(element_ref_m);
......@@ -113,13 +114,10 @@ void CollimatorPhysics::apply(PartBunch &bunch) {
double scalefactor;
// int flagcoll; !! carefull not initialized !!
// Vector_t tmploss;
double angpar, rpar;
double deltatbunch = bunch.getdT();
double deltat = deltatbunch;
double rr = 0.0, rr1 = 0.0, rr2 = 0.0;
if(collshape_m == "CCollimator") {
scalefactor = 1;
} else {
......@@ -140,10 +138,11 @@ void CollimatorPhysics::apply(PartBunch &bunch) {
Vector_t &P = Pincol_m[i];
if(collshape_m == "CCollimator") {
angpar = calculateAngle(R(0), R(1));
rpar = sqrt(R(0) * R(0) + R(1) * R(1));
if(angpar > angstart_m && angpar < angend_m && rpar > rstart_m && rpar < rend_m && fabs(R(2)) > a_m) incoll_m = true;
if( checkPoint(R(0), R(1)) == 1 )
incoll_m = true;
else
incoll_m = false;
} else {
if(collshape_m == "Slit") {
......@@ -236,11 +235,11 @@ void CollimatorPhysics::apply(PartBunch &bunch) {
Vector_t R = bunch.R[i];
Vector_t P = bunch.P[i];
if(collshape_m == "CCollimator") {
angpar = calculateAngle(R(0), R(1));
rpar = sqrt(R(0) * R(0) + R(1) * R(1));
if(angpar > angstart_m && angpar < angend_m && rpar > rstart_m && rpar < rend_m && fabs(R(2)) > a_m) incoll_m = true;
if(collshape_m == "CCollimator") {
if( checkPoint(R(0), R(1)) == 1 )
incoll_m = true;
else
incoll_m = false;
} else {
if(collshape_m == "Slit") {
if(a_m > 0) {
......@@ -339,9 +338,10 @@ void CollimatorPhysics::apply(PartBunch &bunch) {
Vector_t &P = Pincol_m[i];
if(collshape_m == "CCollimator") {
angpar = calculateAngle(R(0), R(1));
rpar = sqrt(R(0) * R(0) + R(1) * R(1));
if(angpar > angstart_m && angpar < angend_m && rpar > rstart_m && rpar < rend_m && fabs(R(2)) > a_m) incoll_m = true;
if( checkPoint(R(0), R(1)) == 1 )
incoll_m = true;
else
incoll_m = false;
} else {
if(collshape_m == "Slit") {
if(a_m > 0) {
......@@ -584,16 +584,47 @@ void CollimatorPhysics::CoulombScat(Vector_t &R, Vector_t &P, double &deltat, d
}
}
void CollimatorPhysics::setCColimatorGeom() {
double slope;
if (xend_m == xstart_m)
slope = 1.0e12;
else
slope = (yend_m - ystart_m) / (xend_m - xstart_m);
double coeff2 = sqrt(1 + slope * slope);
double coeff1 = slope / coeff2;
double halfdist = width_m / 2.0;
geom_m[0].x = xstart_m - halfdist * coeff1;
geom_m[0].y = ystart_m + halfdist / coeff2;
geom_m[1].x = xstart_m + halfdist * coeff1;
geom_m[1].y = ystart_m - halfdist / coeff2;
geom_m[2].x = xend_m + halfdist * coeff1;
geom_m[2].y = yend_m - halfdist / coeff2;
geom_m[3].x = xend_m - halfdist * coeff1;
geom_m[3].y = yend_m + halfdist / coeff2;
geom_m[4].x = geom_m[0].x;
geom_m[4].y = geom_m[0].y;
}
double CollimatorPhysics::calculateAngle(double x, double y) {
double thetaXY = atan2(y, x);
int CollimatorPhysics::checkPoint(const double &x, const double &y) {
int cn = 0;
// if(x < 0) thetaXY = pi + atan(y / x);
// else if((x > 0) && (y >= 0)) thetaXY = atan(y / x);
// else if((x > 0) && (y < 0)) thetaXY = 2.0 * pi + atan(y / x);
// else if((x == 0) && (y > 0)) thetaXY = pi / 2.0;
// else if((x == 0) && (y < 0)) thetaXY = 3.0 / 2.0 * pi;
for(int i = 0; i < 4; i++) {
if(((geom_m[i].y <= y) && (geom_m[i+1].y > y))
|| ((geom_m[i].y > y) && (geom_m[i+1].y <= y))) {
return thetaXY >= 0 ? thetaXY: thetaXY + Physics::two_pi;
float vt = (float)(y - geom_m[i].y) / (geom_m[i+1].y - geom_m[i].y);
if(x < geom_m[i].x + vt * (geom_m[i+1].x - geom_m[i].x))
++cn;
}
}
return (cn & 1); // 0 if even (out), and 1 if odd (in)
}
......@@ -6,13 +6,14 @@
// Class category:
// ------------------------------------------------------------------------
// $Date: 2009/07/20 09:32:31 $
// $Author: bi $
// $Author: Bi, Yang $
//-------------------------------------------------------------------------
#include <vector>
#include "Solvers/SurfacePhysicsHandler.hh"
//#include "Algorithms/PBunchDefs.h"
#include "Algorithms/Vektor.h"
#include "AbsBeamline/Component.h"
class RANLIB_class;
class ElementBase;
......@@ -35,17 +36,17 @@ public:
void EnergyLoss(double &Eng, bool &pdead, double &deltat);
void Rot(Vector_t &P, Vector_t Prot, double normP);
double Rot(double &p1, double &p2, double &scatang);
double calculateAngle(double x, double y);
private:
RANLIB_class *rGen_m;
double a_m;
double b_m;
double xp_m;
double yp_m;
double angstart_m;
double angend_m;
double rstart_m;
double rend_m;
double xstart_m;
double xend_m;
double ystart_m;
double yend_m;
double width_m;
double Begin_m;
......@@ -61,6 +62,7 @@ private:
double n_m;
bool incoll_m;
Point geom_m[5];
int index_m;
std::vector<unsigned> label_m;
......@@ -75,7 +77,9 @@ private:
std::vector<Vector_t> Efincol_m;
std::vector<double> time_m;
std::vector<int> steps_m;
void setCColimatorGeom();
int checkPoint( const double & x, const double & y );
LossDataSink *lossDs_m;
};
......
......@@ -338,25 +338,25 @@ void ParallelCyclotronTracker::visitBeamBeam(const BeamBeam &) {
*/
void ParallelCyclotronTracker::visitCollimator(const Collimator &coll) {
Inform msg("visitCCollimator ");
*gmsg << "* --------- Collimator -----------------------------" << endl;
Collimator* elptr = dynamic_cast<Collimator *>(coll.clone());
myElements.push_back(elptr);
double angstart = elptr->getAngStart();
msg << "AngStart= " << angstart << " [rad]" << endl;
double xstart = elptr->getXStart();
*gmsg << "Xstart= " << xstart << " [mm]" << endl;
double angend = elptr->getAngEnd();
msg << "AngEnd= " << angend << " [rad]" << endl;
double xend = elptr->getXEnd();
*gmsg << "Xend= " << xend << " [mm]" << endl;
double rstart = elptr->getRStart();
msg << "RStart= " << rstart << " [mm]" << endl;
double ystart = elptr->getYStart();
*gmsg << "Ystart= " << ystart << " [mm]" << endl;
double rend = elptr->getREnd();
msg << "REnd= " << rend << " [mm]" << endl;