From e012e8ebf1da22f4034fec9a5c026e1c567df4ea Mon Sep 17 00:00:00 2001 From: shishlo Date: Tue, 5 Mar 2024 14:30:09 -0500 Subject: [PATCH 1/3] Changes were made to fix usage of k1(n) strength of magnet multipoles in the TEAPOT and linac tracking functions and methods. k1 for quad is dB/dr/(p0/charge), but the implementation was correct only for protons and and H- ions. The fix allows to use TEAPOT tracking for all heavy ions. --- py/orbit/py_linac/lattice/LinacAccNodes.py | 38 +++++---- .../lattice/LinacFieldOverlappingNodes.py | 37 +++++---- src/linac/tracking/linac_tracking.cc | 51 ++++++------ src/linac/tracking/linac_tracking.hh | 2 +- src/linac/tracking/wrap_linac_tracking.cc | 8 +- src/teapot/teapotbase.cc | 77 +++++++++++++++---- src/teapot/wrap_teapotbase.cc | 2 +- 7 files changed, 134 insertions(+), 81 deletions(-) diff --git a/py/orbit/py_linac/lattice/LinacAccNodes.py b/py/orbit/py_linac/lattice/LinacAccNodes.py index d1119d8f..400fbb59 100755 --- a/py/orbit/py_linac/lattice/LinacAccNodes.py +++ b/py/orbit/py_linac/lattice/LinacAccNodes.py @@ -301,9 +301,11 @@ def fringeIN(node,paramsDict): 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")/(3.335640952*momentum/abs(charge)) poleArr = node.getParam("poles") klArr = node.getParam("kls") skewArr = node.getParam("skews") @@ -322,8 +324,10 @@ def fringeOUT(node,paramsDict): 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")/(3.335640952*momentum/abs(charge)) poleArr = node.getParam("poles") klArr = node.getParam("kls") skewArr = node.getParam("skews") @@ -400,9 +404,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 without charge sign transformation as kq. + #---- The charge sign will be accounted for inside tracking module + #---- functions. + kq = self.getParam("dB/dr")/(3.335640952*momentum/abs(charge)) nParts = self.getnParts() index = self.getActivePartIndex() length = self.getLength(index) @@ -646,8 +655,7 @@ def track(self, paramsDict): TPB.bend2(bunch, length) TPB.bend1(bunch, length, theta/2.0) return - - + class DCorrectorH(LinacMagnetNode): """ The Horizontal Dipole Corrector. @@ -684,11 +692,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): @@ -727,11 +735,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): @@ -785,7 +793,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") @@ -795,9 +804,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) diff --git a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py index 49159b7f..bb51beae 100755 --- a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py +++ b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py @@ -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() @@ -308,10 +309,10 @@ 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/(3.335640952*momentum/abs(charge)) #------- track through a quad step = part_length/2 self.tracking_module.quad1(bunch,step/4.0, kq) @@ -319,9 +320,8 @@ def track(self, paramsDict): 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 @@ -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/(3.335640952*momentum/abs(charge)) 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 @@ -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/(3.335640952*momentum/abs(charge)) if(abs(kq) == 0.): self.tracking_module.drift(bunch,z_step) continue @@ -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): diff --git a/src/linac/tracking/linac_tracking.cc b/src/linac/tracking/linac_tracking.cc index 5de037d7..08e21972 100644 --- a/src/linac/tracking/linac_tracking.cc +++ b/src/linac/tracking/linac_tracking.cc @@ -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 @@ -127,9 +127,8 @@ namespace linac_tracking { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); - - double kqc = kq * charge; - if(kqc == 0.) + + if(kq == 0. || charge == 0.) { linac_drift(bunch,length); return; @@ -154,8 +153,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*3.335640952*syncPart->getMomentum()/charge; double mass = syncPart->getMass(); double Ekin_s = syncPart->getEnergy(); @@ -175,6 +174,8 @@ namespace linac_tracking double** arr = bunch->coordArr(); int nParts = bunch->getSize(); + double kqc = 0.; + for(int i = 0; i < nParts; i++) { @@ -183,7 +184,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/fabs(charge)); beta_z = p/Etotal; sqrt_kq = sqrt(fabs(kqc)); @@ -263,34 +264,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 = 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(); @@ -325,12 +326,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(); @@ -345,7 +340,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++) { @@ -354,10 +349,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++) { @@ -366,7 +361,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.) diff --git a/src/linac/tracking/linac_tracking.hh b/src/linac/tracking/linac_tracking.hh index aba2d24d..415225e9 100644 --- a/src/linac/tracking/linac_tracking.hh +++ b/src/linac/tracking/linac_tracking.hh @@ -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); /** diff --git a/src/linac/tracking/wrap_linac_tracking.cc b/src/linac/tracking/wrap_linac_tracking.cc index 2a17c93e..b6bb255d 100644 --- a/src/linac/tracking/wrap_linac_tracking.cc +++ b/src/linac/tracking/wrap_linac_tracking.cc @@ -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; } diff --git a/src/teapot/teapotbase.cc b/src/teapot/teapotbase.cc index cffa337a..76e6d25c 100644 --- a/src/teapot/teapotbase.cc +++ b/src/teapot/teapotbase.cc @@ -231,32 +231,28 @@ void wrapbunch(Bunch* bunch, double length) 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; - double kEc = kE * charge; + //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++) { - arr[i][1] += kxc; + arr[i][1] += kx; } } - if(kyc != 0.) + if(ky != 0.) { for(int i = 0; i < bunch->getSize(); i++) { - arr[i][3] += kyc; + arr[i][3] += ky; } } - if(kEc != 0.) + if(kE != 0.) { for(int i = 0; i < bunch->getSize(); i++) { - arr[i][5] += kEc; + arr[i][5] += kE; } } } @@ -284,6 +280,13 @@ void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + + if(charge == 0.){ + return; + } + + charge = fabs(charge)/charge; + double klc = kl * charge; std::complex z, zn; double kl1; @@ -337,6 +340,13 @@ void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + + if(charge == 0.){ + return; + } + + charge = fabs(charge)/charge; + double klc = kl * charge; std::complex z, zn; double kl1; @@ -393,6 +403,13 @@ void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + + if(charge == 0.){ + return; + } + + charge = fabs(charge)/charge; + double klc = kl * charge; std::complex rootm1 = std::complex(0.0, 1.0); @@ -503,6 +520,13 @@ void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + + if(charge == 0.){ + return; + } + + charge = fabs(charge)/charge; + double klc = kl * charge; std::complex rootm1 = std::complex(0.0, 1.0); @@ -612,12 +636,15 @@ void 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. || charge == 0.) { drift(bunch,length); return; } + + double kqc = kq * charge/fabs(charge); + double x_init, xp_init, y_init, yp_init; double sqrt_kq, kqlength; double cx, sx, cy, sy; @@ -775,6 +802,13 @@ void quadfringeIN(Bunch* bunch, double kq, int useCharge) { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + + if(charge == 0.){ + return; + } + + charge = fabs(charge)/charge; + double kqc = kq * charge; double KNL, x_init, xp_init, y_init, yp_init, detM; @@ -842,6 +876,13 @@ void quadfringeOUT(Bunch* bunch, double kq, int useCharge) { double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + + if(charge == 0.){ + return; + } + + charge = fabs(charge)/charge; + double kqc = kq * charge; double KNL, x_init, xp_init, y_init, yp_init, detM; @@ -1359,6 +1400,7 @@ void soln(Bunch* bunch, double length, double B, int useCharge) double charge = +1.0; if(useCharge == 1) charge = bunch->getCharge(); + double Bc = B * charge; double KNL, phase, cs, sn; double cu, cpu, u_init, pu_init, u, pu, phifac; @@ -1526,6 +1568,8 @@ void RingRF(Bunch* bunch, double ring_length, int harmonic_numb, } SyncPart* syncPart = bunch->getSyncPart(); + + double p_synch_in = syncPart->getMomentum(); if(phase_s != 0.) { @@ -1533,6 +1577,10 @@ void RingRF(Bunch* bunch, double ring_length, int harmonic_numb, kin_e += coeff * voltage * sin(phase_s); syncPart->setMomentum(syncPart->energyToMomentum(kin_e)); } + + double p_synch_out = syncPart->getMomentum(); + + double xp_yp_coeff = p_synch_in/p_synch_out; //coordinate array [part. index][x,xp,y,yp,z,dE] double** arr = bunch->coordArr(); @@ -1541,6 +1589,9 @@ void RingRF(Bunch* bunch, double ring_length, int harmonic_numb, { deltaV = voltage * ( sin(harmonic_numb*Factor*arr[i][4] + phase_s)); arr[i][5] += coeff * deltaV; + + arr[i][1] *= xp_yp_coeff; + arr[i][3] *= xp_yp_coeff; } } diff --git a/src/teapot/wrap_teapotbase.cc b/src/teapot/wrap_teapotbase.cc index fe8a29ca..b353130d 100644 --- a/src/teapot/wrap_teapotbase.cc +++ b/src/teapot/wrap_teapotbase.cc @@ -488,7 +488,7 @@ extern "C" {"rotatexy", wrap_rotatexy, METH_VARARGS, "Rotates bunch around z axis "}, {"drift", wrap_drift, METH_VARARGS, "Tracking a bunch through a drift "}, {"drifti", wrap_drifti, METH_VARARGS, "Drifts one macroparticle in the bunch"}, - {"wrapbunch", wrap_wrapbunch, METH_VARARGS, "Tracking a bunch through a wrapbunch routine"}, + {"wrapbunch", wrap_wrapbunch, METH_VARARGS, "Tracking a bunch through a wrapbunch routine"}, {"multp", wrap_multp, METH_VARARGS, "Tracking a bunch through a multipole "}, {"multpfringeIN", wrap_multpfringeIN, METH_VARARGS, "Tracking a bunch through an IN edge of a multipole "}, {"multpfringeOUT", wrap_multpfringeOUT, METH_VARARGS, "Tracking a bunch through an OUT edge of a multipole"}, From 9f616367b3a2ed532c7989bd0b3585b5e5b0a60a Mon Sep 17 00:00:00 2001 From: shishlo Date: Thu, 14 Mar 2024 17:41:27 -0400 Subject: [PATCH 2/3] Changes were made to fix usage of k1(n) strength of magnet multipoles in the TEAPOT and linac tracking functions and methods. k1 for quad is dB/dr/(p0/charge), but the implementation was correct only for protons and and H- ions. The fix allows to use TEAPOT tracking for all heavy ions. Now the Bunch class has method returning B_Rho = B*Rho parameter in [T*m]. --- py/orbit/py_linac/lattice/LinacAccNodes.py | 11 +-- .../lattice/LinacFieldOverlappingNodes.py | 6 +- src/linac/tracking/linac_tracking.cc | 11 +-- src/orbit/Bunch.cc | 15 +++- src/orbit/Bunch.hh | 7 ++ src/orbit/wrap_bunch.cc | 8 ++ src/teapot/teapotbase.cc | 85 ++++++------------- 7 files changed, 70 insertions(+), 73 deletions(-) diff --git a/py/orbit/py_linac/lattice/LinacAccNodes.py b/py/orbit/py_linac/lattice/LinacAccNodes.py index 400fbb59..edeadb8f 100755 --- a/py/orbit/py_linac/lattice/LinacAccNodes.py +++ b/py/orbit/py_linac/lattice/LinacAccNodes.py @@ -297,7 +297,7 @@ 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 @@ -305,7 +305,7 @@ def fringeIN(node,paramsDict): charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() #---- The charge sign will be accounted for inside tracking module functions. - kq = node.getParam("dB/dr")/(3.335640952*momentum/abs(charge)) + kq = node.getParam("dB/dr")/bunch.B_Rho() poleArr = node.getParam("poles") klArr = node.getParam("kls") skewArr = node.getParam("skews") @@ -320,6 +320,7 @@ 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 @@ -327,7 +328,7 @@ def fringeOUT(node,paramsDict): charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() #---- The charge sign will be accounted for inside tracking module functions - kq = node.getParam("dB/dr")/(3.335640952*momentum/abs(charge)) + kq = node.getParam("dB/dr")/bunch.B_Rho() poleArr = node.getParam("poles") klArr = node.getParam("kls") skewArr = node.getParam("skews") @@ -408,10 +409,10 @@ def track(self, paramsDict): charge = bunch.charge() momentum = bunch.getSyncParticle().momentum() #---- The sign of dB/dr will be delivered to tracking module - #---- functions without charge sign transformation as kq. + #---- functions as kq. #---- The charge sign will be accounted for inside tracking module #---- functions. - kq = self.getParam("dB/dr")/(3.335640952*momentum/abs(charge)) + kq = self.getParam("dB/dr")/bunch.B_Rho() nParts = self.getnParts() index = self.getActivePartIndex() length = self.getLength(index) diff --git a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py index bb51beae..63435c9f 100755 --- a/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py +++ b/py/orbit/py_linac/lattice/LinacFieldOverlappingNodes.py @@ -312,7 +312,7 @@ def track(self, paramsDict): dB_dz = 0. if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative((zm+z0)/2) if(abs(G) != 0.): - kq = G/(3.335640952*momentum/abs(charge)) + kq = G/bunch.B_Rho() #------- track through a quad step = part_length/2 self.tracking_module.quad1(bunch,step/4.0, kq) @@ -343,7 +343,7 @@ def track(self, paramsDict): dB_dz = 0. if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative((z0+zp)/2) if(abs(G) != 0.): - kq = G/(3.335640952*momentum/abs(charge)) + 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) @@ -623,7 +623,7 @@ def track(self, paramsDict): G = self.getTotalField(z) dB_dz = 0. if(self.useLongField == True): dB_dz = self.getTotalFieldDerivative(z) - kq = G/(3.335640952*momentum/abs(charge)) + kq = G/bunch.B_Rho() if(abs(kq) == 0.): self.tracking_module.drift(bunch,z_step) continue diff --git a/src/linac/tracking/linac_tracking.cc b/src/linac/tracking/linac_tracking.cc index 08e21972..3ed38749 100644 --- a/src/linac/tracking/linac_tracking.cc +++ b/src/linac/tracking/linac_tracking.cc @@ -125,10 +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(); - - if(kq == 0. || charge == 0.) + if(kq == 0. || bunch->getCharge() == 0.) { linac_drift(bunch,length); return; @@ -154,7 +151,7 @@ namespace linac_tracking } // ==== B*rho = 3.335640952*momentum/charge [T*m] if momentum in GeV/c === - double dB_dr = kq*3.335640952*syncPart->getMomentum()/charge; + double dB_dr = kq*bunch->getB_Rho(); double mass = syncPart->getMass(); double Ekin_s = syncPart->getEnergy(); @@ -184,7 +181,7 @@ namespace linac_tracking Etotal = Ekin + mass; p2 = Ekin*(Ekin+2.0*mass); p = sqrt(p2); - kqc = dB_dr/(3.335640952*p/fabs(charge)); + kqc = dB_dr/(3.335640952*p/bunch->getCharge()); beta_z = p/Etotal; sqrt_kq = sqrt(fabs(kqc)); @@ -281,7 +278,7 @@ namespace linac_tracking void linac_quad3(Bunch* bunch, double length, double dB_dz) { - double charge = charge = bunch->getCharge(); + double charge = bunch->getCharge(); if(dB_dz == 0. || charge == 0.) { diff --git a/src/orbit/Bunch.cc b/src/orbit/Bunch.cc index 6f635f12..f97ccef2 100644 --- a/src/orbit/Bunch.cc +++ b/src/orbit/Bunch.cc @@ -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; diff --git a/src/orbit/Bunch.hh b/src/orbit/Bunch.hh index 2db1bba3..2d49e802 100644 --- a/src/orbit/Bunch.hh +++ b/src/orbit/Bunch.hh @@ -111,6 +111,13 @@ public: double getClassicalRadius(); // m double getCharge(); // sign and value in abs(e-charge) only double getMacroSize(); + + /** + 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 getB_Rho(); double setMass(double mass); // GeV double setCharge(double chrg); // sign and value in abs(e-charge) only diff --git a/src/orbit/wrap_bunch.cc b/src/orbit/wrap_bunch.cc index 8392e59b..4e5d4f2d 100644 --- a/src/orbit/wrap_bunch.cc +++ b/src/orbit/wrap_bunch.cc @@ -504,6 +504,13 @@ namespace wrap_orbit_bunch{ double val = cpp_bunch->getClassicalRadius(); return Py_BuildValue("d",val); } + + //Returns B_Rho of the particle in [Tesla*meter]. Parameter is used in TEAPOT + static PyObject* Bunch_B_Rho(PyObject *self, PyObject *args){ + Bunch* cpp_bunch = (Bunch*) ((pyORBIT_Object *) self)->cpp_obj; + double val = cpp_bunch->getB_Rho(); + return Py_BuildValue("d",val); + } //Sets or returns charge of the macro-particle in e-charge // the action is depended on the number of arguments @@ -1189,6 +1196,7 @@ namespace wrap_orbit_bunch{ { "ringwrap", Bunch_ringwrap ,METH_VARARGS,"Perform the ring wrap. Usage: ringwrap(ring_length)"}, { "mass", Bunch_mass ,METH_VARARGS,"Set mass(value) or get mass() the mass of particle in MeV"}, { "classicalRadius", Bunch_classicalRadius ,METH_VARARGS,"Returns a classical radius of particle in [m]"}, + { "B_Rho", Bunch_B_Rho ,METH_VARARGS,"Returns B*Rho parameter of particle in [T*m]"}, { "charge", Bunch_charge ,METH_VARARGS,"Set charge(value) or get charge() the charge of particle in e-charge"}, { "macroSize", Bunch_macroSize ,METH_VARARGS,"Set macroSize(value) or get macroSize() the charge of particle in e-charge"}, { "initBunchAttr", Bunch_initBunchAttr ,METH_VARARGS,"Reads and initilizes bunch attributes from a bunch file"}, diff --git a/src/teapot/teapotbase.cc b/src/teapot/teapotbase.cc index 76e6d25c..fecd46bd 100644 --- a/src/teapot/teapotbase.cc +++ b/src/teapot/teapotbase.cc @@ -269,6 +269,7 @@ void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge) // i = particle index // pole = multipole number // kl = integrated strength of the kick [m^(-pole)] +// kl already has information about the charge of particle. // skew = 0 - normal, 1 - skew // // RETURNS @@ -278,16 +279,11 @@ void kick(Bunch* bunch, double kx, double ky, double kE, int useCharge) void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(charge == 0.){ + if(bunch->getCharge() == 0.){ return; } - charge = fabs(charge)/charge; - - double klc = kl * charge; + double klc = kl; std::complex z, zn; double kl1; @@ -329,6 +325,7 @@ void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) // pole = multipole number // pole = 0 for dipole, pole = 1 for quad, pole = 2 for sextupole, pole = 3 for octupole // kl = integrated strength of the kick [m^(-pole)] +// kl already has information about the charge of particle. // skew = 0 - normal, 1 - skew // // RETURNS @@ -338,16 +335,11 @@ void multpi(Bunch* bunch, int i, int pole, double kl, int skew, int useCharge) void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(charge == 0.){ + if(bunch->getCharge() == 0.){ return; - } - - charge = fabs(charge)/charge; + } - double klc = kl * charge; + double klc = kl; std::complex z, zn; double kl1; @@ -392,6 +384,7 @@ void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) // bunch = reference to the macro-particle bunch // pole = multipole number // kl = multipole strength +// kl already has information about the charge of particle. // skew = multipole skew // // RETURNS @@ -401,16 +394,11 @@ void multp(Bunch* bunch, int pole, double kl, int skew, int useCharge) void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(charge == 0.){ + if(bunch->getCharge() == 0.){ return; - } - - charge = fabs(charge)/charge; + } - double klc = kl * charge; + double klc = kl; std::complex rootm1 = std::complex(0.0, 1.0); SyncPart* syncPart = bunch->getSyncPart(); @@ -509,6 +497,7 @@ void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) // bunch = reference to the macro-particle bunch // pole = multipole number // kl = multipole strength +// kl already has information about the charge of particle. // skew = multipole skew // // RETURNS @@ -518,16 +507,11 @@ void multpfringeIN(Bunch* bunch, int pole, double kl, int skew, int useCharge) void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(charge == 0.){ + if(bunch->getCharge() == 0.){ return; - } + } - charge = fabs(charge)/charge; - - double klc = kl * charge; + double klc = kl; std::complex rootm1 = std::complex(0.0, 1.0); SyncPart* syncPart = bunch->getSyncPart(); @@ -626,6 +610,7 @@ void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) // bunch = reference to the macro-particle bunch // length = length of transport // kq = quadrupole field strength [m^(-2)] +// kq already has information about the charge of particle. // // RETURNS // Nothing @@ -634,16 +619,13 @@ void multpfringeOUT(Bunch* bunch, int pole, double kl, int skew, int useCharge) void quad1(Bunch* bunch, double length, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(kq == 0. || charge == 0.) + if(kq == 0. || bunch->getCharge() == 0.) { drift(bunch,length); return; } - double kqc = kq * charge/fabs(charge); + double kqc = kq; double x_init, xp_init, y_init, yp_init; double sqrt_kq, kqlength; @@ -792,6 +774,7 @@ void quad3(Bunch* bunch, double length, double kq, int useCharge) // PARAMETERS // bunch = reference to the macro-particle bunch // kq = strength of quad +// kq already has information about the charge of particle. // // RETURNS // Nothing @@ -800,16 +783,11 @@ void quad3(Bunch* bunch, double length, double kq, int useCharge) void quadfringeIN(Bunch* bunch, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(charge == 0.){ + if(bunch->getCharge() == 0.){ return; - } - - charge = fabs(charge)/charge; + } - double kqc = kq * charge; + double kqc = kq; double KNL, x_init, xp_init, y_init, yp_init, detM; SyncPart* syncPart = bunch->getSyncPart(); @@ -866,6 +844,7 @@ void quadfringeIN(Bunch* bunch, double kq, int useCharge) // PARAMETERS // bunch = reference to the macro-particle bunch // kq = strength of quad +// kq already has information about the charge of particle. // // RETURNS // Nothing @@ -874,16 +853,11 @@ void quadfringeIN(Bunch* bunch, double kq, int useCharge) void quadfringeOUT(Bunch* bunch, double kq, int useCharge) { - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - - if(charge == 0.){ + if(bunch->getCharge() == 0.){ return; - } + } - charge = fabs(charge)/charge; - - double kqc = kq * charge; + double kqc = kq; double KNL, x_init, xp_init, y_init, yp_init, detM; SyncPart* syncPart = bunch->getSyncPart(); @@ -1393,15 +1367,12 @@ void bendfringeOUT(Bunch* bunch, double rho) void soln(Bunch* bunch, double length, double B, int useCharge) { //if solenoid field in [T] is zero we have just a drift - if(abs(B) < 1.0e-100){ + if(abs(B) < 1.0e-100 || bunch->getCharge() == 0){ drift(bunch,length); return; } - - double charge = +1.0; - if(useCharge == 1) charge = bunch->getCharge(); - double Bc = B * charge; + double Bc = B * bunch->getCharge(); double KNL, phase, cs, sn; double cu, cpu, u_init, pu_init, u, pu, phifac; From 35e0ae5bfd2079887ffe837fc9a676e6709188a4 Mon Sep 17 00:00:00 2001 From: shishlo Date: Fri, 15 Mar 2024 13:53:18 -0400 Subject: [PATCH 3/3] Update teapotbase.cc Typo. The charge is a double. --- src/teapot/teapotbase.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/teapot/teapotbase.cc b/src/teapot/teapotbase.cc index fecd46bd..0c2b354a 100644 --- a/src/teapot/teapotbase.cc +++ b/src/teapot/teapotbase.cc @@ -1367,7 +1367,7 @@ void bendfringeOUT(Bunch* bunch, double rho) void soln(Bunch* bunch, double length, double B, int useCharge) { //if solenoid field in [T] is zero we have just a drift - if(abs(B) < 1.0e-100 || bunch->getCharge() == 0){ + if(abs(B) < 1.0e-100 || bunch->getCharge() == 0.){ drift(bunch,length); return; }