include/Geant4/G4TwistTrapAlphaSide.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 // G4TwistTrapAlphaSide
0027 //
0028 // Class description:
0029 //
0030 // Class describing a twisted boundary surface for a trapezoid.
0031
0032 // Author: 27-Oct-2004 - O.Link (Oliver.Link@cern.ch)
0033 // --------------------------------------------------------------------
0034 #ifndef G4TWISTTRAPALPHASIDE_HH
0035 #define G4TWISTTRAPALPHASIDE_HH
0036
0037 #include "G4VTwistSurface.hh"
0038
0039 #include <vector>
0040
0041 class G4TwistTrapAlphaSide : public G4VTwistSurface
0042 {
0043   public:
0044
0045     G4TwistTrapAlphaSide(const G4String& name,
0046                                G4double  PhiTwist, // twist angle
0047                                G4double  pDz,      // half z lenght
0048                                G4double  pTheta, // direction between end planes
0049                                G4double  pPhi,   // by polar and azimutal angles
0050                                G4double  pDy1,     // half y length at -pDz
0051                                G4double  pDx1,     // half x length at -pDz,-pDy
0052                                G4double  pDx2,     // half x length at -pDz,+pDy
0053                                G4double  pDy2,     // half y length at +pDz
0054                                G4double  pDx3,     // half x length at +pDz,-pDy
0055                                G4double  pDx4,     // half x length at +pDz,+pDy
0056                                G4double  pAlph,    // tilt angle at +pDz
0057                                G4double  AngleSide // parity
0058                          );
0059
0060     ~G4TwistTrapAlphaSide() override;
0061
0062     G4ThreeVector  GetNormal(const G4ThreeVector& xx,
0063                                            G4bool isGlobal = false) override ;
0064
0065     G4int DistanceToSurface(const G4ThreeVector& gp,
0066                             const G4ThreeVector& gv,
0067                                   G4ThreeVector  gxx[],
0068                                   G4double  distance[],
0069                                   G4int     areacode[],
0070                                   G4bool    isvalid[],
0071                             EValidate validate = kValidateWithTol) override;
0072
0073     G4int DistanceToSurface(const G4ThreeVector& gp,
0074                                   G4ThreeVector  gxx[],
0075                                   G4double       distance[],
0076                                   G4int          areacode[]) override;
0077
0078     G4TwistTrapAlphaSide(__void__&);
0079       // Fake default constructor for usage restricted to direct object
0080       // persistency for clients requiring preallocation of memory for
0081       // persistifiable objects.
0082
0083   private:
0084
0085     G4int GetAreaCode(const G4ThreeVector& xx,
0086                             G4bool withTol = true) override;
0087     void SetCorners() override;
0088     void SetBoundaries() override;
0089
0090     void GetPhiUAtX(const G4ThreeVector& p, G4double& phi, G4double& u);
0091     G4ThreeVector ProjectPoint(const G4ThreeVector& p,
0092                                      G4bool isglobal = false);
0093
0094     inline G4ThreeVector SurfacePoint(G4double phi, G4double u,
0095                                       G4bool isGlobal = false ) override;
0096     inline G4double GetBoundaryMin(G4double phi) override;
0097     inline G4double GetBoundaryMax(G4double phi) override;
0098     inline G4double GetSurfaceArea() override;
0099     void GetFacets( G4int m, G4int n, G4double xyz[][3],
0100                     G4int faces[][4], G4int iside ) override;
0101
0102     inline G4ThreeVector NormAng(G4double phi, G4double u);
0103     inline G4double GetValueA(G4double phi);
0104     inline G4double GetValueB(G4double phi);
0105     inline G4double GetValueD(G4double phi);
0106     inline G4double Xcoef(G4double u,G4double phi);
0107       // To calculate the w(u) function
0108
0109   private:
0110
0111     G4double fTheta;
0112     G4double fPhi ;
0113
0114     G4double fDy1;
0115     G4double fDx1;
0116     G4double fDx2;
0117
0118     G4double fDy2;
0119     G4double fDx3;
0120     G4double fDx4;
0121
0122     G4double fDz;         // Half-length along the z axis
0123
0124     G4double fAlph;
0125     G4double fTAlph;      // std::tan(fAlph)
0126
0127     G4double fPhiTwist;   // twist angle (dphi in surface equation)
0128
0129     G4double fAngleSide;
0130
0131     G4double fDx4plus2;  // fDx4 + fDx2  == a2/2 + a1/2
0132     G4double fDx4minus2; // fDx4 - fDx2          -
0133     G4double fDx3plus1;  // fDx3 + fDx1  == d2/2 + d1/2
0134     G4double fDx3minus1; // fDx3 - fDx1          -
0135     G4double fDy2plus1;  // fDy2 + fDy1  == b2/2 + b1/2
0136     G4double fDy2minus1; // fDy2 - fDy1          -
0137     G4double fa1md1;     // 2 fDx2 - 2 fDx1  == a1 - d1
0138     G4double fa2md2;     // 2 fDx4 - 2 fDx3
0139
0140     G4double fdeltaX;
0141     G4double fdeltaY;
0142 };
0143
0144 //========================================================
0145 // inline functions
0146 //========================================================
0147
0148 inline
0149 G4double G4TwistTrapAlphaSide::GetValueA(G4double phi)
0150 {
0151   return ( fDx4plus2 + fDx4minus2 * ( 2 * phi ) / fPhiTwist  ) ;
0152 }
0153
0154 inline
0155 G4double G4TwistTrapAlphaSide::GetValueD(G4double phi)
0156 {
0157   return ( fDx3plus1 + fDx3minus1 * ( 2 * phi) / fPhiTwist  ) ;
0158 }
0159
0160 inline
0161 G4double G4TwistTrapAlphaSide::GetValueB(G4double phi)
0162 {
0163   return ( fDy2plus1 + fDy2minus1 * ( 2 * phi ) / fPhiTwist ) ;
0164 }
0165
0166
0167 inline
0168 G4double G4TwistTrapAlphaSide::Xcoef(G4double u, G4double phi)
0169 {
0170
0171   return GetValueA(phi)/2. + (GetValueD(phi)-GetValueA(phi))/4.
0172     - u*( ( GetValueD(phi)-GetValueA(phi) )/( 2 * GetValueB(phi) ) - fTAlph );
0173
0174 }
0175
0176 inline G4ThreeVector
0177 G4TwistTrapAlphaSide::SurfacePoint(G4double phi, G4double u , G4bool isGlobal)
0178 {
0179   // function to calculate a point on the surface, given by parameters phi,u
0180
0181   G4ThreeVector SurfPoint ( Xcoef(u,phi) * std::cos(phi)
0182                           - u * std::sin(phi) + fdeltaX*phi/fPhiTwist,
0183                             Xcoef(u,phi) * std::sin(phi)
0184                           + u * std::cos(phi) + fdeltaY*phi/fPhiTwist,
0185                             2*fDz*phi/fPhiTwist  );
0186   if (isGlobal) { return (fRot * SurfPoint + fTrans); }
0187   return SurfPoint;
0188 }
0189
0190 inline
0191 G4double G4TwistTrapAlphaSide::GetBoundaryMin(G4double phi)
0192 {
0193   return -0.5*GetValueB(phi) ;
0194 }
0195
0196 inline
0197 G4double G4TwistTrapAlphaSide::GetBoundaryMax(G4double phi)
0198 {
0199   return 0.5*GetValueB(phi) ;
0200 }
0201
0202 inline
0203 G4double G4TwistTrapAlphaSide::GetSurfaceArea()
0204 {
0205   return (fDz*(std::sqrt(16*fDy1*fDy1
0206              + (fa1md1 + 4*fDy1*fTAlph)*(fa1md1 + 4*fDy1*fTAlph))
0207              + std::sqrt(16*fDy2*fDy2 + (fa2md2 + 4*fDy2*fTAlph)
0208                                       * (fa2md2 + 4*fDy2*fTAlph))))/2. ;
0209 }
0210
0211 inline
0212 G4ThreeVector G4TwistTrapAlphaSide::NormAng( G4double phi, G4double u )
0213 {
0214   // function to calculate the norm at a given point on the surface
0215   // replace a1-d1
0216
0217   G4ThreeVector nvec ( fDy1* fDz*(4*fDy1*std::cos(phi)
0218                      + (fa1md1 + 4*fDy1*fTAlph)*std::sin(phi)),
0219                        -(fDy1* fDz*((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi)
0220                      - 4*fDy1*std::sin(phi))),
0221                        (fDy1*(-8*(fDx3minus1 + fDx4minus2)*fDy1
0222                                 + fa1md1*(fDx2 + fDx3plus1 + fDx4)*fPhiTwist
0223                                 + 4*(fDx2 + fDx3plus1 + fDx4)*fDy1*fPhiTwist
0224                                 *fTAlph + 2*(fDx3minus1 + fDx4minus2)
0225                                 *(fa1md1 + 4*fDy1*fTAlph)*phi)
0226                                 + fPhiTwist*(16*fDy1*fDy1
0227                                 + (fa1md1 + 4*fDy1*fTAlph)
0228                                 *(fa1md1 + 4*fDy1*fTAlph))*u
0229                                 + 4*fDy1*(fa1md1*fdeltaY - 4*fdeltaX*fDy1
0230                                 + 4*fdeltaY*fDy1*fTAlph)* std::cos(phi)
0231                                 - 4*fDy1*(fa1md1*fdeltaX + 4*fDy1*(fdeltaY
0232                                 + fdeltaX*fTAlph))*std::sin(phi))/ 8. ) ;
0233   return nvec.unit();
0234 }
0235
0236 #endif