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:
Double_t dxdt(double t,double vx, double vy, Double_t vz){
Double_t Energy = GetEnergy(vx, vy, vz);
Double_t st = StoppingPower(Energy);
Double_t Ex,Ey,Ez,f1; //Electric field in V/m.
Double_t Bx,By,Bz; // magnetic field in Tesla
Double_t rr,az,po;
Double_t q = 1.6022*TMath::Power(10,-19); //charge of the particle(proton) in (C)
Double_t m = 1.6726*TMath::Power(10,-27); // mass of the particle in kg
Double_t B= -3.0; // Applied magnetic field (T).
Double_t E=TMath::Cos((q*B)/m) * 500; // Applied Electric field(V/m).
Bx = 0;
By = 0; // magnetic field in x and y direction in Tesla.
Bz = B; // magnetic field in the z direction in Tesla.
Ex = 0;
Ey = 0; // Electric field in the x and direction in V/m.
Ez = -E; // Electric field in the z direction.
rr = TMath::Sqrt(TMath::Power(vx,2)+TMath::Power(vy,2)+TMath::Power(vz,2));
az = TMath::ATan2(vy,vx);
po = TMath::ACos(vz/rr);
f1 = q/m * (Ex + vy*Bz-vz*By) - st*TMath::Sin(po)*TMath::Cos(az); //-s*TMath::Sin(po)*TMath::Cos(az) ; //dxdt with energyloss compensation.
double bro = Bz * rr / TMath::Sin(po)/ 1000.0;
//cout<<Energy<<endl;
//cout<<bro<<endl;
//std::cout<<Ex <<" "<< vy<<" "<<vz << " " <<By<< " "<< f1<<std::endl;
return f1;
}
double dydt(double t,double vx, double vy, Double_t vz){
Double_t Energy = GetEnergy(vx, vy, vz);
Double_t st = StoppingPower(Energy);
Double_t Ex,Ey,Ez,f2; //Electric field in V/m.
Double_t Bx,By,Bz; // magnetic field in Tesla
Double_t q = 1.6022*TMath::Power(10,-19); //charge of the particle in C
Double_t B= -3.0; // Applied magnetic field.
Double_t rr,az,po;
Double_t m = 1.6726*TMath::Power(10,-27); // mass of the particle in kg
Double_t E = TMath::Cos((q*B)/m) * 500 ;
Bx = 0; // magnetic field in x and y direction in Tesla.
By = 0;
Bz = B; // magnetic field in the z direction in Tesla.
Ex = 0; // Electric field in the x and direction in V/m.
Ey = 0;
Ez = -E; // Electric field in the z direction.
rr = TMath::Sqrt(TMath::Power(vx,2)+TMath::Power(vy,2)+TMath::Power(vz,2));
az = TMath::ATan(vy/vx);
po = TMath::ACos(vz/rr);
f2 = q/m * (Ey + vz*Bx - vx*Bz) - st*TMath::Sin(po)*TMath::Sin(az);
double bro = Bz * rr / TMath::Sin(po)/ 1000.0;
//cout<<bro<<endl;
//cout<<Energy<<endl;
//std::cout<<f2<<std::endl;
return f2;
}
double dzdt(double t,double vx, double vy,Double_t vz){
Double_t Energy = GetEnergy(vx, vy, vz);
Double_t st = StoppingPower(Energy);
Double_t Ex,Ey,Ez,f3; //Electric field in V/m.
Double_t Bx,By,Bz; // magnetic field in Tesla
Double_t q = 1.6022*pow(10,-19); //charge of the particle in eV
Double_t B= -3.0; // Applied magnetic field.
Double_t rr,az,po;
Double_t m = 1.6726*TMath::Power(10,-27); // mass of the particle in kg
Double_t E = TMath::Cos((q*B)/m) * 500 ;
Bx = 0; // magnetic field in x and y direction in Tesla.
By = 0;
Bz = B; // magnetic field in the z direction in Tesla.
Ex = 0; // Electric field in the x and direction in V/m.
Ey = 0;
Ez = -E; // Electric field in the z direction.
rr = TMath::Sqrt(TMath::Power(vx,2)+TMath::Power(vy,2)+TMath::Power(vz,2));
az = TMath::ATan(vy/vx);
//std::cout<<az<<endl;
po = TMath::ACos(vz/rr);
f3 = q/m * (Ez + vx*By - vy*Bx)- st*TMath::Cos(po);
double bro = Bz * rr / TMath::Sin(po)/ 1000.0;
//cout<<bro<<endl;
//cout<<Energy<<endl<<endl;
//std::cout <<r<< " " << TMath::Power(vy,2) + TMath::Power(vx,2) + TMath::Power(vz,2)<< " " << TMath::Power(vz,2) << std::endl;
return f3;
}
// Functions for positions.
Double_t fx(Double_t t, Double_t vx){
Double_t f1x = vx;
return f1x;
}
Double_t fy(Double_t t,Double_t vy){
Double_t f2y = vy;
return f2y;
}
Double_t fz(Double_t t,Double_t vz){
Double_t f3z = vz;
return f3z;
}
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::cout<<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
}
}
