Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/g3tog4/src/G3Division.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /g3tog4/src/G3Division.cc (Version 11.3.0) and /g3tog4/src/G3Division.cc (Version 9.6.p1)


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