18#ifndef PART_BUNCH_BASE_HPP
19#define PART_BUNCH_BASE_HPP
40template <
class T,
unsigned Dim>
50 stateOfLastBoundP_(unitless),
54 globalMeanR_m(
Vector_t(0.0, 0.0, 0.0)),
55 globalToLocalQuaternion_m(
Quaternion_t(1.0, 0.0, 0.0, 0.0)),
62 couplingConstant_m(0.0),
64 massPerParticle_m(0.0),
69 binemitted_m(nullptr),
77 globalPartPerNode_m(nullptr),
90template <
class T,
unsigned Dim>
92 if (dist_m !=
nullptr) {
93 size_t isBeamEmitted = dist_m->getIfDistEmitting();
95 if (isBeamEmitted > 0)
104template <
class T,
unsigned Dim>
109 int lastEmittedBin = dist_m->getLastEmittedEnergyBin();
111 return lastEmittedBin;
115template <
class T,
unsigned Dim>
117 return dist_m->getNumberOfEmissionSteps();
121template <
class T,
unsigned Dim>
123 return dist_m->getNumberOfEnergyBins();
127template <
class T,
unsigned Dim>
130 size_t isBeamEmitting = dist_m->getIfDistEmitting();
132 if (isBeamEmitting > 0) {
133 *
gmsg <<
"*****************************************************" <<
endl
134 <<
"Warning: attempted to rebin, but not all distribution" <<
endl
135 <<
"particles have been emitted. Rebin failed." <<
endl
136 <<
"*****************************************************" <<
endl;
144template <
class T,
unsigned Dim>
146 bingamma_m = std::unique_ptr<double[]>(
new double[numberOfEnergyBins]);
147 binemitted_m = std::unique_ptr<size_t[]>(
new size_t[numberOfEnergyBins]);
148 for (
int i = 0; i < numberOfEnergyBins; i++)
153template <
class T,
unsigned Dim>
156 if (dist_m !=
nullptr)
157 return dist_m->getNumberOfEnergyBins() > 0;
163template <
class T,
unsigned Dim>
166 if (unit_state_ == unitless)
168 "Cannot make a unitless PartBunch unitless");
170 bool hasToReset =
false;
171 if (!
R.isDirty()) hasToReset =
true;
173 for (
size_t i = 0; i < getLocalNum(); i++) {
175 if (use_dt_per_particle)
181 unit_state_ = unitless;
183 if (hasToReset)
R.resetDirtyFlag();
187template <
class T,
unsigned Dim>
190 if (unit_state_ == units)
192 "Cannot apply units twice to PartBunch");
194 bool hasToReset =
false;
195 if (!
R.isDirty()) hasToReset =
true;
197 for (
size_t i = 0; i < getLocalNum(); i++) {
199 if (use_dt_per_particle)
207 if (hasToReset)
R.resetDirtyFlag();
211template <
class T,
unsigned Dim>
217template <
class T,
unsigned Dim>
219 std::vector<Distribution*> addedDistributions,
221 Inform m(
"setDistribution " );
231template <
class T,
unsigned Dim>
233 size_t numberOfParticles,
241template <
class T,
unsigned Dim>
247template <
class T,
unsigned Dim>
249 return (pbin_m !=
nullptr);
253template <
class T,
unsigned Dim>
259template <
class T,
unsigned Dim>
265template <
class T,
unsigned Dim>
267 if (dist_m !=
nullptr)
268 return dist_m->getIfDistEmitting();
274template <
class T,
unsigned Dim>
276 if (pbin_m !=
nullptr)
277 return pbin_m->weHaveBins();
283template <
class T,
unsigned Dim>
287 setEnergyBins(pbin_m->getNBins());
291template <
class T,
unsigned Dim>
294 setEnergyBins(pbin_m->getNBins());
298template <
class T,
unsigned Dim>
300 return dist_m->emitParticles(
this, eZ);
304template <
class T,
unsigned Dim>
306 size_t numParticles = getLocalNum();
308 setTotalNum(numParticles);
312template <
class T,
unsigned Dim>
321template <
class T,
unsigned Dim>
323 if (pbin_m !=
nullptr)
324 return pbin_m->getLastemittedBin();
330template <
class T,
unsigned Dim>
332 if (pbin_m !=
nullptr) {
333 pbin_m->setPartNum(bin, num);
338template <
class T,
unsigned Dim>
341 const int emittedBins = dist_m->getNumberOfEnergyBins();
345 for (
int i = 0; i < emittedBins; i++)
348 for (
unsigned int n = 0;
n < getLocalNum();
n++)
351 std::unique_ptr<size_t[]> particlesInBin(
new size_t[emittedBins]);
352 reduce(bingamma_m.get(), bingamma_m.get() + emittedBins, bingamma_m.get(),
OpAddAssign());
353 reduce(binemitted_m.get(), binemitted_m.get() + emittedBins, particlesInBin.get(),
OpAddAssign());
354 for (
int i = 0; i < emittedBins; i++) {
355 size_t &pInBin = particlesInBin[i];
357 bingamma_m[i] /= pInBin;
359 <<
" gamma = " << std::setw(8) << std::scientific
360 << std::setprecision(5) << bingamma_m[i]
361 <<
"; NpInBin= " << std::setw(8)
362 << std::setfill(
' ') << pInBin <<
endl);
369 particlesInBin.reset();
373 ERRORMSG(
"sum(Bins)= " << s <<
" != sum(R)= " << getTotalNum() <<
endl;);
375 if (emittedBins >= 2) {
376 for (
int i = 1; i < emittedBins; i++) {
377 if (binemitted_m[i - 1] != 0 && binemitted_m[i] != 0)
379 <<
"between bin " << i - 1 <<
" and " << i <<
endl);
385template <
class T,
unsigned Dim>
388 const int emittedBins = pbin_m->getLastemittedBin();
390 for (
int i = 0; i < emittedBins; i++)
392 for (
unsigned int n = 0;
n < getLocalNum();
n++) {
393 if ( this->Bin[
n] > -1 ) {
394 bingamma_m[this->Bin[
n]] +=
std::sqrt(1.0 +
dot(this->P[
n], this->P[
n]));
398 allreduce(*bingamma_m.get(), emittedBins, std::plus<double>());
400 for (
int i = 0; i < emittedBins; i++) {
401 if (pbin_m->getTotalNumPerBin(i) > 0) {
402 bingamma_m[i] /= pbin_m->getTotalNumPerBin(i);
406 INFOMSG(
"Bin " << i <<
" : particle number = " << pbin_m->getTotalNumPerBin(i)
407 <<
" gamma = " << bingamma_m[i] <<
endl);
412template <
class T,
unsigned Dim>
414 momentsComputer_m.computeDebyeLength(*
this, rmsDensity_m);
419template <
class T,
unsigned Dim>
421 return bingamma_m[bin];
425template <
class T,
unsigned Dim>
431template <
class T,
unsigned Dim>
433 this->Q =
where(
eq(this->Bin, bin), this->qi_m, 0.0);
437template <
class T,
unsigned Dim>
440 std::size_t localnum = 0;
443 for (
unsigned long k = 0; k < getLocalNum(); ++ k)
451 gather(&localnum, &globalPartPerNode_m[0], 1);
453 size_t npOutside = std::accumulate(globalPartPerNode_m.get(),
455 std::plus<size_t>());
470template <
class T,
unsigned Dim>
472 std::vector<double>& lineDensity,
473 std::pair<double, double>& meshInfo) {
475 get_bounds(rmin, rmax);
479 this->updateDomainLength(grid);
483 double length = rmax(2) - rmin(2);
484 double zmin = rmin(2) - dh_m * length, zmax = rmax(2) + dh_m * length;
485 double hz = (zmax - zmin) / (nBins - 2);
486 double perMeter = 1.0 / hz;
489 lineDensity.resize(nBins, 0.0);
490 std::fill(lineDensity.begin(), lineDensity.end(), 0.0);
492 const unsigned int lN = getLocalNum();
493 for (
unsigned int i = 0; i < lN; ++ i) {
494 const double z =
R[i](2) - 0.5 * hz;
495 unsigned int idx = (z - zmin) / hz;
496 double tau = (z - zmin) / hz - idx;
498 lineDensity[idx] += Q[i] * (1.0 - tau) * perMeter;
499 lineDensity[idx + 1] += Q[i] * tau * perMeter;
502 reduce(&(lineDensity[0]), &(lineDensity[0]) + nBins, &(lineDensity[0]),
OpAddAssign());
504 meshInfo.first = zmin;
505 meshInfo.second = hz;
509template <
class T,
unsigned Dim>
517 if ( !(
R.isDirty() || ID.isDirty() ) && stateOfLastBoundP_ == unit_state_)
return;
519 stateOfLastBoundP_ = unit_state_;
521 if (!isGridFixed()) {
522 const int dimIdx = (dcBeam_m? 2: 3);
530 this->updateDomainLength(nr_m);
532 get_bounds(rmin_m, rmax_m);
537 for (
int i = 0; i < dimIdx; i++) {
539 if (length < 1
e-10) {
543 rmax_m[i] += dh_m * length;
544 rmin_m[i] -= dh_m * length;
546 hr_m[i] = (rmax_m[i] - rmin_m[i]) / (nr_m[i] - 1);
549 rmax_m[2] = rmin_m[2] + periodLength_m;
550 hr_m[2] = periodLength_m / (nr_m[2] - 1);
552 for (
int i = 0; i < dimIdx; ++ i) {
553 volume *=
std::abs(rmax_m[i] - rmin_m[i]);
556 if (getIfBeamEmitting() && dist_m !=
nullptr) {
558 double percent =
std::max(1.0 / (nr_m[2] - 1), dist_m->getPercentageEmitted());
559 double length =
std::abs(rmax_m[2] - rmin_m[2]) / (1.0 + 2 * dh_m);
560 if (percent < 1.0 && percent > 0.0) {
561 rmax_m[2] -= dh_m * length;
562 rmin_m[2] = rmax_m[2] - length / percent;
566 rmax_m[2] += dh_m * length;
567 rmin_m[2] -= dh_m * length;
569 hr_m[2] = (rmax_m[2] - rmin_m[2]) / (nr_m[2] - 1);
573 if (volume < 1
e-21 && getTotalNum() > 1 &&
std::abs(
sum(Q)) > 0.0) {
578 if (hr_m[0] * hr_m[1] * hr_m[2] <= 0) {
582 Vector_t origin = rmin_m -
Vector_t(hr_m[0] / 2.0, hr_m[1] / 2.0, hr_m[2] / 2.0);
583 this->updateFields(hr_m, origin);
590 double gammaz =
sum(this->P)[2] / getTotalNum();
610template <
class T,
unsigned Dim>
614 const int dimIdx = 3;
617 std::unique_ptr<size_t[]> countLost;
619 const int tempN = pbin_m->getLastemittedBin();
620 countLost = std::unique_ptr<size_t[]>(
new size_t[tempN]);
621 for (
int ii = 0; ii < tempN; ii++) countLost[ii] = 0;
624 this->updateDomainLength(nr_m);
627 get_bounds(rmin_m, rmax_m);
630 len = rmax_m - rmin_m;
632 calcBeamParameters();
635 if (checkfactor != 0) {
638 Vector_t rmean = momentsComputer_m.getMeanPosition();
639 Vector_t rrms = momentsComputer_m.getStandardDeviationPosition();
640 if(checkfactor < 0) {
642 if (len[0] > checkfactor * rrms[0] ||
643 len[1] > checkfactor * rrms[1] ||
644 len[2] > checkfactor * rrms[2])
646 for(
unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
650 if (
std::abs(
R[ii](0) - rmean(0)) > checkfactor * rrms[0] ||
651 std::abs(
R[ii](1) - rmean(1)) > checkfactor * rrms[1] ||
652 std::abs(
R[ii](2) - rmean(2)) > checkfactor * rrms[2])
658 countLost[Bin[ii]] += 1 ;
667 if (len[0] > checkfactor * rrms[0] ||
668 len[2] > checkfactor * rrms[2])
670 for(
unsigned int ii = 0; ii < this->getLocalNum(); ii++) {
674 if (
std::abs(
R[ii](0) - rmean(0)) > checkfactor * rrms[0] ||
675 std::abs(
R[ii](2) - rmean(2)) > checkfactor * rrms[2])
681 countLost[Bin[ii]] += 1 ;
691 for (
int i = 0; i < dimIdx; i++) {
692 double length =
std::abs(rmax_m[i] - rmin_m[i]);
693 rmax_m[i] += dh_m * length;
694 rmin_m[i] -= dh_m * length;
695 hr_m[i] = (rmax_m[i] - rmin_m[i]) / (nr_m[i] - 1);
699 this->updateFields(hr_m, rmin_m);
702 pbin_m->updatePartInBin_cyc(countLost.get());
708 countTotalNumPerBunch();
718template <
class T,
unsigned Dim>
721 this->updateDomainLength(nr_m);
723 std::vector<size_t> tmpbinemitted;
728 const size_t localNum = getLocalNum();
730 double rzmean = momentsComputer_m.getMeanPosition()(2);
731 double rzrms = momentsComputer_m.getStandardDeviationPosition()(2);
732 const bool haveEnergyBins = weHaveEnergyBins();
733 if (haveEnergyBins) {
734 tmpbinemitted.resize(getNumberOfEnergyBins(), 0.0);
736 for (
unsigned int i = 0; i < localNum; i++) {
740 }
else if (haveEnergyBins) {
741 tmpbinemitted[Bin[i]]++;
747 calcBeamParameters();
748 gatherLoadBalanceStatistics();
750 if (haveEnergyBins) {
751 const int lastBin = dist_m->getLastEmittedEnergyBin();
752 for (
int i = 0; i < lastBin; i++) {
753 binemitted_m[i] = tmpbinemitted[i];
761template <
class T,
unsigned Dim>
764 std::unique_ptr<size_t[]> tmpbinemitted;
766 const size_t localNum = getLocalNum();
767 const size_t totalNum = getTotalNum();
770 if (weHaveEnergyBins()) {
771 tmpbinemitted = std::unique_ptr<size_t[]>(
new size_t[getNumberOfEnergyBins()]);
772 for (
int i = 0; i < getNumberOfEnergyBins(); i++) {
773 tmpbinemitted[i] = 0;
775 for (
unsigned int i = 0; i < localNum; i++) {
780 tmpbinemitted[Bin[i]]++;
784 for (
size_t i = 0; i < localNum; i++) {
793 performDestroy(
true);
796 calcBeamParameters();
797 gatherLoadBalanceStatistics();
799 if (weHaveEnergyBins()) {
800 const int lastBin = dist_m->getLastEmittedEnergyBin();
801 for (
int i = 0; i < lastBin; i++) {
802 binemitted_m[i] = tmpbinemitted[i];
805 size_t newTotalNum = getLocalNum();
808 setTotalNum(newTotalNum);
810 return totalNum - newTotalNum;
814template <
class T,
unsigned Dim>
819template <
class T,
unsigned Dim>
824template <
class T,
unsigned Dim>
829template <
class T,
unsigned Dim>
834template <
class T,
unsigned Dim>
840template <
class T,
unsigned Dim>
842 return this->
R[i](0);
846template <
class T,
unsigned Dim>
848 return this->
R[i](1);
852template <
class T,
unsigned Dim>
854 return this->
R[i](2);
858template <
class T,
unsigned Dim>
864template <
class T,
unsigned Dim>
870template <
class T,
unsigned Dim>
875template <
class T,
unsigned Dim>
878 this->getLocalBounds(rmin, rmax);
882 for (
unsigned int i = 0; i <
Dim; ++i) {
884 min[2*i + 1] = -rmax[i];
889 for (
unsigned int i = 0; i <
Dim; ++i) {
891 rmax[i] = -
min[2*i + 1];
896template <
class T,
unsigned Dim>
898 const size_t localNum = getLocalNum();
900 double maxValue = 1e8;
901 rmin =
Vector_t(maxValue, maxValue, maxValue);
902 rmax =
Vector_t(-maxValue, -maxValue, -maxValue);
908 for (
size_t i = 1; i < localNum; ++ i) {
909 for (
unsigned short d = 0; d < 3u; ++ d) {
910 if (rmin(d) >
R[i](d)) rmin(d) =
R[i](d);
911 else if (rmax(d) <
R[i](d)) rmax(d) =
R[i](d);
917template <
class T,
unsigned Dim>
920 get_bounds(rmin, rmax);
922 std::pair<Vector_t, double> sphere;
923 sphere.first = 0.5 * (rmin + rmax);
924 sphere.second =
std::sqrt(
dot(rmax - sphere.first, rmax - sphere.first));
930template <
class T,
unsigned Dim>
933 getLocalBounds(rmin, rmax);
935 std::pair<Vector_t, double> sphere;
936 sphere.first = 0.5 * (rmin + rmax);
937 sphere.second =
std::sqrt(
dot(rmax - sphere.first, rmax - sphere.first));
943template <
class T,
unsigned Dim>
947 size_t i = getLocalNum();
950 R[i] = particle.
getR();
951 P[i] = particle.
getP();
954 msg <<
"Created one particle i= " << i <<
endl;
958template <
class T,
unsigned Dim>
969template <
class T,
unsigned Dim>
971 R[ii] = particle.
getR();
972 P[ii] = particle.
getP();
976template <
class T,
unsigned Dim>
980 R[ii](2), Q[ii], M[ii]);
985template <
class T,
unsigned Dim>
987 double& axmax,
double& aymax) {
991 for (
unsigned int ii = 0; ii < getLocalNum(); ii++) {
993 particle = getParticle(ii);
1018template <
class T,
unsigned Dim>
1024template <
class T,
unsigned Dim>
1030template <
class T,
unsigned Dim>
1036template <
class T,
unsigned Dim>
1042template <
class T,
unsigned Dim>
1048template <
class T,
unsigned Dim>
1054template <
class T,
unsigned Dim>
1060template <
class T,
unsigned Dim>
1062 return momentsComputer_m.getMeanGamma();
1066template <
class T,
unsigned Dim>
1068 return momentsComputer_m.getMeanKineticEnergy();
1072template <
class T,
unsigned Dim>
1074 return momentsComputer_m.getTemperature();
1077template <
class T,
unsigned Dim>
1079 return momentsComputer_m.getDebyeLength();
1082template <
class T,
unsigned Dim>
1084 return momentsComputer_m.getPlasmaParameter();
1087template <
class T,
unsigned Dim>
1089 return rmsDensity_m;
1092template <
class T,
unsigned Dim>
1098template <
class T,
unsigned Dim>
1104template <
class T,
unsigned Dim>
1106 return momentsComputer_m.getMeanPosition();
1110template <
class T,
unsigned Dim>
1112 return momentsComputer_m.getStandardDeviationPosition();
1116template <
class T,
unsigned Dim>
1118 return momentsComputer_m.getStandardDeviationRP();
1122template <
class T,
unsigned Dim>
1124 return momentsComputer_m.getMeanPosition();
1128template <
class T,
unsigned Dim>
1130 return momentsComputer_m.getStandardDeviationMomentum();
1134template <
class T,
unsigned Dim>
1136 return momentsComputer_m.getMeanMomentum();
1140template <
class T,
unsigned Dim>
1143 return dist_m->get_pmean();
1145 double gamma = 0.1 / getM() + 1;
1150template <
class T,
unsigned Dim>
1152 return momentsComputer_m.getGeometricEmittance();
1156template <
class T,
unsigned Dim>
1158 return momentsComputer_m.getNormalizedEmittance();
1162template <
class T,
unsigned Dim>
1164 return momentsComputer_m.getHalo();
1167template <
class T,
unsigned Dim>
1169 return momentsComputer_m.get68Percentile();
1172template <
class T,
unsigned Dim>
1174 return momentsComputer_m.get95Percentile();
1177template <
class T,
unsigned Dim>
1179 return momentsComputer_m.get99Percentile();
1182template <
class T,
unsigned Dim>
1184 return momentsComputer_m.get99_99Percentile();
1187template <
class T,
unsigned Dim>
1189 return momentsComputer_m.getNormalizedEmittance68Percentile();
1192template <
class T,
unsigned Dim>
1194 return momentsComputer_m.getNormalizedEmittance95Percentile();
1197template <
class T,
unsigned Dim>
1199 return momentsComputer_m.getNormalizedEmittance99Percentile();
1202template <
class T,
unsigned Dim>
1204 return momentsComputer_m.getNormalizedEmittance99_99Percentile();
1207template <
class T,
unsigned Dim>
1209 return momentsComputer_m.getDx();
1213template <
class T,
unsigned Dim>
1215 return momentsComputer_m.getDy();
1219template <
class T,
unsigned Dim>
1221 return momentsComputer_m.getDDx();
1225template <
class T,
unsigned Dim>
1227 return momentsComputer_m.getDDy();
1231template <
class T,
unsigned Dim>
1237template <
class T,
unsigned Dim>
1243template <
class T,
unsigned Dim>
1247 globalPartPerNode_m[i] = 0;
1249 std::size_t localnum = getLocalNum();
1250 gather(&localnum, &globalPartPerNode_m[0], 1);
1254template <
class T,
unsigned Dim>
1256 return globalPartPerNode_m[p];
1260template <
class T,
unsigned Dim>
1266template <
class T,
unsigned Dim>
1270 get_bounds(rmin_m, rmax_m);
1271 momentsComputer_m.compute(*
this);
1275template <
class T,
unsigned Dim>
1277 return couplingConstant_m;
1281template <
class T,
unsigned Dim>
1283 couplingConstant_m =
c;
1287template <
class T,
unsigned Dim>
1290 if (getTotalNum() != 0)
1293 WARNMSG(
"Could not set total charge in PartBunch::setCharge based on getTotalNum" <<
endl);
1297template <
class T,
unsigned Dim>
1303template <
class T,
unsigned Dim>
1305 massPerParticle_m = mass;
1306 if (getTotalNum() != 0)
1310template <
class T,
unsigned Dim>
1312 massPerParticle_m = mass;
1316template <
class T,
unsigned Dim>
1322template <
class T,
unsigned Dim>
1327template <
class T,
unsigned Dim>
1329 return massPerParticle_m;
1332template <
class T,
unsigned Dim>
1335 fs_m->initSolver(
this);
1342 this->initialize(fs_m->getFieldLayout());
1350template <
class T,
unsigned Dim>
1353 return fs_m->hasValidSolver();
1360template <
class T,
unsigned Dim>
1363 return fs_m->getFieldSolverType();
1370template <
class T,
unsigned Dim>
1376template <
class T,
unsigned Dim>
1378 return stepsPerTurn_m;
1382template <
class T,
unsigned Dim>
1384 globalTrackStep_m =
n;
1388template <
class T,
unsigned Dim>
1390 return globalTrackStep_m;
1394template <
class T,
unsigned Dim>
1396 localTrackStep_m =
n;
1400template <
class T,
unsigned Dim>
1402 globalTrackStep_m++; localTrackStep_m++;
1406template <
class T,
unsigned Dim>
1408 return localTrackStep_m;
1412template <
class T,
unsigned Dim>
1415 bunchTotalNum_m.resize(
n);
1416 bunchLocalNum_m.resize(
n);
1420template <
class T,
unsigned Dim>
1426template <
class T,
unsigned Dim>
1429 bunchTotalNum_m[
n] = totalnum;
1433template <
class T,
unsigned Dim>
1436 return bunchTotalNum_m[
n];
1440template <
class T,
unsigned Dim>
1443 bunchLocalNum_m[
n] = localnum;
1447template <
class T,
unsigned Dim>
1450 return bunchLocalNum_m[
n];
1454template <
class T,
unsigned Dim>
1456 bunchTotalNum_m.clear();
1457 bunchTotalNum_m.resize(numBunch_m);
1458 bunchLocalNum_m.clear();
1459 bunchLocalNum_m.resize(numBunch_m);
1461 for (
size_t i = 0; i < this->getLocalNum(); ++i) {
1463 ++bunchLocalNum_m[this->bunchNum[i]];
1466 allreduce(bunchLocalNum_m.data(), bunchTotalNum_m.data(),
1467 bunchLocalNum_m.size(), std::plus<size_t>());
1469 size_t totalnum = std::accumulate(bunchTotalNum_m.begin(),
1470 bunchTotalNum_m.end(), 0);
1472 if ( totalnum != this->getTotalNum() )
1473 throw OpalException(
"PartBunchBase::countTotalNumPerBunch()",
1474 "Sum of total number of particles per bunch (" +
1475 std::to_string(totalnum) +
") != total number of particles (" +
1476 std::to_string(this->getTotalNum()) +
")");
1480template <
class T,
unsigned Dim>
1482 globalMeanR_m = globalMeanR;
1486template <
class T,
unsigned Dim>
1488 return globalMeanR_m;
1492template <
class T,
unsigned Dim>
1495 globalToLocalQuaternion_m = globalToLocalQuaternion;
1499template <
class T,
unsigned Dim>
1501 return globalToLocalQuaternion_m;
1505template <
class T,
unsigned Dim>
1507 SteptoLastInj_m =
n;
1511template <
class T,
unsigned Dim>
1513 return SteptoLastInj_m;
1517template <
class T,
unsigned Dim>
1520 const int emittedBins = pbin_m->getLastemittedBin();
1521 double phi[emittedBins];
1522 double px[emittedBins];
1523 double py[emittedBins];
1524 double meanPhi = 0.0;
1526 for (
int ii = 0; ii < emittedBins; ii++) {
1532 for (
unsigned int ii = 0; ii < getLocalNum(); ii++) {
1533 px[Bin[ii]] += P[ii](0);
1534 py[Bin[ii]] += P[ii](1);
1539 for (
int ii = 0; ii < emittedBins; ii++) {
1540 phi[ii] = calculateAngle(px[ii], py[ii]);
1545 meanPhi /= emittedBins;
1556template <
class T,
unsigned Dim>
1562 int maxbin = pbin_m->getNBins();
1563 size_t partInBin[maxbin];
1564 for (
int ii = 0; ii < maxbin; ii++) partInBin[ii] = 0;
1566 double pMin0 = 1.0e9;
1568 double maxbinIndex = 0;
1570 for (
unsigned long int n = 0;
n < getLocalNum();
n++) {
1572 if (pMin0 > temp_betagamma)
1573 pMin0 = temp_betagamma;
1578 double asinh0 = std::asinh(pMin);
1579 for (
unsigned long int n = 0;
n < getLocalNum();
n++) {
1582 int itsBinID =
std::floor((std::asinh(temp_betagamma) - asinh0) / eta + 1.0E-6);
1584 if (maxbinIndex < itsBinID) {
1585 maxbinIndex = itsBinID;
1588 if (itsBinID >= maxbin) {
1589 ERRORMSG(
"The bin number limit is " << maxbin <<
", please increase the energy interval and try again" <<
endl);
1592 partInBin[itsBinID]++;
1597 pbin_m->resetPartInBin_cyc(partInBin, maxbinIndex);
1607template <
class T,
unsigned Dim>
1609 int maxbin = pbin_m->getNBins();
1610 std::size_t partInBin[maxbin];
1611 for (
int i = 0; i < maxbin; ++i) {
1615 for (std::size_t i = 0; i < getLocalNum(); ++i) {
1616 partInBin[Bin[i]]++;
1620 pbin_m->resetPartInBin_cyc(partInBin, numBunch_m - 1);
1629template <
class T,
unsigned Dim>
1631 return reference->getQ();
1635template <
class T,
unsigned Dim>
1637 return reference->getM();
1641template <
class T,
unsigned Dim>
1643 return reference->getP();
1647template <
class T,
unsigned Dim>
1649 return reference->getE();
1653template <
class T,
unsigned Dim>
1655 return refPOrigin_m;
1658template <
class T,
unsigned Dim>
1660 refPOrigin_m = origin;
1664template <
class T,
unsigned Dim>
1669template <
class T,
unsigned Dim>
1675template <
class T,
unsigned Dim>
1677 return dist_m->getType();
1681template <
class T,
unsigned Dim>
1683 const_cast<PartData *
>(reference)->setQ(q);
1687template <
class T,
unsigned Dim>
1689 const_cast<PartData *
>(reference)->setM(m);
1693template <
class T,
unsigned Dim>
1695 return momentsComputer_m.getStdKineticEnergy();
1699template <
class T,
unsigned Dim>
1701 return reference->getBeta();
1705template <
class T,
unsigned Dim>
1707 return reference->getGamma();
1711template <
class T,
unsigned Dim>
1717template <
class T,
unsigned Dim>
1723template <
class T,
unsigned Dim>
1729template <
class T,
unsigned Dim>
1735template <
class T,
unsigned Dim>
1737 return dist_m->getEmissionDeltaT();
1741template <
class T,
unsigned Dim>
1743 binemitted_m[binNumber]++;
1747template <
class T,
unsigned Dim>
1749 momentsComputer_m.computeMeanKineticEnergy(*
this);
1752template <
class T,
unsigned Dim>
1755 if (getTotalNum() != 0) {
1758 double pathLength = get_sPos();
1760 os << std::scientific;
1762 os <<
"* ************** B U N C H ********************************************************* \n";
1763 os <<
"* NP = " << getTotalNum() <<
"\n";
1770 if (getTotalNum() >= 2) {
1772 os <<
"* rms momenta = " << std::setw(12) << std::setprecision(5) << get_prms() <<
" [beta gamma]\n";
1774 os <<
"* mean momenta = " << std::setw(12) << std::setprecision(5) << get_pmean() <<
" [beta gamma]\n";
1775 os <<
"* rms emittance = " << std::setw(12) << std::setprecision(5) << get_emit() <<
" (not normalized)\n";
1776 os <<
"* rms correlation = " << std::setw(12) << std::setprecision(5) << get_rprms() <<
"\n";
1779 os <<
"* dh = " << std::setw(13) << std::setprecision(5) << dh_m * 100 <<
" [%]\n";
1783 os <<
"* ********************************************************************************** " <<
endl;
1790template <
class T,
unsigned Dim>
1792 double thetaXY =
atan2(y, x);
1798template <
class T,
unsigned Dim>
1808template <
class T,
unsigned Dim>
1810 throw OpalException(
"PartBunchBase<T, Dim>::runTests() ",
"No test supported.");
1814template <
class T,
unsigned Dim>
1819template <
class T,
unsigned Dim>
1821 if (i >= getLocalNum() || j >= getLocalNum() || i == j)
return;
1823 std::swap(
R[i],
R[j]);
1824 std::swap(P[i], P[j]);
1825 std::swap(Q[i], Q[j]);
1826 std::swap(M[i], M[j]);
1827 std::swap(Phi[i], Phi[j]);
1828 std::swap(Ef[i], Ef[j]);
1829 std::swap(Eftmp[i], Eftmp[j]);
1830 std::swap(Bf[i], Bf[j]);
1831 std::swap(Bin[i], Bin[j]);
1832 std::swap(dt[i], dt[j]);
1833 std::swap(PType[i], PType[j]);
1834 std::swap(POrigin[i], POrigin[j]);
1835 std::swap(TriID[i], TriID[j]);
1836 std::swap(cavityGapCrossed[i], cavityGapCrossed[j]);
1837 std::swap(bunchNum[i], bunchNum[j]);
1841template <
class T,
unsigned Dim>
1843 throw OpalException(
"PartBunchBase<T, Dim>::setBCAllPeriodic() ",
"Not supported BC.");
1847template <
class T,
unsigned Dim>
1849 throw OpalException(
"PartBunchBase<T, Dim>::setBCAllOpen() ",
"Not supported BC.");
1853template <
class T,
unsigned Dim>
1855 throw OpalException(
"PartBunchBase<T, Dim>::setBCForDCBeam() ",
"Not supported BC.");
1859template <
class T,
unsigned Dim>
1864template <
class T,
unsigned Dim>
1893 globalPartPerNode_m = std::unique_ptr<size_t[]>(
new size_t[
Ippl::getNodes()]);
1899template <
class T,
unsigned Dim>
1901 return pbase_m->getTotalNum();
1904template <
class T,
unsigned Dim>
1906 return pbase_m->getLocalNum();
1910template <
class T,
unsigned Dim>
1912 return pbase_m->getDestroyNum();
1915template <
class T,
unsigned Dim>
1917 return pbase_m->getGhostNum();
1920template <
class T,
unsigned Dim>
1922 pbase_m->setTotalNum(
n);
1925template <
class T,
unsigned Dim>
1927 pbase_m->setLocalNum(
n);
1930template <
class T,
unsigned Dim>
1932 return pbase_m->getLayout();
1935template <
class T,
unsigned Dim>
1937 return pbase_m->getLayout();
1940template <
class T,
unsigned Dim>
1942 return pbase_m->getUpdateFlag(f);
1945template <
class T,
unsigned Dim>
1947 pbase_m->setUpdateFlag(f, val);
1950template <
class T,
unsigned Dim>
1952 return pbase_m->singleInitNode();
1955template <
class T,
unsigned Dim>
1960template <
class T,
unsigned Dim>
1969template <
class T,
unsigned Dim>
1972 pbase_m->update(canSwap);
1978template <
class T,
unsigned Dim>
1980 pbase_m->createWithID(
id);
1983template <
class T,
unsigned Dim>
1988template <
class T,
unsigned Dim>
1990 pbase_m->globalCreate(np);
1993template <
class T,
unsigned Dim>
1995 pbase_m->destroy(M, I, doNow);
1998template <
class T,
unsigned Dim>
2000 pbase_m->performDestroy(updateLocalNum);
2003template <
class T,
unsigned Dim>
2005 pbase_m->ghostDestroy(M, I);
2008template <
class T,
unsigned Dim>
2013template <
class T,
unsigned Dim>
2018 for (
unsigned int i = 0; i <
Dim; i++) {
2019 rpmean(2*i)= get_rmean()(i);
2020 rpmean((2*i)+1)= get_pmean()(i);
2024 for (
unsigned int i = 0; i < 2 *
Dim; i++) {
2025 for (
unsigned int j = 0; j <= i; j++) {
2026 sigmaMatrix[i][j] = momentsComputer_m.getMoments6x6()[i][j] - rpmean(i) * rpmean(j);
2027 sigmaMatrix[j][i] = sigmaMatrix[i][j];
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
Inform & operator<<(Inform &os, PartBunchBase< T, Dim > &p)
Tps< T > pow(const Tps< T > &x, int y)
Integer power.
Tps< T > sqrt(const Tps< T > &x)
Square root.
void allreduce(const T *input, T *output, int count, Op op)
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
void gather(const T *input, T *output, int count, int root=0)
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
T::PETE_Expr_t::PETE_Return_t min(const PETE_Expr< T > &expr, NDIndex< D > &loc)
PETE_TBTree< OpNE, Index::PETE_Expr_t, PETE_Scalar< double > > ne(const Index &idx, double x)
T::PETE_Expr_t::PETE_Return_t sum(const PETE_Expr< T > &expr)
PETE_TUTree< FnFloor, typename T::PETE_Expr_t > floor(const PETE_Expr< T > &l)
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
PETE_TTTree< OpWhere, typename Cond_t::PETE_Expr_t, typename True_t::PETE_Expr_t, PETE_Scalar< Vektor< T, Dim > > > where(const PETE_Expr< Cond_t > &c, const PETE_Expr< True_t > &t, const Vektor< T, Dim > &f)
void bounds(const PETE_Expr< T1 > &expr, Vektor< T2, D > &minval, Vektor< T2, D > &maxval)
PETE_TBTree< FnArcTan2, PETE_Scalar< Vektor< T1, Dim > >, typename T2::PETE_Expr_t > atan2(const Vektor< T1, Dim > &l, const PETE_Expr< T2 > &r)
Inform & level2(Inform &inf)
Inform & endl(Inform &inf)
Inform & level1(Inform &inf)
constexpr double two_pi
The value of.
constexpr double e
The value of.
constexpr double c
The velocity of light in m/s.
std::string getChargeString(double charge, unsigned int precision=3)
std::string getEnergyString(double energyInMeV, unsigned int precision=3)
std::string getTimeString(double time, unsigned int precision=3)
std::string getLengthString(double spos, unsigned int precision=3)
bool eq(double x, double y)
boost::function< boost::tuple< double, bool >(arguments_t)> type
ParticleLayout< T, Dim > & getLayout()
void setEnergyBins(int numberOfEnergyBins)
double get_meanKineticEnergy() const
void setDistribution(Distribution *d, std::vector< Distribution * > addedDistributions, size_t &np)
virtual void resetInterpolationCache(bool clearCache=false)
void setPType(const std::string &type)
const PartData * getReference() const
ParticleOrigin getPOrigin() const
void setMass(double mass)
int getSteptoLastInj() const
Vector_t get_99Percentile() const
void boundp_destroyCycl()
void get_bounds(Vector_t &rmin, Vector_t &rmax) const
bool resetPartBinID2(const double eta)
reset Bin[] for each particle according to the method given in paper PAST-AB(064402) by G....
void setNumBunch(short n)
double getMassPerParticle() const
double getQ() const
Access to reference data.
double getCouplingConstant() const
Vector_t get_normalizedEps_99Percentile() const
double getChargePerParticle() const
get the macro particle charge
virtual void set_meshEnlargement(double dh)
virtual double getBeta(int i)
virtual double getPx0(int i)
size_t getLocalNum() const
FMatrix< double, 2 *Dim, 2 *Dim > getSigmaMatrix() const
void setLocalTrackStep(long long n)
step in a TRACK command
void setLocalBinCount(size_t num, int bin)
double get_debyeLength() const
bool getUpdateFlag(UpdateFlags_t f) const
void setParticle(FVector< double, 6 > z, int ii)
size_t emitParticles(double eZ)
Emit particles in the given bin i.e. copy the particles from the bin structure into the particle cont...
void maximumAmplitudes(const FMatrix< double, 6, 6 > &D, double &axmax, double &aymax)
Return maximum amplitudes.
double get_plasmaParameter() const
void setMassZeroPart(double mass)
size_t getTotalNum() const
virtual void updateFields(const Vector_t &hr, const Vector_t &origin)
Vector_t get_95Percentile() const
virtual void setZ(int i, double zcoo)
double calcMeanPhi()
calculate average angle of longitudinal direction of bins
Quaternion_t getGlobalToLocalQuaternion()
void setBeamFrequency(double v)
Vector_t get_normalizedEps_68Percentile() const
void switchToUnitlessPositions(bool use_dt_per_particle=false)
size_t calcNumPartsOutside(Vector_t x)
returns the number of particles outside of a box defined by x
Inform & print(Inform &os)
FieldSolverType getFieldSolverType() const
Return the fieldsolver type if we have a fieldsolver.
double getEmissionDeltaT()
void setLocalNumPerBunch(size_t numpart, short n)
size_t getLocalNumPerBunch(short n) const
void setGlobalToLocalQuaternion(Quaternion_t globalToLocalQuaternion)
Vector_t get_rrms() const
double get_rmsDensity() const
virtual double getPy(int i)
double getBinGamma(int bin)
Get gamma of one bin.
double getInitialBeta() const
int getLastEmittedEnergyBin()
void calcBeamParameters()
void getLocalBounds(Vector_t &rmin, Vector_t &rmax) const
void setSteptoLastInj(int n)
virtual double getX0(int i)
double getInitialGamma() const
virtual void setBCAllPeriodic()
void setChargeZeroPart(double q)
virtual double getY0(int i)
ParticleType getPType() const
Vector_t get_origin() const
Vector_t get_prms() const
void createWithID(unsigned id)
PartBunchBase(AbstractParticle< T, Dim > *pb, const PartData *ref)
double getCharge() const
get the total charge per simulation particle
virtual double getPy0(int i)
double get_temperature() const
size_t getTotalNumPerBunch(short n) const
virtual double getPx(int i)
Vector_t get_pmean_Distribution() const
size_t getNumberOfEmissionSteps()
void calcLineDensity(unsigned int nBins, std::vector< double > &lineDensity, std::pair< double, double > &meshInfo)
calculates the 1d line density (not normalized) and append it to a file.
Vector_t get_68Percentile() const
void setup(AbstractParticle< T, Dim > *pb)
virtual void setBCAllOpen()
void setTotalNum(size_t n)
long long getLocalTrackStep() const
void setGlobalMeanR(Vector_t globalMeanR)
void setCouplingConstant(double c)
short getNumBunch() const
virtual void do_binaryRepart()
Vector_t get_norm_emit() const
void calcGammas()
Compute the gammas of all bins.
int getStepsPerTurn() const
Vector_t get_rprms() const
void setPOrigin(ParticleOrigin)
Vector_t get_maxExtent() const
Vector_t get_normalizedEps_95Percentile() const
void gatherLoadBalanceStatistics()
void countTotalNumPerBunch()
DistributionType getDistType() const
void get_PBounds(Vector_t &min, Vector_t &max) const
virtual void setBinCharge(int bin, double q)
Set the charge of one bin to the value of q and all other to zero.
Vector_t get_halo() const
void performDestroy(bool updateLocalNum=false)
Vector_t getGlobalMeanR()
void ghostDestroy(size_t M, size_t I)
int getNumberOfEnergyBins()
void setGlobalTrackStep(long long n)
step in multiple TRACK commands
std::pair< Vector_t, double > getBoundingSphere()
Vector_t get_centroid() const
Vector_t get_99_99Percentile() const
void push_back(OpalParticle const &p)
virtual double getX(int i)
virtual double getPz(int i)
void setLocalNum(size_t n)
Vector_t get_pmean() const
void setTotalNumPerBunch(size_t numpart, short n)
virtual void setBCForDCBeam()
void destroy(size_t M, size_t I, bool doNow=false)
void globalCreate(size_t np)
Vector_t get_normalizedEps_99_99Percentile() const
OpalParticle getParticle(int ii)
virtual void swap(unsigned int i, unsigned int j)
void setPBins(PartBins *pbin)
size_t getDestroyNum() const
virtual double getGamma(int i)
virtual void setSolver(FieldSolver *fs)
long long getGlobalTrackStep() const
size_t getGhostNum() const
void setTEmission(double t)
size_t getLoadBalance(int p) const
std::pair< Vector_t, double > getLocalBoundingSphere()
virtual Vector_t get_hr() const
Vector_t get_emit() const
double calculateAngle(double x, double y)
angle range [0~2PI) degree
void calcDebyeLength()
Compute the (global) Debye length for the beam.
void iterateEmittedBin(int binNumber)
bool singleInitNode() const
void setUpdateFlag(UpdateFlags_t f, bool val)
virtual double getY(int i)
void switchOffUnitlessPositions(bool use_dt_per_particle=false)
void setStepsPerTurn(int n)
virtual double getZ(int i)
Vector_t get_rmean() const
static OpalData * getInstance()
double getPy() const
Get vertical momentum (no dimension).
const Vector_t & getR() const
Get position in m.
const Vector_t & getP() const
Get momentum.
double getPz() const
Get relative momentum error (no dimension).
double getY() const
Get vertical displacement in m.
double getZ() const
Get longitudinal displacement c*t in m.
double getPx() const
Get horizontal momentum (no dimension).
double getX() const
Get horizontal position in m.
An abstract sequence of beam line components.
A templated representation for vectors.
static ParticleType getParticleType(const std::string &str)
void createOpalCycl(PartBunchBase< double, 3 > *beam, size_t numberOfParticles, double current, const Beamline &bl)
void createOpalT(PartBunchBase< double, 3 > *beam, std::vector< Distribution * > addedDistributions, size_t &numberOfParticles)
The base class for all OPAL exceptions.
virtual void addAttribute(ParticleAttribBase &pa)=0
void setCacheDimension(int d, T length)
std::ios_base::fmtflags FmtFlags_t
virtual const char * what() const
virtual const std::string & where() const
static TimerRef getTimer(const char *nm)
static void stopTimer(TimerRef t)
static void startTimer(TimerRef t)
Vektor< double, 3 > Vector_t