OPAL (Object Oriented Parallel Accelerator Library)  2024.1
OPAL
Stripper.cpp
Go to the documentation of this file.
1 //
2 // Class Stripper
3 // The Stripper element defines the interface for a stripping foil
4 //
5 // Copyright (c) 2011, Jianjun Yang, Paul Scherrer Institut, Villigen PSI, Switzerland
6 // Copyright (c) 2017-2020, Paul Scherrer Institut, Villigen PSI, Switzerland
7 // All rights reserved
8 //
9 // Implemented as part of the PhD thesis
10 // "Beam dynamics in high intensity cyclotrons including neighboring bunch effects"
11 // and the paper
12 // "Beam dynamics in high intensity cyclotrons including neighboring bunch effects:
13 // Model, implementation, and application"
14 // (https://journals.aps.org/prab/pdf/10.1103/PhysRevSTAB.13.064201)
15 //
16 // This file is part of OPAL.
17 //
18 // OPAL is free software: you can redistribute it and/or modify
19 // it under the terms of the GNU General Public License as published by
20 // the Free Software Foundation, either version 3 of the License, or
21 // (at your option) any later version.
22 //
23 // You should have received a copy of the GNU General Public License
24 // along with OPAL. If not, see <https://www.gnu.org/licenses/>.
25 //
26 #include "AbsBeamline/Stripper.h"
27 
30 #include "Physics/Physics.h"
31 #include "Physics/Units.h"
32 #include "Structure/LossDataSink.h"
33 
34 extern Inform *gmsg;
35 extern Inform *gmsgALL;
36 
38 {}
39 
40 Stripper::Stripper(const std::string &name):
41  PluginElement(name),
42  opcharge_m(0.0),
43  opmass_m(0.0),
44  opyield_m(1.0),
45  stop_m(true)
46 {}
47 
49  PluginElement(right),
50  opcharge_m(right.opcharge_m),
51  opmass_m(right.opmass_m),
52  opyield_m(right.opyield_m),
53  stop_m(right.stop_m)
54 {}
55 
57 
58 void Stripper::accept(BeamlineVisitor &visitor) const {
59  visitor.visitStripper(*this);
60 }
61 
63  *gmsg << "* Finalize stripper " << getName() << endl;
64 }
65 
66 void Stripper::setOPCharge(double charge) {
67  opcharge_m = charge;
68 }
69 
70 void Stripper::setOPMass(double mass) {
71  opmass_m = mass;
72 }
73 
74 void Stripper::setOPYield(double yield) {
75  opyield_m = yield;
76 }
77 
78 void Stripper::setStop(bool stopflag) {
79  stop_m = stopflag;
80 }
81 
82 double Stripper::getOPCharge() const {
83  return opcharge_m;
84 }
85 
86 double Stripper::getOPMass() const {
87  return opmass_m;
88 }
89 
90 double Stripper::getOPYield() const {
91  return opyield_m;
92 }
93 
94 bool Stripper::getStop () const {
95  return stop_m;
96 }
97 
99  Vector_t rmin, rmax, strippoint;
100  bunch->get_bounds(rmin, rmax);
101  // interested in absolute maximum
102  double xmax = std::max(std::abs(rmin(0)), std::abs(rmax(0)));
103  double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
104  double rbunch_max = std::hypot(xmax, ymax);
105 
106  if (rbunch_max > rmin_m - 1e-2) {
107  return true;
108  }
109  return false;
110 }
111 
112 //change the stripped particles to outcome particles
113 bool Stripper::doCheck(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
114  bool flagNeedUpdate = false;
115  Vector_t strippoint;
116 
117  size_t count = 0;
118  size_t tempnum = bunch->getLocalNum();
119 
120  for (unsigned int i = 0; i < tempnum; ++i) {
121  if (bunch->POrigin[i] != ParticleOrigin::REGULAR) continue;
122 
123  double tangle = calculateIncidentAngle(bunch->P[i](0), bunch->P[i](1));
124  changeWidth(bunch, i, tstep, tangle);
125  int pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
126  if (pflag == 0) continue;
127 
128  // dist1 > 0, right hand, dt > 0; dist1 < 0, left hand, dt < 0
129  double dist1 = (A_m * bunch->R[i](0) + B_m * bunch->R[i](1) + C_m) / R_m; // [m]
130  double dist2 = dist1 * std::sqrt(1.0 + 1.0 / tangle / tangle);
131  double dt = dist2 / (std::sqrt(1.0 - 1.0 / (1.0 + dot(bunch->P[i], bunch->P[i]))) * Physics::c);
132  strippoint(0) = (B_m * B_m * bunch->R[i](0) - A_m * B_m* bunch->R[i](1) - A_m * C_m) / (R_m * R_m);
133  strippoint(1) = (A_m * A_m * bunch->R[i](1) - A_m * B_m* bunch->R[i](0) - B_m * C_m) / (R_m * R_m);
134  strippoint(2) = bunch->R[i](2);
135  lossDs_m->addParticle(OpalParticle(bunch->ID[i],
136  strippoint, bunch->P[i],
137  t+dt, bunch->Q[i], bunch->M[i]),
138  std::make_pair(turnnumber, bunch->bunchNum[i]));
139 
140  flagNeedUpdate = true;
141  if (stop_m) {
142  bunch->Bin[i] = -1;
143  *gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i] << " is deleted by stripper " << getName() << endl;
144  } else {
145  *gmsgALL << level4 << getName() << ": Particle " << bunch->ID[i] << " collide in stripper " << getName() << endl;
146  // change charge and mass of PartData when the reference particle hits the stripper.
147  if (bunch->ID[i] == 0)
149 
150  // change the mass and charge
151  bunch->M[i] = opmass_m;
152  bunch->Q[i] = opcharge_m * Physics::q_e;
153  bunch->POrigin[i] = ParticleOrigin::STRIPPED;
154 
155  int j = 1;
156  //create new particles
157  while (j < opyield_m){
158  bunch->create(1);
159  size_t index = tempnum + count;
160  bunch->R[index] = bunch->R[i];
161  bunch->P[index] = bunch->P[i];
162  bunch->Q[index] = bunch->Q[i];
163  bunch->M[index] = bunch->M[i];
164  // once the particle is stripped, change POrigin from 0 to 1 as a flag so as to avoid repetitive stripping.
165  bunch->POrigin[index] = ParticleOrigin::STRIPPED;
166  if (bunch->weHaveBins())
167  bunch->Bin[index] = bunch->Bin[i];
168 
169  count++;
170  j++;
171  }
172  }
173  }
174  return flagNeedUpdate;
175 }
176 
177 bool Stripper::doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpdate) {
178  reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
179 
180  if (!stop_m){
181  // change charge and mass of PartData when the reference particle hits the stripper.
182  if (bunch->getPOrigin() == ParticleOrigin::STRIPPED) {
183  bunch->resetM(opmass_m * Units::GeV2eV);
184  bunch->resetQ(opcharge_m); // elementary charge
185  }
186  }
187 
188  return flagNeedUpdate;
189 }
190 
192  return ElementType::STRIPPER;
193 }
Tps< T > sqrt(const Tps< T > &x)
Square root.
Definition: TpsMath.h:91
constexpr double c
The velocity of light in m/s.
Definition: Physics.h:45
bool getStop() const
Definition: Stripper.cpp:94
double getOPYield() const
Definition: Stripper.cpp:90
ParticleAttrib< ParticleOrigin > POrigin
bool stop_m
Flag if particles should be stripped or stopped.
Definition: Stripper.h:74
double getOPCharge() const
Definition: Stripper.cpp:82
PETE_TUTree< FnAbs, typename T::PETE_Expr_t > abs(const PETE_Expr< T > &l)
ParticleAttrib< Vector_t > P
ParticleAttrib< short > bunchNum
void get_bounds(Vector_t &rmin, Vector_t &rmax) const
ParticleAttrib< double > M
Stripper()
Definition: Stripper.cpp:37
void setStop(bool stopflag)
Definition: Stripper.cpp:78
virtual ElementType getType() const override
Override implementation of PluginElement.
Definition: Stripper.cpp:191
void resetM(double m)
void setOPCharge(double charge)
Member variable access.
Definition: Stripper.cpp:66
virtual const std::string & getName() const
Get element name.
void setOPMass(double mass)
Definition: Stripper.cpp:70
Inform & endl(Inform &inf)
Definition: Inform.cpp:42
std::unique_ptr< LossDataSink > lossDs_m
Pointer to Loss instance.
double opcharge_m
Charge number of the out-coming particle.
Definition: Stripper.h:71
T::PETE_Expr_t::PETE_Return_t max(const PETE_Expr< T > &expr, NDIndex< D > &loc)
Definition: ReductionLoc.h:84
virtual void doFinalise() override
Virtual hook for finalise.
Definition: Stripper.cpp:62
double getOPMass() const
Definition: Stripper.cpp:86
bool weHaveBins() const
double rmin_m
radius closest to the origin
ParticleOrigin getPOrigin() const
double opmass_m
Mass of the out-coming particle.
Definition: Stripper.h:72
int checkPoint(const double &x, const double &y) const
Checks if coordinate is within element.
Inform & level4(Inform &inf)
Definition: Inform.cpp:48
size_t getLocalNum() const
Definition: Inform.h:42
virtual void visitStripper(const Stripper &)=0
Apply the algorithm to a particle stripper.
ElementType
Definition: ElementBase.h:88
Inform * gmsgALL
Definition: Main.cpp:71
void setPOrigin(ParticleOrigin)
void setOPYield(double yield)
Definition: Stripper.cpp:74
ParticleAttrib< double > Q
virtual bool doFinaliseCheck(PartBunchBase< double, 3 > *bunch, bool flagNeedUpdate) override
Virtual hook for finaliseCheck.
Definition: Stripper.cpp:177
constexpr double q_e
The elementary charge in As.
Definition: Physics.h:69
double calculateIncidentAngle(double xp, double yp) const
Calculate angle of particle/bunch wrt to element.
void changeWidth(PartBunchBase< double, 3 > *bunch, int i, const double tstep, const double tangle)
Change probe width depending on step size and angle of particle.
virtual bool doPreCheck(PartBunchBase< double, 3 > *) override
Virtual hook for preCheck.
Definition: Stripper.cpp:98
const std::string name
void resetQ(double q)
ParticlePos_t & R
ParticleIndex_t & ID
double C_m
Geometric lengths used in calculations.
constexpr double GeV2eV
Definition: Units.h:68
constexpr double e
The value of .
Definition: Physics.h:39
void create(size_t M)
ParticleAttrib< int > Bin
virtual void accept(BeamlineVisitor &) const override
Apply visitor to Stripper.
Definition: Stripper.cpp:58
virtual bool doCheck(PartBunchBase< double, 3 > *bunch, const int turnnumber, const double t, const double tstep) override
Record hits when bunch particles pass.
Definition: Stripper.cpp:113
virtual ~Stripper()
Definition: Stripper.cpp:56
double opyield_m
Yield of the out-coming particle.
Definition: Stripper.h:73
Inform * gmsg
Definition: Main.cpp:70
double dot(const Vector3D &lhs, const Vector3D &rhs)
Vector dot product.
Definition: Vector3D.cpp:118
bool reduce(Communicate &, InputIterator, InputIterator, OutputIterator, const ReduceOp &, bool *IncludeVal=0)
Definition: GlobalComm.hpp:55