Skip to content

Commit

Permalink
Merge pull request #88 from shishlo/master
Browse files Browse the repository at this point in the history 
Changes were made to fix usage of k1(n) strength of magnet multipoles…
  • Loading branch information
@azukov
azukov authored Mar 15, 2024
2 parents 8e54b39 + 35e0ae5 commit 437306a
Show file tree
Hide file tree
Showing 10 changed files with 160 additions and 110 deletions.
41 changes: 25 additions & 16 deletions py/orbit/py_linac/lattice/LinacAccNodes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -297,13 +297,15 @@ def __init__(self, name = "quad"):
self.setType("linacQuad")

def fringeIN(node,paramsDict):
# B*rho = 3.335640952*momentum [T*m] if momentum in GeV/c
# B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c
usageIN = node.getUsage()
if(not usageIN):
return
bunch = paramsDict["bunch"]
bunch = paramsDict["bunch"]
charge = bunch.charge()
momentum = bunch.getSyncParticle().momentum()
kq = node.getParam("dB/dr")/(3.335640952*momentum)
#---- The charge sign will be accounted for inside tracking module functions.
kq = node.getParam("dB/dr")/bunch.B_Rho()
poleArr = node.getParam("poles")
klArr = node.getParam("kls")
skewArr = node.getParam("skews")
Expand All @@ -318,12 +320,15 @@ def fringeIN(node,paramsDict):
TPB.multpfringeIN(bunch,pole,k,skew)

def fringeOUT(node,paramsDict):
# B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c
usageOUT = node.getUsage()
if(not usageOUT):
return
bunch = paramsDict["bunch"]
charge = bunch.charge()
momentum = bunch.getSyncParticle().momentum()
kq = node.getParam("dB/dr")/(3.335640952*momentum)
#---- The charge sign will be accounted for inside tracking module functions
kq = node.getParam("dB/dr")/bunch.B_Rho()
poleArr = node.getParam("poles")
klArr = node.getParam("kls")
skewArr = node.getParam("skews")
Expand Down Expand Up @@ -400,9 +405,14 @@ def track(self, paramsDict):
The Quad Combined Function TEAPOT class implementation
of the AccNode class track(probe) method.
"""
bunch = paramsDict["bunch"]
bunch = paramsDict["bunch"]
charge = bunch.charge()
momentum = bunch.getSyncParticle().momentum()
kq = self.getParam("dB/dr")/(3.335640952*momentum)
#---- The sign of dB/dr will be delivered to tracking module
#---- functions as kq.
#---- The charge sign will be accounted for inside tracking module
#---- functions.
kq = self.getParam("dB/dr")/bunch.B_Rho()
nParts = self.getnParts()
index = self.getActivePartIndex()
length = self.getLength(index)
Expand Down Expand Up @@ -646,8 +656,7 @@ def track(self, paramsDict):
TPB.bend2(bunch, length)
TPB.bend1(bunch, length, theta/2.0)
return



class DCorrectorH(LinacMagnetNode):
"""
The Horizontal Dipole Corrector.
Expand Down Expand Up @@ -684,11 +693,11 @@ def track(self, paramsDict):
length = self.getParam("effLength")/nParts
field = self.getParam("B")
bunch = paramsDict["bunch"]
charge = bunch.charge()
syncPart = bunch.getSyncParticle()
momentum = syncPart.momentum()
# dp/p = Q*c*B*L/p p in GeV/c c = 2.99792*10^8/10^9
# Q is used inside kick-method
kick = field*length*0.299792/momentum
kick = -field*charge*length*0.299792/momentum
self.tracking_module.kick(bunch,kick,0.,0.)

class DCorrectorV(LinacMagnetNode):
Expand Down Expand Up @@ -727,11 +736,11 @@ def track(self, paramsDict):
length = self.getParam("effLength")/nParts
field = self.getParam("B")
bunch = paramsDict["bunch"]
charge = bunch.charge()
syncPart = bunch.getSyncParticle()
momentum = syncPart.momentum()
# dp/p = Q*c*B*L/p p in GeV/c, c = 2.99792*10^8/10^9
# Q is used inside kick-method
kick = field*length*0.299792/momentum
kick = field*charge*length*0.299792/momentum
self.tracking_module.kick(bunch,0,kick,0.)

class ThickKick(LinacMagnetNode):
Expand Down Expand Up @@ -785,7 +794,8 @@ def track(self, paramsDict):
"""
The Thick Kick class implementation of the AccNode class track(probe) method.
"""
bunch = paramsDict["bunch"]
bunch = paramsDict["bunch"]
charge = bunch.charge()
momentum = bunch.getSyncParticle().momentum()
Bx = self.getParam("Bx")
By = self.getParam("By")
Expand All @@ -795,9 +805,8 @@ def track(self, paramsDict):
#print "debug name =",self.getName()," Bx=",Bx," By=",By," L=",self.getLength(index)," index=",index
#==========================================
# dp/p = Q*c*B*L/p p in GeV/c, c = 2.99792*10^8/10^9
# Q is used inside kick-method
kickY = Bx*length*0.299792/momentum
kickX = By*length*0.299792/momentum
kickY = +Bx*charge*length*0.299792/momentum
kickX = -By*charge*length*0.299792/momentum
self.tracking_module.drift(bunch, length/2.0)
self.tracking_module.kick(bunch,kickX,kickY,0.)
self.tracking_module.drift(bunch, length/2.0)
Expand Down
37 changes: 18 additions & 19 deletions py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -283,7 +283,8 @@ def track(self, paramsDict):
nParts = self.getnParts()
index = self.getActivePartIndex()
part_length = self.getLength(index)
bunch = paramsDict["bunch"]
bunch = paramsDict["bunch"]
charge = bunch.charge()
syncPart = bunch.getSyncParticle()
eKin_in = syncPart.kinEnergy()
momentum = syncPart.momentum()
Expand All @@ -308,20 +309,19 @@ def track(self, paramsDict):
Ep = self.getEzFiledInternal(zp,rfCavity,E0L,rf_ampl)
#------- track through a quad
G = self.getTotalField((zm+z0)/2)
GP = 0.
if(self.useLongField == True): GP = self.getTotalFieldDerivative((zm+z0)/2)
dB_dz = 0.
if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative((zm+z0)/2)
if(abs(G) != 0.):
kq = G/(3.335640952*momentum)
kq = G/bunch.B_Rho()
#------- track through a quad
step = part_length/2
self.tracking_module.quad1(bunch,step/4.0, kq)
self.tracking_module.quad2(bunch,step/2.0)
self.tracking_module.quad1(bunch,step/2.0, kq)
self.tracking_module.quad2(bunch,step/2.0)
self.tracking_module.quad1(bunch,step/4.0, kq)
if(abs(GP) != 0.):
kqP = GP/(3.335640952*momentum)
self.tracking_module.quad3(bunch,step, kqP)
if(abs(dB_dz) != 0.):
self.tracking_module.quad3(bunch,step,dB_dz)
else:
self.tracking_module.drift(bunch,part_length/2)
self.part_pos += part_length/2
Expand All @@ -340,19 +340,18 @@ def track(self, paramsDict):
self.cppGapModel.trackBunch(bunch,part_length/2,Em,E0,Ep,frequency,phase+delta_phase+modePhase)
#------- track through a quad
G = self.getTotalField((z0+zp)/2)
GP = 0.
if(self.useLongField == True): GP = self.getTotalFieldDerivative((z0+zp)/2)
dB_dz = 0.
if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative((z0+zp)/2)
if(abs(G) != 0.):
kq = G/(3.335640952*momentum)
kq = G/bunch.B_Rho()
step = part_length/2
self.tracking_module.quad1(bunch,step/4.0, kq)
self.tracking_module.quad2(bunch,step/2.0)
self.tracking_module.quad1(bunch,step/2.0, kq)
self.tracking_module.quad2(bunch,step/2.0)
self.tracking_module.quad1(bunch,step/4.0, kq)
if(abs(GP) != 0.):
kqP = GP/(3.335640952*momentum)
self.tracking_module.quad3(bunch,step, kqP)
if(abs(dB_dz) != 0.):
self.tracking_module.quad3(bunch,step,dB_dz)
else:
self.tracking_module.drift(bunch,part_length/2)
#---- advance the particle position
Expand Down Expand Up @@ -615,15 +614,16 @@ def track(self, paramsDict):
length = self.getLength(index)
if(index == 0): self.z_value = - self.getLength()/2
bunch = paramsDict["bunch"]
charge = bunch.charge()
momentum = bunch.getSyncParticle().momentum()
n_steps = int(length/self.z_step)+1
z_step = length/n_steps
for z_ind in range(n_steps):
z = self.z_value + z_step*(z_ind+0.5)
G = self.getTotalField(z)
GP = 0.
if(self.useLongField == True): GP = self.getTotalFieldDerivative(z)
kq = G/(3.335640952*momentum)
dB_dz = 0.
if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative(z)
kq = G/bunch.B_Rho()
if(abs(kq) == 0.):
self.tracking_module.drift(bunch,z_step)
continue
Expand All @@ -633,9 +633,8 @@ def track(self, paramsDict):
self.tracking_module.quad1(bunch,z_step/2.0, kq)
self.tracking_module.quad2(bunch,z_step/2.0)
self.tracking_module.quad1(bunch,z_step/4.0, kq)
if(abs(GP) != 0.):
kqP = GP/(3.335640952*momentum)
self.tracking_module.quad3(bunch,z_step, kqP)
if(abs(dB_dz) != 0.):
self.tracking_module.quad3(bunch,z_step, dB_dz)
self.z_value += length

def getTotalField(self,z_from_center):
Expand Down
52 changes: 22 additions & 30 deletions src/linac/tracking/linac_tracking.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ namespace linac_tracking
// bunch = reference to the macro-particle bunch
// length = length of transport
// kq = quadrupole field strength [m^(-2)]
// kq = G/(B*rho) B*rho = 3.33564*P0 [T*m] where P0 in GeV/c
// kq = G/(B*rho) B*rho = 3.33564*P0/charge [T*m] where P0 in GeV/c
// G = dB/dr [T/m] quad gradient
//
// RETURNS
Expand All @@ -125,11 +125,7 @@ namespace linac_tracking

void linac_quad1(Bunch* bunch, double length, double kq, int useCharge)
{
double charge = +1.0;
if(useCharge == 1) charge = bunch->getCharge();

double kqc = kq * charge;
if(kqc == 0.)
if(kq == 0. || bunch->getCharge() == 0.)
{
linac_drift(bunch,length);
return;
Expand All @@ -154,8 +150,8 @@ namespace linac_tracking
syncPart->setTime(syncPart->getTime() + delta_t);
}

// ==== B*rho = 3.335640952*momentum [T*m] if momentum in GeV/c ===
double dB_dr = kqc*3.335640952*syncPart->getMomentum();
// ==== B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c ===
double dB_dr = kq*bunch->getB_Rho();

double mass = syncPart->getMass();
double Ekin_s = syncPart->getEnergy();
Expand All @@ -175,6 +171,8 @@ namespace linac_tracking
double** arr = bunch->coordArr();
int nParts = bunch->getSize();

double kqc = 0.;

for(int i = 0; i < nParts; i++)
{

Expand All @@ -183,7 +181,7 @@ namespace linac_tracking
Etotal = Ekin + mass;
p2 = Ekin*(Ekin+2.0*mass);
p = sqrt(p2);
kqc = dB_dr/(3.335640952*p);
kqc = dB_dr/(3.335640952*p/bunch->getCharge());
beta_z = p/Etotal;

sqrt_kq = sqrt(fabs(kqc));
Expand Down Expand Up @@ -263,34 +261,34 @@ namespace linac_tracking
//
// DESCRIPTION
// Quadrupole element 3: longitudinal field component of the quad field
// Bz(z) = x*y*dG(z)/dz
// Bz(x,y,z) = x*y*dG(z)/dz
// This function performs the transverse kicks correction, so the length
// is just a parameter.
//
// PARAMETERS
// bunch = reference to the macro-particle bunch
// length = length of transport
// kq = quadrupole field strength derivative [m^(-2)]
// kq = (dG/dz)/(B*rho) B*rho = 3.33564*P0 [T*m] where P0 in GeV/c
// G = dB/dr [T/m] quad gradient
// 0.299792458 = (1/v_light) * 10^9
//
// RETURNS
// Nothing
//
///////////////////////////////////////////////////////////////////////////

void linac_quad3(Bunch* bunch, double length, double kq, int useCharge)
void linac_quad3(Bunch* bunch, double length, double dB_dz)
{
double charge = +1.0;
if(useCharge == 1) charge = bunch->getCharge();

double kqc = kq * charge;
if(kqc == 0.)
double charge = bunch->getCharge();

if(dB_dz == 0. || charge == 0.)
{
return;
}

double kick_coeff = kqc*length;

SyncPart* syncPart = bunch->getSyncPart();
//momentum in GeV/c
double momentum = syncPart->getMomentum();
double kick_coeff = 0.299792458*charge*dB_dz*length/momentum;

//coordinate array [part. index][x,xp,y,yp,z,dE]
double** arr = bunch->coordArr();
Expand Down Expand Up @@ -325,12 +323,6 @@ namespace linac_tracking

void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge)
{
double charge = +1.0;
if(useCharge == 1) charge = bunch->getCharge();

double kxc = kx * charge;
double kyc = ky * charge;

SyncPart* syncPart = bunch->getSyncPart();

double mass = syncPart->getMass();
Expand All @@ -345,7 +337,7 @@ namespace linac_tracking

//coordinate array [part. index][x,xp,y,yp,z,dE]
double** arr = bunch->coordArr();
if(kxc != 0.)
if(kx != 0.)
{
for(int i = 0; i < bunch->getSize(); i++)
{
Expand All @@ -354,10 +346,10 @@ namespace linac_tracking
p2 = Ekin*(Ekin+2.0*mass);
p_z = sqrt(p2 - (arr[i][1]*arr[i][1] + arr[i][3]*arr[i][3])*p2_s);
coeff = p_s/p_z;
arr[i][1] += kxc*coeff;
arr[i][1] += kx*coeff;
}
}
if(kyc != 0.)
if(ky != 0.)
{
for(int i = 0; i < bunch->getSize(); i++)
{
Expand All @@ -366,7 +358,7 @@ namespace linac_tracking
p2 = Ekin*(Ekin+2.0*mass);
p_z = sqrt(p2 - (arr[i][1]*arr[i][1] + arr[i][3]*arr[i][3])*p2_s);
coeff = p_s/p_z;
arr[i][3] += kyc*coeff;
arr[i][3] += ky*coeff;
}
}
if(kE != 0.)
Expand Down
2 changes: 1 addition & 1 deletion src/linac/tracking/linac_tracking.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,7 +46,7 @@ namespace linac_tracking
This function performs the transverse kicks correction, so the length
is just a parameter.
*/
void linac_quad3(Bunch* bunch, double length, double kq, int useCharge);
void linac_quad3(Bunch* bunch, double length, double dB_dz);


/**
Expand Down
8 changes: 4 additions & 4 deletions src/linac/tracking/wrap_linac_tracking.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,15 +75,15 @@ extern "C"
{
PyObject* pyBunch;
double length;
double kq;
double dB_dz = 0.;
int useCharge = 1;
if(!PyArg_ParseTuple( args, "Odd:linac_quad3",
&pyBunch, &length,&kq))
&pyBunch, &length,&dB_dz))
{
error("linac tracking - linac_quad3(pyBunch,length,kq) - cannot parse arguments!");
error("linac tracking - linac_quad3(pyBunch,length,dB_dz) - cannot parse arguments!");
}
Bunch* cpp_bunch = (Bunch*) ((pyORBIT_Object *) pyBunch)->cpp_obj;
linac_tracking::linac_quad3(cpp_bunch, length, kq, useCharge);
linac_tracking::linac_quad3(cpp_bunch, length, dB_dz);
Py_INCREF(Py_None);
return Py_None;
}
Expand Down
15 changes: 14 additions & 1 deletion src/orbit/Bunch.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -771,7 +771,20 @@ void Bunch::compress()
double Bunch::getMass(){ return mass;}
double Bunch::getCharge(){ return charge;}
double Bunch::getClassicalRadius(){ return classicalRadius;}
double Bunch::getMacroSize(){ return macroSizeForAll;}
double Bunch::getMacroSize(){ return macroSizeForAll;}

/**
Return the B*Rho parameter of the particle
B*Rho = momentum/charge - [T*m]
or B*Rho = 3.335640952*momentum[GeV/c]/charge[electron charges]
*/
double Bunch::getB_Rho(){
double B_Rho = 1.0e+36;
if(charge != 0. && syncPart != NULL){
B_Rho = 3.335640952*syncPart->getMomentum()/charge;
}
return B_Rho;
}

double Bunch::setMass(double val){
mass = val;
Expand Down
Loading

0 comments on commit 437306a

Please sign in to comment.