Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // G4ParameterisationTrd[X/Y/Z] implementation 26 // G4ParameterisationTrd[X/Y/Z] implementation 27 // 27 // 28 // 26.05.03 - P.Arce, Initial version 28 // 26.05.03 - P.Arce, Initial version 29 // 08.04.04 - I.Hrivnacova, Implemented reflec 29 // 08.04.04 - I.Hrivnacova, Implemented reflection 30 // 21.04.10 - M.Asai, Added gaps 30 // 21.04.10 - M.Asai, Added gaps 31 // ------------------------------------------- 31 // -------------------------------------------------------------------- 32 32 33 #include "G4ParameterisationTrd.hh" 33 #include "G4ParameterisationTrd.hh" 34 34 35 #include <iomanip> 35 #include <iomanip> 36 #include "G4ThreeVector.hh" 36 #include "G4ThreeVector.hh" 37 #include "G4RotationMatrix.hh" 37 #include "G4RotationMatrix.hh" 38 #include "G4VPhysicalVolume.hh" 38 #include "G4VPhysicalVolume.hh" 39 #include "G4LogicalVolume.hh" 39 #include "G4LogicalVolume.hh" 40 #include "G4ReflectedSolid.hh" 40 #include "G4ReflectedSolid.hh" 41 #include "G4Trd.hh" 41 #include "G4Trd.hh" 42 #include "G4Trap.hh" 42 #include "G4Trap.hh" 43 43 44 //-------------------------------------------- 44 //-------------------------------------------------------------------------- 45 G4VParameterisationTrd:: 45 G4VParameterisationTrd:: 46 G4VParameterisationTrd( EAxis axis, G4int nDiv 46 G4VParameterisationTrd( EAxis axis, G4int nDiv, G4double width, 47 G4double offset, G4VSo 47 G4double offset, G4VSolid* msolid, 48 DivisionType divType ) 48 DivisionType divType ) 49 : G4VDivisionParameterisation( axis, nDiv, 49 : G4VDivisionParameterisation( axis, nDiv, width, offset, divType, msolid ) 50 { 50 { 51 auto msol = (G4Trd*)(msolid); << 51 G4Trd* msol = (G4Trd*)(msolid); 52 if (msolid->GetEntityType() == "G4ReflectedS 52 if (msolid->GetEntityType() == "G4ReflectedSolid") 53 { 53 { 54 // Get constituent solid 54 // Get constituent solid 55 G4VSolid* mConstituentSolid 55 G4VSolid* mConstituentSolid 56 = ((G4ReflectedSolid*)msolid)->GetConst 56 = ((G4ReflectedSolid*)msolid)->GetConstituentMovedSolid(); 57 msol = (G4Trd*)(mConstituentSolid); 57 msol = (G4Trd*)(mConstituentSolid); 58 58 59 // Create a new solid with inversed parame 59 // Create a new solid with inversed parameters 60 auto newSolid << 60 G4Trd* newSolid 61 = new G4Trd(msol->GetName(), 61 = new G4Trd(msol->GetName(), 62 msol->GetXHalfLength2(), mso 62 msol->GetXHalfLength2(), msol->GetXHalfLength1(), 63 msol->GetYHalfLength2(), mso 63 msol->GetYHalfLength2(), msol->GetYHalfLength1(), 64 msol->GetZHalfLength()); 64 msol->GetZHalfLength()); 65 msol = newSolid; 65 msol = newSolid; 66 fmotherSolid = newSolid; 66 fmotherSolid = newSolid; 67 fReflectedSolid = true; 67 fReflectedSolid = true; 68 fDeleteSolid = true; 68 fDeleteSolid = true; 69 } 69 } 70 } 70 } 71 71 72 //-------------------------------------------- 72 //------------------------------------------------------------------------ 73 G4VParameterisationTrd::~G4VParameterisationTr << 73 G4VParameterisationTrd::~G4VParameterisationTrd() >> 74 { >> 75 } 74 76 75 //-------------------------------------------- 77 //------------------------------------------------------------------------ 76 G4ParameterisationTrdX:: 78 G4ParameterisationTrdX:: 77 G4ParameterisationTrdX( EAxis axis, G4int nDiv 79 G4ParameterisationTrdX( EAxis axis, G4int nDiv, 78 G4double width, G4doub 80 G4double width, G4double offset, 79 G4VSolid* msolid, Divi 81 G4VSolid* msolid, DivisionType divType ) 80 : G4VParameterisationTrd( axis, nDiv, width 82 : G4VParameterisationTrd( axis, nDiv, width, offset, msolid, divType ) 81 { 83 { 82 CheckParametersValidity(); 84 CheckParametersValidity(); 83 SetType( "DivisionTrdX" ); 85 SetType( "DivisionTrdX" ); 84 86 85 auto msol = (G4Trd*)(fmotherSolid); << 87 G4Trd* msol = (G4Trd*)(fmotherSolid); 86 if( divType == DivWIDTH ) 88 if( divType == DivWIDTH ) 87 { 89 { 88 fnDiv = CalculateNDiv( msol->GetXHalfLengt 90 fnDiv = CalculateNDiv( msol->GetXHalfLength1()+msol->GetXHalfLength2(), 89 width, offset ); 91 width, offset ); 90 } 92 } 91 else if( divType == DivNDIV ) 93 else if( divType == DivNDIV ) 92 { 94 { 93 fwidth = CalculateWidth( msol->GetXHalfLen 95 fwidth = CalculateWidth( msol->GetXHalfLength1()+msol->GetXHalfLength2(), 94 nDiv, offset ); 96 nDiv, offset ); 95 } 97 } 96 98 97 #ifdef G4DIVDEBUG 99 #ifdef G4DIVDEBUG 98 if( verbose >= 1 ) 100 if( verbose >= 1 ) 99 { 101 { 100 G4cout << " G4ParameterisationTrdX - ## di 102 G4cout << " G4ParameterisationTrdX - ## divisions " << fnDiv << " = " 101 << nDiv << G4endl 103 << nDiv << G4endl 102 << " Offset " << foffset << " = " < 104 << " Offset " << foffset << " = " << offset << G4endl 103 << " Width " << fwidth << " = " << 105 << " Width " << fwidth << " = " << width << G4endl; 104 } 106 } 105 #endif 107 #endif 106 108 107 G4double mpDx1 = msol->GetXHalfLength1(); 109 G4double mpDx1 = msol->GetXHalfLength1(); 108 G4double mpDx2 = msol->GetXHalfLength2(); 110 G4double mpDx2 = msol->GetXHalfLength2(); 109 if( std::fabs(mpDx1 - mpDx2) > kCarTolerance 111 if( std::fabs(mpDx1 - mpDx2) > kCarTolerance ) 110 { 112 { 111 bDivInTrap = true; 113 bDivInTrap = true; 112 } 114 } 113 } 115 } 114 116 115 //-------------------------------------------- 117 //------------------------------------------------------------------------ 116 G4ParameterisationTrdX::~G4ParameterisationTrd 118 G4ParameterisationTrdX::~G4ParameterisationTrdX() 117 = default; << 119 { >> 120 } 118 121 119 //-------------------------------------------- 122 //------------------------------------------------------------------------ 120 G4double G4ParameterisationTrdX::GetMaxParamet 123 G4double G4ParameterisationTrdX::GetMaxParameter() const 121 { 124 { 122 auto msol = (G4Trd*)(fmotherSolid); << 125 G4Trd* msol = (G4Trd*)(fmotherSolid); 123 return (msol->GetXHalfLength1()+msol->GetXHa 126 return (msol->GetXHalfLength1()+msol->GetXHalfLength2()); 124 } 127 } 125 128 126 //-------------------------------------------- 129 //------------------------------------------------------------------------ 127 void 130 void 128 G4ParameterisationTrdX:: 131 G4ParameterisationTrdX:: 129 ComputeTransformation( const G4int copyNo, 132 ComputeTransformation( const G4int copyNo, 130 G4VPhysicalVolume *phys 133 G4VPhysicalVolume *physVol ) const 131 { 134 { 132 auto msol = (G4Trd*)(fmotherSolid ); << 135 G4Trd* msol = (G4Trd*)(fmotherSolid ); 133 G4double mdx = ( msol->GetXHalfLength1() + m 136 G4double mdx = ( msol->GetXHalfLength1() + msol->GetXHalfLength2() ) / 2.; 134 //----- translation 137 //----- translation 135 G4ThreeVector origin(0.,0.,0.); 138 G4ThreeVector origin(0.,0.,0.); 136 G4double posi; 139 G4double posi; 137 posi = -mdx + foffset + (copyNo+0.5)*fwidth; 140 posi = -mdx + foffset + (copyNo+0.5)*fwidth; 138 if( faxis == kXAxis ) 141 if( faxis == kXAxis ) 139 { 142 { 140 origin.setX( posi ); 143 origin.setX( posi ); 141 } 144 } 142 else 145 else 143 { 146 { 144 std::ostringstream message; 147 std::ostringstream message; 145 message << "Only axes along X are allowed 148 message << "Only axes along X are allowed ! Axis: " << faxis; 146 G4Exception("G4ParameterisationTrdX::Compu 149 G4Exception("G4ParameterisationTrdX::ComputeTransformation()", 147 "GeomDiv0002", FatalException, 150 "GeomDiv0002", FatalException, message); 148 } 151 } 149 152 150 #ifdef G4DIVDEBUG 153 #ifdef G4DIVDEBUG 151 if( verbose >= 2 ) 154 if( verbose >= 2 ) 152 { 155 { 153 G4cout << std::setprecision(8) 156 G4cout << std::setprecision(8) 154 << " G4ParameterisationTrdX::Comput 157 << " G4ParameterisationTrdX::ComputeTransformation() " 155 << copyNo << G4endl 158 << copyNo << G4endl 156 << " Position: " << origin << " - A 159 << " Position: " << origin << " - Axis: " << faxis << G4endl; 157 } 160 } 158 #endif 161 #endif 159 162 160 //----- set translation 163 //----- set translation 161 physVol->SetTranslation( origin ); 164 physVol->SetTranslation( origin ); 162 } 165 } 163 166 164 //-------------------------------------------- 167 //-------------------------------------------------------------------------- 165 void 168 void 166 G4ParameterisationTrdX:: 169 G4ParameterisationTrdX:: 167 ComputeDimensions( G4Trd& trd, [[maybe_unused] << 170 ComputeDimensions( G4Trd& trd, const G4int /* copyNo */, >> 171 const G4VPhysicalVolume* ) const 168 { 172 { 169 auto msol = (G4Trd*)(fmotherSolid); << 173 G4Trd* msol = (G4Trd*)(fmotherSolid); 170 G4double pDy1 = msol->GetYHalfLength1(); 174 G4double pDy1 = msol->GetYHalfLength1(); 171 G4double pDy2 = msol->GetYHalfLength2(); 175 G4double pDy2 = msol->GetYHalfLength2(); 172 G4double pDz = msol->GetZHalfLength(); 176 G4double pDz = msol->GetZHalfLength(); 173 G4double pDx = fwidth/2. - fhgap; 177 G4double pDx = fwidth/2. - fhgap; 174 178 175 trd.SetAllParameters ( pDx, pDx, pDy1, pDy2, 179 trd.SetAllParameters ( pDx, pDx, pDy1, pDy2, pDz ); 176 180 177 #ifdef G4DIVDEBUG 181 #ifdef G4DIVDEBUG 178 if( verbose >= 2 ) 182 if( verbose >= 2 ) 179 { 183 { 180 G4cout << " G4ParameterisationTrdX::Comput 184 G4cout << " G4ParameterisationTrdX::ComputeDimensions():" 181 << copyNo << G4endl; << 185 // << copyNo >> 186 << G4endl; 182 trd.DumpInfo(); 187 trd.DumpInfo(); 183 } 188 } 184 #endif 189 #endif 185 } 190 } 186 191 187 //-------------------------------------------- 192 //-------------------------------------------------------------------------- 188 void 193 void 189 G4ParameterisationTrdX::ComputeDimensions( G4T 194 G4ParameterisationTrdX::ComputeDimensions( G4Trap& trap, const G4int copyNo, 190 con 195 const G4VPhysicalVolume* ) const 191 { 196 { 192 auto msol = (G4Trd*)(fmotherSolid); << 197 G4Trd* msol = (G4Trd*)(fmotherSolid); 193 G4double pDy1 = msol->GetYHalfLength1(); 198 G4double pDy1 = msol->GetYHalfLength1(); 194 G4double pDy2 = msol->GetYHalfLength2(); 199 G4double pDy2 = msol->GetYHalfLength2(); 195 G4double pDz = msol->GetZHalfLength(); 200 G4double pDz = msol->GetZHalfLength(); 196 //fwidth is at Z=0 201 //fwidth is at Z=0 197 G4double pDx1 = msol->GetXHalfLength1(); 202 G4double pDx1 = msol->GetXHalfLength1(); 198 G4double pDx2 = msol->GetXHalfLength2(); 203 G4double pDx2 = msol->GetXHalfLength2(); 199 // G4double pDxAVE = (pDx1+pDx2)/2.; 204 // G4double pDxAVE = (pDx1+pDx2)/2.; 200 G4double xChangeRatio = (pDx2-pDx1) / (pDx2+ 205 G4double xChangeRatio = (pDx2-pDx1) / (pDx2+pDx1); 201 G4double fWidChange = xChangeRatio*fwidth; 206 G4double fWidChange = xChangeRatio*fwidth; 202 G4double fWid1 = fwidth - fWidChange; 207 G4double fWid1 = fwidth - fWidChange; 203 G4double fWid2 = fwidth + fWidChange; 208 G4double fWid2 = fwidth + fWidChange; 204 G4double fOffsetChange = xChangeRatio*foffse 209 G4double fOffsetChange = xChangeRatio*foffset/2.; 205 G4double fOffset1 = foffset - fOffsetChange; 210 G4double fOffset1 = foffset - fOffsetChange; 206 G4double fOffset2 = foffset + fOffsetChange; 211 G4double fOffset2 = foffset + fOffsetChange; 207 G4double cxy1 = -pDx1+fOffset1 + (copyNo+0.5 212 G4double cxy1 = -pDx1+fOffset1 + (copyNo+0.5)*fWid1;// centre of the side at y=-pDy1 208 G4double cxy2 = -pDx2+fOffset2 + (copyNo+0.5 213 G4double cxy2 = -pDx2+fOffset2 + (copyNo+0.5)*fWid2;// centre of the side at y=+pDy1 209 G4double alp = std::atan( (cxy2-cxy1)/(pDz*2 214 G4double alp = std::atan( (cxy2-cxy1)/(pDz*2.) ); 210 215 211 pDx1 = fwidth/2. - fWidChange/2.; 216 pDx1 = fwidth/2. - fWidChange/2.; 212 pDx2 = fwidth/2. + fWidChange/2.; 217 pDx2 = fwidth/2. + fWidChange/2.; 213 218 214 219 215 trap.SetAllParameters ( pDz, 220 trap.SetAllParameters ( pDz, 216 alp, 221 alp, 217 0., 222 0., 218 pDy1, 223 pDy1, 219 pDx1, 224 pDx1, 220 pDx1, 225 pDx1, 221 0., 226 0., 222 pDy2, 227 pDy2, 223 pDx2 - fhgap, 228 pDx2 - fhgap, 224 pDx2 - fhgap * pDx2/pDx1, 229 pDx2 - fhgap * pDx2/pDx1, 225 0.); 230 0.); 226 231 227 #ifdef G4DIVDEBUG 232 #ifdef G4DIVDEBUG 228 if( verbose >= 2 ) 233 if( verbose >= 2 ) 229 { 234 { 230 G4cout << " G4ParameterisationTrdX::Comput 235 G4cout << " G4ParameterisationTrdX::ComputeDimensions():" 231 << copyNo << G4endl; 236 << copyNo << G4endl; 232 trap.DumpInfo(); 237 trap.DumpInfo(); 233 } 238 } 234 #endif 239 #endif 235 } 240 } 236 241 237 242 238 //-------------------------------------------- 243 //-------------------------------------------------------------------------- 239 G4ParameterisationTrdY:: 244 G4ParameterisationTrdY:: 240 G4ParameterisationTrdY( EAxis axis, G4int nDiv 245 G4ParameterisationTrdY( EAxis axis, G4int nDiv, 241 G4double width, G4doub 246 G4double width, G4double offset, 242 G4VSolid* msolid, Divi 247 G4VSolid* msolid, DivisionType divType) 243 : G4VParameterisationTrd( axis, nDiv, width, 248 : G4VParameterisationTrd( axis, nDiv, width, offset, msolid, divType ) 244 { 249 { 245 CheckParametersValidity(); 250 CheckParametersValidity(); 246 SetType( "DivisionTrdY" ); 251 SetType( "DivisionTrdY" ); 247 252 248 auto msol = (G4Trd*)(fmotherSolid); << 253 G4Trd* msol = (G4Trd*)(fmotherSolid); 249 if( divType == DivWIDTH ) 254 if( divType == DivWIDTH ) 250 { 255 { 251 fnDiv = CalculateNDiv( 2*msol->GetYHalfLen 256 fnDiv = CalculateNDiv( 2*msol->GetYHalfLength1(), 252 width, offset ); 257 width, offset ); 253 } 258 } 254 else if( divType == DivNDIV ) 259 else if( divType == DivNDIV ) 255 { 260 { 256 fwidth = CalculateWidth( 2*msol->GetYHalfL 261 fwidth = CalculateWidth( 2*msol->GetYHalfLength1(), 257 nDiv, offset ); 262 nDiv, offset ); 258 } 263 } 259 264 260 #ifdef G4DIVDEBUG 265 #ifdef G4DIVDEBUG 261 if( verbose >= 1 ) 266 if( verbose >= 1 ) 262 { 267 { 263 G4cout << " G4ParameterisationTrdY no div 268 G4cout << " G4ParameterisationTrdY no divisions " << fnDiv 264 << " = " << nDiv << G4endl 269 << " = " << nDiv << G4endl 265 << " Offset " << foffset << " = " 270 << " Offset " << foffset << " = " << offset << G4endl 266 << " width " << fwidth << " = " << 271 << " width " << fwidth << " = " << width << G4endl; 267 } 272 } 268 #endif 273 #endif 269 } 274 } 270 275 271 //-------------------------------------------- 276 //------------------------------------------------------------------------ 272 G4ParameterisationTrdY::~G4ParameterisationTrd << 277 G4ParameterisationTrdY::~G4ParameterisationTrdY() >> 278 { >> 279 } 273 280 274 //-------------------------------------------- 281 //------------------------------------------------------------------------ 275 G4double G4ParameterisationTrdY::GetMaxParamet 282 G4double G4ParameterisationTrdY::GetMaxParameter() const 276 { 283 { 277 auto msol = (G4Trd*)(fmotherSolid); << 284 G4Trd* msol = (G4Trd*)(fmotherSolid); 278 return (msol->GetYHalfLength1()+msol->GetYHa 285 return (msol->GetYHalfLength1()+msol->GetYHalfLength2()); 279 } 286 } 280 287 281 //-------------------------------------------- 288 //-------------------------------------------------------------------------- 282 void 289 void 283 G4ParameterisationTrdY:: 290 G4ParameterisationTrdY:: 284 ComputeTransformation( const G4int copyNo, G4V 291 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume* physVol ) const 285 { 292 { 286 auto msol = (G4Trd*)(fmotherSolid ); << 293 G4Trd* msol = (G4Trd*)(fmotherSolid ); 287 G4double mdy = ( msol->GetYHalfLength1() + 294 G4double mdy = ( msol->GetYHalfLength1() + msol->GetYHalfLength2() ) / 2.; 288 295 289 //----- translation 296 //----- translation 290 G4ThreeVector origin(0.,0.,0.); 297 G4ThreeVector origin(0.,0.,0.); 291 G4double posi = -mdy + foffset + (copyNo+0.5 298 G4double posi = -mdy + foffset + (copyNo+0.5)*fwidth; 292 299 293 if( faxis == kYAxis ) 300 if( faxis == kYAxis ) 294 { 301 { 295 origin.setY( posi ); 302 origin.setY( posi ); 296 } 303 } 297 else 304 else 298 { 305 { 299 std::ostringstream message; 306 std::ostringstream message; 300 message << "Only axes along Y are allowed 307 message << "Only axes along Y are allowed ! Axis: " << faxis; 301 G4Exception("G4ParameterisationTrdY::Compu 308 G4Exception("G4ParameterisationTrdY::ComputeTransformation()", 302 "GeomDiv0002", FatalException, 309 "GeomDiv0002", FatalException, message); 303 } 310 } 304 311 305 #ifdef G4DIVDEBUG 312 #ifdef G4DIVDEBUG 306 if( verbose >= 2 ) 313 if( verbose >= 2 ) 307 { 314 { 308 G4cout << std::setprecision(8) 315 G4cout << std::setprecision(8) 309 << " G4ParameterisationTrdY::Comput 316 << " G4ParameterisationTrdY::ComputeTransformation " << copyNo 310 << " pos " << origin << " rot mat " 317 << " pos " << origin << " rot mat " << " axis " << faxis << G4endl; 311 } 318 } 312 #endif 319 #endif 313 320 314 //----- set translation 321 //----- set translation 315 physVol->SetTranslation( origin ); 322 physVol->SetTranslation( origin ); 316 } 323 } 317 324 318 //-------------------------------------------- 325 //-------------------------------------------------------------------------- 319 void 326 void 320 G4ParameterisationTrdY:: 327 G4ParameterisationTrdY:: 321 ComputeDimensions(G4Trd& trd, const G4int, con 328 ComputeDimensions(G4Trd& trd, const G4int, const G4VPhysicalVolume*) const 322 { 329 { 323 //---- The division along Y of a Trd will re 330 //---- The division along Y of a Trd will result a Trd, only 324 //--- if Y at -Z and +Z are equal, else use 331 //--- if Y at -Z and +Z are equal, else use the G4Trap version 325 auto msol = (G4Trd*)(fmotherSolid); << 332 G4Trd* msol = (G4Trd*)(fmotherSolid); 326 333 327 G4double pDx1 = msol->GetXHalfLength1(); 334 G4double pDx1 = msol->GetXHalfLength1(); 328 G4double pDx2 = msol->GetXHalfLength2(); 335 G4double pDx2 = msol->GetXHalfLength2(); 329 G4double pDz = msol->GetZHalfLength(); 336 G4double pDz = msol->GetZHalfLength(); 330 G4double pDy = fwidth/2. - fhgap; 337 G4double pDy = fwidth/2. - fhgap; 331 338 332 trd.SetAllParameters ( pDx1, pDx2, pDy, pDy, 339 trd.SetAllParameters ( pDx1, pDx2, pDy, pDy, pDz ); 333 340 334 #ifdef G4DIVDEBUG 341 #ifdef G4DIVDEBUG 335 if( verbose >= 2 ) 342 if( verbose >= 2 ) 336 { 343 { 337 G4cout << " G4ParameterisationTrdY::Comput 344 G4cout << " G4ParameterisationTrdY::ComputeDimensions():" << G4endl; 338 trd.DumpInfo(); 345 trd.DumpInfo(); 339 } 346 } 340 #endif 347 #endif 341 } 348 } 342 349 343 //-------------------------------------------- 350 //-------------------------------------------------------------------------- 344 void 351 void 345 G4ParameterisationTrdY::ComputeDimensions( G4T 352 G4ParameterisationTrdY::ComputeDimensions( G4Trap& trap, const G4int copyNo, 346 con 353 const G4VPhysicalVolume* ) const 347 { 354 { 348 auto msol = (G4Trd*)(fmotherSolid); << 355 G4Trd* msol = (G4Trd*)(fmotherSolid); 349 G4double pDx1 = msol->GetXHalfLength1(); 356 G4double pDx1 = msol->GetXHalfLength1(); 350 G4double pDx2 = msol->GetXHalfLength2(); 357 G4double pDx2 = msol->GetXHalfLength2(); 351 G4double pDz = msol->GetZHalfLength(); 358 G4double pDz = msol->GetZHalfLength(); 352 //fwidth is at Z=0 359 //fwidth is at Z=0 353 G4double pDy1 = msol->GetYHalfLength1(); 360 G4double pDy1 = msol->GetYHalfLength1(); 354 G4double pDy2 = msol->GetYHalfLength2(); 361 G4double pDy2 = msol->GetYHalfLength2(); 355 // G4double pDxAVE = (pDy1+pDy2)/2.; 362 // G4double pDxAVE = (pDy1+pDy2)/2.; 356 G4double yChangeRatio = (pDy2-pDy1) / (pDy2+ 363 G4double yChangeRatio = (pDy2-pDy1) / (pDy2+pDy1); 357 G4double fWidChange = yChangeRatio*fwidth; 364 G4double fWidChange = yChangeRatio*fwidth; 358 G4double fWid1 = fwidth - fWidChange; 365 G4double fWid1 = fwidth - fWidChange; 359 G4double fWid2 = fwidth + fWidChange; 366 G4double fWid2 = fwidth + fWidChange; 360 G4double fOffsetChange = yChangeRatio*foffse 367 G4double fOffsetChange = yChangeRatio*foffset/2.; 361 G4double fOffset1 = foffset - fOffsetChange; 368 G4double fOffset1 = foffset - fOffsetChange; 362 G4double fOffset2 = foffset + fOffsetChange; 369 G4double fOffset2 = foffset + fOffsetChange; 363 G4double cyx1 = -pDy1+fOffset1 + (copyNo+0.5 370 G4double cyx1 = -pDy1+fOffset1 + (copyNo+0.5)*fWid1;// centre of the side at y=-pDy1 364 G4double cyx2 = -pDy2+fOffset2 + (copyNo+0.5 371 G4double cyx2 = -pDy2+fOffset2 + (copyNo+0.5)*fWid2;// centre of the side at y=+pDy1 365 G4double alp = std::atan( (cyx2-cyx1)/(pDz*2 372 G4double alp = std::atan( (cyx2-cyx1)/(pDz*2.) ); 366 373 367 pDy1 = fwidth/2. - fWidChange/2.; 374 pDy1 = fwidth/2. - fWidChange/2.; 368 pDy2 = fwidth/2. + fWidChange/2.; 375 pDy2 = fwidth/2. + fWidChange/2.; 369 376 370 377 371 trap.SetAllParameters ( pDz, 378 trap.SetAllParameters ( pDz, 372 alp, 379 alp, 373 90*CLHEP::deg, 380 90*CLHEP::deg, 374 pDy1, 381 pDy1, 375 pDx1, 382 pDx1, 376 pDx1, 383 pDx1, 377 0., 384 0., 378 pDy2, 385 pDy2, 379 pDx2 - fhgap, 386 pDx2 - fhgap, 380 pDx2 - fhgap * pDx2/pDx1, 387 pDx2 - fhgap * pDx2/pDx1, 381 0.); 388 0.); 382 389 383 #ifdef G4DIVDEBUG 390 #ifdef G4DIVDEBUG 384 if( verbose >= 2 ) 391 if( verbose >= 2 ) 385 { 392 { 386 G4cout << " G4ParameterisationTrdY::Comput 393 G4cout << " G4ParameterisationTrdY::ComputeDimensions():" 387 << copyNo << G4endl; 394 << copyNo << G4endl; 388 trap.DumpInfo(); 395 trap.DumpInfo(); 389 } 396 } 390 #endif 397 #endif 391 } 398 } 392 399 393 //-------------------------------------------- 400 //-------------------------------------------------------------------------- 394 G4ParameterisationTrdZ:: 401 G4ParameterisationTrdZ:: 395 G4ParameterisationTrdZ( EAxis axis, G4int nDiv 402 G4ParameterisationTrdZ( EAxis axis, G4int nDiv, 396 G4double width, G4doub 403 G4double width, G4double offset, 397 G4VSolid* msolid, Divi 404 G4VSolid* msolid, DivisionType divType ) 398 : G4VParameterisationTrd( axis, nDiv, width, 405 : G4VParameterisationTrd( axis, nDiv, width, offset, msolid, divType ) 399 { 406 { 400 CheckParametersValidity(); 407 CheckParametersValidity(); 401 SetType( "DivTrdZ" ); 408 SetType( "DivTrdZ" ); 402 409 403 auto msol = (G4Trd*)(fmotherSolid); << 410 G4Trd* msol = (G4Trd*)(fmotherSolid); 404 if( divType == DivWIDTH ) 411 if( divType == DivWIDTH ) 405 { 412 { 406 fnDiv = CalculateNDiv( 2*msol->GetZHalfLen 413 fnDiv = CalculateNDiv( 2*msol->GetZHalfLength(), 407 width, offset ); 414 width, offset ); 408 } 415 } 409 else if( divType == DivNDIV ) 416 else if( divType == DivNDIV ) 410 { 417 { 411 fwidth = CalculateWidth( 2*msol->GetZHalfL 418 fwidth = CalculateWidth( 2*msol->GetZHalfLength(), 412 nDiv, offset ); 419 nDiv, offset ); 413 } 420 } 414 421 415 #ifdef G4DIVDEBUG 422 #ifdef G4DIVDEBUG 416 if( verbose >= 1 ) 423 if( verbose >= 1 ) 417 { 424 { 418 G4cout << " G4ParameterisationTrdZ no divi 425 G4cout << " G4ParameterisationTrdZ no divisions " << fnDiv 419 << " = " << nDiv << G4endl 426 << " = " << nDiv << G4endl 420 << " Offset " << foffset << " = " < 427 << " Offset " << foffset << " = " << offset << G4endl 421 << " Width " << fwidth << " = " << 428 << " Width " << fwidth << " = " << width << G4endl; 422 } 429 } 423 #endif 430 #endif 424 } 431 } 425 432 426 //-------------------------------------------- 433 //------------------------------------------------------------------------ 427 G4ParameterisationTrdZ::~G4ParameterisationTrd << 434 G4ParameterisationTrdZ::~G4ParameterisationTrdZ() >> 435 { >> 436 } 428 437 429 //-------------------------------------------- 438 //------------------------------------------------------------------------ 430 G4double G4ParameterisationTrdZ::GetMaxParamet 439 G4double G4ParameterisationTrdZ::GetMaxParameter() const 431 { 440 { 432 auto msol = (G4Trd*)(fmotherSolid); << 441 G4Trd* msol = (G4Trd*)(fmotherSolid); 433 return 2*msol->GetZHalfLength(); 442 return 2*msol->GetZHalfLength(); 434 } 443 } 435 444 436 //-------------------------------------------- 445 //-------------------------------------------------------------------------- 437 void 446 void 438 G4ParameterisationTrdZ:: 447 G4ParameterisationTrdZ:: 439 ComputeTransformation(const G4int copyNo, G4VP 448 ComputeTransformation(const G4int copyNo, G4VPhysicalVolume *physVol) const 440 { 449 { 441 auto msol = (G4Trd*)(fmotherSolid ); << 450 G4Trd* msol = (G4Trd*)(fmotherSolid ); 442 G4double mdz = msol->GetZHalfLength(); 451 G4double mdz = msol->GetZHalfLength(); 443 452 444 //----- translation 453 //----- translation 445 G4ThreeVector origin(0.,0.,0.); 454 G4ThreeVector origin(0.,0.,0.); 446 G4double posi = -mdz + OffsetZ() + (copyNo+0 455 G4double posi = -mdz + OffsetZ() + (copyNo+0.5)*fwidth; 447 if( faxis == kZAxis ) 456 if( faxis == kZAxis ) 448 { 457 { 449 origin.setZ( posi ); 458 origin.setZ( posi ); 450 } 459 } 451 else 460 else 452 { 461 { 453 std::ostringstream message; 462 std::ostringstream message; 454 message << "Only axes along Z are allowed 463 message << "Only axes along Z are allowed ! Axis: " << faxis; 455 G4Exception("G4ParameterisationTrdZ::Compu 464 G4Exception("G4ParameterisationTrdZ::ComputeTransformation()", 456 "GeomDiv0002", FatalException, 465 "GeomDiv0002", FatalException, message); 457 } 466 } 458 467 459 #ifdef G4DIVDEBUG 468 #ifdef G4DIVDEBUG 460 if( verbose >= 1 ) 469 if( verbose >= 1 ) 461 { 470 { 462 G4cout << std::setprecision(8) << " G4Para 471 G4cout << std::setprecision(8) << " G4ParameterisationTrdZ: " 463 << copyNo << G4endl 472 << copyNo << G4endl 464 << " Position: " << origin << " - O 473 << " Position: " << origin << " - Offset: " << foffset 465 << " - Width: " << fwidth << " Axis 474 << " - Width: " << fwidth << " Axis " << faxis << G4endl; 466 } 475 } 467 #endif 476 #endif 468 477 469 //----- set translation 478 //----- set translation 470 physVol->SetTranslation( origin ); 479 physVol->SetTranslation( origin ); 471 } 480 } 472 481 473 //-------------------------------------------- 482 //-------------------------------------------------------------------------- 474 void 483 void 475 G4ParameterisationTrdZ:: 484 G4ParameterisationTrdZ:: 476 ComputeDimensions(G4Trd& trd, const G4int copy 485 ComputeDimensions(G4Trd& trd, const G4int copyNo, 477 const G4VPhysicalVolume*) co 486 const G4VPhysicalVolume*) const 478 { 487 { 479 //---- The division along Z of a Trd will re 488 //---- The division along Z of a Trd will result a Trd 480 auto msol = (G4Trd*)(fmotherSolid); << 489 G4Trd* msol = (G4Trd*)(fmotherSolid); 481 490 482 G4double pDx1 = msol->GetXHalfLength1(); 491 G4double pDx1 = msol->GetXHalfLength1(); 483 G4double DDx = (msol->GetXHalfLength2() - ms 492 G4double DDx = (msol->GetXHalfLength2() - msol->GetXHalfLength1() ); 484 G4double pDy1 = msol->GetYHalfLength1(); 493 G4double pDy1 = msol->GetYHalfLength1(); 485 G4double DDy = (msol->GetYHalfLength2() - ms 494 G4double DDy = (msol->GetYHalfLength2() - msol->GetYHalfLength1() ); 486 G4double pDz = fwidth/2. - fhgap; 495 G4double pDz = fwidth/2. - fhgap; 487 G4double zLength = 2*msol->GetZHalfLength(); 496 G4double zLength = 2*msol->GetZHalfLength(); 488 497 489 trd.SetAllParameters( pDx1+DDx*(OffsetZ()+co 498 trd.SetAllParameters( pDx1+DDx*(OffsetZ()+copyNo*fwidth+fhgap)/zLength, 490 pDx1+DDx*(OffsetZ()+(c 499 pDx1+DDx*(OffsetZ()+(copyNo+1)*fwidth-fhgap)/zLength, 491 pDy1+DDy*(OffsetZ()+co 500 pDy1+DDy*(OffsetZ()+copyNo*fwidth+fhgap)/zLength, 492 pDy1+DDy*(OffsetZ()+(c 501 pDy1+DDy*(OffsetZ()+(copyNo+1)*fwidth-fhgap)/zLength, 493 pDz ); 502 pDz ); 494 503 495 #ifdef G4DIVDEBUG 504 #ifdef G4DIVDEBUG 496 if( verbose >= 1 ) 505 if( verbose >= 1 ) 497 { 506 { 498 G4cout << " G4ParameterisationTrdZ::Comput 507 G4cout << " G4ParameterisationTrdZ::ComputeDimensions()" 499 << " - Mother TRD " << G4endl; 508 << " - Mother TRD " << G4endl; 500 msol->DumpInfo(); 509 msol->DumpInfo(); 501 G4cout << " - Parameterised TRD: " 510 G4cout << " - Parameterised TRD: " 502 << copyNo << G4endl; 511 << copyNo << G4endl; 503 trd.DumpInfo(); 512 trd.DumpInfo(); 504 } 513 } 505 #endif 514 #endif 506 } 515 } 507 516