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 // G4ParameterisationPolycone[Rho/Phi/Z] imple << 26 // >> 27 // $Id: G4ParameterisationPolycone.cc,v 1.17 2009-05-18 19:30:29 ivana Exp $ >> 28 // GEANT4 tag $Name: not supported by cvs2svn $ >> 29 // >> 30 // class G4ParameterisationPolycone Implementation file 27 // 31 // 28 // 26.05.03 - P.Arce, Initial version 32 // 26.05.03 - P.Arce, Initial version 29 // 08.04.04 - I.Hrivnacova, Implemented reflec 33 // 08.04.04 - I.Hrivnacova, Implemented reflection 30 //-------------------------------------------- 34 //--------------------------------------------------------------------- 31 35 32 #include "G4ParameterisationPolycone.hh" 36 #include "G4ParameterisationPolycone.hh" 33 37 34 #include <iomanip> 38 #include <iomanip> 35 #include "G4ThreeVector.hh" 39 #include "G4ThreeVector.hh" 36 #include "G4RotationMatrix.hh" 40 #include "G4RotationMatrix.hh" 37 #include "G4VPhysicalVolume.hh" 41 #include "G4VPhysicalVolume.hh" 38 #include "G4LogicalVolume.hh" 42 #include "G4LogicalVolume.hh" 39 #include "G4ReflectedSolid.hh" 43 #include "G4ReflectedSolid.hh" 40 44 41 //-------------------------------------------- 45 //----------------------------------------------------------------------- 42 G4VParameterisationPolycone:: 46 G4VParameterisationPolycone:: 43 G4VParameterisationPolycone( EAxis axis, G4int 47 G4VParameterisationPolycone( EAxis axis, G4int nDiv, G4double width, 44 G4double offset, 48 G4double offset, G4VSolid* msolid, 45 DivisionType divT 49 DivisionType divType ) 46 : G4VDivisionParameterisation( axis, nDiv, 50 : G4VDivisionParameterisation( axis, nDiv, width, offset, divType, msolid ) 47 { 51 { 48 /*#ifdef G4MULTITHREADED << 52 G4Polycone* msol = (G4Polycone*)(msolid); 49 std::ostringstream message; << 53 if ((msolid->GetEntityType() != "G4ReflectedSolid") && (msol->IsGeneric())) 50 message << "Divisions for G4Polycone curren << 54 { 51 << G4endl << 55 G4String message = 52 << "Sorry! Solid: " << msolid->GetN << 56 "Sorry, generic construct for G4Polycone NOT supported.\n Solid: " 53 G4Exception("G4VParameterisationPolycone::G << 57 + msol->GetName(); 54 "GeomDiv0001", FatalException, << 58 G4Exception("G4VParameterisationPolycone::G4VParameterisationPolycone()", 55 #endif */ << 59 "NotSupported", FatalException, message); 56 auto msol = (G4Polycone*)(msolid); << 60 } 57 if (msolid->GetEntityType() == "G4ReflectedS 61 if (msolid->GetEntityType() == "G4ReflectedSolid") 58 { 62 { 59 // Get constituent solid << 63 // Get constituent solid 60 // << 61 G4VSolid* mConstituentSolid 64 G4VSolid* mConstituentSolid 62 = ((G4ReflectedSolid*)msolid)->GetConst 65 = ((G4ReflectedSolid*)msolid)->GetConstituentMovedSolid(); 63 msol = (G4Polycone*)(mConstituentSolid); 66 msol = (G4Polycone*)(mConstituentSolid); 64 67 65 // Get parameters 68 // Get parameters 66 // << 67 G4int nofZplanes = msol->GetOriginalPara 69 G4int nofZplanes = msol->GetOriginalParameters()->Num_z_planes; 68 G4double* zValues = msol->GetOriginalPara 70 G4double* zValues = msol->GetOriginalParameters()->Z_values; 69 G4double* rminValues = msol->GetOriginalP 71 G4double* rminValues = msol->GetOriginalParameters()->Rmin; 70 G4double* rmaxValues = msol->GetOriginalP 72 G4double* rmaxValues = msol->GetOriginalParameters()->Rmax; 71 73 72 // Invert z values 74 // Invert z values 73 // << 75 G4double* zValuesRefl = new double[nofZplanes]; 74 auto zValuesRefl = new G4double[nofZplanes << 76 for (G4int i=0; i<nofZplanes; i++) zValuesRefl[i] = - zValues[i]; 75 for (G4int i=0; i<nofZplanes; ++i) { zVal << 76 77 77 auto newSolid << 78 G4Polycone* newSolid 78 = new G4Polycone(msol->GetName(), 79 = new G4Polycone(msol->GetName(), 79 msol->GetStartPhi(), 80 msol->GetStartPhi(), 80 msol->GetEndPhi() - mso 81 msol->GetEndPhi() - msol->GetStartPhi(), 81 nofZplanes, zValuesRefl 82 nofZplanes, zValuesRefl, rminValues, rmaxValues); 82 83 83 delete [] zValuesRefl; 84 delete [] zValuesRefl; 84 85 85 msol = newSolid; 86 msol = newSolid; 86 fmotherSolid = newSolid; 87 fmotherSolid = newSolid; 87 fReflectedSolid = true; 88 fReflectedSolid = true; 88 fDeleteSolid = true; 89 fDeleteSolid = true; 89 } 90 } 90 } 91 } 91 92 92 //-------------------------------------------- 93 //--------------------------------------------------------------------- 93 G4VParameterisationPolycone::~G4VParameterisat << 94 G4VParameterisationPolycone::~G4VParameterisationPolycone() >> 95 { >> 96 } 94 97 95 //-------------------------------------------- 98 //--------------------------------------------------------------------- 96 G4ParameterisationPolyconeRho:: 99 G4ParameterisationPolyconeRho:: 97 G4ParameterisationPolyconeRho( EAxis axis, G4i 100 G4ParameterisationPolyconeRho( EAxis axis, G4int nDiv, 98 G4double width, 101 G4double width, G4double offset, 99 G4VSolid* msoli 102 G4VSolid* msolid, DivisionType divType ) 100 : G4VParameterisationPolycone( axis, nDiv, 103 : G4VParameterisationPolycone( axis, nDiv, width, offset, msolid, divType ) 101 { 104 { 102 CheckParametersValidity(); 105 CheckParametersValidity(); 103 SetType( "DivisionPolyconeRho" ); 106 SetType( "DivisionPolyconeRho" ); 104 107 105 auto msol = (G4Polycone*)(fmotherSolid); << 108 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 106 G4PolyconeHistorical* origparamMother = msol 109 G4PolyconeHistorical* origparamMother = msol->GetOriginalParameters(); 107 110 108 if( divType == DivWIDTH ) 111 if( divType == DivWIDTH ) 109 { 112 { 110 fnDiv = CalculateNDiv( origparamMother->Rm 113 fnDiv = CalculateNDiv( origparamMother->Rmax[0] 111 - origparamMother->Rm 114 - origparamMother->Rmin[0], width, offset ); 112 } 115 } 113 else if( divType == DivNDIV ) 116 else if( divType == DivNDIV ) 114 { 117 { 115 fwidth = CalculateWidth( origparamMother-> 118 fwidth = CalculateWidth( origparamMother->Rmax[0] 116 - origparamMother-> 119 - origparamMother->Rmin[0], nDiv, offset ); 117 } 120 } 118 121 119 #ifdef G4DIVDEBUG 122 #ifdef G4DIVDEBUG 120 if( verbose >= -1 ) 123 if( verbose >= -1 ) 121 { 124 { 122 G4cout << " G4ParameterisationPolyconeRho 125 G4cout << " G4ParameterisationPolyconeRho - # divisions " << fnDiv 123 << " = " << nDiv << G4endl 126 << " = " << nDiv << G4endl 124 << " Offset " << foffset << " = " < 127 << " Offset " << foffset << " = " << offset << G4endl 125 << " Width " << fwidth << " = " << 128 << " Width " << fwidth << " = " << width << G4endl; 126 } 129 } 127 #endif 130 #endif 128 } 131 } 129 132 130 //-------------------------------------------- 133 //--------------------------------------------------------------------- 131 G4ParameterisationPolyconeRho::~G4Parameterisa << 134 G4ParameterisationPolyconeRho::~G4ParameterisationPolyconeRho() >> 135 { >> 136 } 132 137 133 //-------------------------------------------- 138 //--------------------------------------------------------------------- 134 void G4ParameterisationPolyconeRho::CheckParam 139 void G4ParameterisationPolyconeRho::CheckParametersValidity() 135 { 140 { 136 G4VDivisionParameterisation::CheckParameters 141 G4VDivisionParameterisation::CheckParametersValidity(); 137 142 138 auto msol = (G4Polycone*)(fmotherSolid); << 143 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 139 144 140 if( fDivisionType == DivNDIVandWIDTH || fDiv 145 if( fDivisionType == DivNDIVandWIDTH || fDivisionType == DivWIDTH ) 141 { 146 { 142 std::ostringstream message; << 147 G4cerr << "WARNING - " 143 message << "In solid " << msol->GetName() << 148 << "G4ParameterisationPolyconeRho::CheckParametersValidity()" 144 << "Division along R will be done << 149 << G4endl 145 << "different for each solid secti << 150 << " Solid " << msol->GetName() << G4endl 146 << "WIDTH will not be used !"; << 151 << " Division along R will be done with a width " 147 G4Exception("G4VParameterisationPolycone:: << 152 << "different for each solid section." << G4endl 148 "GeomDiv1001", JustWarning, me << 153 << " WIDTH will not be used !" << G4endl; 149 } 154 } 150 if( foffset != 0. ) 155 if( foffset != 0. ) 151 { 156 { 152 std::ostringstream message; << 157 G4cerr << "WARNING - " 153 message << "In solid " << msol->GetName() << 158 << "G4ParameterisationPolyconeRho::CheckParametersValidity()" 154 << "Division along R will be done << 159 << G4endl 155 << "different for each solid secti << 160 << " Solid " << msol->GetName() << G4endl 156 << "OFFSET will not be used !"; << 161 << " Division along R will be done with a width " 157 G4Exception("G4VParameterisationPolycone:: << 162 << "different for each solid section." << G4endl 158 "GeomDiv1001", JustWarning, me << 163 << " OFFSET will not be used !" << G4endl; 159 } 164 } >> 165 160 } 166 } 161 167 162 //-------------------------------------------- 168 //------------------------------------------------------------------------ 163 G4double G4ParameterisationPolyconeRho::GetMax 169 G4double G4ParameterisationPolyconeRho::GetMaxParameter() const 164 { 170 { 165 auto msol = (G4Polycone*)(fmotherSolid); << 171 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 166 G4PolyconeHistorical* original_pars = msol-> 172 G4PolyconeHistorical* original_pars = msol->GetOriginalParameters(); 167 return original_pars->Rmax[0] - original_par 173 return original_pars->Rmax[0] - original_pars->Rmin[0]; 168 } 174 } 169 175 170 176 171 //-------------------------------------------- 177 //--------------------------------------------------------------------- 172 void 178 void 173 G4ParameterisationPolyconeRho:: 179 G4ParameterisationPolyconeRho:: 174 ComputeTransformation( const G4int, G4VPhysica 180 ComputeTransformation( const G4int, G4VPhysicalVolume* physVol ) const 175 { 181 { 176 //----- translation 182 //----- translation 177 G4ThreeVector origin(0.,0.,0.); 183 G4ThreeVector origin(0.,0.,0.); 178 //----- set translation 184 //----- set translation 179 physVol->SetTranslation( origin ); 185 physVol->SetTranslation( origin ); 180 186 181 //----- calculate rotation matrix: unit 187 //----- calculate rotation matrix: unit 182 188 183 #ifdef G4DIVDEBUG 189 #ifdef G4DIVDEBUG 184 if( verbose >= 2 ) 190 if( verbose >= 2 ) 185 { 191 { 186 G4cout << " G4ParameterisationPolyconeRho 192 G4cout << " G4ParameterisationPolyconeRho " << G4endl 187 << " foffset: " << foffset 193 << " foffset: " << foffset 188 << " - fwidth: " << fwidth << G4end 194 << " - fwidth: " << fwidth << G4endl; 189 } 195 } 190 #endif 196 #endif 191 197 192 ChangeRotMatrix( physVol ); 198 ChangeRotMatrix( physVol ); 193 199 194 #ifdef G4DIVDEBUG 200 #ifdef G4DIVDEBUG 195 if( verbose >= 2 ) 201 if( verbose >= 2 ) 196 { 202 { 197 G4cout << std::setprecision(8) << " G4Para 203 G4cout << std::setprecision(8) << " G4ParameterisationPolyconeRho " 198 << G4endl 204 << G4endl 199 << " Position: (0,0,0)" << 205 << " Position: " << origin/mm 200 << " - Width: " << fwidth/CLHEP::de << 206 << " - Width: " << fwidth/deg 201 << " - Axis: " << faxis << G4endl; 207 << " - Axis: " << faxis << G4endl; 202 } 208 } 203 #endif 209 #endif 204 } 210 } 205 211 206 //-------------------------------------------- 212 //--------------------------------------------------------------------- 207 void 213 void 208 G4ParameterisationPolyconeRho:: 214 G4ParameterisationPolyconeRho:: 209 ComputeDimensions( G4Polycone& pcone, const G4 215 ComputeDimensions( G4Polycone& pcone, const G4int copyNo, 210 const G4VPhysicalVolume* ) 216 const G4VPhysicalVolume* ) const 211 { 217 { 212 auto msol = (G4Polycone*)(fmotherSolid); << 218 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 213 219 214 G4PolyconeHistorical* origparamMother = msol 220 G4PolyconeHistorical* origparamMother = msol->GetOriginalParameters(); 215 G4PolyconeHistorical origparam( *origparamMo 221 G4PolyconeHistorical origparam( *origparamMother ); 216 G4int nZplanes = origparamMother->Num_z_plan 222 G4int nZplanes = origparamMother->Num_z_planes; 217 223 218 G4double width = 0.; 224 G4double width = 0.; 219 for( G4int ii = 0; ii < nZplanes; ++ii ) << 225 for( G4int ii = 0; ii < nZplanes; ii++ ) 220 { 226 { 221 width = CalculateWidth( origparamMother->R 227 width = CalculateWidth( origparamMother->Rmax[ii] 222 - origparamMother->R 228 - origparamMother->Rmin[ii], fnDiv, foffset ); 223 origparam.Rmin[ii] = origparamMother->Rmin 229 origparam.Rmin[ii] = origparamMother->Rmin[ii]+foffset+width*copyNo; 224 origparam.Rmax[ii] = origparamMother->Rmin 230 origparam.Rmax[ii] = origparamMother->Rmin[ii]+foffset+width*(copyNo+1); 225 } 231 } 226 232 227 pcone.SetOriginalParameters(&origparam); // 233 pcone.SetOriginalParameters(&origparam); // copy values & transfer pointers 228 pcone.Reset(); // 234 pcone.Reset(); // reset to new solid parameters 229 235 230 #ifdef G4DIVDEBUG 236 #ifdef G4DIVDEBUG 231 if( verbose >= -2 ) 237 if( verbose >= -2 ) 232 { 238 { 233 G4cout << "G4ParameterisationPolyconeRho:: 239 G4cout << "G4ParameterisationPolyconeRho::ComputeDimensions()" << G4endl 234 << "-- Parametrised pcone copy-numb 240 << "-- Parametrised pcone copy-number: " << copyNo << G4endl; 235 pcone.DumpInfo(); 241 pcone.DumpInfo(); 236 } 242 } 237 #endif 243 #endif 238 } 244 } 239 245 240 //-------------------------------------------- 246 //--------------------------------------------------------------------- 241 G4ParameterisationPolyconePhi:: 247 G4ParameterisationPolyconePhi:: 242 G4ParameterisationPolyconePhi( EAxis axis, G4i 248 G4ParameterisationPolyconePhi( EAxis axis, G4int nDiv, 243 G4double width, 249 G4double width, G4double offset, 244 G4VSolid* msoli 250 G4VSolid* msolid, DivisionType divType ) 245 : G4VParameterisationPolycone( axis, nDiv, 251 : G4VParameterisationPolycone( axis, nDiv, width, offset, msolid, divType ) 246 { 252 { 247 CheckParametersValidity(); 253 CheckParametersValidity(); 248 SetType( "DivisionPolyconePhi" ); 254 SetType( "DivisionPolyconePhi" ); 249 255 250 auto msol = (G4Polycone*)(fmotherSolid); << 256 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 251 G4double deltaPhi = msol->GetEndPhi() - msol 257 G4double deltaPhi = msol->GetEndPhi() - msol->GetStartPhi(); 252 258 253 if( divType == DivWIDTH ) 259 if( divType == DivWIDTH ) 254 { 260 { 255 fnDiv = CalculateNDiv( deltaPhi, width, of 261 fnDiv = CalculateNDiv( deltaPhi, width, offset ); 256 } 262 } 257 else if( divType == DivNDIV ) 263 else if( divType == DivNDIV ) 258 { 264 { 259 fwidth = CalculateWidth( deltaPhi, nDiv, o 265 fwidth = CalculateWidth( deltaPhi, nDiv, offset ); 260 } 266 } 261 267 262 #ifdef G4DIVDEBUG 268 #ifdef G4DIVDEBUG 263 if( verbose >= 1 ) 269 if( verbose >= 1 ) 264 { 270 { 265 G4cout << " G4ParameterisationPolyconePhi 271 G4cout << " G4ParameterisationPolyconePhi - # divisions " << fnDiv 266 << " = " << nDiv << G4endl 272 << " = " << nDiv << G4endl 267 << " Offset " << foffset/CLHEP::deg << 273 << " Offset " << foffset/deg << " = " << offset/deg << G4endl 268 << " Width " << fwidth/CLHEP::deg < << 274 << " Width " << fwidth/deg << " = " << width/deg << G4endl; 269 } 275 } 270 #endif 276 #endif 271 } 277 } 272 278 273 //-------------------------------------------- 279 //--------------------------------------------------------------------- 274 G4ParameterisationPolyconePhi::~G4Parameterisa << 280 G4ParameterisationPolyconePhi::~G4ParameterisationPolyconePhi() >> 281 { >> 282 } 275 283 276 //-------------------------------------------- 284 //------------------------------------------------------------------------ 277 G4double G4ParameterisationPolyconePhi::GetMax 285 G4double G4ParameterisationPolyconePhi::GetMaxParameter() const 278 { 286 { 279 auto msol = (G4Polycone*)(fmotherSolid); << 287 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 280 return msol->GetEndPhi() - msol->GetStartPhi 288 return msol->GetEndPhi() - msol->GetStartPhi(); 281 } 289 } 282 290 283 //-------------------------------------------- 291 //--------------------------------------------------------------------- 284 void 292 void 285 G4ParameterisationPolyconePhi:: 293 G4ParameterisationPolyconePhi:: 286 ComputeTransformation( const G4int copyNo, G4V << 294 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume *physVol ) const 287 { 295 { 288 //----- translation 296 //----- translation 289 G4ThreeVector origin(0.,0.,0.); 297 G4ThreeVector origin(0.,0.,0.); 290 //----- set translation 298 //----- set translation 291 physVol->SetTranslation( origin ); 299 physVol->SetTranslation( origin ); 292 300 293 //----- calculate rotation matrix (so that a 301 //----- calculate rotation matrix (so that all volumes point to the centre) 294 G4double posi = foffset + copyNo*fwidth; 302 G4double posi = foffset + copyNo*fwidth; 295 303 296 #ifdef G4DIVDEBUG 304 #ifdef G4DIVDEBUG 297 if( verbose >= 2 ) 305 if( verbose >= 2 ) 298 { 306 { 299 G4cout << " G4ParameterisationPolyconePhi << 307 G4cout << " G4ParameterisationPolyconePhi - position: " << posi/deg 300 << G4endl 308 << G4endl 301 << " copyNo: " << copyNo << " - fof << 309 << " copyNo: " << copyNo << " - foffset: " << foffset/deg 302 << " - fwidth: " << fwidth/CLHEP::d << 310 << " - fwidth: " << fwidth/deg << G4endl; 303 } 311 } 304 #endif 312 #endif 305 313 306 ChangeRotMatrix( physVol, -posi ); 314 ChangeRotMatrix( physVol, -posi ); 307 315 308 #ifdef G4DIVDEBUG 316 #ifdef G4DIVDEBUG 309 if( verbose >= 2 ) 317 if( verbose >= 2 ) 310 { 318 { 311 G4cout << std::setprecision(8) << " G4Para 319 G4cout << std::setprecision(8) << " G4ParameterisationPolyconePhi " 312 << copyNo << G4endl 320 << copyNo << G4endl 313 << " Position: (0,0,0) - Width: " << fwid << 321 << " Position: " << origin << " - Width: " << fwidth 314 << " - Axis: " << faxis << G4endl; 322 << " - Axis: " << faxis << G4endl; 315 } 323 } 316 #endif 324 #endif 317 } 325 } 318 326 319 //-------------------------------------------- 327 //--------------------------------------------------------------------- 320 void 328 void 321 G4ParameterisationPolyconePhi:: 329 G4ParameterisationPolyconePhi:: 322 ComputeDimensions( G4Polycone& pcone, const G4 330 ComputeDimensions( G4Polycone& pcone, const G4int, 323 const G4VPhysicalVolume* ) 331 const G4VPhysicalVolume* ) const 324 { 332 { 325 auto msol = (G4Polycone*)(fmotherSolid); << 333 G4Polycone* msol = (G4Polycone*)(fmotherSolid); 326 334 327 G4PolyconeHistorical* origparamMother = msol 335 G4PolyconeHistorical* origparamMother = msol->GetOriginalParameters(); 328 G4PolyconeHistorical origparam( *origparamMo 336 G4PolyconeHistorical origparam( *origparamMother ); 329 origparam.Start_angle = origparamMother->Sta 337 origparam.Start_angle = origparamMother->Start_angle; 330 origparam.Opening_angle = fwidth; 338 origparam.Opening_angle = fwidth; 331 339 332 pcone.SetOriginalParameters(&origparam); // 340 pcone.SetOriginalParameters(&origparam); // copy values & transfer pointers 333 pcone.Reset(); // 341 pcone.Reset(); // reset to new solid parameters 334 342 335 #ifdef G4DIVDEBUG 343 #ifdef G4DIVDEBUG 336 if( verbose >= 2 ) 344 if( verbose >= 2 ) 337 { 345 { 338 G4cout << "G4ParameterisationPolyconePhi:: 346 G4cout << "G4ParameterisationPolyconePhi::ComputeDimensions():" << G4endl; 339 pcone.DumpInfo(); 347 pcone.DumpInfo(); 340 } 348 } 341 #endif 349 #endif 342 } 350 } 343 351 344 //-------------------------------------------- 352 //--------------------------------------------------------------------- 345 G4ParameterisationPolyconeZ:: 353 G4ParameterisationPolyconeZ:: 346 G4ParameterisationPolyconeZ( EAxis axis, G4int 354 G4ParameterisationPolyconeZ( EAxis axis, G4int nDiv, 347 G4double width, G 355 G4double width, G4double offset, 348 G4VSolid* msolid, 356 G4VSolid* msolid, DivisionType divType) 349 : G4VParameterisationPolycone( axis, nDiv, w 357 : G4VParameterisationPolycone( axis, nDiv, width, offset, msolid, divType ), >> 358 fNSegment(0), 350 fOrigParamMother(((G4Polycone*)fmotherSoli 359 fOrigParamMother(((G4Polycone*)fmotherSolid)->GetOriginalParameters()) 351 { 360 { 352 361 353 CheckParametersValidity(); 362 CheckParametersValidity(); 354 SetType( "DivisionPolyconeZ" ); 363 SetType( "DivisionPolyconeZ" ); 355 364 356 if( divType == DivWIDTH ) 365 if( divType == DivWIDTH ) 357 { 366 { 358 fnDiv = 367 fnDiv = 359 CalculateNDiv( fOrigParamMother->Z_value 368 CalculateNDiv( fOrigParamMother->Z_values[fOrigParamMother->Num_z_planes-1] 360 - fOrigParamMother->Z_val 369 - fOrigParamMother->Z_values[0] , width, offset ); 361 } 370 } 362 else if( divType == DivNDIV ) 371 else if( divType == DivNDIV ) 363 { 372 { 364 fwidth = 373 fwidth = 365 CalculateNDiv( fOrigParamMother->Z_value 374 CalculateNDiv( fOrigParamMother->Z_values[fOrigParamMother->Num_z_planes-1] 366 - fOrigParamMother->Z_val 375 - fOrigParamMother->Z_values[0] , nDiv, offset ); 367 } 376 } 368 377 369 #ifdef G4DIVDEBUG 378 #ifdef G4DIVDEBUG 370 if( verbose >= 1 ) 379 if( verbose >= 1 ) 371 { 380 { 372 G4cout << " G4ParameterisationPolyconeZ - 381 G4cout << " G4ParameterisationPolyconeZ - # divisions " << fnDiv << " = " 373 << nDiv << G4endl 382 << nDiv << G4endl 374 << " Offset " << foffset << " = " < 383 << " Offset " << foffset << " = " << offset << G4endl 375 << " Width " << fwidth << " = " << 384 << " Width " << fwidth << " = " << width << G4endl; 376 } 385 } 377 #endif 386 #endif 378 } 387 } 379 388 380 //-------------------------------------------- 389 //--------------------------------------------------------------------- 381 G4ParameterisationPolyconeZ::~G4Parameterisati << 390 G4ParameterisationPolyconeZ::~G4ParameterisationPolyconeZ() >> 391 { >> 392 } 382 393 383 //-------------------------------------------- 394 //------------------------------------------------------------------------ 384 G4double G4ParameterisationPolyconeZ::GetR(G4d 395 G4double G4ParameterisationPolyconeZ::GetR(G4double z, 385 G4d 396 G4double z1, G4double r1, 386 G4d 397 G4double z2, G4double r2) const 387 { 398 { 388 // Linear parameterisation: 399 // Linear parameterisation: 389 // r = az + b 400 // r = az + b 390 // a = (r1 - r2)/(z1-z2) 401 // a = (r1 - r2)/(z1-z2) 391 // b = r1 - a*z1 402 // b = r1 - a*z1 392 403 393 return (r1-r2)/(z1-z2)*z + ( r1 - (r1-r2)/(z 404 return (r1-r2)/(z1-z2)*z + ( r1 - (r1-r2)/(z1-z2)*z1 ) ; 394 } 405 } 395 406 396 //-------------------------------------------- 407 //------------------------------------------------------------------------ 397 G4double G4ParameterisationPolyconeZ::GetRmin( 408 G4double G4ParameterisationPolyconeZ::GetRmin(G4double z, G4int nseg) const 398 { 409 { 399 // Get Rmin in the given z position for the gi 410 // Get Rmin in the given z position for the given polycone segment 400 411 401 return GetR(z, 412 return GetR(z, 402 fOrigParamMother->Z_values[nseg] 413 fOrigParamMother->Z_values[nseg], 403 fOrigParamMother->Rmin[nseg], 414 fOrigParamMother->Rmin[nseg], 404 fOrigParamMother->Z_values[nseg+ 415 fOrigParamMother->Z_values[nseg+1], 405 fOrigParamMother->Rmin[nseg+1]); 416 fOrigParamMother->Rmin[nseg+1]); 406 } 417 } 407 418 408 //-------------------------------------------- 419 //------------------------------------------------------------------------ 409 G4double G4ParameterisationPolyconeZ::GetRmax( 420 G4double G4ParameterisationPolyconeZ::GetRmax(G4double z, G4int nseg) const 410 { 421 { 411 // Get Rmax in the given z position for the gi 422 // Get Rmax in the given z position for the given polycone segment 412 423 413 return GetR(z, 424 return GetR(z, 414 fOrigParamMother->Z_values[nseg] 425 fOrigParamMother->Z_values[nseg], 415 fOrigParamMother->Rmax[nseg], 426 fOrigParamMother->Rmax[nseg], 416 fOrigParamMother->Z_values[nseg+ 427 fOrigParamMother->Z_values[nseg+1], 417 fOrigParamMother->Rmax[nseg+1]); 428 fOrigParamMother->Rmax[nseg+1]); 418 } 429 } 419 430 420 //-------------------------------------------- 431 //------------------------------------------------------------------------ 421 G4double G4ParameterisationPolyconeZ::GetMaxPa 432 G4double G4ParameterisationPolyconeZ::GetMaxParameter() const 422 { 433 { 423 return std::abs (fOrigParamMother->Z_values[ 434 return std::abs (fOrigParamMother->Z_values[fOrigParamMother->Num_z_planes-1] 424 -fOrigParamMother->Z_values[ << 435 -fOrigParamMother->Z_values[0]); 425 } 436 } 426 437 427 //-------------------------------------------- 438 //--------------------------------------------------------------------- 428 void G4ParameterisationPolyconeZ::CheckParamet 439 void G4ParameterisationPolyconeZ::CheckParametersValidity() 429 { 440 { 430 G4VDivisionParameterisation::CheckParameters 441 G4VDivisionParameterisation::CheckParametersValidity(); 431 442 432 // Division will be following the mother pol 443 // Division will be following the mother polycone segments 433 // << 444 if( fDivisionType == DivNDIV ) { 434 if( fDivisionType == DivNDIV ) << 445 if( fnDiv > fOrigParamMother->Num_z_planes-1 ) { 435 { << 446 G4cerr << "ERROR - " 436 if( fnDiv > fOrigParamMother->Num_z_planes << 447 << "G4ParameterisationPolyconeZ::CheckParametersValidity()" 437 { << 448 << G4endl 438 std::ostringstream error; << 449 << " Division along Z will be done splitting in the defined" 439 error << "Configuration not supported." << 440 << "Division along Z will be done << 441 << G4endl 450 << G4endl 442 << "Z planes, i.e, the number of << 451 << " z_planes, i.e, the number of division would be :" 443 << fOrigParamMother->Num_z_planes << 452 << " " << fOrigParamMother->Num_z_planes-1 444 << ", instead of: " << fnDiv << " << 453 << " instead of " << fnDiv << " !" >> 454 << G4endl; 445 G4Exception("G4ParameterisationPolyconeZ 455 G4Exception("G4ParameterisationPolyconeZ::CheckParametersValidity()", 446 "GeomDiv0001", FatalExceptio << 456 "IllegalConstruct", FatalException, >> 457 "Not supported configuration."); 447 } 458 } 448 } 459 } 449 460 450 // Division will be done within one polycone 461 // Division will be done within one polycone segment 451 // with applying given width and offset 462 // with applying given width and offset 452 // << 463 if( fDivisionType == DivNDIVandWIDTH || fDivisionType == DivWIDTH ) { 453 if( fDivisionType == DivNDIVandWIDTH || fDiv << 454 { << 455 // Check if divided region does not span o 464 // Check if divided region does not span over more 456 // than one z segment 465 // than one z segment 457 466 458 G4int isegstart = -1; // number of the se 467 G4int isegstart = -1; // number of the segment containing start position 459 G4int isegend = -1; // number of the se 468 G4int isegend = -1; // number of the segment containing end position 460 469 461 if ( !fReflectedSolid ) << 470 if ( ! fReflectedSolid ) { 462 { << 463 // The start/end position of the divided 471 // The start/end position of the divided region 464 // << 465 G4double zstart 472 G4double zstart 466 = fOrigParamMother->Z_values[0] + foff 473 = fOrigParamMother->Z_values[0] + foffset; 467 G4double zend 474 G4double zend 468 = fOrigParamMother->Z_values[0] + foff 475 = fOrigParamMother->Z_values[0] + foffset + fnDiv* fwidth; 469 476 470 G4int counter = 0; 477 G4int counter = 0; 471 while ( isegend < 0 && counter < fOrigPa << 478 while ( isegend < 0 && counter < fOrigParamMother->Num_z_planes-1 ) { 472 { << 473 // first segment 479 // first segment 474 if ( zstart >= fOrigParamMother->Z_val 480 if ( zstart >= fOrigParamMother->Z_values[counter] && 475 zstart < fOrigParamMother->Z_val << 481 zstart < fOrigParamMother->Z_values[counter+1] ) { 476 { << 477 isegstart = counter; 482 isegstart = counter; 478 } 483 } 479 // last segment 484 // last segment 480 if ( zend > fOrigParamMother->Z_value 485 if ( zend > fOrigParamMother->Z_values[counter] && 481 zend <= fOrigParamMother->Z_value << 486 zend <= fOrigParamMother->Z_values[counter+1] ) { 482 { << 483 isegend = counter; 487 isegend = counter; 484 } 488 } 485 ++counter; 489 ++counter; 486 } // Loop checking, 06.08.2015, G.Cosmo << 490 } 487 } 491 } 488 else << 492 else { 489 { << 490 // The start/end position of the divided 493 // The start/end position of the divided region 491 // << 492 G4double zstart 494 G4double zstart 493 = fOrigParamMother->Z_values[0] - foff 495 = fOrigParamMother->Z_values[0] - foffset; 494 G4double zend 496 G4double zend 495 = fOrigParamMother->Z_values[0] - ( fo 497 = fOrigParamMother->Z_values[0] - ( foffset + fnDiv* fwidth); 496 498 497 G4int counter = 0; 499 G4int counter = 0; 498 while ( isegend < 0 && counter < fOrigPa << 500 while ( isegend < 0 && counter < fOrigParamMother->Num_z_planes-1 ) { 499 { << 500 // first segment 501 // first segment 501 if ( zstart <= fOrigParamMother->Z_val 502 if ( zstart <= fOrigParamMother->Z_values[counter] && 502 zstart > fOrigParamMother->Z_val << 503 zstart > fOrigParamMother->Z_values[counter+1] ) { 503 { << 504 isegstart = counter; 504 isegstart = counter; 505 } 505 } 506 // last segment 506 // last segment 507 if ( zend < fOrigParamMother->Z_value 507 if ( zend < fOrigParamMother->Z_values[counter] && 508 zend >= fOrigParamMother->Z_value << 508 zend >= fOrigParamMother->Z_values[counter+1] ) { 509 { << 510 isegend = counter; 509 isegend = counter; 511 } 510 } 512 ++counter; 511 ++counter; 513 } // Loop checking, 06.08.2015, G.Cosmo << 512 } 514 } 513 } 515 514 516 515 517 if ( isegstart != isegend ) << 516 if ( isegstart != isegend ) { 518 { << 517 G4cerr << "WARNING - " 519 std::ostringstream message; << 518 << "G4ParameterisationPolyconeZ::CheckParametersValidity()" 520 message << "Condiguration not supported. << 519 << G4endl 521 << "Division with user defined w << 520 << " Division with user defined width." << G4endl 522 << "Solid " << fmotherSolid->Get << 521 << " Solid " << fmotherSolid->GetName() << G4endl 523 << "Divided region is not betwee << 522 << " Divided region is not between two z planes." >> 523 << G4endl; >> 524 524 G4Exception("G4ParameterisationPolyconeZ 525 G4Exception("G4ParameterisationPolyconeZ::CheckParametersValidity()", 525 "GeomDiv0001", FatalExceptio << 526 "IllegalConstruct", FatalException, >> 527 "Not supported configuration."); 526 } 528 } 527 529 528 fNSegment = isegstart; 530 fNSegment = isegstart; 529 } 531 } 530 } 532 } 531 533 532 //-------------------------------------------- 534 //--------------------------------------------------------------------- 533 void 535 void 534 G4ParameterisationPolyconeZ:: 536 G4ParameterisationPolyconeZ:: 535 ComputeTransformation( const G4int copyNo, G4V 537 ComputeTransformation( const G4int copyNo, G4VPhysicalVolume* physVol) const 536 { 538 { 537 G4double posi = 0.; << 539 if ( fDivisionType == DivNDIV ) { 538 if ( fDivisionType == DivNDIV ) << 539 { << 540 // The position of the centre of copyNo-th 540 // The position of the centre of copyNo-th mother polycone segment 541 // << 541 G4double posi 542 posi = ( fOrigParamMother->Z_values[copyNo << 542 = ( fOrigParamMother->Z_values[copyNo] 543 + fOrigParamMother->Z_valu << 543 + fOrigParamMother->Z_values[copyNo+1])/2; 544 physVol->SetTranslation( G4ThreeVector(0, 544 physVol->SetTranslation( G4ThreeVector(0, 0, posi) ); 545 } 545 } 546 546 547 if ( fDivisionType == DivNDIVandWIDTH || fDi << 547 if ( fDivisionType == DivNDIVandWIDTH || fDivisionType == DivWIDTH ) { 548 { << 549 // The position of the centre of copyNo-th 548 // The position of the centre of copyNo-th division 550 // << 549 G4double posi = fOrigParamMother->Z_values[0]; 551 posi = fOrigParamMother->Z_values[0]; << 552 550 553 if ( !fReflectedSolid ) << 551 if ( ! fReflectedSolid ) 554 posi += foffset + (2*copyNo + 1) * fwidt 552 posi += foffset + (2*copyNo + 1) * fwidth/2.; 555 else 553 else 556 posi -= foffset + (2*copyNo + 1) * fwidt 554 posi -= foffset + (2*copyNo + 1) * fwidth/2.; 557 555 558 physVol->SetTranslation( G4ThreeVector(0, 556 physVol->SetTranslation( G4ThreeVector(0, 0, posi) ); 559 } 557 } 560 558 561 //----- calculate rotation matrix: unit 559 //----- calculate rotation matrix: unit 562 560 563 #ifdef G4DIVDEBUG 561 #ifdef G4DIVDEBUG 564 if( verbose >= 2 ) 562 if( verbose >= 2 ) 565 { 563 { 566 G4cout << " G4ParameterisationPolyconeZ - 564 G4cout << " G4ParameterisationPolyconeZ - position: " << posi << G4endl 567 << " copyNo: " << copyNo << " - fof << 565 << " copyNo: " << copyNo << " - foffset: " << foffset/deg 568 << " - fwidth: " << fwidth/CLHEP::d << 566 << " - fwidth: " << fwidth/deg << G4endl; 569 } 567 } 570 #endif 568 #endif 571 569 572 ChangeRotMatrix( physVol ); 570 ChangeRotMatrix( physVol ); 573 571 574 #ifdef G4DIVDEBUG 572 #ifdef G4DIVDEBUG 575 if( verbose >= 2 ) 573 if( verbose >= 2 ) 576 { 574 { 577 G4cout << std::setprecision(8) << " G4Para 575 G4cout << std::setprecision(8) << " G4ParameterisationPolyconeZ " 578 << copyNo << G4endl 576 << copyNo << G4endl 579 << " Position: (0,0,0) - Width: " < << 577 << " Position: " << origin << " - Width: " << fwidth 580 << " - Axis: " << faxis << G4endl; 578 << " - Axis: " << faxis << G4endl; 581 } 579 } 582 #endif 580 #endif 583 } 581 } 584 582 585 //-------------------------------------------- 583 //--------------------------------------------------------------------- 586 void 584 void 587 G4ParameterisationPolyconeZ:: 585 G4ParameterisationPolyconeZ:: 588 ComputeDimensions( G4Polycone& pcone, const G4 586 ComputeDimensions( G4Polycone& pcone, const G4int copyNo, 589 const G4VPhysicalVolume* ) 587 const G4VPhysicalVolume* ) const 590 { 588 { 591 589 592 // Define division solid 590 // Define division solid 593 // << 594 G4PolyconeHistorical origparam; 591 G4PolyconeHistorical origparam; 595 G4int nz = 2; 592 G4int nz = 2; 596 origparam.Num_z_planes = nz; 593 origparam.Num_z_planes = nz; 597 origparam.Start_angle = fOrigParamMother->St 594 origparam.Start_angle = fOrigParamMother->Start_angle; 598 origparam.Opening_angle = fOrigParamMother-> 595 origparam.Opening_angle = fOrigParamMother->Opening_angle; 599 596 600 // Define division solid z sections 597 // Define division solid z sections 601 // << 602 origparam.Z_values = new G4double[nz]; 598 origparam.Z_values = new G4double[nz]; 603 origparam.Rmin = new G4double[nz]; 599 origparam.Rmin = new G4double[nz]; 604 origparam.Rmax = new G4double[nz]; 600 origparam.Rmax = new G4double[nz]; 605 601 606 if ( fDivisionType == DivNDIV ) << 602 if ( fDivisionType == DivNDIV ) { 607 { << 608 // The position of the centre of copyNo-th 603 // The position of the centre of copyNo-th mother polycone segment 609 G4double posi = (fOrigParamMother->Z_value 604 G4double posi = (fOrigParamMother->Z_values[copyNo] 610 + fOrigParamMother->Z_value 605 + fOrigParamMother->Z_values[copyNo+1])/2; 611 606 612 origparam.Z_values[0] = fOrigParamMother-> 607 origparam.Z_values[0] = fOrigParamMother->Z_values[copyNo] - posi; 613 origparam.Z_values[1] = fOrigParamMother-> 608 origparam.Z_values[1] = fOrigParamMother->Z_values[copyNo+1] - posi; 614 origparam.Rmin[0] = fOrigParamMother->Rmin 609 origparam.Rmin[0] = fOrigParamMother->Rmin[copyNo]; 615 origparam.Rmin[1] = fOrigParamMother->Rmin 610 origparam.Rmin[1] = fOrigParamMother->Rmin[copyNo+1]; 616 origparam.Rmax[0] = fOrigParamMother->Rmax 611 origparam.Rmax[0] = fOrigParamMother->Rmax[copyNo]; 617 origparam.Rmax[1] = fOrigParamMother->Rmax 612 origparam.Rmax[1] = fOrigParamMother->Rmax[copyNo+1]; 618 } 613 } 619 614 620 if ( fDivisionType == DivNDIVandWIDTH || fDi << 615 if ( fDivisionType == DivNDIVandWIDTH || fDivisionType == DivWIDTH ) { 621 { << 616 if ( ! fReflectedSolid ) { 622 if ( !fReflectedSolid ) << 623 { << 624 origparam.Z_values[0] = - fwidth/2.; 617 origparam.Z_values[0] = - fwidth/2.; 625 origparam.Z_values[1] = fwidth/2.; 618 origparam.Z_values[1] = fwidth/2.; 626 619 627 // The position of the centre of copyNo- 620 // The position of the centre of copyNo-th division 628 // << 621 G4double posi 629 G4double posi = fOrigParamMother->Z_valu << 622 = fOrigParamMother->Z_values[0] + foffset + (2*copyNo + 1) * fwidth/2.; 630 + foffset + (2*copyNo + 1) << 631 623 632 // The first and last z sides z values 624 // The first and last z sides z values 633 // << 634 G4double zstart = posi - fwidth/2.; 625 G4double zstart = posi - fwidth/2.; 635 G4double zend = posi + fwidth/2.; 626 G4double zend = posi + fwidth/2.; 636 origparam.Rmin[0] = GetRmin(zstart, fNSe 627 origparam.Rmin[0] = GetRmin(zstart, fNSegment); 637 origparam.Rmax[0] = GetRmax(zstart, fNSe 628 origparam.Rmax[0] = GetRmax(zstart, fNSegment); 638 origparam.Rmin[1] = GetRmin(zend, fNSegm 629 origparam.Rmin[1] = GetRmin(zend, fNSegment); 639 origparam.Rmax[1] = GetRmax(zend, fNSegm 630 origparam.Rmax[1] = GetRmax(zend, fNSegment); 640 } 631 } 641 else << 632 else { 642 { << 643 origparam.Z_values[0] = fwidth/2.; 633 origparam.Z_values[0] = fwidth/2.; 644 origparam.Z_values[1] = - fwidth/2.; 634 origparam.Z_values[1] = - fwidth/2.; 645 635 646 // The position of the centre of copyNo- 636 // The position of the centre of copyNo-th division 647 // << 637 G4double posi 648 G4double posi = fOrigParamMother->Z_valu << 638 = fOrigParamMother->Z_values[0] - ( foffset + (2*copyNo + 1) * fwidth/2.); 649 - ( foffset + (2*copyNo + << 650 639 651 // The first and last z sides z values 640 // The first and last z sides z values 652 // << 653 G4double zstart = posi + fwidth/2.; 641 G4double zstart = posi + fwidth/2.; 654 G4double zend = posi - fwidth/2.; 642 G4double zend = posi - fwidth/2.; 655 origparam.Rmin[0] = GetRmin(zstart, fNSe 643 origparam.Rmin[0] = GetRmin(zstart, fNSegment); 656 origparam.Rmax[0] = GetRmax(zstart, fNSe 644 origparam.Rmax[0] = GetRmax(zstart, fNSegment); 657 origparam.Rmin[1] = GetRmin(zend, fNSegm 645 origparam.Rmin[1] = GetRmin(zend, fNSegment); 658 origparam.Rmax[1] = GetRmax(zend, fNSegm 646 origparam.Rmax[1] = GetRmax(zend, fNSegment); 659 } 647 } 660 648 661 // It can happen due to rounding errors 649 // It can happen due to rounding errors 662 // << 663 if ( origparam.Rmin[0] < 0.0 ) origpara 650 if ( origparam.Rmin[0] < 0.0 ) origparam.Rmin[0] = 0.0; 664 if ( origparam.Rmin[nz-1] < 0.0 ) origpara 651 if ( origparam.Rmin[nz-1] < 0.0 ) origparam.Rmin[1] = 0.0; 665 } 652 } 666 653 667 pcone.SetOriginalParameters(&origparam); // 654 pcone.SetOriginalParameters(&origparam); // copy values & transfer pointers 668 pcone.Reset(); // 655 pcone.Reset(); // reset to new solid parameters 669 656 670 #ifdef G4DIVDEBUG 657 #ifdef G4DIVDEBUG 671 if( verbose >= 2 ) 658 if( verbose >= 2 ) 672 { 659 { 673 G4cout << "G4ParameterisationPolyconeZ::Co 660 G4cout << "G4ParameterisationPolyconeZ::ComputeDimensions()" << G4endl 674 << "-- Parametrised pcone copy-numb 661 << "-- Parametrised pcone copy-number: " << copyNo << G4endl; 675 pcone.DumpInfo(); 662 pcone.DumpInfo(); 676 } 663 } 677 #endif 664 #endif 678 } 665 } 679 666