diff --git a/src/musrPrimaryGeneratorAction.cc b/src/musrPrimaryGeneratorAction.cc index a144e6a..2baf7d9 100644 --- a/src/musrPrimaryGeneratorAction.cc +++ b/src/musrPrimaryGeneratorAction.cc @@ -271,8 +271,10 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) G4double px, py, pz; px = p*sin(xangle); py = p*sin(yangle); - // works for beam along +-z only for now - pz = zDirection * std::sqrt(p*p - px*px - py*py); + pz = std::sqrt(p*p - px*px - py*py); + + + // Now rotate so that beam propagates along /gun/direction // Proper rotation of beam direction as follows: // hat{n} is beam direction unit vecrot and vec{p} is momentum vector @@ -280,7 +282,25 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) // where R is a rotation matrix around (hat{z} x hat{n}) - cross product // with an angle of acos(hat{z} dot hat{n}) - dot product // i.e. R is a rotation around (-n_y, n_x, 0) with angle acos(n_z). - // G4RotationMatrix* R = new G4RotationMatrix(G4ThreeVector(-nx,ny,0),acos(nz)*deg); + + if (xDirection == 0 & yDirection == 0) { + // Rotation does not work for beam direction along z. + pz = zDirection * pz; + // No change to the beam spot... + } else { + G4ThreeVector* PVec = new G4ThreeVector(px,py,pz); + PVec->rotate(std::acos(zDirection),G4ThreeVector(-yDirection,xDirection,0)); + px = PVec->x(); + py = PVec->y(); + pz = PVec->z(); + + // Rotate also beam spot + G4ThreeVector* RVec = new G4ThreeVector(x-x0,y-y0,z-z0); + RVec->rotate(std::acos(zDirection),G4ThreeVector(-yDirection,xDirection,0)); + x = x0+RVec->x(); + y = y0+RVec->y(); + z = z0+RVec->z(); + } // Assign spin G4double xpolaris=0, ypolaris=0, zpolaris=0;