| 1 | //
|
|---|
| 2 | // ********************************************************************
|
|---|
| 3 | // * License and Disclaimer *
|
|---|
| 4 | // * *
|
|---|
| 5 | // * The Geant4 software is copyright of the Copyright Holders of *
|
|---|
| 6 | // * the Geant4 Collaboration. It is provided under the terms and *
|
|---|
| 7 | // * conditions of the Geant4 Software License, included in the file *
|
|---|
| 8 | // * LICENSE and available at http://cern.ch/geant4/license . These *
|
|---|
| 9 | // * include a list of copyright holders. *
|
|---|
| 10 | // * *
|
|---|
| 11 | // * Neither the authors of this software system, nor their employing *
|
|---|
| 12 | // * institutes,nor the agencies providing financial support for this *
|
|---|
| 13 | // * work make any representation or warranty, express or implied, *
|
|---|
| 14 | // * regarding this software system or assume any liability for its *
|
|---|
| 15 | // * use. Please see the license in the file LICENSE and URL above *
|
|---|
| 16 | // * for the full disclaimer and the limitation of liability. *
|
|---|
| 17 | // * *
|
|---|
| 18 | // * This code implementation is the result of the scientific and *
|
|---|
| 19 | // * technical work of the GEANT4 collaboration. *
|
|---|
| 20 | // * By using, copying, modifying or distributing the software (or *
|
|---|
| 21 | // * any work based on the software) you agree to acknowledge its *
|
|---|
| 22 | // * use in resulting scientific publications, and indicate your *
|
|---|
| 23 | // * acceptance of all terms of the Geant4 Software license. *
|
|---|
| 24 | // ********************************************************************
|
|---|
| 25 | //
|
|---|
| 26 | //
|
|---|
| 27 | // $Id: G4TwistTrapAlphaSide.hh,v 1.5 2006/06/29 18:47:46 gunter Exp $
|
|---|
| 28 | //
|
|---|
| 29 | // --------------------------------------------------------------------
|
|---|
| 30 | // GEANT 4 class header file
|
|---|
| 31 | //
|
|---|
| 32 | //
|
|---|
| 33 | // G4TwistTrapAlphaSide
|
|---|
| 34 | //
|
|---|
| 35 | // Class description:
|
|---|
| 36 | //
|
|---|
| 37 | // Class describing a twisted boundary surface for a trapezoid.
|
|---|
| 38 |
|
|---|
| 39 | // Author:
|
|---|
| 40 | //
|
|---|
| 41 | // 27-Oct-2004 - O.Link (Oliver.Link@cern.ch)
|
|---|
| 42 | //
|
|---|
| 43 | // --------------------------------------------------------------------
|
|---|
| 44 | #ifndef __G4TWISTTRAPALPHASIDE__
|
|---|
| 45 | #define __G4TWISTTRAPALPHASIDE__
|
|---|
| 46 |
|
|---|
| 47 | #include "G4VTwistSurface.hh"
|
|---|
| 48 |
|
|---|
| 49 | #include <vector>
|
|---|
| 50 |
|
|---|
| 51 | class G4TwistTrapAlphaSide : public G4VTwistSurface
|
|---|
| 52 | {
|
|---|
| 53 | public: // with description
|
|---|
| 54 |
|
|---|
| 55 | G4TwistTrapAlphaSide(const G4String &name,
|
|---|
| 56 | G4double PhiTwist, // twist angle
|
|---|
| 57 | G4double pDz, // half z lenght
|
|---|
| 58 | G4double pTheta, // direction between end planes
|
|---|
| 59 | G4double pPhi, // by polar and azimutal angles
|
|---|
| 60 | G4double pDy1, // half y length at -pDz
|
|---|
| 61 | G4double pDx1, // half x length at -pDz,-pDy
|
|---|
| 62 | G4double pDx2, // half x length at -pDz,+pDy
|
|---|
| 63 | G4double pDy2, // half y length at +pDz
|
|---|
| 64 | G4double pDx3, // half x length at +pDz,-pDy
|
|---|
| 65 | G4double pDx4, // half x length at +pDz,+pDy
|
|---|
| 66 | G4double pAlph, // tilt angle at +pDz
|
|---|
| 67 | G4double AngleSide // parity
|
|---|
| 68 | );
|
|---|
| 69 |
|
|---|
| 70 | virtual ~G4TwistTrapAlphaSide();
|
|---|
| 71 |
|
|---|
| 72 | virtual G4ThreeVector GetNormal(const G4ThreeVector &xx,
|
|---|
| 73 | G4bool isGlobal = false) ;
|
|---|
| 74 |
|
|---|
| 75 | virtual G4int DistanceToSurface(const G4ThreeVector &gp,
|
|---|
| 76 | const G4ThreeVector &gv,
|
|---|
| 77 | G4ThreeVector gxx[],
|
|---|
| 78 | G4double distance[],
|
|---|
| 79 | G4int areacode[],
|
|---|
| 80 | G4bool isvalid[],
|
|---|
| 81 | EValidate validate = kValidateWithTol);
|
|---|
| 82 |
|
|---|
| 83 | virtual G4int DistanceToSurface(const G4ThreeVector &gp,
|
|---|
| 84 | G4ThreeVector gxx[],
|
|---|
| 85 | G4double distance[],
|
|---|
| 86 | G4int areacode[]);
|
|---|
| 87 |
|
|---|
| 88 |
|
|---|
| 89 | public: // without description
|
|---|
| 90 |
|
|---|
| 91 | G4TwistTrapAlphaSide(__void__&);
|
|---|
| 92 | // Fake default constructor for usage restricted to direct object
|
|---|
| 93 | // persistency for clients requiring preallocation of memory for
|
|---|
| 94 | // persistifiable objects.
|
|---|
| 95 |
|
|---|
| 96 | private:
|
|---|
| 97 |
|
|---|
| 98 | virtual G4int GetAreaCode(const G4ThreeVector &xx,
|
|---|
| 99 | G4bool withTol = true);
|
|---|
| 100 | virtual void SetCorners();
|
|---|
| 101 | virtual void SetBoundaries();
|
|---|
| 102 |
|
|---|
| 103 | void GetPhiUAtX(G4ThreeVector p, G4double &phi, G4double &u);
|
|---|
| 104 | G4ThreeVector ProjectPoint(const G4ThreeVector &p,
|
|---|
| 105 | G4bool isglobal = false);
|
|---|
| 106 |
|
|---|
| 107 | virtual G4ThreeVector SurfacePoint(G4double phi, G4double u,
|
|---|
| 108 | G4bool isGlobal = false );
|
|---|
| 109 | virtual G4double GetBoundaryMin(G4double phi);
|
|---|
| 110 | virtual G4double GetBoundaryMax(G4double phi);
|
|---|
| 111 | virtual G4double GetSurfaceArea();
|
|---|
| 112 | virtual void GetFacets( G4int m, G4int n, G4double xyz[][3],
|
|---|
| 113 | G4int faces[][4], G4int iside );
|
|---|
| 114 |
|
|---|
| 115 | inline G4ThreeVector NormAng(G4double phi, G4double u);
|
|---|
| 116 | inline G4double GetValueA(G4double phi);
|
|---|
| 117 | inline G4double GetValueB(G4double phi);
|
|---|
| 118 | inline G4double GetValueD(G4double phi);
|
|---|
| 119 | inline G4double Xcoef(G4double u,G4double phi);
|
|---|
| 120 | // To calculate the w(u) function
|
|---|
| 121 |
|
|---|
| 122 | private:
|
|---|
| 123 |
|
|---|
| 124 | G4double fTheta;
|
|---|
| 125 | G4double fPhi ;
|
|---|
| 126 |
|
|---|
| 127 | G4double fDy1;
|
|---|
| 128 | G4double fDx1;
|
|---|
| 129 | G4double fDx2;
|
|---|
| 130 |
|
|---|
| 131 | G4double fDy2;
|
|---|
| 132 | G4double fDx3;
|
|---|
| 133 | G4double fDx4;
|
|---|
| 134 |
|
|---|
| 135 | G4double fDz; // Half-length along the z axis
|
|---|
| 136 |
|
|---|
| 137 | G4double fAlph;
|
|---|
| 138 | G4double fTAlph; // std::tan(fAlph)
|
|---|
| 139 |
|
|---|
| 140 | G4double fPhiTwist; // twist angle (dphi in surface equation)
|
|---|
| 141 |
|
|---|
| 142 | G4double fAngleSide;
|
|---|
| 143 |
|
|---|
| 144 | G4double fDx4plus2; // fDx4 + fDx2 == a2/2 + a1/2
|
|---|
| 145 | G4double fDx4minus2; // fDx4 - fDx2 -
|
|---|
| 146 | G4double fDx3plus1; // fDx3 + fDx1 == d2/2 + d1/2
|
|---|
| 147 | G4double fDx3minus1; // fDx3 - fDx1 -
|
|---|
| 148 | G4double fDy2plus1; // fDy2 + fDy1 == b2/2 + b1/2
|
|---|
| 149 | G4double fDy2minus1; // fDy2 - fDy1 -
|
|---|
| 150 | G4double fa1md1; // 2 fDx2 - 2 fDx1 == a1 - d1
|
|---|
| 151 | G4double fa2md2; // 2 fDx4 - 2 fDx3
|
|---|
| 152 |
|
|---|
| 153 | G4double fdeltaX;
|
|---|
| 154 | G4double fdeltaY;
|
|---|
| 155 | };
|
|---|
| 156 |
|
|---|
| 157 | //========================================================
|
|---|
| 158 | // inline functions
|
|---|
| 159 | //========================================================
|
|---|
| 160 |
|
|---|
| 161 | inline
|
|---|
| 162 | G4double G4TwistTrapAlphaSide::GetValueA(G4double phi)
|
|---|
| 163 | {
|
|---|
| 164 | return ( fDx4plus2 + fDx4minus2 * ( 2 * phi ) / fPhiTwist ) ;
|
|---|
| 165 | }
|
|---|
| 166 |
|
|---|
| 167 | inline
|
|---|
| 168 | G4double G4TwistTrapAlphaSide::GetValueD(G4double phi)
|
|---|
| 169 | {
|
|---|
| 170 | return ( fDx3plus1 + fDx3minus1 * ( 2 * phi) / fPhiTwist ) ;
|
|---|
| 171 | }
|
|---|
| 172 |
|
|---|
| 173 | inline
|
|---|
| 174 | G4double G4TwistTrapAlphaSide::GetValueB(G4double phi)
|
|---|
| 175 | {
|
|---|
| 176 | return ( fDy2plus1 + fDy2minus1 * ( 2 * phi ) / fPhiTwist ) ;
|
|---|
| 177 | }
|
|---|
| 178 |
|
|---|
| 179 |
|
|---|
| 180 | inline
|
|---|
| 181 | G4double G4TwistTrapAlphaSide::Xcoef(G4double u, G4double phi)
|
|---|
| 182 | {
|
|---|
| 183 |
|
|---|
| 184 | return GetValueA(phi)/2. + (GetValueD(phi)-GetValueA(phi))/4.
|
|---|
| 185 | - u*( ( GetValueD(phi)-GetValueA(phi) )/( 2 * GetValueB(phi) ) - fTAlph );
|
|---|
| 186 |
|
|---|
| 187 | }
|
|---|
| 188 |
|
|---|
| 189 | inline G4ThreeVector
|
|---|
| 190 | G4TwistTrapAlphaSide::SurfacePoint(G4double phi, G4double u , G4bool isGlobal)
|
|---|
| 191 | {
|
|---|
| 192 | // function to calculate a point on the surface, given by parameters phi,u
|
|---|
| 193 |
|
|---|
| 194 | G4ThreeVector SurfPoint ( Xcoef(u,phi) * std::cos(phi)
|
|---|
| 195 | - u * std::sin(phi) + fdeltaX*phi/fPhiTwist,
|
|---|
| 196 | Xcoef(u,phi) * std::sin(phi)
|
|---|
| 197 | + u * std::cos(phi) + fdeltaY*phi/fPhiTwist,
|
|---|
| 198 | 2*fDz*phi/fPhiTwist );
|
|---|
| 199 | if (isGlobal) { return (fRot * SurfPoint + fTrans); }
|
|---|
| 200 | return SurfPoint;
|
|---|
| 201 | }
|
|---|
| 202 |
|
|---|
| 203 | inline
|
|---|
| 204 | G4double G4TwistTrapAlphaSide::GetBoundaryMin(G4double phi)
|
|---|
| 205 | {
|
|---|
| 206 | return -0.5*GetValueB(phi) ;
|
|---|
| 207 | }
|
|---|
| 208 |
|
|---|
| 209 | inline
|
|---|
| 210 | G4double G4TwistTrapAlphaSide::GetBoundaryMax(G4double phi)
|
|---|
| 211 | {
|
|---|
| 212 | return 0.5*GetValueB(phi) ;
|
|---|
| 213 | }
|
|---|
| 214 |
|
|---|
| 215 | inline
|
|---|
| 216 | G4double G4TwistTrapAlphaSide::GetSurfaceArea()
|
|---|
| 217 | {
|
|---|
| 218 | return (fDz*(std::sqrt(16*fDy1*fDy1
|
|---|
| 219 | + (fa1md1 + 4*fDy1*fTAlph)*(fa1md1 + 4*fDy1*fTAlph))
|
|---|
| 220 | + std::sqrt(16*fDy2*fDy2 + (fa2md2 + 4*fDy2*fTAlph)
|
|---|
| 221 | * (fa2md2 + 4*fDy2*fTAlph))))/2. ;
|
|---|
| 222 | }
|
|---|
| 223 |
|
|---|
| 224 | inline
|
|---|
| 225 | G4ThreeVector G4TwistTrapAlphaSide::NormAng( G4double phi, G4double u )
|
|---|
| 226 | {
|
|---|
| 227 | // function to calculate the norm at a given point on the surface
|
|---|
| 228 | // replace a1-d1
|
|---|
| 229 |
|
|---|
| 230 | G4ThreeVector nvec ( fDy1* fDz*(4*fDy1*std::cos(phi)
|
|---|
| 231 | + (fa1md1 + 4*fDy1*fTAlph)*std::sin(phi)),
|
|---|
| 232 | -(fDy1* fDz*((fa1md1 + 4*fDy1*fTAlph)*std::cos(phi)
|
|---|
| 233 | - 4*fDy1*std::sin(phi))),
|
|---|
| 234 | (fDy1*(-8*(fDx3minus1 + fDx4minus2)*fDy1
|
|---|
| 235 | + fa1md1*(fDx2 + fDx3plus1 + fDx4)*fPhiTwist
|
|---|
| 236 | + 4*(fDx2 + fDx3plus1 + fDx4)*fDy1*fPhiTwist
|
|---|
| 237 | *fTAlph + 2*(fDx3minus1 + fDx4minus2)
|
|---|
| 238 | *(fa1md1 + 4*fDy1*fTAlph)*phi)
|
|---|
| 239 | + fPhiTwist*(16*fDy1*fDy1
|
|---|
| 240 | + (fa1md1 + 4*fDy1*fTAlph)
|
|---|
| 241 | *(fa1md1 + 4*fDy1*fTAlph))*u
|
|---|
| 242 | + 4*fDy1*(fa1md1*fdeltaY - 4*fdeltaX*fDy1
|
|---|
| 243 | + 4*fdeltaY*fDy1*fTAlph)* std::cos(phi)
|
|---|
| 244 | - 4*fDy1*(fa1md1*fdeltaX + 4*fDy1*(fdeltaY
|
|---|
| 245 | + fdeltaX*fTAlph))*std::sin(phi))/ 8. ) ;
|
|---|
| 246 | return nvec.unit();
|
|---|
| 247 | }
|
|---|
| 248 |
|
|---|
| 249 | #endif
|
|---|