1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G3Division.cc,v 1.17 2006/06/29 18:12:53 gunter Exp $ |
---|
28 | // GEANT4 tag $Name: $ |
---|
29 | // |
---|
30 | // by I.Hrivnacova, V.Berejnoi 13.10.99 |
---|
31 | |
---|
32 | #include "G3Division.hh" |
---|
33 | #include "G3VolTableEntry.hh" |
---|
34 | #include "G3toG4MakeSolid.hh" |
---|
35 | #include "G4Para.hh" |
---|
36 | #include "G3Pos.hh" |
---|
37 | #include "G4LogicalVolume.hh" |
---|
38 | #include "G4VPhysicalVolume.hh" |
---|
39 | #include "G4PVPlacement.hh" |
---|
40 | #include "G4PVReplica.hh" |
---|
41 | #ifndef G3G4_NO_REFLECTION |
---|
42 | #include "G4ReflectionFactory.hh" |
---|
43 | #endif |
---|
44 | |
---|
45 | G3VolTableEntry* G4CreateVTE(G4String vname, G4String shape, G4int nmed, |
---|
46 | G4double Rpar[], G4int npar); |
---|
47 | |
---|
48 | G3Division::G3Division(G3DivType type, G3VolTableEntry* vte, |
---|
49 | G3VolTableEntry* mvte, G4int nofDivisions, |
---|
50 | G4int iaxis, G4int nmed, G4double c0, G4double step) |
---|
51 | : fType(type), |
---|
52 | fVTE(vte), |
---|
53 | fMVTE(mvte), |
---|
54 | fNofDivisions(nofDivisions), |
---|
55 | fIAxis(iaxis), |
---|
56 | fNmed(nmed), |
---|
57 | fC0(c0), |
---|
58 | fStep(step), |
---|
59 | fLowRange(0.), |
---|
60 | fHighRange(0.), |
---|
61 | fWidth(0.), |
---|
62 | fOffset(0.), |
---|
63 | fAxis(kXAxis) |
---|
64 | { |
---|
65 | fVTE->SetHasNegPars(true); |
---|
66 | } |
---|
67 | |
---|
68 | G3Division::G3Division(G3VolTableEntry* vte, G3VolTableEntry* mvte, |
---|
69 | const G3Division& division) |
---|
70 | : fVTE(vte), |
---|
71 | fMVTE(mvte) |
---|
72 | { |
---|
73 | // only "input" parameters are copied from division |
---|
74 | fType = division.fType; |
---|
75 | fNofDivisions = division.fNofDivisions; |
---|
76 | fIAxis = division.fIAxis; |
---|
77 | fNmed = division.fNmed; |
---|
78 | fC0 = division.fC0; |
---|
79 | fStep = division.fStep; |
---|
80 | |
---|
81 | // other parameters are set as in standard constructor |
---|
82 | fLowRange = 0.; |
---|
83 | fHighRange = 0.; |
---|
84 | fWidth = 0.; |
---|
85 | fOffset = 0.; |
---|
86 | fAxis = kXAxis; |
---|
87 | fVTE->SetHasNegPars(true); |
---|
88 | } |
---|
89 | |
---|
90 | G3Division::~G3Division() |
---|
91 | {} |
---|
92 | |
---|
93 | // public methods |
---|
94 | |
---|
95 | void G3Division::UpdateVTE() |
---|
96 | { |
---|
97 | if (fVTE->HasNegPars() && !(fMVTE->HasNegPars())) { |
---|
98 | |
---|
99 | // set nmed from mother |
---|
100 | if (fNmed == 0) fNmed = fMVTE->GetNmed(); |
---|
101 | fVTE->SetNmed(fNmed); |
---|
102 | |
---|
103 | SetRangeAndAxis(); |
---|
104 | |
---|
105 | // create envelope (if necessary) |
---|
106 | // and solid |
---|
107 | G3VolTableEntry* envVTE = 0; |
---|
108 | if (fType == kDvn) envVTE = Dvn(); |
---|
109 | else if (fType == kDvn2) envVTE = Dvn2(); |
---|
110 | else if (fType == kDvt) envVTE = Dvt(); |
---|
111 | else if (fType == kDvt2) envVTE = Dvt2(); |
---|
112 | |
---|
113 | if (envVTE) { |
---|
114 | // reset mother <-> daughter |
---|
115 | fMVTE->ReplaceDaughter(fVTE, envVTE); |
---|
116 | fVTE->ReplaceMother(fMVTE, envVTE); |
---|
117 | envVTE->AddDaughter(fVTE); |
---|
118 | envVTE->AddMother(fMVTE); |
---|
119 | |
---|
120 | // replace mother with envelope |
---|
121 | fMVTE = envVTE; |
---|
122 | } |
---|
123 | } |
---|
124 | } |
---|
125 | |
---|
126 | void G3Division::CreatePVReplica() |
---|
127 | { |
---|
128 | G4String name = fVTE->GetName(); |
---|
129 | G4LogicalVolume* lv = fVTE->GetLV(); |
---|
130 | G4LogicalVolume* mlv = fMVTE->GetLV(); |
---|
131 | |
---|
132 | G4String shape = fMVTE->GetShape(); |
---|
133 | if (shape == "PARA") { |
---|
134 | // The para volume cannot be replicated using G4PVReplica. |
---|
135 | // (Replicating a volume along a cartesian axis means "slicing" it |
---|
136 | // with slices -perpendicular- to that axis.) |
---|
137 | |
---|
138 | // position the replicated elements |
---|
139 | for (G4int i=0; i<fNofDivisions; i++) { |
---|
140 | G4ThreeVector position = G4ThreeVector(); |
---|
141 | position[fIAxis-1] = fLowRange + fWidth/2. + i*fWidth; |
---|
142 | if (position.y()!=0.) |
---|
143 | position.setX(position.y()*((G4Para*)lv->GetSolid())->GetTanAlpha()); |
---|
144 | |
---|
145 | #ifndef G3G4_NO_REFLECTION |
---|
146 | G4ReflectionFactory::Instance() |
---|
147 | ->Place(G4Translate3D(position), name, lv, mlv, 0, i); |
---|
148 | |
---|
149 | #else |
---|
150 | new G4PVPlacement(0, position, lv, name, mlv, 0, i); |
---|
151 | |
---|
152 | #endif |
---|
153 | } |
---|
154 | |
---|
155 | // G4PVReplica cannot be created |
---|
156 | return; |
---|
157 | } |
---|
158 | |
---|
159 | #ifdef G3G4DEBUG |
---|
160 | G4cout << "Create G4PVReplica name " << name << " logical volume name " |
---|
161 | << lv->GetName() << " mother logical volme name " |
---|
162 | << mlv->GetName() << " axis " << fAxis << " ndivisions " |
---|
163 | << fNofDivisions << " width " << fWidth << " Offset " |
---|
164 | << fOffset << G4endl; |
---|
165 | #endif |
---|
166 | |
---|
167 | #ifndef G3G4_NO_REFLECTION |
---|
168 | G4ReflectionFactory::Instance() |
---|
169 | ->Replicate(name, lv, mlv, fAxis, fNofDivisions, fWidth, fOffset); |
---|
170 | |
---|
171 | #else |
---|
172 | new G4PVReplica(name, lv, mlv, fAxis, fNofDivisions, fWidth, fOffset); |
---|
173 | |
---|
174 | #endif |
---|
175 | } |
---|
176 | |
---|
177 | // private methods |
---|
178 | |
---|
179 | void G3Division::Exception(G4String where, G4String what) { |
---|
180 | G4Exception("G3Division::" + where + " for " + what + " is not implemented"); |
---|
181 | } |
---|
182 | |
---|
183 | void G3Division::SetRangeAndAxis() |
---|
184 | // set fHighRange, fLowRange, fAxis |
---|
185 | { |
---|
186 | G4String shape = fMVTE->GetShape(); |
---|
187 | G4double *Rpar = fMVTE->GetRpar(); |
---|
188 | |
---|
189 | switch (fIAxis) { |
---|
190 | case 1: fAxis = kXAxis; |
---|
191 | break; |
---|
192 | case 2: fAxis = kYAxis; |
---|
193 | break; |
---|
194 | case 3: fAxis = kZAxis; |
---|
195 | break; |
---|
196 | default: G4Exception("G3Division: wrong iaxis defenition"); |
---|
197 | } |
---|
198 | |
---|
199 | if ( shape == "BOX" ) { |
---|
200 | fHighRange = Rpar[fIAxis-1]*cm; |
---|
201 | fLowRange = -fHighRange; |
---|
202 | } |
---|
203 | else if ( shape == "TRD1" ) { |
---|
204 | if (fIAxis == 1){ |
---|
205 | fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm); |
---|
206 | } |
---|
207 | else if( fIAxis == 2) { |
---|
208 | fHighRange = Rpar[2]*cm; |
---|
209 | } |
---|
210 | else if( fIAxis == 3) { |
---|
211 | fHighRange = Rpar[3]*cm; |
---|
212 | } |
---|
213 | fLowRange = - fHighRange; |
---|
214 | } |
---|
215 | else if ( shape == "TRD2" ) { |
---|
216 | if (fIAxis == 1){ |
---|
217 | fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm); |
---|
218 | } |
---|
219 | else if( fIAxis == 2) { |
---|
220 | fHighRange = std::max(Rpar[2]*cm, Rpar[3]*cm); |
---|
221 | } |
---|
222 | else if( fIAxis == 3) { |
---|
223 | fHighRange = Rpar[4]*cm; |
---|
224 | } |
---|
225 | } |
---|
226 | else if ( shape == "TRAP" ) { |
---|
227 | if ( fIAxis == 3 ) fHighRange = Rpar[0]*cm; |
---|
228 | else fHighRange = 0.; |
---|
229 | fLowRange = -fHighRange; |
---|
230 | } |
---|
231 | else if ( shape == "TUBE" ) { |
---|
232 | if (fIAxis == 1){ |
---|
233 | fHighRange = Rpar[1]*cm; |
---|
234 | fLowRange = Rpar[0]*cm; |
---|
235 | fAxis = kRho; |
---|
236 | } |
---|
237 | else if( fIAxis == 2) { |
---|
238 | fHighRange = 360.*deg; |
---|
239 | fLowRange = 0.; |
---|
240 | fAxis = kPhi; |
---|
241 | } |
---|
242 | else if( fIAxis == 3) { |
---|
243 | fHighRange = Rpar[2]*cm; |
---|
244 | fLowRange = -fHighRange; |
---|
245 | } |
---|
246 | } |
---|
247 | else if ( shape == "TUBS" ) { |
---|
248 | if (fIAxis == 1){ |
---|
249 | fHighRange = Rpar[1]*cm; |
---|
250 | fLowRange = Rpar[0]*cm; |
---|
251 | fAxis = kRho; |
---|
252 | } |
---|
253 | else if( fIAxis == 2) { |
---|
254 | |
---|
255 | fLowRange = Rpar[3]*deg; |
---|
256 | fHighRange = Rpar[4]*deg - fLowRange; |
---|
257 | if ( Rpar[4]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg; |
---|
258 | fHighRange = fHighRange + fLowRange; |
---|
259 | fAxis = kPhi; |
---|
260 | } |
---|
261 | else if( fIAxis == 3) { |
---|
262 | fHighRange = Rpar[2]*cm; |
---|
263 | fLowRange = -fHighRange; |
---|
264 | } |
---|
265 | } |
---|
266 | else if ( shape == "CONE" ) { |
---|
267 | if (fIAxis == 1){ |
---|
268 | fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm); |
---|
269 | fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm); |
---|
270 | fAxis = kRho; |
---|
271 | } |
---|
272 | else if( fIAxis == 2) { |
---|
273 | |
---|
274 | fLowRange = 0.; |
---|
275 | fHighRange = 360.*deg; |
---|
276 | fAxis = kPhi; |
---|
277 | } |
---|
278 | else if( fIAxis == 3) { |
---|
279 | fHighRange = Rpar[0]*cm; |
---|
280 | fLowRange = -fHighRange; |
---|
281 | } |
---|
282 | } |
---|
283 | else if ( shape == "CONS" ) { |
---|
284 | if (fIAxis == 1){ |
---|
285 | fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm); |
---|
286 | fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm); |
---|
287 | fAxis = kRho; |
---|
288 | } |
---|
289 | else if( fIAxis == 2) { |
---|
290 | |
---|
291 | fLowRange = Rpar[5]*deg; |
---|
292 | fHighRange = Rpar[6]*deg - fLowRange; |
---|
293 | if ( Rpar[6]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg; |
---|
294 | fHighRange = fHighRange + fLowRange; |
---|
295 | fAxis = kPhi; |
---|
296 | } |
---|
297 | else if( fIAxis == 3) { |
---|
298 | fHighRange = Rpar[2]*cm; |
---|
299 | fLowRange = -fHighRange; |
---|
300 | } |
---|
301 | } |
---|
302 | else if ( shape == "SPHE" ) { |
---|
303 | if (fIAxis == 1){ |
---|
304 | fHighRange = Rpar[1]*cm; |
---|
305 | fLowRange = Rpar[0]*cm; |
---|
306 | fAxis = kRho; |
---|
307 | } |
---|
308 | else if( fIAxis == 2) { |
---|
309 | fLowRange = std::min(Rpar[2]*deg,Rpar[3]*deg); |
---|
310 | fHighRange = std::max(Rpar[2]*deg,Rpar[3]*deg); |
---|
311 | fAxis = kPhi; |
---|
312 | } |
---|
313 | else if( fIAxis == 3) { |
---|
314 | fLowRange = std::min(Rpar[4]*deg,Rpar[5]*deg); |
---|
315 | fHighRange = std::max(Rpar[4]*deg,Rpar[5]*deg); |
---|
316 | fAxis = kPhi; // ?????? |
---|
317 | } |
---|
318 | } |
---|
319 | else if ( shape == "PARA" ) { |
---|
320 | fHighRange = Rpar[fIAxis-1]*cm; |
---|
321 | fLowRange = -fHighRange; |
---|
322 | } |
---|
323 | else if ( shape == "PGON" ) { |
---|
324 | G4int i; |
---|
325 | G4int nz = G4int(Rpar[3]); |
---|
326 | |
---|
327 | G4double pPhi1 = Rpar[0]*deg; |
---|
328 | G4double dPhi = Rpar[1]*deg; |
---|
329 | |
---|
330 | G4double *DzArray = new G4double[nz]; |
---|
331 | G4double *Rmax = new G4double[nz]; |
---|
332 | G4double *Rmin = new G4double[nz]; |
---|
333 | G4double rangehi[3], rangelo[3]; |
---|
334 | rangehi[0] = -kInfinity ; |
---|
335 | rangelo[0] = kInfinity ; |
---|
336 | rangehi[2] = -kInfinity ; |
---|
337 | rangelo[2] = kInfinity ; |
---|
338 | |
---|
339 | for(i=0; i<nz; i++) |
---|
340 | { |
---|
341 | G4int i4=3*i+4; |
---|
342 | G4int i5=i4+1; |
---|
343 | G4int i6=i4+2; |
---|
344 | |
---|
345 | DzArray[i] = Rpar[i4]*cm; |
---|
346 | Rmin[i] = Rpar[i5]*cm; |
---|
347 | Rmax[i] = Rpar[i6]*cm; |
---|
348 | rangelo[0] = std::min(rangelo[0], Rmin[i]); |
---|
349 | rangehi[0] = std::max(rangehi[0], Rmax[i]); |
---|
350 | rangelo[2] = std::min(rangelo[2], DzArray[i]); |
---|
351 | rangehi[2] = std::max(rangehi[2], DzArray[i]); |
---|
352 | } |
---|
353 | for (i=0;i<nz;i++){ |
---|
354 | assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]); |
---|
355 | } |
---|
356 | rangehi[1] = pPhi1 + dPhi; |
---|
357 | rangelo[1] = pPhi1; |
---|
358 | fHighRange = rangehi[fIAxis-1]; |
---|
359 | fLowRange = rangelo[fIAxis-1]; |
---|
360 | if (fIAxis == 1)fAxis = kRho; |
---|
361 | else if (fIAxis == 2)fAxis = kPhi; |
---|
362 | else if (fIAxis == 3)fAxis = kZAxis; |
---|
363 | |
---|
364 | delete [] DzArray; |
---|
365 | delete [] Rmin; |
---|
366 | delete [] Rmax; |
---|
367 | |
---|
368 | } |
---|
369 | else if ( shape == "PCON" ) { |
---|
370 | |
---|
371 | G4int i; |
---|
372 | G4double pPhi1 = Rpar[0]*deg; |
---|
373 | G4double dPhi = Rpar[1]*deg; |
---|
374 | G4int nz = G4int(Rpar[2]); |
---|
375 | |
---|
376 | G4double *DzArray = new G4double[nz]; |
---|
377 | G4double *Rmax = new G4double[nz]; |
---|
378 | G4double *Rmin = new G4double[nz]; |
---|
379 | G4double rangehi[3],rangelo[3]; |
---|
380 | |
---|
381 | rangehi[0] = -kInfinity ; |
---|
382 | rangelo[0] = kInfinity ; |
---|
383 | rangehi[2] = -kInfinity ; |
---|
384 | rangelo[2] = kInfinity ; |
---|
385 | |
---|
386 | for(i=0; i<nz; i++){ |
---|
387 | G4int i4=3*i+3; |
---|
388 | G4int i5=i4+1; |
---|
389 | G4int i6=i4+2; |
---|
390 | |
---|
391 | DzArray[i] = Rpar[i4]*cm; |
---|
392 | Rmin[i] = Rpar[i5]*cm; |
---|
393 | Rmax[i] = Rpar[i6]*cm; |
---|
394 | rangelo[0] = std::min(rangelo[0], Rmin[i]); |
---|
395 | rangehi[0] = std::max(rangehi[0], Rmax[i]); |
---|
396 | rangelo[2] = std::min(rangelo[2], DzArray[i]); |
---|
397 | rangehi[2] = std::max(rangehi[2], DzArray[i]); |
---|
398 | } |
---|
399 | for (i=0;i<nz;i++){ |
---|
400 | assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]); |
---|
401 | } |
---|
402 | rangehi[1] = pPhi1 + dPhi; |
---|
403 | rangelo[1] = pPhi1; |
---|
404 | fHighRange = rangehi[fIAxis-1]; |
---|
405 | fLowRange = rangelo[fIAxis-1]; |
---|
406 | if (fIAxis == 1)fAxis = kRho; |
---|
407 | else if (fIAxis == 2)fAxis = kPhi; |
---|
408 | else if (fIAxis == 3)fAxis = kZAxis; |
---|
409 | |
---|
410 | |
---|
411 | delete [] DzArray; |
---|
412 | delete [] Rmin; |
---|
413 | delete [] Rmax; |
---|
414 | } |
---|
415 | else if ( shape == "ELTU" || shape == "HYPE" || shape == "GTRA" || |
---|
416 | shape == "CTUB") { |
---|
417 | Exception("SetRangeAndAxis", shape); |
---|
418 | } |
---|
419 | else { |
---|
420 | Exception("SetRangeAndAxis", "Unknown shape" + shape); |
---|
421 | } |
---|
422 | |
---|
423 | // verbose |
---|
424 | #ifdef G3G4DEBUG |
---|
425 | G4cout << "Shape " << shape << " SetRangeAndAxis: " |
---|
426 | << fLowRange << " " << fHighRange << " " << fAxis << G4endl; |
---|
427 | #endif |
---|
428 | } |
---|
429 | |
---|
430 | G3VolTableEntry* G3Division::CreateEnvelope(G4String shape, G4double hi, |
---|
431 | G4double lo, G4double par[], G4int npar) |
---|
432 | // create new VTE with G3Pos corresponding to the |
---|
433 | // envelope of divided volume |
---|
434 | { |
---|
435 | // verbose |
---|
436 | // G4cout << " G3Division::CreateEnvelope " << "fIAaxis= " << fIAxis |
---|
437 | // << " hi= " << hi |
---|
438 | // << " lo= " << lo |
---|
439 | // << G4endl; |
---|
440 | |
---|
441 | G4double *Rpar = new G4double[npar+2]; |
---|
442 | for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];} |
---|
443 | G4double pos[3] = {0.,0.,0.}; |
---|
444 | |
---|
445 | if ( shape == "BOX" ) { |
---|
446 | Rpar[fIAxis-1] = (hi - lo)/2./cm; |
---|
447 | pos [fIAxis-1] = (hi + lo)/2.; |
---|
448 | } |
---|
449 | else if ( shape == "TRD1" ) { |
---|
450 | if ( fIAxis == 1 || fIAxis == 2 ) { |
---|
451 | Exception("CreateEnvelope","TRD1-x,y"); |
---|
452 | } |
---|
453 | else if ( fIAxis == 3 ) { |
---|
454 | // x = x1 + (c-z1)(x2 -x1)/(z2-z1) |
---|
455 | G4double tn, x1, z1; |
---|
456 | tn = (Rpar[1] - Rpar[0])/(2.* Rpar[3]); |
---|
457 | x1 = Rpar[0]; z1 = -Rpar[3]; |
---|
458 | Rpar[0] = x1 + tn * (lo/cm - z1); |
---|
459 | Rpar[1] = x1 + tn * (hi/cm - z1); |
---|
460 | Rpar[3] = (hi - lo)/2./cm; |
---|
461 | pos[2] = (hi + lo)/2.; |
---|
462 | } |
---|
463 | } |
---|
464 | else if ( shape == "TRD2" ) { |
---|
465 | if ( fIAxis == 1 || fIAxis == 2) { |
---|
466 | Exception("CreateEnvelope","TRD2-x,y"); |
---|
467 | } |
---|
468 | else if ( fIAxis == 3 ) { |
---|
469 | // x = x1 + (c-z1)(x2 -x1)/(z2-z1) |
---|
470 | // y = y1 + (c-z1)(y2 -y1)/(z2-z1) |
---|
471 | G4double tn1, tn2, x1, y1, z1; |
---|
472 | tn1 = (Rpar[1] - Rpar[0])/(2.* Rpar[4]); |
---|
473 | tn2 = (Rpar[3] - Rpar[2])/(2.* Rpar[4]); |
---|
474 | x1 = Rpar[0]; y1 = Rpar[2]; z1 = -Rpar[3]; |
---|
475 | Rpar[0] = x1 + tn1 * (lo/cm - z1); |
---|
476 | Rpar[1] = x1 + tn1 * (hi/cm - z1); |
---|
477 | Rpar[2] = y1 + tn2 * (lo/cm - z1); |
---|
478 | Rpar[3] = y1 + tn2 * (hi/cm - z1); |
---|
479 | Rpar[4] = (hi - lo)/2./cm; |
---|
480 | pos[2] = (hi + lo)/2.; |
---|
481 | } |
---|
482 | } |
---|
483 | else if ( shape == "TRAP" ) { |
---|
484 | Exception("CreateEnvelope","TRAP-x,y,z"); |
---|
485 | } |
---|
486 | else if ( shape == "TUBE" ) { |
---|
487 | if ( fIAxis == 1 ) { |
---|
488 | Rpar[0] = lo/cm; |
---|
489 | Rpar[1] = hi/cm; |
---|
490 | } |
---|
491 | else if ( fIAxis == 2 ) { |
---|
492 | Rpar[3] = lo/deg; |
---|
493 | Rpar[4] = hi/deg; |
---|
494 | npar = npar + 2; |
---|
495 | shape = "TUBS"; |
---|
496 | } |
---|
497 | else if ( fIAxis == 3 ) { |
---|
498 | Rpar[2] = (hi - lo)/2./cm; |
---|
499 | pos [2] = (hi + lo)/2.; |
---|
500 | } |
---|
501 | } |
---|
502 | else if ( shape == "TUBS" ) { |
---|
503 | if ( fIAxis == 1 ) { |
---|
504 | Rpar[0] = lo/cm; |
---|
505 | Rpar[1] = hi/cm; |
---|
506 | } |
---|
507 | else if ( fIAxis == 2 ) { |
---|
508 | Rpar[3] = lo/deg; |
---|
509 | Rpar[4] = hi/deg; |
---|
510 | } |
---|
511 | else if ( fIAxis == 3 ) { |
---|
512 | Rpar[2] = (hi - lo)/2./cm; |
---|
513 | pos [2] = (hi + lo)/2.; |
---|
514 | } |
---|
515 | } |
---|
516 | else if ( shape == "CONE" ) { |
---|
517 | if ( fIAxis == 1) { |
---|
518 | Exception("CreateEnvelope","CONE-x,z"); |
---|
519 | } |
---|
520 | else if ( fIAxis == 2 ) { |
---|
521 | Rpar[5] = lo/deg; |
---|
522 | Rpar[6] = hi/deg; |
---|
523 | npar = npar + 2; |
---|
524 | shape = "CONS"; |
---|
525 | } |
---|
526 | else if ( fIAxis == 3 ) { |
---|
527 | G4double tn1, tn2, rmin, rmax, z1; |
---|
528 | tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); |
---|
529 | tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); |
---|
530 | rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0]; |
---|
531 | Rpar[1] = rmin + tn1 * (lo/cm - z1); |
---|
532 | Rpar[3] = rmin + tn1 * (hi/cm - z1); |
---|
533 | Rpar[2] = rmax + tn2 * (lo/cm - z1); |
---|
534 | Rpar[4] = rmax + tn2 * (hi/cm - z1); |
---|
535 | Rpar[0] = (hi - lo)/2./cm; |
---|
536 | pos[2] = (hi + lo)/2.; |
---|
537 | } |
---|
538 | } |
---|
539 | else if ( shape == "CONS" ) { |
---|
540 | if ( fIAxis == 1 ) { |
---|
541 | Exception("CreateEnvelope","CONS-x"); |
---|
542 | } |
---|
543 | else if ( fIAxis == 2 ) { |
---|
544 | Rpar[5] = lo/deg; |
---|
545 | Rpar[6] = hi/deg; |
---|
546 | } |
---|
547 | else if ( fIAxis == 3 ) { |
---|
548 | G4double tn1, tn2, rmin, rmax, z1; |
---|
549 | tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); |
---|
550 | tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); |
---|
551 | rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0]; |
---|
552 | Rpar[1] = rmin + tn1 * (lo/cm - z1); |
---|
553 | Rpar[3] = rmin + tn1 * (hi/cm - z1); |
---|
554 | Rpar[2] = rmax + tn2 * (lo/cm - z1); |
---|
555 | Rpar[4] = rmax + tn2 * (hi/cm - z1); |
---|
556 | Rpar[0] = (hi - lo)/2./cm; |
---|
557 | pos[2] = (hi + lo)/2.; |
---|
558 | } |
---|
559 | } |
---|
560 | else if ( shape == "SPHE" ) { |
---|
561 | Exception("CreateEnvelope","SPHE-x,y,z"); |
---|
562 | } |
---|
563 | else if ( shape == "PARA" ) { |
---|
564 | Exception("CreateEnvelope","PARA-x,y,z"); |
---|
565 | } |
---|
566 | else if ( shape == "PGON" ) { |
---|
567 | if ( fIAxis == 2) { |
---|
568 | Rpar[0] = lo/deg; |
---|
569 | Rpar[1] = hi/deg; |
---|
570 | // rotm = ??? |
---|
571 | } |
---|
572 | else { |
---|
573 | Exception("CreateEnvelope","PGON-x,z"); |
---|
574 | } |
---|
575 | } |
---|
576 | else if ( shape == "PCON" ) { |
---|
577 | if ( fIAxis == 2) { |
---|
578 | Rpar[0] = lo/deg; |
---|
579 | Rpar[1] = hi/deg; |
---|
580 | // rotm = ??? |
---|
581 | } |
---|
582 | else { |
---|
583 | Exception("CreateEnvelope","PCON-x,z"); |
---|
584 | } |
---|
585 | } |
---|
586 | else { |
---|
587 | Exception("CreateEnvelope", "Unknown shape" + shape); |
---|
588 | } |
---|
589 | |
---|
590 | // create new VTE corresponding to envelope |
---|
591 | G4String envName = fVTE->GetName() + "_ENV"; |
---|
592 | G3VolTableEntry* envVTE |
---|
593 | = G4CreateVTE(envName, shape, fNmed, Rpar, npar); |
---|
594 | |
---|
595 | // create a G3Pos object and add it to envVTE |
---|
596 | G4String motherName = fMVTE->GetMasterClone()->GetName(); |
---|
597 | G4ThreeVector* offset = new G4ThreeVector(pos[0],pos[1],pos[2]); |
---|
598 | G4String only = "ONLY"; |
---|
599 | G3Pos* aG3Pos = new G3Pos(motherName, 1, offset, 0, only); |
---|
600 | envVTE->AddG3Pos(aG3Pos); |
---|
601 | |
---|
602 | delete [] Rpar; |
---|
603 | |
---|
604 | return envVTE; |
---|
605 | } |
---|
606 | |
---|
607 | void G3Division::CreateSolid(G4String shape, G4double par[], G4int npar) |
---|
608 | // create the solid corresponding to divided volume |
---|
609 | // and set the fOffset for replica |
---|
610 | { |
---|
611 | G4double *Rpar = new G4double[npar+2]; |
---|
612 | for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];} |
---|
613 | |
---|
614 | // verbose |
---|
615 | // G4cout << "G3Division::CreateSolid volume before: " |
---|
616 | // << fVTE->GetName() << " " << shape << G4endl; |
---|
617 | // G4cout << " npar,Rpar: " << npar; |
---|
618 | // for (G4int ii = 0; ii < npar; ++ii) G4cout << " " << Rpar[ii]; |
---|
619 | // G4cout << G4endl; |
---|
620 | |
---|
621 | if ( shape == "BOX" ) { |
---|
622 | if ( fIAxis == 1 ) Rpar[0] = fWidth/2./cm; |
---|
623 | else if ( fIAxis == 2 ) Rpar[1] = fWidth/2./cm; |
---|
624 | else if ( fIAxis == 3 ) Rpar[2] = fWidth/2./cm; |
---|
625 | } |
---|
626 | else if ( shape == "TRD1" ) { |
---|
627 | if ( fIAxis == 1 || fIAxis == 2 ) { |
---|
628 | Exception("CreateSolid", "TRD1-x,y"); |
---|
629 | } |
---|
630 | else if ( fIAxis == 3 ) { |
---|
631 | Rpar[3] = fWidth/2./cm; |
---|
632 | } |
---|
633 | } |
---|
634 | else if ( shape == "TRD2" ) { |
---|
635 | if ( fIAxis == 1 || fIAxis == 2 ) { |
---|
636 | Exception("CreateSolid", "TRD2-x,y"); |
---|
637 | } |
---|
638 | else if ( fIAxis == 3 ) { |
---|
639 | Rpar[4] = fWidth/2./cm; |
---|
640 | } |
---|
641 | } |
---|
642 | else if ( shape == "TRAP" ) { |
---|
643 | if ( fIAxis == 1 || fIAxis == 2) { |
---|
644 | Exception("CreateSolid", "TRAP-x,y"); |
---|
645 | } |
---|
646 | else if ( fIAxis == 3 ) { |
---|
647 | Rpar[0] = fWidth/2./cm; |
---|
648 | } |
---|
649 | } |
---|
650 | else if ( shape == "TUBE" ) { |
---|
651 | if ( fIAxis == 1 ) { |
---|
652 | Rpar[1] = Rpar[0] + fWidth/cm; |
---|
653 | fOffset = Rpar[0]*cm; |
---|
654 | } |
---|
655 | else if ( fIAxis == 2 ) { |
---|
656 | Rpar[3] = 0.; |
---|
657 | Rpar[4] = fWidth/deg; |
---|
658 | shape = "TUBS"; |
---|
659 | npar = npar + 2; |
---|
660 | } |
---|
661 | else if ( fIAxis == 3 ) { |
---|
662 | Rpar[2] = fWidth/2./cm; |
---|
663 | } |
---|
664 | } |
---|
665 | else if ( shape == "TUBS" ) { |
---|
666 | if ( fIAxis == 1 ) { |
---|
667 | Rpar[1] = Rpar[0] + fWidth/cm; |
---|
668 | fOffset = Rpar[0]*cm; |
---|
669 | } |
---|
670 | else if ( fIAxis == 2 ) { |
---|
671 | fOffset = Rpar[3]*deg; |
---|
672 | Rpar[3] = 0.; |
---|
673 | Rpar[4] = fWidth/deg; |
---|
674 | } |
---|
675 | else if ( fIAxis == 3 ) { |
---|
676 | Rpar[2] = fWidth/2./cm; |
---|
677 | } |
---|
678 | } |
---|
679 | else if ( shape == "CONE" ) { |
---|
680 | if ( fIAxis == 1 ) { |
---|
681 | Exception("CreateSolid", "CONE-x"); |
---|
682 | } |
---|
683 | else if ( fIAxis == 2 ) { |
---|
684 | Rpar[5] = 0.; |
---|
685 | Rpar[6] = fWidth/deg; |
---|
686 | shape = "CONS"; |
---|
687 | npar = npar + 2; |
---|
688 | } |
---|
689 | else if ( fIAxis == 3 ) { |
---|
690 | Rpar[0] = fWidth/2./cm; |
---|
691 | } |
---|
692 | } |
---|
693 | else if ( shape == "CONS" ) { |
---|
694 | if ( fIAxis == 1 ) { |
---|
695 | Exception("CreateSolid", "CONS-x"); |
---|
696 | } |
---|
697 | else if ( fIAxis == 2 ) { |
---|
698 | fOffset = Rpar[5]*deg; |
---|
699 | Rpar[5] = 0.; |
---|
700 | Rpar[6] = fWidth/deg; |
---|
701 | } |
---|
702 | else if ( fIAxis == 3 ) { |
---|
703 | Rpar[0] = fWidth/2./cm; |
---|
704 | } |
---|
705 | } |
---|
706 | else if (shape == "PARA") { |
---|
707 | if ( fIAxis == 1 ) { |
---|
708 | Rpar[0] = fWidth/2./cm; |
---|
709 | } |
---|
710 | else if ( Rpar[4] == 0. && Rpar[5] == 0. ) { |
---|
711 | // only special case for axis 2,3 is supported |
---|
712 | if ( fIAxis == 2 ) { |
---|
713 | Rpar[1] = fWidth/2./cm; |
---|
714 | } |
---|
715 | else if ( fIAxis == 3) { |
---|
716 | Rpar[2] = fWidth/2./cm; |
---|
717 | } |
---|
718 | } |
---|
719 | else |
---|
720 | Exception("CreateSolid", shape); |
---|
721 | } |
---|
722 | else if (shape == "SPHE") { |
---|
723 | Exception("CreateSolid", shape); |
---|
724 | } |
---|
725 | else if ( shape == "PGON" ) { |
---|
726 | if ( fIAxis == 2 ) { |
---|
727 | fOffset = Rpar[0]*deg; |
---|
728 | Rpar[0] = 0.; |
---|
729 | Rpar[1] = fWidth/deg; |
---|
730 | Rpar[2] = 1.; |
---|
731 | } |
---|
732 | else |
---|
733 | Exception("CreateSolid", shape); |
---|
734 | } |
---|
735 | else if ( shape == "PCON" ) { |
---|
736 | if ( fIAxis == 2 ) { |
---|
737 | fOffset = Rpar[0]*deg; |
---|
738 | Rpar[0] = 0.; |
---|
739 | Rpar[1] = fWidth/deg; |
---|
740 | } |
---|
741 | else { |
---|
742 | Exception("CreateSolid", shape); |
---|
743 | } |
---|
744 | } |
---|
745 | else { |
---|
746 | Exception("CreateSolid", "Unknown shape" + shape); |
---|
747 | } |
---|
748 | |
---|
749 | // create solid and set it to fVTE |
---|
750 | G4bool hasNegPars; |
---|
751 | G4bool deferred; |
---|
752 | G4bool okAxis[3]; |
---|
753 | G4VSolid* solid |
---|
754 | = G3toG4MakeSolid(fVTE->GetName(), shape, Rpar, npar, hasNegPars, deferred, okAxis); |
---|
755 | |
---|
756 | if (hasNegPars) { |
---|
757 | G4String name = fVTE->GetName(); |
---|
758 | G4Exception("CreateSolid VTE " + name + " has negative parameters."); |
---|
759 | } |
---|
760 | |
---|
761 | // update vte |
---|
762 | fVTE->SetSolid(solid); |
---|
763 | fVTE->SetNRpar(npar, Rpar); |
---|
764 | fVTE->SetHasNegPars(hasNegPars); |
---|
765 | |
---|
766 | // verbose |
---|
767 | // G4cout << "G3Division::CreateSolid volume after: " |
---|
768 | // << fVTE->GetName() << " " << shape << G4endl; |
---|
769 | // G4cout << " npar,Rpar: " << npar; |
---|
770 | // for (G4int iii = 0; iii < npar; ++iii) G4cout << " " << Rpar[iii]; |
---|
771 | // G4cout << G4endl; |
---|
772 | } |
---|
773 | |
---|
774 | |
---|
775 | G3VolTableEntry* G3Division::Dvn() |
---|
776 | { |
---|
777 | // no envelope need to be created |
---|
778 | |
---|
779 | // get parameters from mother |
---|
780 | G4String shape = fMVTE->GetShape(); |
---|
781 | G4double* Rpar = fMVTE->GetRpar(); |
---|
782 | G4int npar = fMVTE->GetNpar(); |
---|
783 | |
---|
784 | // set width for replica and create solid |
---|
785 | fWidth = (fHighRange - fLowRange)/fNofDivisions; |
---|
786 | CreateSolid(shape, Rpar, npar); |
---|
787 | |
---|
788 | return 0; |
---|
789 | } |
---|
790 | |
---|
791 | G3VolTableEntry* G3Division::Dvn2() |
---|
792 | { |
---|
793 | // to be defined as const of this class |
---|
794 | G4double Rmin = 0.0001*cm; |
---|
795 | |
---|
796 | G4String shape = fMVTE->GetShape(); |
---|
797 | G4double* Rpar = fMVTE->GetRpar(); |
---|
798 | G4int npar = fMVTE->GetNpar(); |
---|
799 | |
---|
800 | G4double c0 = fC0; |
---|
801 | if (fAxis == kPhi) c0 = c0*deg; |
---|
802 | else c0 = c0*cm; |
---|
803 | |
---|
804 | // create envelope (if needed) |
---|
805 | G3VolTableEntry* envVTE = 0; |
---|
806 | if( std::abs(c0 - fLowRange) > Rmin) { |
---|
807 | envVTE = CreateEnvelope(shape, fHighRange, c0, Rpar, npar); |
---|
808 | Rpar = envVTE->GetRpar(); |
---|
809 | npar = envVTE->GetNpar(); |
---|
810 | } |
---|
811 | |
---|
812 | // set width for replica and create solid |
---|
813 | fWidth = (fHighRange - c0)/fNofDivisions; |
---|
814 | CreateSolid(shape, Rpar, npar); |
---|
815 | |
---|
816 | return envVTE; |
---|
817 | } |
---|
818 | |
---|
819 | G3VolTableEntry* G3Division::Dvt() |
---|
820 | { |
---|
821 | // to be defined as const of this class |
---|
822 | G4double Rmin = 0.0001*cm; |
---|
823 | |
---|
824 | // get parameters from mother |
---|
825 | G4String shape = fMVTE->GetShape(); |
---|
826 | G4double* Rpar = fMVTE->GetRpar(); |
---|
827 | G4int npar = fMVTE->GetNpar(); |
---|
828 | |
---|
829 | // calculate the number of divisions |
---|
830 | G4int ndvmx = fNofDivisions; |
---|
831 | G4double step = fStep; |
---|
832 | |
---|
833 | if (fAxis == kPhi) step = step*deg; |
---|
834 | else step = step*cm; |
---|
835 | |
---|
836 | G4int ndiv = G4int((fHighRange - fLowRange + Rmin)/step); |
---|
837 | // to be added warning |
---|
838 | if (ndvmx > 255) ndvmx = 255; |
---|
839 | if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx; |
---|
840 | |
---|
841 | // create envVTE (if needed) |
---|
842 | G3VolTableEntry* envVTE = 0; |
---|
843 | G4double delta = std::abs((fHighRange - fLowRange) - ndiv*step); |
---|
844 | if (delta > Rmin) { |
---|
845 | envVTE |
---|
846 | = CreateEnvelope(shape, fHighRange-delta/2., fLowRange+delta/2., |
---|
847 | Rpar, npar); |
---|
848 | Rpar = envVTE->GetRpar(); |
---|
849 | npar = envVTE->GetNpar(); |
---|
850 | } |
---|
851 | |
---|
852 | // set width for replica and create solid |
---|
853 | fWidth = step; |
---|
854 | fNofDivisions = ndiv; |
---|
855 | CreateSolid(shape, Rpar, npar); |
---|
856 | |
---|
857 | return envVTE; |
---|
858 | } |
---|
859 | |
---|
860 | G3VolTableEntry* G3Division::Dvt2() |
---|
861 | { |
---|
862 | // to be defined as const of this class |
---|
863 | G4double Rmin = 0.0001*cm; |
---|
864 | |
---|
865 | // get parameters from mother |
---|
866 | G4String shape = fMVTE->GetShape(); |
---|
867 | G4double* Rpar = fMVTE->GetRpar(); |
---|
868 | G4int npar = fMVTE->GetNpar(); |
---|
869 | |
---|
870 | // calculate the number of divisions |
---|
871 | G4int ndvmx = fNofDivisions; |
---|
872 | G4double step = fStep; |
---|
873 | G4double c0 = fC0; |
---|
874 | |
---|
875 | if(fAxis == kPhi){ |
---|
876 | step = step*deg; |
---|
877 | c0 = c0*deg; |
---|
878 | } |
---|
879 | else { |
---|
880 | step = step*cm; |
---|
881 | c0 = c0*cm; |
---|
882 | } |
---|
883 | |
---|
884 | G4int ndiv = G4int((fHighRange - c0 + Rmin)/step); |
---|
885 | // to be added warning |
---|
886 | if (ndvmx > 255) ndvmx = 255; |
---|
887 | if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx; |
---|
888 | |
---|
889 | // create envelope (if needed) |
---|
890 | G3VolTableEntry* envVTE = 0; |
---|
891 | G4double delta = std::abs((fHighRange - c0) - ndiv*step); |
---|
892 | if (std::abs(c0 - fLowRange) > Rmin) { |
---|
893 | envVTE |
---|
894 | = CreateEnvelope(shape, fHighRange-delta/2., c0+delta/2., Rpar, npar); |
---|
895 | Rpar = envVTE->GetRpar(); |
---|
896 | npar = envVTE->GetNpar(); |
---|
897 | } |
---|
898 | |
---|
899 | // set with for replica and create solid |
---|
900 | fWidth = step; |
---|
901 | fNofDivisions = ndiv; |
---|
902 | CreateSolid(shape, Rpar, npar); |
---|
903 | |
---|
904 | return envVTE; |
---|
905 | } |
---|