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