src/QweakSimTrap.cc

Go to the documentation of this file.

#include "QweakSimTrap.hh"
#include "G4SystemOfUnits.hh"

G4double QweakSimTrap::FixPlane(G4ThreeVector pt[8], G4int iface[4])
{
  G4ThreeVector& p1 = pt[iface[0]];
  G4ThreeVector& p2 = pt[iface[1]];
  G4ThreeVector& p3 = pt[iface[2]];
  G4ThreeVector& p4 = pt[iface[3]];

  struct TrapSidePlane
  {
    G4double a,b,c,d;    // Normal unit vector (a,b,c)  and offset (d)
  } plane;

  G4ThreeVector normal = ((p4 - p2).cross(p3 - p1)).unit();
  if (std::abs(normal.x()) < DBL_EPSILON) normal.setX(0);
  if (std::abs(normal.y()) < DBL_EPSILON) normal.setY(0);
  if (std::abs(normal.z()) < DBL_EPSILON) normal.setZ(0);
  normal = normal.unit();

  G4ThreeVector centre = (p1 + p2 + p3 + p4)*0.25;
  plane.a =  normal.x();
  plane.b =  normal.y();
  plane.c =  normal.z();
  plane.d = -normal.dot(centre);

  // compute distances and check planarity
  G4double sd1 = (normal.dot(p1) + plane.d);
  G4double sd2 = (normal.dot(p2) + plane.d);
  G4double sd3 = (normal.dot(p3) + plane.d);
  G4double sd4 = (normal.dot(p4) + plane.d);
  if (std::sqrt(sd1*sd1+sd2*sd2+sd3*sd3+sd4*sd4)/4.0 > 1*mm) {
    G4cerr << "Warning! Not even QweakSimTrap is going to help here." << G4endl;
    G4cerr << "There is something wrong with your QweakSimTrap parameters." << G4endl;
  }

  // fix along xhat
  G4ThreeVector xhat(1.0,0.0,0.0);
  G4double xd = (normal.dot(centre + xhat) + plane.d);
  p1 -= sd1/xd * xhat;
  p2 -= sd2/xd * xhat;
  p3 -= sd3/xd * xhat;
  p4 -= sd4/xd * xhat;

  // compute distances and check planarity
  G4double sd1c = (normal.dot(p1) + plane.d);
  G4double sd2c = (normal.dot(p2) + plane.d);
  G4double sd3c = (normal.dot(p3) + plane.d);
  G4double sd4c = (normal.dot(p4) + plane.d);

  return std::sqrt(sd1c*sd1c+sd2c*sd2c+sd3c*sd3c+sd4c*sd4c)/4.0;
}

QweakSimTrap::QweakSimTrap(const G4String& pName,
  G4double pDz, G4double pTheta, G4double pPhi,
  G4double pDy1, G4double pDx1, G4double pDx2, G4double pAlp1,
  G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlp2)
: G4Trap(pName)
{
  G4cout << "You are using a QweakSim trapezoid." << G4endl;

  G4double fDz = pDz;
  G4double fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
  G4double fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);

  G4double fDy1, fDx1, fDx2, fTalpha1;
  G4double fDy2, fDx3, fDx4, fTalpha2;
  fDy1 = pDy1; fDx1 = pDx1; fDx2 = pDx2; fTalpha1 = std::tan(pAlp1);
  fDy2 = pDy2; fDx3 = pDx3; fDx4 = pDx4; fTalpha2 = std::tan(pAlp2);

  G4double DzTthetaCphi = fDz*fTthetaCphi;
  G4double DzTthetaSphi = fDz*fTthetaSphi;
  G4double Dy1Talpha1   = fDy1*fTalpha1;
  G4double Dy2Talpha2   = fDy2*fTalpha2;

  G4ThreeVector pt[8] =
  {
    G4ThreeVector(-DzTthetaCphi-Dy1Talpha1-fDx1,-DzTthetaSphi-fDy1,-fDz),
    G4ThreeVector(-DzTthetaCphi-Dy1Talpha1+fDx1,-DzTthetaSphi-fDy1,-fDz),
    G4ThreeVector(-DzTthetaCphi+Dy1Talpha1-fDx2,-DzTthetaSphi+fDy1,-fDz),
    G4ThreeVector(-DzTthetaCphi+Dy1Talpha1+fDx2,-DzTthetaSphi+fDy1,-fDz),
    G4ThreeVector( DzTthetaCphi-Dy2Talpha2-fDx3, DzTthetaSphi-fDy2, fDz),
    G4ThreeVector( DzTthetaCphi-Dy2Talpha2+fDx3, DzTthetaSphi-fDy2, fDz),
    G4ThreeVector( DzTthetaCphi+Dy2Talpha2-fDx4, DzTthetaSphi+fDy2, fDz),
    G4ThreeVector( DzTthetaCphi+Dy2Talpha2+fDx4, DzTthetaSphi+fDy2, fDz)
  };

  G4int iface[2][4] = { {0,2,6,4}, {1,3,5,7} }; // ~-X, ~+X
  for (G4int i=0; i<2; ++i)
  {
    G4double d = FixPlane(pt,iface[i]);
  }

  fDx1 = 0.5*(pt[1]-pt[0]).x();
  fDx2 = 0.5*(pt[3]-pt[2]).x();
  fDx3 = 0.5*(pt[5]-pt[4]).x();
  fDx4 = 0.5*(pt[7]-pt[6]).x();

  // Pass to G4Trap
  SetAllParameters(pDz, pTheta, pPhi, pDy1, fDx1, fDx2, pAlp1, pDy2, fDx3, fDx4, pAlp2);
}