harryharns
- 6
- 0
- TL;DR
- I have a simulated data of charged particles in a magnetic field. I have selected clusters, each cluster contains a set of points(x,z) and I want to perform RK4 between the first and second clusters and fill the positions in a histogram.
I have a simulated data of charged particles in a magnetic field. I have selected clusters, each cluster contains a set of points(x,z) and I want to perform RK4 between the first and second clusters and fill the positions in a histogram.
I have selected the clusters with the initial conditions(positions and velocity). but the thing is I have performed the rk4 but it is not behaving like I want and I am not able to perform the rk4 between the two clusters. Here is my code:
<Moderator's note: Please use code tags when posting code>
I have selected the clusters with the initial conditions(positions and velocity). but the thing is I have performed the rk4 but it is not behaving like I want and I am not able to perform the rk4 between the two clusters. Here is my code:
C++:
firstCluster = hitClusterArray->at(0);
secondCluster = hitClusterArray->at(1);
firstPosition = firstCluster.GetPosition();
secondPosition = secondCluster.GetPosition();
Double_t iniPx, iniPy, iniPz;
Double_t iniVx1, iniVy1, iniVz1;
Double_t secPosX,secPosY,secPosZ;
Double_t iniPosX,iniPosY,iniPosZ;
secPosX = secondPosition.X();
secPosY = secondPosition.Y();
secPosZ = secondPosition.Z();
iniPosX = firstPosition.X();
iniPosY = firstPosition.Y();
iniPosZ = firstPosition.Z();
phi = TMath::ATan2(secPosY- iniPosY, secPosX - iniPosX);
Double_t phiDeg;
if (phi < 0){
phiDeg = 360+ phi*TMath::RadToDeg();
} else {
phiDeg = phi*TMath::RadToDeg();
}
iniPx = p * TMath::Cos(phiDeg*TMath::DegToRad()) * TMath::Sin(theta); // in MeV/c
iniPy = p * TMath::Sin(phiDeg*TMath::DegToRad()) * TMath::Sin(theta); // in MeV/c
iniPz = p * TMath::Cos(theta); // in MeV/c iniVx1 = (iniPx * 5.344286e-22)/m; //meter per second [m/s]
iniVy1 = (iniPy * 5.344286e-22)/m;
iniVz1 = (iniPz * 5.344286e-22)/m;
int clusterCount = 0;
for (auto iCluster = 0; iCluster < hitClusterArray->size(); ++iCluster) {
auto cluster = hitClusterArray->at(iCluster);
auto pos = cluster.GetPosition();
Double_t clusterPosX = pos.X();
Double_t clusterPosY = pos.Y();
Double_t clusterPosZ = pos.Z();
//std::count<<iniPosX << " " << iniPosY << " "<< iniPosZ <<endl << endl;
for (pp = 0; pp< 100; pp++){
l1x[pp] = fx(t1,iniPx1);
l1y[pp] = fy(t1,iniPy1);
l1z[pp] = fz(t1,iniPz1);
l1vx[pp] = dxdt(t1,iniPx1,iniPy1,iniPz1);
l1vy[pp] = dydt(t1,iniPx1,iniPy1,iniPz1);
l1vz[pp] = dzdt(t1,iniPx1,iniPy1,iniPz1);
l2x[pp] = fx(t1+h*0.5,iniPx1+0.5*h*l1x[pp]);
l2y[pp] = fy(t1+h*0.5,iniPy1+0.5*h*l1y[pp]);
l2z[pp] = fz(t1+h*0.5,iniPz1+0.5*h*l1z[pp]);
l2vx[pp] =dxdt(t1+h*0.5,iniPx1+l1vx[pp]*0.5*h,iniPy1+l1vy[pp]*h*0.5,iniPz1+l1vz[pp]*h*0.5);
l2vy[pp] =dydt(t1+h*0.5,iniPx1+l1vx[pp]*h*0.5,iniPy1+l1vy[pp]*h*0.5,iniPz1+l1vz[pp]*h*0.5);
l2vz[pp] =dzdt(t1+h*0.5,iniPx1+l1vx[pp]*h*0.5,iniPy1+l1vy[pp]*h*0.5,iniPz1+l1vz[pp]*h*0.5);
l3x[pp] = fx(t1+h*0.5,iniPx1+0.5*h*l2x[pp]);
l3y[pp] = fy(t1+h*0.5,iniPy1+0.5*h*l2y[pp]);
l3z[pp] = fz(t1+h*0.5,iniPz1+0.5*h*l2z[pp]);
l3vx[pp] =dxdt(t1+h*0.5,iniPx1+l2vx[pp]*0.5*h,iniPy1+l2vy[pp]*h*0.5,iniPz1+l2vz[pp]*h*0.5);
l3vy[pp] =dydt(t1+h*0.5,iniPx1+l2vx[pp]*h*0.5,iniPy1+l2vy[pp]*h*0.5,iniPz1+l2vz[pp]*h*0.5);
l3vz[pp] =dzdt(t1+h*0.5,iniPx1+l2vx[pp]*h*0.5,iniPy1+l2vy[pp]*h*0.5,iniPz1+l2vz[pp]*h*0.5);
l4x[pp] = fx(t1+h,iniPx1+l3x[pp]*h);
l4y[pp] = fy(t1+h,iniPy1+l3y[pp]*h);
l4z[pp] = fz(t1+h,iniPz1+l3z[pp]*h);
l4vx[pp] = dxdt(t1+h,iniPx1+l3vx[pp]*h,iniPy1+l3vy[pp]*h,iniPz1+l3vz[pp]*h);
l4vy[pp] = dydt(t1+h,iniPx1+l3vx[pp]*h,iniPy1+l3vy[pp]*h,iniPz1+l3vz[pp]*h);
l4vz[pp] = dzdt(t1+h,iniPx1+l3vx[pp]*h,iniPy1+l3vy[pp]*h,iniPz1+l3vz[pp]*h); iniPx1 = iniPx1 + h/6 *( l1vx[pp] + 2*l2vx[pp] + 2*l3vx[pp] + l4vx[pp]);
iniPy1 = iniPy1 + h/6 *(l1vy[pp] + 2*l2vy[pp] + 2*l3vy[pp] + l4vy[pp]);
iniPz1 = iniPz1 + h/6 *(l1vz[pp] + 2*l2vz[pp] + 2*l3vz[pp] + l4vz[pp]);
iniPosX = iniPosX + h/6 *( l1x[pp] + 2*l2x[pp] + 2*l3x[pp] + l4x[pp]);
iniPosY = iniPosY + h/6 *( l1y[pp] + 2*l2y[pp] + 2*l3y[pp] + l4y[pp]);
iniPosZ = iniPosZ + h/6 *( l1z[pp] + 2*l2z[pp] + 2*l3z[pp] + l4z[pp]);
kx_vs_ky->Fill(iniPosX, iniPosY);
}
//Update the initial conditions for the next cluster
iniPosX = clusterPosX;
iniPosY = clusterPosY;
iniPosZ = clusterPosZ;
++clusterCount;
if (clusterCount > 2) {
break; // stop iterating over clusters after the second one
}
}
Last edited by a moderator:
