114 bool revBeam,
bool revTrack,
115 int maxSTEPS,
int timeIntegrator,
118 const double& mbPara,
119 const std::string& mbMode,
120 const std::string& mbBinning)
121 :
Tracker(beamline, bunch, reference, revBeam, revTrack)
123 , maxSteps_m(maxSTEPS)
124 , lastDumpedStep_m(0)
125 , myNode_m(
Ippl::myNode())
126 , initialLocalNum_m(bunch->getLocalNum())
127 , initialTotalNum_m(bunch->getTotalNum())
128 , opalRing_m(nullptr)
129 , itsStepper_mp(nullptr)
134 if ( numBunch > 1 ) {
137 mbPara, mbMode, mbBinning)
158 if ( timeIntegrator == 0 ) {
160 }
else if ( timeIntegrator == 1) {
162 }
else if ( timeIntegrator == 2) {
189 Inform msg(
"bgf_main_collision_test ");
219 const short& bunchNr)
221 if ( prevAzimuth < 0.0 ) {
226 azimuth = theta + plus;
228 double dtheta = theta - prevAzimuth;
232 if ( dtheta > 180 ) {
243 return 1000.0 *
std::sqrt(meanR(0) * meanR(0) + meanR(1) * meanR(1));
251 std::vector<double> dotP(dl.size());
258 allreduce(dotP.data(), dotP.size(), std::plus<double>());
261 double sum = std::accumulate(dotP.begin(), dotP.end() - 1, 0);
265 for (
short b = 0; b < (short)dotP.size() - 1; ++b) {
277 for (
size_t i = 0; i < dotP.size(); ++i) {
278 double const gamma =
std::sqrt(1.0 + dotP[i]);
279 double const beta =
std::sqrt(dotP[i]) / gamma;
280 dl[i] =
c_mmtns * dt * 1.0e-3 * beta;
293 for (
unsigned int i=0; i<numFiles; i++) {
294 std::string SfileName2 = SfileName;
296 SfileName2 += std::string(
"-Angle" + std::to_string(
int(i)) +
".dat");
300 SfileName2 += std::string(
"-afterEachTurn.dat");
303 outfTheta_m.emplace_back(
new std::ofstream(SfileName2.c_str()));
305 outfTheta_m.back()->setf(std::ios::scientific, std::ios::floatfield);
307 <<
"theta [deg] beta_theta*gamma "
329 *gmsg <<
"* ----------------------------- Adding Ring ------------------------------ *" <<
endl;
343 if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
344 throw OpalException(
"Error in ParallelCyclotronTracker::visitRing",
345 "PHIINIT is out of [-180, 180)!");
360 double BcParameter[8] = {};
365 *gmsg <<
"* Initial beam radius = " <<
referenceR <<
" [mm] " <<
endl;
368 *gmsg <<
"* Total reference momentum = " <<
referencePtot * 1000.0
370 *gmsg <<
"* Reference azimuthal momentum = " <<
referencePt * 1000.0
372 *gmsg <<
"* Reference radial momentum = " << referencePr * 1000.0
387 *gmsg <<
"* -------------------------- Adding Cyclotron ---------------------------- *" <<
endl;
398 *gmsg << endl <<
"* This is a Spiral Inflector Simulation! This means the following:" <<
endl;
399 *gmsg <<
"* 1.) It is up to the user to provide appropriate geometry, electric and magnetic fields!" <<
endl;
400 *gmsg <<
"* (Use BANDRF type cyclotron and use RFMAPFN to load both magnetic" <<
endl;
401 *gmsg <<
"* and electric fields, setting SUPERPOSE to an array of TRUE values.)" <<
endl;
402 *gmsg <<
"* 2.) For high currents it is strongly recommended to use the SAAMG fieldsolver," <<
endl;
403 *gmsg <<
"* FFT does not give the correct results (boundary conditions are missing)." <<
endl;
404 *gmsg <<
"* 3.) The whole geometry will be meshed and used for the fieldsolver." <<
endl;
405 *gmsg <<
"* There will be no transformations of the bunch into a local frame und consequently," <<
endl;
406 *gmsg <<
"* the problem will be treated non-relativistically!" <<
endl;
407 *gmsg <<
"* (This is not an issue for spiral inflectors as they are typically < 100 keV/amu.)" <<
endl;
408 *gmsg << endl <<
"* Note: For now, multi-bunch mode (MBM) needs to be de-activated for spiral inflector" <<
endl;
409 *gmsg <<
"* and space charge needs to be solved every time-step. numBunch_m and scSolveFreq are reset." <<
endl;
426 if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
427 throw OpalException(
"Error in ParallelCyclotronTracker::visitCyclotron",
428 "PHIINIT is out of [-180, 180)!");
439 if(insqrt > -1.0
e-10) {
445 throw OpalException(
"Error in ParallelCyclotronTracker::visitCyclotron",
465 throw OpalException(
"Error in ParallelCyclotronTracker::visitCyclotron",
466 "You are trying a local restart from a global h5 file!");
472 throw OpalException(
"Error in ParallelCyclotronTracker::visitCyclotron",
473 "You are trying a global restart from a local h5 file!");
481 if(referenceTheta <= -180.0 || referenceTheta > 180.0) {
482 throw OpalException(
"Error in ParallelCyclotronTracker::visitCyclotron",
483 "PHIINIT is out of [-180, 180)!");
491 *gmsg <<
"* Bunch global starting position:" <<
endl;
496 *gmsg <<
"* Bunch global starting momenta:" <<
endl;
499 *gmsg <<
"* Reference total momentum (beta * gamma) = " <<
referencePtot * 1000.0 <<
" [MCU]" <<
endl;
500 *gmsg <<
"* Reference azimuthal momentum (Pt) = " <<
referencePt * 1000.0 <<
" [MCU]" <<
endl;
501 *gmsg <<
"* Reference radial momentum (Pr) = " <<
referencePr * 1000.0 <<
" [MCU]" <<
endl;
502 *gmsg <<
"* Reference axial momentum (Pz) = " <<
referencePz * 1000.0 <<
" [MCU]" <<
endl;
506 *gmsg <<
"* " << sym <<
"-fold field symmetry " <<
endl;
513 *gmsg <<
"* Field map file name = " << fmfn <<
" " <<
endl;
516 *gmsg <<
"* Type of cyclotron = " << type <<
" " <<
endl;
520 *gmsg <<
"* Radial aperture = " << rmin <<
" ... " << rmax<<
" [m] "<<
endl;
524 *gmsg <<
"* Vertical aperture = " << zmin <<
" ... " << zmax<<
" [m]"<<
endl;
528 *gmsg <<
"* Harmonic number h = " << h <<
" " <<
endl;
546 double BcParameter[8] = {};
560 *gmsg <<
"In BeamBeam tracker is missing " <<
endl;
565 *gmsg <<
"* ------------------------------ Beam Stripping ------------------------------" <<
endl;
570 double BcParameter[8] = {};
574 *gmsg <<
"* Pressure = " << pressure <<
" [mbar]" <<
endl;
575 BcParameter[0] = pressure;
579 *gmsg <<
"* Temperature = " << temperature <<
" [K]" <<
endl;
582 *gmsg <<
"* Residual gas = " << gas <<
endl;
585 *gmsg << std::boolalpha <<
"* Particles stripped will be deleted after interaction -> " << stop <<
endl;
589 BcParameter[1] = temperature;
590 BcParameter[2] = stop;
603 *gmsg <<
"* ------------------------------ Collimator ------------------------------" <<
endl;
609 *gmsg <<
"* Xstart = " << xstart <<
" [mm]" <<
endl;
611 double xend = elptr->
getXEnd();
612 *gmsg <<
"* Xend = " << xend <<
" [mm]" <<
endl;
615 *gmsg <<
"* Ystart = " << ystart <<
" [mm]" <<
endl;
617 double yend = elptr->
getYEnd();
618 *gmsg <<
"* Yend = " << yend <<
" [mm]" <<
endl;
621 *gmsg <<
"* Zstart = " << zstart <<
" [mm]" <<
endl;
623 double zend = elptr->
getZEnd();
624 *gmsg <<
"* Zend = " << zend <<
" [mm]" <<
endl;
627 *gmsg <<
"* Width = " << width <<
" [mm]" <<
endl;
631 double BcParameter[8] = {};
633 BcParameter[0] = 0.001 * xstart ;
634 BcParameter[1] = 0.001 * xend;
635 BcParameter[2] = 0.001 * ystart ;
636 BcParameter[3] = 0.001 * yend;
637 BcParameter[4] = 0.001 * width ;
706 "ParallelCylcotronTracker::visitOffset",
707 "Attempt to place an offset when Ring not defined");
742 *gmsg <<
"In Multipole; L= " << mult.
getElementLength() <<
" however the element is missing " <<
endl;
752 *gmsg <<
"Adding MultipoleT" <<
endl;
756 throw OpalException(
"ParallelCyclotronTracker::visitMultipoleT",
757 "Need to define a RINGDEFINITION to use MultipoleT element");
768 *gmsg <<
"Adding MultipoleTStraight" <<
endl;
772 throw OpalException(
"ParallelCyclotronTracker::visitMultipoleTStraight",
773 "Need to define a RINGDEFINITION to use MultipoleTStraight element");
775 myElements.push_back(dynamic_cast<MultipoleTStraight *>(multTstraight.
clone()));
784 *gmsg <<
"Adding MultipoleTCurvedConstRadius" <<
endl;
788 throw OpalException(
"ParallelCyclotronTracker::visitMultipoleTCurvedConstRadius",
789 "Need to define a RINGDEFINITION to use MultipoleTCurvedConstRadius element");
791 myElements.push_back(dynamic_cast<MultipoleTCurvedConstRadius *>(multTccurv.
clone()));
800 *gmsg <<
"Adding MultipoleTCurvedVarRadius" <<
endl;
804 throw OpalException(
"ParallelCyclotronTracker::visitMultipoleTCurvedVarRadius",
805 "Need to define a RINGDEFINITION to use MultipoleTCurvedVarRadius element");
807 myElements.push_back(dynamic_cast<MultipoleTCurvedVarRadius *>(multTvcurv.
clone()));
816 *gmsg <<
"* ------------------------------ Probe ------------------------------" <<
endl;
823 *gmsg <<
"* XStart = " << xstart <<
" [mm]" <<
endl;
825 double xend = elptr->
getXEnd();
826 *gmsg <<
"* XEnd = " << xend <<
" [mm]" <<
endl;
829 *gmsg <<
"* YStart = " << ystart <<
" [mm]" <<
endl;
831 double yend = elptr->
getYEnd();
832 *gmsg <<
"* YEnd = " << yend <<
" [mm]" <<
endl;
837 double BcParameter[8] = {};
839 BcParameter[0] = 0.001 * xstart ;
840 BcParameter[1] = 0.001 * xend;
841 BcParameter[2] = 0.001 * ystart ;
842 BcParameter[3] = 0.001 * yend;
843 BcParameter[4] = 0.001 ;
855 *gmsg <<
"In RBend; L= " << bend.
getElementLength() <<
" however the element is missing " <<
endl;
860 *gmsg <<
"Adding SBend3D" <<
endl;
864 throw OpalException(
"ParallelCyclotronTracker::visitSBend3D",
865 "Need to define a RINGDEFINITION to use SBend3D element");
869 *gmsg <<
"Adding ScalingFFAMagnet" <<
endl;
873 throw OpalException(
"ParallelCyclotronTracker::visitScalingFFAMagnet",
874 "Need to define a RINGDEFINITION to use ScalingFFAMagnet element");
879 *gmsg <<
"Adding Variable RF Cavity" <<
endl;
883 throw OpalException(
"ParallelCyclotronTracker::visitVariableRFCavity",
884 "Need to define a RINGDEFINITION to use VariableRFCavity element");
889 *gmsg <<
"Adding Variable RF Cavity with Fringe Field" <<
endl;
890 if (opalRing_m != NULL)
891 opalRing_m->appendElement(cav);
893 throw OpalException(
"ParallelCyclotronTracker::visitVariableRFCavityFringeField",
894 "Need to define a RINGDEFINITION to use VariableRFCavity element");
898 *gmsg <<
"Adding Vertical FFA Magnet" <<
endl;
902 throw OpalException(
"ParallelCyclotronTracker::visitVerticalFFAMagnet",
903 "Need to define a RINGDEFINITION to use VerticalFFAMagnet element");
913 *gmsg <<
"* ------------------------------ RFCavity ------------------------------" <<
endl;
920 throw OpalException(
"ParallelCyclotronTracker::visitRFCavity",
921 "The ParallelCyclotronTracker can only play with cyclotron type RF system currently ...");
924 double rmin = elptr->
getRmin();
925 *gmsg <<
"* Minimal radius of cavity= " << rmin <<
" [mm]" <<
endl;
927 double rmax = elptr->
getRmax();
928 *gmsg <<
"* Maximal radius of cavity= " << rmax <<
" [mm]" <<
endl;
931 *gmsg <<
"* RF frequency (2*pi*f)= " << rff <<
" [rad/s]" <<
endl;
934 *gmsg <<
"* RF Field map file name= " << fmfn <<
endl;
937 *gmsg <<
"* Cavity azimuth position= " << angle <<
" [deg] " <<
endl;
940 *gmsg <<
"* Cavity gap width= " << gap <<
" [mm] " <<
endl;
943 *gmsg <<
"* Cavity Shift distance= " << pdis <<
" [mm] " <<
endl;
945 double phi0 = elptr->
getPhi0();
946 *gmsg <<
"* Initial RF phase (t=0)= " << phi0 <<
" [deg] " <<
endl;
953 std::shared_ptr<AbstractTimeDependence> freq_atd =
nullptr;
954 std::shared_ptr<AbstractTimeDependence> ampl_atd =
nullptr;
955 std::shared_ptr<AbstractTimeDependence> phase_atd =
nullptr;
958 unityVec.push_back(1.);
959 unityVec.push_back(0.);
960 unityVec.push_back(0.);
961 unityVec.push_back(0.);
987 double BcParameter[8] = {};
989 BcParameter[0] = 0.001 * rmin;
990 BcParameter[1] = 0.001 * rmax;
991 BcParameter[2] = 0.001 * pdis;
992 BcParameter[3] = angle;
1003 *gmsg <<
"In RFQuadrupole; L = " << rfq.
getElementLength() <<
" however the element is missing " <<
endl;
1013 *gmsg <<
"In SBend; L = " << bend.
getElementLength() <<
" however the element is missing " <<
endl;
1023 *gmsg <<
"In Separator L= " << sep.
getElementLength() <<
" however the element is missing " <<
endl;
1034 *gmsg <<
endl <<
"* ------------------------------ Septum ------------------------------- *" <<
endl;
1040 *gmsg <<
"* XStart = " << xstart <<
" [mm]" <<
endl;
1042 double xend = elptr->
getXEnd();
1043 *gmsg <<
"* XEnd = " << xend <<
" [mm]" <<
endl;
1046 *gmsg <<
"* YStart = " << ystart <<
" [mm]" <<
endl;
1048 double yend = elptr->
getYEnd();
1049 *gmsg <<
"* YEnd = " << yend <<
" [mm]" <<
endl;
1052 *gmsg <<
"* Width = " << width <<
" [mm]" <<
endl;
1058 double BcParameter[8] = {};
1060 BcParameter[0] = 0.001 * xstart ;
1061 BcParameter[1] = 0.001 * xend;
1062 BcParameter[2] = 0.001 * ystart ;
1063 BcParameter[3] = 0.001 * yend;
1064 BcParameter[4] = 0.001 * width ;
1079 *gmsg <<
"Solenoid: no position of the element given!" <<
endl;
1114 *gmsg <<
"* ------------------------------ Stripper ------------------------------" <<
endl;
1122 *gmsg <<
"* XStart = " << xstart <<
" [mm]" <<
endl;
1124 double xend = elptr->
getXEnd();
1125 *gmsg <<
"* XEnd = " << xend <<
" [mm]" <<
endl;
1128 *gmsg <<
"* YStart = " << ystart <<
" [mm]" <<
endl;
1130 double yend = elptr->
getYEnd();
1131 *gmsg <<
"* YEnd = " << yend <<
" [mm]" <<
endl;
1134 *gmsg <<
"* Charge of outcoming particle = +e * " << opcharge <<
endl;
1137 *gmsg <<
"* Mass of the outcoming particle = " << opmass <<
" [GeV/c^2]" <<
endl;
1140 *gmsg << std::boolalpha <<
"* Particles stripped will be deleted after interaction -> " << stop <<
endl;
1144 double BcParameter[8] = {};
1146 BcParameter[0] = 0.001 * xstart ;
1147 BcParameter[1] = 0.001 * xend;
1148 BcParameter[2] = 0.001 * ystart ;
1149 BcParameter[3] = 0.001 * yend;
1150 BcParameter[4] = 0.001;
1151 BcParameter[5] = opcharge;
1152 BcParameter[6] = opmass;
1168 localpair->first = elementType;
1170 for(
int i = 0; i < 8; i++)
1171 *(((localpair->second).first) + i) = *(BcParameter + i);
1173 (localpair->second).second = elptr;
1198 int maxnlp = 111111;
1201 *gmsg << s <<
" min local particle number " << minnlp <<
" max local particle number: " << maxnlp <<
endl;
1233 *gmsg <<
"* ---------------------------------------------------" <<
endl;
1234 *gmsg <<
"* The selected Beam line elements are :" <<
endl;
1252 *gmsg <<
"* ---------------------------------------------------" <<
endl;
1260 for(
int k = 0; k < 2; k++)
1270 std::placeholders::_1,
1271 std::placeholders::_2,
1272 std::placeholders::_3,
1273 std::placeholders::_4);
1278 *gmsg <<
"* 4th order Runge-Kutta integrator" <<
endl;
1282 *gmsg <<
"* 2nd order Leap-Frog integrator" <<
endl;
1286 *gmsg <<
"* Multiple time stepping (MTS) integrator" <<
endl;
1292 "Invalid name of TIMEINTEGRATOR in Track command");
1301 *gmsg <<
"* ----------------------------------------------- *" <<
endl;
1302 *gmsg <<
"* Finalizing i.e. write data and close files :" <<
endl;
1303 for(
auto fd : FieldDimensions) {
1304 ((fd->second).second)->finalise();
1306 *gmsg <<
"* ----------------------------------------------- *" <<
endl;
1321 double t = 0, dt = 0, oldReferenceTheta = 0;
1326 *gmsg <<
"MTS: Number of substeps per step is " << numSubsteps <<
endl;
1328 double const dt_inner = dt / double(numSubsteps);
1330 *gmsg <<
"MTS: The inner time step is therefore " << dt_inner <<
" [ns]" <<
endl;
1334 bool flagTransition =
false;
1337 *gmsg <<
"* ---------------------------- Start tracking ----------------------------" <<
endl;
1344 bool finishedTurn =
false;
1358 for(
int n = 0;
n < numSubsteps; ++
n)
1401 temp_meanTheta, finishedTurn);
1404 oldReferenceTheta, temp_meanTheta);
1406 oldReferenceTheta = temp_meanTheta;
1413 finishedTurn =
true;
1417 <<
", Total number of live particles: "
1435 *gmsg <<
"* ---------------------------- DONE TRACKING PARTICLES -------------------------------- * " <<
endl;
1470 double t = 0, dt = 0, oldReferenceTheta = 0;
1489 *gmsg <<
"* --------------------------------- Start tracking ------------------------------------ *" <<
endl;
1493 bool finishedTurn =
false;
1499 seoMode_m(t, dt, finishedTurn, Ttime, Tdeltr, Tdeltz, TturnNumber);
1515 throw OpalException(
"ParallelCyclotronTracker::GenericTracker()",
1516 "No such tracking mode.");
1525 *gmsg <<
"* ---------------------------- DONE TRACKING PARTICLES -------------------------------- * " <<
endl;
1569 double distNew = (Rnew[0] * sinx - Rnew[1] * cosx) - PerpenDistance;
1570 double distOld = (Rold[0] * sinx - Rold[1] * cosx) - PerpenDistance;
1571 if(distOld > 0.0 && distNew <= 0.0) flag =
true;
1573 Dold = 1.0e3 * distOld;
1581 double rmin = rfcavity->
getRmin();
1582 double rmax = rfcavity->
getRmax();
1583 double nomalRadius = (radius - rmin) / (rmax - rmin);
1585 if(nomalRadius <= 1.0 && nomalRadius >= 0.0) {
1587 for(
int j = 0; j < 3; j++)
1593 for(
int j = 0; j < 3; j++)
1601 struct adder :
public std::unary_function<double, void> {
1620 int lastTurn,
double &nur,
double &nuz) {
1623 int Ndat = t.size();
1630 for(
int i = 0; i < Ndat; i++)
1635 for(
int i = 0; i < Ndat; i++)
1637 double ti = *(t.begin());
1638 double tf = t[t.size()-1];
1639 double T = (tf - ti);
1642 for(
int i = 0; i < Ndat; i++) {
1649 *gmsg <<
"* ************************************* nuR ******************************************* *" <<
endl;
1650 *gmsg << endl <<
"* ===> " << Ndat <<
" data points Ti=" << ti <<
" Tf= " << tf <<
" -> T= " << T <<
endl;
1656 stat = tune->
lombAnalysis(t, r, nhis_lomb, T / lastTurn);
1658 *gmsg <<
"* TUNE: Lomb analysis failed" <<
endl;
1659 *gmsg <<
"* ************************************************************************************* *" <<
endl;
1669 *gmsg <<
"* ************************************* nuZ ******************************************* *" <<
endl;
1670 *gmsg << endl <<
"* ===> " << Ndat <<
" data points Ti=" << ti <<
" Tf= " << tf <<
" -> T= " << T <<
endl;
1674 stat = tune->
lombAnalysis(t, z, nhis_lomb, T / lastTurn);
1676 *gmsg <<
"* TUNE: Lomb analysis failed" <<
endl;
1677 *gmsg <<
"* ************************************************************************************* *" <<
endl;
1692 throw OpalException(
"ParallelCyclotronTracker::getHarmonicNumber()",
1693 std::string(
"The first item in the FieldDimensions list does not ")
1694 +std::string(
"seem to be an Ring or a Cyclotron element"));
1706 for(
int d = 0; d < 3; ++d) {
1725 for(
int d = 0; d < 3; ++d) {
1744 double phi,
Vector_t const translationToGlobal) {
1746 particleVectors -= translationToGlobal;
1754 particleVectors[i] =
dot(rotation, particleVectors[i]);
1760 double phi,
Vector_t const translationToGlobal) {
1768 particleVectors[i] =
dot(rotation, particleVectors[i]);
1771 particleVectors += translationToGlobal;
1782 particleVectors -= meanR;
1795 reverseQuaternion(0) *= -1.0;
1801 particleVectors += meanR;
1814 particleVectors -= meanR;
1864 particleVectors += meanR;
1891 Vector_t const quaternionVectorComponent =
Vector_t(quaternion(1), quaternion(2), quaternion(3));
1892 double const quaternionScalarComponent = quaternion(0);
1896 particleVectors[i] = 2.0f *
dot(quaternionVectorComponent, particleVectors[i]) * quaternionVectorComponent +
1897 (quaternionScalarComponent * quaternionScalarComponent -
1898 dot(quaternionVectorComponent, quaternionVectorComponent)) * particleVectors[i] + 2.0f *
1899 quaternionScalarComponent *
cross(quaternionVectorComponent, particleVectors[i]);
1905 double tolerance = 1.0e-10;
1906 double length2 =
dot(quaternion, quaternion);
1908 if (
fabs(length2) > tolerance &&
fabs(length2 - 1.0f) > tolerance) {
1910 double length =
sqrt(length2);
1911 quaternion /= length;
1917 double tolerance = 1.0e-10;
1918 double length2 =
dot(vector, vector);
1920 if (
fabs(length2) > tolerance &&
fabs(length2 - 1.0f) > tolerance) {
1922 double length =
sqrt(length2);
1936 particleVectors[i] =
dot(rotation, particleVectors[i]);
1947 myVector =
dot(rotation, myVector);
1959 particleVectors[i] =
dot(rotation, particleVectors[i]);
1970 myVector =
dot(rotation, myVector);
1979 double k_cos_theta =
dot(u, v);
1981 double tolerance1 = 1.0e-5;
1982 double tolerance2 = 1.0e-8;
1985 if (
fabs(k_cos_theta / k + 1.0) < tolerance1) {
1992 if (
dot(resultVectorComponent, resultVectorComponent) < tolerance2) {
1997 double halfAngle = 0.5 *
pi;
1998 double sinHalfAngle =
sin(halfAngle);
2000 resultVectorComponent *= sinHalfAngle;
2003 k_cos_theta =
cos(halfAngle);
2007 resultVectorComponent =
cross(u, v);
2011 quaternion(0) = k_cos_theta + k;
2012 quaternion(1) = resultVectorComponent(0);
2013 quaternion(2) = resultVectorComponent(1);
2014 quaternion(3) = resultVectorComponent(2);
2030 bool flagNeedUpdate =
false;
2038 double const distNew = (
itsBunch_m->
R[i][0] * ccd.sinAzimuth -
itsBunch_m->
R[i][1] * ccd.cosAzimuth) - ccd.perpenDistance;
2039 bool tagCrossing =
false;
2041 if(distNew <= 0.0) {
2042 distOld = (oldR[0] * ccd.sinAzimuth - oldR[1] * ccd.cosAzimuth) - ccd.perpenDistance;
2043 if(distOld > 0.0) tagCrossing =
true;
2046 double const dt1 = distOld /
sqrt(
dot(v, v));
2047 double const dt2 = h - dt1;
2064 return flagNeedUpdate;
2071 bool flagNeedUpdate =
false;
2084 return flagNeedUpdate;
2092 bool flagNeedUpdate =
push(0.5 * h);
2107 flagNeedUpdate |=
kick(h);
2110 flagNeedUpdate |=
push(0.5 * h);
2142 *gmsg <<
"* Total number of particles after PluginElement= "
2161 allreduce(flagNeedUpdate, 1, std::logical_or<bool>());
2163 if(flagNeedUpdate) {
2165 std::vector<size_t> locLostParticleNum(bunchCount, 0);
2168 std::unique_ptr<size_t[]> localBinCount;
2171 localBinCount = std::unique_ptr<size_t[]>(
new size_t[leb]);
2172 for (
int i = 0; i < leb; ++i)
2173 localBinCount[i] = 0;
2190 for (
int i = 0; i < leb; ++i) {
2195 std::vector<size_t> localnum(bunchCount + 1);
2196 for (
size_t i = 0; i < localnum.size() - 1; ++i) {
2212 std::vector<size_t> totalnum(bunchCount + 1);
2215 allreduce(localnum.data(), totalnum.data(), localnum.size(), std::plus<size_t>());
2218 for (
short i = 0; i < bunchCount; ++i) {
2222 size_t sum = std::accumulate(totalnum.begin(),
2223 totalnum.end() - 1, 0);
2225 if ( sum != totalnum[bunchCount] ) {
2226 throw OpalException(
"ParallelCyclotronTracker::deleteParticle()",
2227 "Total number of particles " + std::to_string(totalnum[bunchCount]) +
2228 " != " + std::to_string(sum) +
" (sum over all bunches)");
2231 size_t globLostParticleNum = 0;
2232 size_t locNumLost = std::accumulate(locLostParticleNum.begin(),
2233 locLostParticleNum.end(), 0);
2235 reduce(locNumLost, globLostParticleNum, 1, std::plus<size_t>());
2237 if ( globLostParticleNum > 0 ) {
2238 *gmsg <<
"At step " <<
step_m
2239 <<
", lost " << globLostParticleNum <<
" particles "
2240 <<
"on stripper, collimator, septum, geometry, "
2241 <<
"out of cyclotron aperture or beam stripping"
2245 if (totalnum[bunchCount] == 0) {
2247 return flagNeedUpdate;
2257 double const psi = 0.5 *
pi -
acos(meanP(2) /
sqrt(
dot(meanP, meanP)));
2284 return flagNeedUpdate;
2367 *gmsg <<
"* Restart in the local frame" <<
endl;
2389 *gmsg <<
"* Restart in the global frame" <<
endl;
2407 double const psi = 0.5 *
pi -
acos(meanP(2) /
sqrt(
dot(meanP, meanP)));
2409 double radius =
std::sqrt(meanR[0] * meanR[0] + meanR[1] * meanR[1]);
2428 *gmsg <<
endl <<
"* *********************** Bunch information in local frame: ************************";
2456 *gmsg << endl <<
"* *********************** Bunch information in global frame: ***********************";
2476 int found[2] = {-1, -1};
2492 int numberOfPart = 0;
2494 while(notReceived > 0) {
2500 ERRORMSG(
"Could not receive from client nodes in main." <<
endl);
2504 rmsg->
get(&numberOfPart);
2506 for(
int i = 0; i < numberOfPart; ++i) {
2509 for(
int ii = 0; ii < 6; ii++) {
2517 for(
int i = 0; i < counter; ++i) {
2521 for(
int j = 0; j < 3; ++j) {
2530 for(
auto tmpid : tmpi) {
2542 for(
int ii = 0; ii < 6; ii++) {
2555 for(
int i = 0; i < counter; i++) {
2559 for(
int j = 0; j < 3; j++) {
2567 ERRORMSG(
"Ippl::Comm->send(smsg, 0, tag) failed " <<
endl);
2599 double phi = 0.0, psi = 0.0;
2674 (((*DumpSindex)->second).second)->apply(meanR, meanP, temp_t,
extE_m,
extB_m);
2739 double const betagamma_temp =
sqrt(
dot(meanP, meanP));
2744 double const theta =
atan2(meanR(1), meanR(0));
2750 double const psi = 0.5 *
pi -
acos(meanP(2) /
sqrt(
dot(meanP, meanP)));
2772 (((*DumpSindex)->second).second)->apply(meanR, meanP, temp_t,
extE_m,
extB_m);
2778 bool dumpLocalFrame =
true;
2779 std::string dumpString =
"local";
2781 dumpLocalFrame =
false;
2782 dumpString =
"global";
2785 if (dumpLocalFrame ==
true) {
2814 if (dumpLocalFrame ==
true) {
2825 *gmsg <<
endl <<
"* Integration step " <<
step_m + 1
2826 <<
" (no phase space dump for <= 2 particles)" <<
endl;
2829 <<
" (" << dumpString <<
" frame) at integration step " <<
step_m + 1 <<
endl;
2834 *gmsg <<
"* T = " << temp_t <<
" ns"
2836 *gmsg <<
"* E = " << E <<
" MeV"
2837 <<
", beta * gamma = " << betagamma_temp <<
endl;
2838 *gmsg <<
"* Bunch position: R = " <<
referenceR <<
" mm"
2853 const bool& finishedTurn)
2861 if (!(
step_m + 1 % 1000))
2942 *gmsg <<
"* Beginning of this run is at t = " << t <<
" [ns]" <<
endl;
2943 *gmsg <<
"* The time step is set to dt = " << dt <<
" [ns]" <<
endl;
2947 *gmsg <<
"* MBM: Time interval between neighbour bunches is set to "
2949 *gmsg <<
"* MBM: The particles energy bin reset frequency is set to "
2962 *gmsg <<
"* ------------------------- STATIC EQUILIBRIUM ORBIT MODE ----------------------------- *" << endl
2963 <<
"* Instruction: When the total particle number is equal to 2, SEO mode is triggered *" << endl
2964 <<
"* automatically. This mode does NOT include any RF cavities. The initial distribution *" << endl
2965 <<
"* file must be specified. In the file the first line is for reference particle and the *" << endl
2966 <<
"* second line is for off-center particle. The tune is calculated by FFT routines based *" << endl
2967 <<
"* on these two particles. *" << endl
2968 <<
"* ---------------- NOTE: SEO MODE ONLY WORKS SERIALLY ON SINGLE NODE ------------------ *" <<
endl;
2971 throw OpalException(
"Error in ParallelCyclotronTracker::initializeTracking_m",
2972 "SEO MODE ONLY WORKS SERIALLY ON SINGLE NODE!");
2976 *gmsg <<
"* ------------------------------ SINGLE PARTICLE MODE --------------------------------- *" << endl
2977 <<
"* Instruction: When the total particle number is equal to 1, single particle mode is *" << endl
2978 <<
"* triggered automatically. The initial distribution file must be specified which should *" << endl
2979 <<
"* contain only one line for the single particle *" << endl
2980 <<
"* ---------NOTE: SINGLE PARTICLE MODE ONLY WORKS SERIALLY ON A SINGLE NODE ------------ *" <<
endl;
2983 throw OpalException(
"Error in ParallelCyclotronTracker::initializeTracking_m",
2984 "SINGLE PARTICLE MODE ONLY WORKS SERIALLY ON A SINGLE NODE!");
2993 throw OpalException(
"ParallelCyclotronTracker::initializeTracking_m()",
2994 "No such tracking mode.");
2997 return std::make_tuple(t, dt, oldReferenceTheta);
3008 double FinalEnergy = (
sqrt(1.0 + FinalMomentum2) - 1.0) *
itsBunch_m->
getM() * 1.0e-6;
3009 *gmsg <<
"* Final energy of reference particle = " << FinalEnergy <<
" [MeV]" <<
endl;
3010 *gmsg <<
"* Total phase space dump number(includes the initial distribution) = " <<
lastDumpedStep_m + 1 <<
endl;
3011 *gmsg <<
"* One can restart simulation from the last dump step (--restart " <<
lastDumpedStep_m <<
")" <<
endl;
3023 *gmsg <<
"* **************** The result for tune calulation (NO space charge) ******************* *" << endl
3024 <<
"* Number of tracked turns: " << TturnNumber.back() <<
endl;
3026 getTunes(Ttime, Tdeltr, Tdeltz, TturnNumber.back(), nur, nuz);
3038 *gmsg <<
"*" <<
endl;
3040 *gmsg <<
"* Cave: Turn number is not correct for restart mode"<<
endl;
3050 *gmsg <<
endl <<
"* *********************** Bunch information in global frame: ***********************";
3059 *gmsg << endl <<
"* No Particles left in bunch!" <<
endl;
3060 *gmsg <<
"* **********************************************************************************" <<
endl;
3073 double r_tuning[2], z_tuning[2] ;
3109 Ttime.push_back(t * 1.0
e-9);
3110 Tdeltz.push_back(z_tuning[1]);
3111 Tdeltr.push_back(r_tuning[1] - r_tuning[0]);
3118 bool& finishedTurn,
double& oldReferenceTheta) {
3155 temp_meanTheta, finishedTurn);
3158 oldReferenceTheta, temp_meanTheta);
3160 oldReferenceTheta = temp_meanTheta;
3186 static bool flagTransition =
false;
3266 finishedTurn =
true;
3270 <<
", Total number of live particles: "
3286 bool tag_crossing =
false;
3287 double DistOld = 0.0;
3292 rfcav =
static_cast<RFCavity *
>(((*sindex)->second).second);
3296 if ( tag_crossing ) {
3299 double oldMomentum2 =
dot(Pold, Pold);
3300 double oldBetgam =
sqrt(oldMomentum2);
3301 double oldGamma =
sqrt(1.0 + oldMomentum2);
3302 double oldBeta = oldBetgam / oldGamma;
3303 double dt1 = DistOld / (
Physics::c * oldBeta * 1.0e-6);
3304 double dt2 = dt - dt1;
3311 if (dt / dt1 < 1.0e9)
3315 RFkick(rfcav, t, dt1, i);
3320 if (dt / dt2 < 1.0e9)
3330 const double& oldReferenceTheta,
3331 const double& temp_meanTheta)
3333 for (
unsigned int i=0; i<=2; i++) {
3337 <<
", Time = " << t <<
" [ns]" <<
std::endl
3338 <<
" " << std::hypot(
R(0),
R(1))
3339 <<
" " << P(0) *
std::cos(temp_meanTheta) + P(1) *
std::sin(temp_meanTheta)
3341 <<
" " << - P(0) *
std::sin(temp_meanTheta) + P(1) *
std::cos(temp_meanTheta)
3351 const double& temp_meanTheta,
3358 finishedTurn =
true;
3363 <<
", Time = " << t <<
" [ns]" <<
std::endl
3365 <<
" " << P(0) *
std::cos(temp_meanTheta) +
3368 <<
" " << - P(0) *
std::sin(temp_meanTheta) +
3472 bool outOfBound = (((*sindex)->second).second)->apply(i, t, Efield, Bfield);
3498 if ( flagTransition ) {
3499 *gmsg <<
"* MBM: Saving beam distribution at turn " <<
turnnumber_m <<
endl;
3500 *gmsg <<
"* MBM: After one revolution, Multi-Bunch Mode will be invoked" <<
endl;
3510 throw OpalException(
"ParallelCyclotronTracker::injectBunch()",
3511 "Unknown return value " + std::to_string(result));
3544 std::vector<double> lpaths(1);
3579 for (
short b = 0; b <
mbHandler_m->getNumBunch(); ++b) {
virtual std::string getFieldMapFN() const
virtual double getMaxR() const
virtual double getRmin() const
void updatePathLength(const double &dt)
ParticleAttrib< Vector_t > P
void singleParticleDump()
double getXStart() const
Member variable access.
static Vector_t const yaxis
bool computeExternalFields_m(const size_t &i, const double &t, Vector_t &Efield, Vector_t &Bfield)
double calculateAngle(double x, double y)
virtual void visitOffset(const Offset &)
Apply the algorithm to a Offset.
int lombAnalysis(double *x, double *y, int Ndat, int nhis)
virtual ElementBase * clone() const override
PETE_TUTree< FnArcCos, typename T::PETE_Expr_t > acos(const PETE_Expr< T > &l)
virtual std::string getPressureMapFN() const
constexpr double e
The value of .
void destroy(size_t M, size_t I, bool doNow=false)
void finalizeTracking_m(dvector_t &Ttime, dvector_t &Tdeltr, dvector_t &Tdeltz, ivector_t &TturnNumber)
ParticleAttrib< Vector_t > Ef
virtual ~ParallelCyclotronTracker()
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
void setOpenMode(OPENMODE openMode)
double getBeamRInit() const
void dumpAngle(const double &theta, double &prevAzimuth, double &azimuth, const short &bunchNr=0)
std::tuple< double, double, double > initializeTracking_m()
void operator()(double x)
virtual bool getSpiralFlag() const
std::string getFieldSolverType() const
Return the fieldsolver type if we have a fieldsolver.
std::vector< double > azimuth_angle_m
the different azimuthal angles for the outfTheta_m output files
Interface for septum magnet.
long long getLocalTrackStep() const
virtual double getZinit() const
double getZStart()
Member variable access.
void bgf_main_collision_test()
void writeMultiBunchStatistics(PartBunchBase< double, 3 > *beam, MultiBunchHandler *mbhandler)
Interface for electrostatic separator.
virtual void visitMarker(const Marker &)
Apply the algorithm to a Marker.
Interface for a Cyclotron.
void updateAzimuthAndRadius()
std::list< Component * > myElements
void setLocalNumPerBunch(size_t numpart, short n)
PETE_TUTree< FnFabs, typename T::PETE_Expr_t > fabs(const PETE_Expr< T > &l)
void setGlobalMeanR(Vector_t globalMeanR)
virtual void visitCCollimator(const CCollimator &)
Apply the algorithm to a collimator.
Interface for beam position monitors.
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
virtual double getPHIinit() const
Interface for RF Quadrupole.
double getPressure() const
The base class for all OPAL exceptions.
Tps< T > sin(const Tps< T > &x)
Sine.
virtual void computeSelfFields_cycl(double gamma)=0
double getBeamPRInit() const
virtual bool checkBeamStripping(PartBunchBase< double, 3 > *bunch, Cyclotron *cycl, const int turnnumber, const double t, const double tstep)
void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion)
ParticleAttrib< double > Q
std::vector< PluginElement * > pluginElements_m
virtual void visitDegrader(const Degrader &)
Apply the algorithm to a Degrader.
virtual void visitVerticalFFAMagnet(const VerticalFFAMagnet &bend)
Apply the algorithm to a VerticalFFAMagnet.
virtual void visitRFCavity(const RFCavity &)
Apply the algorithm to a RFCavity.
IpplTimings::TimerRef DumpTimer_m
static Vector_t const zaxis
virtual void do_binaryRepart()
Interface for general corrector.
virtual double getSymmetry() const
std::vector< int > ivector_t
Interface for beam diagnostics.
4-th order Runnge-Kutta stepper
double computeRadius(const Vector_t &meanR) const
double get_meanKineticEnergy() const
void setLocalTrackStep(long long n)
step in a TRACK command
struct ParallelCyclotronTracker::settings setup_m
virtual void visitRFQuadrupole(const RFQuadrupole &)
Apply the algorithm to a RFQuadrupole.
const std::string & getCyclotronType() const
IpplTimings::TimerRef IntegrationTimer_m
virtual void visitVariableRFCavity(const VariableRFCavity &cav)
Apply the algorithm to a VariabelRFCavity.
virtual const std::string & getName() const
Get element name.
ParticleAttrib< short > bunchNum
double getHarmonicNumber() const
int maxSteps_m
The maximal number of steps the system is integrated.
void localToGlobal(ParticleAttrib< Vector_t > &vectorArray, double phi, Vector_t const translationToGlobal=0)
virtual void visitSeptum(const Septum &)
Apply the algorithm to a Septum.
void globalToLocal(ParticleAttrib< Vector_t > &vectorArray, double phi, Vector_t const translationToGlobal=0)
Interface for drift space.
virtual void visitDrift(const Drift &)
Apply the algorithm to a Drift.
double getBeta() const
The relativistic beta per particle.
virtual ElementBase * clone() const override
virtual double getRmin() const
constexpr double rad2deg
The conversion factor from radians to degrees.
Inform & level4(Inform &inf)
void rotateAroundZ(ParticleAttrib< Vector_t > &particleVectors, double const phi)
void update_m(double &t, const double &dt, const bool &finishedTurn)
Update time and path length, write to output files.
void setTotalNumPerBunch(size_t numpart, short n)
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
size_t getTotalNum() const
void getQuaternionTwoVectors(Vector_t u, Vector_t v, Quaternion_t &quaternion)
virtual void visitSolenoid(const Solenoid &)
Apply the algorithm to a Solenoid.
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
Interface for general multipole.
void gapCrossKick_m(size_t i, double t, double dt, const Vector_t &Rold, const Vector_t &Pold)
int next_tag(int t, int s=1000)
Vector_t calcMeanR(short bunchNr=-1) const
virtual double getMinZ() const
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
virtual void visitMultipoleTCurvedVarRadius(const MultipoleTCurvedVarRadius &)
Apply the algorithm to a MultipoleTCurvedVarRadius.
int statDumpFreq
The frequency to dump statistical values, e.e. dump data when stepstatDumpFreq==0.
double bega
The reference variables.
T deg(T x)
Convert radians to degrees.
void bunchMode_m(double &t, const double dt, bool &finishedTurn)
ParallelCyclotronTracker()
virtual ElementBase * clone() const override
std::unique_ptr< Stepper< function_t > > itsStepper_mp
IpplTimings::TimerRef DelParticleTimer_m
virtual double getAzimuth() const
unsigned int delPartFreq
The frequency to delete particles (currently: OPAL-cycl only)
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
virtual double getCyclHarm() const
static OpalData * getInstance()
virtual void visitProbe(const Probe &)
Apply the algorithm to a Probe.
virtual double getMinR() const
static void writeFields(Component *field)
void calcBeamParameters()
Vector_t getNextPosition() const
virtual void execute()
Apply the algorithm to the top-level beamline.
void rotateAroundX(ParticleAttrib< Vector_t > &particleVectors, double const psi)
Template class for beam lines.
double getSymmetry() const
int sptDumpFreq
The frequency to dump single particle trajectory of particles with ID = 0 & 1.
ParticleAttrib< short > PType
virtual double getPhi0() const
virtual double getGapWidth() const
virtual void visitFlexibleCollimator(const FlexibleCollimator &)
Apply the algorithm to a flexible collimator.
void initTrackOrbitFile()
constexpr double pi
The value of .
virtual double getBScale() const
void initDistInGlobalFrame()
virtual void visitRBend(const RBend &)
Apply the algorithm to a RBend.
bool getFieldsAtPoint(const double &t, const size_t &Pindex, Vector_t &Efield, Vector_t &Bfield)
virtual double getPerpenDistance() const
void countTotalNumPerBunch()
virtual void accept(BeamlineVisitor &visitor) const =0
Apply visitor.
virtual double getElementLength() const
Get design length.
bool getTunes(dvector_t &t, dvector_t &r, dvector_t &z, int lastTurn, double &nur, double &nuz)
virtual std::string getComponentType() const override
void injectBunch(bool &flagTransition)
virtual void visitMultipoleTCurvedConstRadius(const MultipoleTCurvedConstRadius &)
Apply the algorithm to a MultipoleTCurvedConstRadius.
int repartFreq
The frequency to do particles repartition for better load balance between nodes.
virtual void visitMultipole(const Multipole &)
Apply the algorithm to a Multipole.
void bunchDumpPhaseSpaceData()
int partInside(const Vector_t &r, const Vector_t &v, const double dt, int Parttype, const double Qloss, Vector_t &intecoords, int &triId)
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField)
void allreduce(const T *input, T *output, int count, Op op)
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
Pure virtual implementation of Component.
constexpr double c
The velocity of light in m/s.
bool isTurnDone()
Check if turn done.
Vector3D cross(const Vector3D &lhs, const Vector3D &rhs)
Vector cross product.
Interface for cyclotron collimator.
size_t getTotalNumPerBunch(short n) const
virtual void visitCyclotronValley(const CyclotronValley &)
Apply the algorithm to a CyclotronValley.it is empty for cyclotrontracker .
static void startTimer(TimerRef t)
static void writeFields(Component *field)
virtual void visitRing(const Ring &ring)
Apply the algorithm to an Ring.
virtual void visitScalingFFAMagnet(const ScalingFFAMagnet &bend)
Apply the algorithm to a ScalingFFAMagnet.
double getChargePerParticle() const
get the macro particle charge
void kick(const Vector_t &R, Vector_t &P, const Vector_t &Ef, const Vector_t &Bf, const double &dt) const
void computePathLengthUpdate(std::vector< double > &dl, const double &dt)
Abstract beam-beam interaction.
virtual bool getStop() const
const size_t initialLocalNum_m
std::string getFrequencyModelName()
virtual void iterate(BeamlineVisitor &, bool r2l) const
Apply visitor to all elements of the line.
virtual double getPZinit() const
void computeSpaceChargeFields_m()
IpplTimings::TimerRef TransformTimer_m
virtual std::string getFieldMapFN() const
std::string getTypeString() const
void dumpAzimuthAngles_m(const double &t, const Vector_t &R, const Vector_t &P, const double &oldReferenceTheta, const double &temp_meanTheta)
double getTemperature() const
bool RFkick(RFCavity *rfcavity, const double t, const double dt, const int Pindex)
virtual ElementBase * clone() const override
Interface for cyclotron valley.
PETE_TBTree< FnArcTan2, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > atan2(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
void dumpH5(PartBunchBase< double, 3 > *beam, Vector_t FDext[]) const
unsigned int getNumberOfTrimcoils() const
IpplTimings::TimerRef BinRepartTimer_m
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
virtual double getPRinit() const
Vektor< double, 3 > Vector_t
std::unique_ptr< MultiBunchHandler > mbHandler_m
bool hasMultiBunch() const
std::vector< std::unique_ptr< std::ofstream > > outfTheta_m
output coordinates at different azimuthal angles and one after every turn
int rebinFreq
The frequency to reset energy bin ID for all particles.
virtual void visitCorrector(const Corrector &)
Apply the algorithm to a Corrector.
std::ofstream outfTrackOrbit_m
stepper::INTEGRATOR stepper_m
An abstract sequence of beam line components.
void setTotalNum(size_t n)
const PartData itsReference
The reference information.
std::pair< ElementBase::ElementType, element_pair > type_pair
void borisExternalFields(double h)
int getFieldFlag(const std::string &type) const
virtual void visitVariableRFCavityFringeField(const VariableRFCavityFringeField &cav)
Apply the algorithm to a VariabelRFCavity.
void openFiles(size_t numFiles, std::string fn)
@ open / close output coordinate files
size_t getLocalNum() const
virtual ElementBase * clone() const =0
Return clone.
Vector_t calcMeanP() const
beamline_list FieldDimensions
size_t getLocalNumPerBunch(short n) const
Message & get(const T &cval)
virtual double getPScale() const
Tps< T > sqrt(const Tps< T > &x)
Square root.
std::function< bool(const double &, const size_t &, Vector_t &, Vector_t &)> function_t
virtual void visitSBend3D(const SBend3D &)
Apply the algorithm to a SBend3D.
virtual void visitCyclotron(const Cyclotron &cycl)
Apply the algorithm to a Cyclotorn.
int psDumpFreq
The frequency to dump the phase space, i.e.dump data when steppsDumpFreq==0.
double getBeamPhiInit() const
void setLocalBinCount(size_t num, int bin)
Message & put(const T &val)
bool psDumpEachTurn
phase space dump flag for OPAL-cycl
virtual void visitMultipoleT(const MultipoleT &)
Apply the algorithm to a MultipoleT.
double getGamma() const
The relativistic gamma per particle.
short getNumBunch() const
DumpFrame psDumpFrame
flag to decide in which coordinate frame the phase space will be dumped for OPAL-cycl ...
virtual double getRmax() const
void performDestroy(bool updateLocalNum=false)
virtual double getCycFrequency() const
bool isMultiBunch() const
void singleMode_m(double &t, const double dt, bool &finishedTurn, double &oldReferenceTheta)
void getMomentaKick(const double normalRadius, double momentum[], const double t, const double dtCorrt, const int PID, const double restMass, const int chargenumber)
used in OPAL-cycl
ParticleAttrib< short > cavityGapCrossed
particle just crossed cavity gap (for ParallelCyclotronTracker)
virtual double getRmax() const
IpplTimings::TimerRef PluginElemTimer_m
static Vector_t const xaxis
The positive axes unit vectors.
virtual void visitBeamStripping(const BeamStripping &)
Apply the algorithm to a Beam Stripping.
virtual double getRinit() const
virtual double getSinAzimuth() const
Tps< T > cos(const Tps< T > &x)
Cosine.
bool deleteParticle(bool=false)
virtual void visitBeamline(const Beamline &)
Apply the algorithm to a beam line.
constexpr double q_e
The elementary charge in As.
double getQ() const
Access to reference data.
std::vector< CavityCrossData > cavCrossDatas_m
void normalizeVector(Vector_t &vector)
void appendElement(const Component &element)
constexpr double deg2rad
The conversion factor from degrees to radians.
virtual void visitBeamBeam(const BeamBeam &)
Apply the algorithm to a BeamBeam.
virtual std::string getResidualGas() const
void setMultiBunchInitialPathLengh(MultiBunchHandler *mbhandler_p)
Vector_t getNextNormal() const
bool applyPluginElements(const double dt)
int scSolveFreq
The frequency to solve space charge fields.
ParticleAttrib< int > Bin
void dumpSDDS(PartBunchBase< double, 3 > *beam, Vector_t FDext[], const double &azimuth=-1) const
virtual void visitMonitor(const Monitor &)
Apply the algorithm to a Monitor.
Ring describes a ring type geometry for tracking.
virtual void visitSBend(const SBend &)
Apply the algorithm to a SBend.
virtual void setBinCharge(int bin, double q)
Set the charge of one bin to the value of q and all other to zero.
ParticleAttrib< Vector_t > Bf
std::string getAmplitudeModelName()
double getHarmonicNumber()
ParticleAttrib< double > M
void updateGeometry(Vector_t startPosition, Vector_t startDirection)
virtual void visitSeparator(const Separator &)
Apply the algorithm to a Separator.
virtual bool hasAttribute(const std::string &aKey) const
Test for existence of an attribute.
std::string::iterator iterator
void rotateWithQuaternion(ParticleAttrib< Vector_t > &vectorArray, Quaternion_t const quaternion)
virtual void visitLambertson(const Lambertson &)
Apply the algorithm to a Lambertson.
virtual double getCosAzimuth() const
virtual double getMaxZ() const
ElementBase * clone() const override
static TimerRef getTimer(const char *nm)
bool for_each(const BareFieldIterator< T, D > &p, SameFieldID s, C)
Interface for a single beam element.
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField) override
Pure virtual implementation of Component.
void seoMode_m(double &t, const double dt, bool &finishedTurn, dvector_t &Ttime, dvector_t &Tdeltr, dvector_t &Tdeltz, ivector_t &TturnNumber)
static void stopTimer(TimerRef t)
virtual void initialise(PartBunchBase< double, 3 > *bunch, double &startField, double &endField)
void dumpThetaEachTurn_m(const double &t, const Vector_t &R, const Vector_t &P, const double &temp_meanTheta, bool &finishedTurn)
std::string getInputBasename()
get input file name without extension
bool checkGapCross(Vector_t Rold, Vector_t Rnew, RFCavity *rfcavity, double &DistOld)
Message * receive_block(int &node, int &tag)
const size_t initialTotalNum_m
void normalizeQuaternion(Quaternion_t &quaternion)
int getStepsPerTurn() const
void updateTime(const double &dt)
static std::shared_ptr< AbstractTimeDependence > getTimeDependence(std::string name)
virtual void visitDiagnostic(const Diagnostic &)
Apply the algorithm to a Diagnostic.
static Communicate * Comm
void checkNumPart(std::string s)
virtual void visitMultipoleTStraight(const MultipoleTStraight &)
Apply the algorithm to a MultipoleTStraight.
void buildupFieldList(double BcParameter[], ElementBase::ElementType elementType, Component *elptr)
double calculateAngle2(double x, double y)
PartBunchBase< double, 3 > * itsBunch_m
The bunch of particles to be tracked.
double getOPCharge() const
std::vector< double > dvector_t
Inform & endl(Inform &inf)
BoundaryGeometry * getGlobalGeometry()
virtual void visitStripper(const Stripper &)
Apply the algorithm to a charge stripper.
std::string getPhaseModelName()
Interface for a Lambertson septum.
virtual void visitParallelPlate(const ParallelPlate &)
Apply the algorithm to a ParallelPlate, it is empty for cyclotrontracker .
std::unique_ptr< LossDataSink > lossDs_m
Track particles or bunches.