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