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