diff --git a/src/detector_sim.cpp b/src/detector_sim.cpp
index 385b9d2279b78ff38e7e5ae588d3c598b93b7a07..4e84b8512000b54c6fc2932be09b8704731a510d 100644
--- a/src/detector_sim.cpp
+++ b/src/detector_sim.cpp
@@ -190,48 +190,49 @@ int main(int argc, char **argv){
//simple gaussian multiple scattering
if(detector_geo->GetMultipleScattStatus()){
-
//get the angle by wich the track is reflected (the full angle between track before and after scattering)
- double theta_ms = sqrt(2.) * detector_geo->GetMultipleScattering()->GetThetaGauss(curr_layer->GetXoverXZero(),
- curr_particle->GetBeta(),
+ double theta_ms = sqrt(2.) * detector_geo->GetMultipleScattering()->GetThetaGauss(curr_layer->GetXoverXZero(),
+ curr_particle->GetBeta(),
curr_particle->GetMomentumMod(),
curr_particle->GetCharge());
//Express the current momentum in the coordinate system where z is along the momentum
double sinalpha = temp_state->GetMomentum()[0]
- / sqrt(pow(temp_state->GetMomentum()[0], 2.) + pow(temp_state->GetMomentum()[2], 2.));
-
- double sinbeta = temp_state->GetMomentum()[1] / sqrt(temp_state->GetMomentum().Norm2Sqr());
+ /sqrt(pow(temp_state->GetMomentum()[0], 2.)+ pow(temp_state->GetMomentum()[2], 2.));
+ double sinbeta = temp_state->GetMomentum()[1]/sqrt(temp_state->GetMomentum().Norm2Sqr());
double cosalpha = sqrt(1. - pow(sinalpha, 2.));
- double cosbeta = sqrt(1.- pow(sinbeta, 2.));
+ double cosbeta = sqrt(1. - pow(sinbeta, 2.));
double rot1[9] = {cosalpha, 0, -sinalpha,
0 , 1, 0,
sinalpha, 0, cosalpha};
-
- TMatrixD Rot1 = TMatrixD(3, 3, rot1);
+
+ TMatrixD Rot1 = TMatrixD(3,3, rot1);
TMatrixD Rot1Inv(TMatrixD::kInverted,Rot1);
double rot2[9] = {1, 0, 0,
0, cosbeta, -sinbeta,
0, sinbeta, cosbeta};
-
- TMatrixD Rot2 = TMatrixD(3, 3, rot2);
+
+ TMatrixD Rot2 = TMatrixD(3,3, rot2);
TMatrixD Rot2Inv(TMatrixD::kInverted,Rot2);
- TVectorD p_rot = Rot2 * Rot1 * temp_state->GetMomentum();
+ //actually not necessary to calculate the momentum vector in the rotated
+ //coordiante system, since it is anyway along the z axis, but let's do it
+ //anyway
+ TVectorD p_rot = Rot2*Rot1*temp_state->GetMomentum();
//Rotate the momentum vector due to multiple scattering
//"directional" angle phi is uniformly distributed
double phi = random_ms->Uniform(-TMath::Pi(), TMath::Pi());
TVectorD p_rot_ms = TVectorD(3);
-
- p_rot_ms[0] = sin(theta_ms) * sin(phi) * sqrt(p_rot.Norm2Sqr());
- p_rot_ms[1] = sin(theta_ms) * cos(phi) * sqrt(p_rot.Norm2Sqr());
- p_rot_ms[2] = cos(theta_ms) * sqrt(p_rot.Norm2Sqr());
+
+ p_rot_ms[0]=sin(theta_ms)*sin(phi)*sqrt(p_rot.Norm2Sqr());
+ p_rot_ms[1]=sin(theta_ms)*cos(phi)*sqrt(p_rot.Norm2Sqr());
+ p_rot_ms[2]=cos(theta_ms)*sqrt(p_rot.Norm2Sqr());
//Rotate this momentum back to the lab frame
- TVectorD p_ms = Rot1Inv * Rot2Inv * p_rot_ms;
+ TVectorD p_ms = Rot1Inv*Rot2Inv*p_rot_ms;
temp_state->UpdateTxTy(p_ms);