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.2.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: G3Division.cc,v 1.17 2006/06/29 18:12:53 gunter Exp $
                                                   >>  28 // GEANT4 tag $Name: geant4-09-02 $
 27 //                                                 29 //
 28 // by I.Hrivnacova, V.Berejnoi 13.10.99            30 // by I.Hrivnacova, V.Berejnoi 13.10.99
 29                                                    31 
 30 #include <assert.h>                            << 
 31                                                << 
 32 #include "G3Division.hh"                           32 #include "G3Division.hh"
 33 #include "G3VolTableEntry.hh"                      33 #include "G3VolTableEntry.hh"
 34 #include "G3toG4MakeSolid.hh"                      34 #include "G3toG4MakeSolid.hh"
 35 #include "G4Para.hh"                               35 #include "G4Para.hh"
 36 #include "G3Pos.hh"                                36 #include "G3Pos.hh"
 37 #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   G4Exception("G3Division::" + where + " for " + what + " is not implemented");
182   G4String err_message = "G3Division::" + wher << 
183                        + what + " is not imple << 
184   G4Exception("G3Division::Exception()", "G3to << 
185               FatalException, err_message);    << 
186   return;                                      << 
187 }                                                 181 }  
188                                                   182 
189 void G3Division::SetRangeAndAxis()                183 void G3Division::SetRangeAndAxis()
190 // set fHighRange, fLowRange, fAxis               184 // set fHighRange, fLowRange, fAxis
191 {                                                 185 {
192     G4String shape = fMVTE->GetShape();           186     G4String shape = fMVTE->GetShape();
193     G4double *Rpar = fMVTE->GetRpar();            187     G4double *Rpar = fMVTE->GetRpar();
194                                                   188     
195     switch (fIAxis) {                             189     switch (fIAxis) {
196       case 1: fAxis = kXAxis;                     190       case 1: fAxis = kXAxis;
197               break;                              191               break;
198       case 2: fAxis = kYAxis;                     192       case 2: fAxis = kYAxis;
199               break;                              193               break;
200       case 3: fAxis = kZAxis;                     194       case 3: fAxis = kZAxis;
201               break;                              195               break;
202       default: G4Exception("G3Division::SetRan << 196       default: G4Exception("G3Division: wrong iaxis defenition");
203                             FatalException, "W << 
204     }                                             197     }
205                                                   198 
206     if ( shape == "BOX" ) {                       199     if ( shape == "BOX" ) {
207       fHighRange = Rpar[fIAxis-1]*cm;             200       fHighRange = Rpar[fIAxis-1]*cm;
208       fLowRange = -fHighRange;                    201       fLowRange = -fHighRange;
209     }                                             202     }
210     else if ( shape == "TRD1" ) {                 203     else if ( shape == "TRD1" ) {
211       if (fIAxis == 1){                           204       if (fIAxis == 1){
212         fHighRange = std::max(Rpar[0]*cm, Rpar    205         fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm);
213       }                                           206       }
214       else if( fIAxis == 2) {                     207       else if( fIAxis == 2) {
215        fHighRange = Rpar[2]*cm;                   208        fHighRange = Rpar[2]*cm;
216       }                                           209       }
217       else if( fIAxis == 3) {                     210       else if( fIAxis == 3) {
218        fHighRange = Rpar[3]*cm;                   211        fHighRange = Rpar[3]*cm;
219       }                                           212       }
220       fLowRange = - fHighRange;                   213       fLowRange = - fHighRange;
221     }                                             214     }
222     else if ( shape == "TRD2" ) {                 215     else if ( shape == "TRD2" ) {
223       if (fIAxis == 1){                           216       if (fIAxis == 1){
224         fHighRange = std::max(Rpar[0]*cm, Rpar    217         fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm);
225       }                                           218       }
226       else if( fIAxis == 2) {                     219       else if( fIAxis == 2) {
227         fHighRange = std::max(Rpar[2]*cm, Rpar    220         fHighRange = std::max(Rpar[2]*cm, Rpar[3]*cm);
228       }                                           221       }
229       else if( fIAxis == 3) {                     222       else if( fIAxis == 3) {
230        fHighRange = Rpar[4]*cm;                   223        fHighRange = Rpar[4]*cm;
231       }                                           224       }
232     }                                             225     }
233     else if ( shape == "TRAP" ) {                 226     else if ( shape == "TRAP" ) {
234       if ( fIAxis == 3 ) fHighRange = Rpar[0]*    227       if ( fIAxis == 3 ) fHighRange = Rpar[0]*cm;
235       else               fHighRange = 0.;         228       else               fHighRange = 0.;
236       fLowRange = -fHighRange;                    229       fLowRange = -fHighRange;
237     }                                             230     }
238     else if ( shape == "TUBE" ) {                 231     else if ( shape == "TUBE" ) {
239       if (fIAxis == 1){                           232       if (fIAxis == 1){
240         fHighRange = Rpar[1]*cm;                  233         fHighRange = Rpar[1]*cm;
241         fLowRange = Rpar[0]*cm;                   234         fLowRange = Rpar[0]*cm;
242         fAxis = kRho;                             235         fAxis = kRho;
243       }                                           236       }
244       else if( fIAxis == 2) {                     237       else if( fIAxis == 2) {
245         fHighRange = 360.*deg;                    238         fHighRange = 360.*deg;
246         fLowRange = 0.;                           239         fLowRange = 0.;
247         fAxis = kPhi;                             240         fAxis = kPhi;
248       }                                           241       }
249       else if( fIAxis == 3) {                     242       else if( fIAxis == 3) {
250        fHighRange = Rpar[2]*cm;                   243        fHighRange = Rpar[2]*cm;
251        fLowRange = -fHighRange;                   244        fLowRange = -fHighRange;
252       }                                           245       }
253     }                                             246     }
254     else if ( shape == "TUBS" ) {                 247     else if ( shape == "TUBS" ) {
255       if (fIAxis == 1){                           248       if (fIAxis == 1){
256         fHighRange = Rpar[1]*cm;                  249         fHighRange = Rpar[1]*cm;
257         fLowRange = Rpar[0]*cm;                   250         fLowRange = Rpar[0]*cm;
258         fAxis = kRho;                             251         fAxis = kRho;
259       }                                           252       }
260       else if( fIAxis == 2) {                     253       else if( fIAxis == 2) {
261                                                   254 
262        fLowRange = Rpar[3]*deg;                   255        fLowRange = Rpar[3]*deg;
263        fHighRange = Rpar[4]*deg - fLowRange;      256        fHighRange = Rpar[4]*deg - fLowRange;
264        if ( Rpar[4]*deg <= fLowRange )fHighRan    257        if ( Rpar[4]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg;
265        fHighRange = fHighRange + fLowRange;       258        fHighRange = fHighRange + fLowRange;
266        fAxis = kPhi;                              259        fAxis = kPhi;
267       }                                           260       }
268       else if( fIAxis == 3) {                     261       else if( fIAxis == 3) {
269        fHighRange = Rpar[2]*cm;                   262        fHighRange = Rpar[2]*cm;
270        fLowRange = -fHighRange;                   263        fLowRange = -fHighRange;
271       }                                           264       }
272     }                                             265     }
273     else if ( shape == "CONE" ) {                 266     else if ( shape == "CONE" ) {
274       if (fIAxis == 1){                           267       if (fIAxis == 1){
275         fHighRange = std::max(Rpar[2]*cm,Rpar[    268         fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm);
276         fLowRange = std::max(Rpar[1]*cm,Rpar[3    269         fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm);
277         fAxis = kRho;                             270         fAxis = kRho;
278       }                                           271       }
279       else if( fIAxis == 2) {                     272       else if( fIAxis == 2) {
280                                                   273 
281        fLowRange = 0.;                            274        fLowRange = 0.;
282        fHighRange = 360.*deg;                     275        fHighRange = 360.*deg;
283        fAxis = kPhi;                              276        fAxis = kPhi;
284       }                                           277       }
285       else if( fIAxis == 3) {                     278       else if( fIAxis == 3) {
286        fHighRange = Rpar[0]*cm;                   279        fHighRange = Rpar[0]*cm;
287        fLowRange = -fHighRange;                   280        fLowRange = -fHighRange;
288       }                                           281       }
289     }                                             282     }
290     else if ( shape == "CONS" ) {                 283     else if ( shape == "CONS" ) {
291       if (fIAxis == 1){                           284       if (fIAxis == 1){
292         fHighRange = std::max(Rpar[2]*cm,Rpar[    285         fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm);
293         fLowRange = std::max(Rpar[1]*cm,Rpar[3    286         fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm);
294         fAxis = kRho;                             287         fAxis = kRho;
295       }                                           288       }
296       else if( fIAxis == 2) {                     289       else if( fIAxis == 2) {
297                                                   290 
298        fLowRange = Rpar[5]*deg;                   291        fLowRange = Rpar[5]*deg;
299        fHighRange = Rpar[6]*deg - fLowRange;      292        fHighRange = Rpar[6]*deg - fLowRange;
300        if ( Rpar[6]*deg <= fLowRange )fHighRan    293        if ( Rpar[6]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg;
301        fHighRange = fHighRange + fLowRange;       294        fHighRange = fHighRange + fLowRange;
302        fAxis = kPhi;                              295        fAxis = kPhi;
303       }                                           296       }
304       else if( fIAxis == 3) {                     297       else if( fIAxis == 3) {
305        fHighRange = Rpar[2]*cm;                   298        fHighRange = Rpar[2]*cm;
306        fLowRange = -fHighRange;                   299        fLowRange = -fHighRange;
307       }                                           300       }
308     }                                             301     }
309     else if ( shape == "SPHE" ) {                 302     else if ( shape == "SPHE" ) {
310       if (fIAxis == 1){                           303       if (fIAxis == 1){
311         fHighRange = Rpar[1]*cm;                  304         fHighRange = Rpar[1]*cm;
312         fLowRange = Rpar[0]*cm;                   305         fLowRange = Rpar[0]*cm;
313         fAxis = kRho;                             306         fAxis = kRho;
314       }                                           307       }
315       else if( fIAxis == 2) {                     308       else if( fIAxis == 2) {
316        fLowRange = std::min(Rpar[2]*deg,Rpar[3    309        fLowRange = std::min(Rpar[2]*deg,Rpar[3]*deg);
317        fHighRange = std::max(Rpar[2]*deg,Rpar[    310        fHighRange = std::max(Rpar[2]*deg,Rpar[3]*deg);
318        fAxis = kPhi;                              311        fAxis = kPhi;
319       }                                           312       }
320       else if( fIAxis == 3) {                     313       else if( fIAxis == 3) {
321        fLowRange = std::min(Rpar[4]*deg,Rpar[5    314        fLowRange = std::min(Rpar[4]*deg,Rpar[5]*deg);
322        fHighRange = std::max(Rpar[4]*deg,Rpar[    315        fHighRange = std::max(Rpar[4]*deg,Rpar[5]*deg);
323        fAxis = kPhi; // ??????                    316        fAxis = kPhi; // ?????? 
324       }                                           317       }
325     }                                             318     }
326     else if ( shape == "PARA" ) {                 319     else if ( shape == "PARA" ) {
327       fHighRange = Rpar[fIAxis-1]*cm;             320       fHighRange = Rpar[fIAxis-1]*cm;
328       fLowRange = -fHighRange;                    321       fLowRange = -fHighRange;
329     }                                             322     }
330     else if ( shape == "PGON" ) {                 323     else if ( shape == "PGON" ) {
331         G4int i;                                  324         G4int i;
332         G4int nz = G4int(Rpar[3]);                325         G4int nz = G4int(Rpar[3]);
333                                                   326 
334         G4double pPhi1 = Rpar[0]*deg;             327         G4double pPhi1 = Rpar[0]*deg;
335         G4double dPhi  = Rpar[1]*deg;             328         G4double dPhi  = Rpar[1]*deg;
336                                                   329     
337         G4double *DzArray = new G4double[nz];     330         G4double *DzArray = new G4double[nz];
338         G4double *Rmax    = new G4double[nz];     331         G4double *Rmax    = new G4double[nz];
339         G4double *Rmin    = new G4double[nz];     332         G4double *Rmin    = new G4double[nz];
340         G4double rangehi[3], rangelo[3];          333         G4double rangehi[3], rangelo[3];
341         rangehi[0] = -kInfinity  ;                334         rangehi[0] = -kInfinity  ;
342         rangelo[0] =  kInfinity ;                 335         rangelo[0] =  kInfinity ;
343         rangehi[2] = -kInfinity ;                 336         rangehi[2] = -kInfinity ;
344         rangelo[2] =  kInfinity ;                 337         rangelo[2] =  kInfinity ;
345                                                   338 
346         for(i=0; i<nz; i++)                       339         for(i=0; i<nz; i++) 
347         {                                         340         {
348             G4int i4=3*i+4;                       341             G4int i4=3*i+4;
349             G4int i5=i4+1;                        342             G4int i5=i4+1;
350             G4int i6=i4+2;                        343             G4int i6=i4+2;
351                                                   344             
352             DzArray[i] = Rpar[i4]*cm;             345             DzArray[i] = Rpar[i4]*cm;
353             Rmin[i] = Rpar[i5]*cm;                346             Rmin[i] = Rpar[i5]*cm;
354             Rmax[i] = Rpar[i6]*cm;                347             Rmax[i] = Rpar[i6]*cm;
355             rangelo[0] = std::min(rangelo[0],     348             rangelo[0] = std::min(rangelo[0], Rmin[i]);
356             rangehi[0] = std::max(rangehi[0],     349             rangehi[0] = std::max(rangehi[0], Rmax[i]);
357             rangelo[2] = std::min(rangelo[2],     350             rangelo[2] = std::min(rangelo[2], DzArray[i]);
358             rangehi[2] = std::max(rangehi[2],     351             rangehi[2] = std::max(rangehi[2], DzArray[i]);
359         }                                         352         }
360         for (i=0;i<nz;i++){                       353         for (i=0;i<nz;i++){
361             assert(Rmin[i]>=0 && Rmax[i]>=Rmin    354             assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]);
362         }                                         355         }
363         rangehi[1] = pPhi1 + dPhi;                356         rangehi[1] = pPhi1 + dPhi;
364         rangelo[1] = pPhi1;                       357         rangelo[1] = pPhi1;
365         fHighRange = rangehi[fIAxis-1];           358         fHighRange = rangehi[fIAxis-1];
366         fLowRange = rangelo[fIAxis-1];            359         fLowRange = rangelo[fIAxis-1];
367         if      (fIAxis == 1)fAxis = kRho;        360         if      (fIAxis == 1)fAxis = kRho;
368         else if (fIAxis == 2)fAxis = kPhi;        361         else if (fIAxis == 2)fAxis = kPhi;
369         else if (fIAxis == 3)fAxis = kZAxis;      362         else if (fIAxis == 3)fAxis = kZAxis;
370                                                   363 
371         delete [] DzArray;                        364         delete [] DzArray;
372         delete [] Rmin;                           365         delete [] Rmin;
373         delete [] Rmax;                           366         delete [] Rmax;
374                                                   367 
375     }                                             368     }
376     else if ( shape == "PCON" ) {                 369     else if ( shape == "PCON" ) {
377                                                   370 
378         G4int i;                                  371         G4int i;
379         G4double pPhi1 = Rpar[0]*deg;             372         G4double pPhi1 = Rpar[0]*deg;
380         G4double dPhi  = Rpar[1]*deg;             373         G4double dPhi  = Rpar[1]*deg;    
381         G4int nz = G4int(Rpar[2]);                374         G4int nz = G4int(Rpar[2]);
382                                                   375     
383         G4double *DzArray = new G4double[nz];     376         G4double *DzArray = new G4double[nz];
384         G4double *Rmax    = new G4double[nz];     377         G4double *Rmax    = new G4double[nz];
385         G4double *Rmin    = new G4double[nz];     378         G4double *Rmin    = new G4double[nz];
386         G4double rangehi[3],rangelo[3];           379         G4double rangehi[3],rangelo[3];
387                                                   380 
388         rangehi[0] = -kInfinity  ;                381         rangehi[0] = -kInfinity  ;
389         rangelo[0] =  kInfinity ;                 382         rangelo[0] =  kInfinity ;
390         rangehi[2] = -kInfinity ;                 383         rangehi[2] = -kInfinity ;
391         rangelo[2] =  kInfinity ;                 384         rangelo[2] =  kInfinity ;
392                                                   385         
393         for(i=0; i<nz; i++){                      386         for(i=0; i<nz; i++){
394             G4int i4=3*i+3;                       387             G4int i4=3*i+3;
395             G4int i5=i4+1;                        388             G4int i5=i4+1;
396             G4int i6=i4+2;                        389             G4int i6=i4+2;
397                                                   390             
398             DzArray[i] = Rpar[i4]*cm;             391             DzArray[i] = Rpar[i4]*cm;
399             Rmin[i] = Rpar[i5]*cm;                392             Rmin[i] = Rpar[i5]*cm;
400             Rmax[i] = Rpar[i6]*cm;                393             Rmax[i] = Rpar[i6]*cm;
401             rangelo[0] = std::min(rangelo[0],     394             rangelo[0] = std::min(rangelo[0], Rmin[i]);
402             rangehi[0] = std::max(rangehi[0],     395             rangehi[0] = std::max(rangehi[0], Rmax[i]);
403             rangelo[2] = std::min(rangelo[2],     396             rangelo[2] = std::min(rangelo[2], DzArray[i]);
404             rangehi[2] = std::max(rangehi[2],     397             rangehi[2] = std::max(rangehi[2], DzArray[i]);
405         }                                         398         }
406         for (i=0;i<nz;i++){                       399         for (i=0;i<nz;i++){
407             assert(Rmin[i]>=0 && Rmax[i]>=Rmin    400             assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]);
408         }                                         401         }
409         rangehi[1] = pPhi1 + dPhi;                402         rangehi[1] = pPhi1 + dPhi;
410         rangelo[1] = pPhi1;                       403         rangelo[1] = pPhi1;
411         fHighRange = rangehi[fIAxis-1];           404         fHighRange = rangehi[fIAxis-1];
412         fLowRange = rangelo[fIAxis-1];            405         fLowRange = rangelo[fIAxis-1];
413         if      (fIAxis == 1)fAxis = kRho;        406         if      (fIAxis == 1)fAxis = kRho;
414         else if (fIAxis == 2)fAxis = kPhi;        407         else if (fIAxis == 2)fAxis = kPhi;
415         else if (fIAxis == 3)fAxis = kZAxis;      408         else if (fIAxis == 3)fAxis = kZAxis;
416                                                   409 
417                                                   410 
418         delete [] DzArray;                        411         delete [] DzArray;
419         delete [] Rmin;                           412         delete [] Rmin;
420         delete [] Rmax;                           413         delete [] Rmax;
421     }                                             414     }
422     else if ( shape == "ELTU" ||  shape == "HY    415     else if ( shape == "ELTU" ||  shape == "HYPE" || shape == "GTRA" ||
423          shape == "CTUB") {                       416          shape == "CTUB") {
424        Exception("SetRangeAndAxis", shape);       417        Exception("SetRangeAndAxis", shape);
425     }                                             418     }
426     else {                                        419     else {
427        Exception("SetRangeAndAxis", "Unknown s    420        Exception("SetRangeAndAxis", "Unknown shape" + shape);
428     }                                             421     }  
429                                                   422 
430     // verbose                                    423     // verbose
431     #ifdef G3G4DEBUG                              424     #ifdef G3G4DEBUG
432       G4cout << "Shape " << shape << " SetRang    425       G4cout << "Shape " << shape << " SetRangeAndAxis: " 
433        << fLowRange << " " << fHighRange << "     426        << fLowRange << " " << fHighRange << " " << fAxis << G4endl;
434     #endif                                        427     #endif
435 }                                                 428 }
436                                                   429 
437 G3VolTableEntry* G3Division::CreateEnvelope(G4    430 G3VolTableEntry* G3Division::CreateEnvelope(G4String shape, G4double hi, 
438                                G4double lo, G4    431                                G4double lo, G4double par[], G4int npar)
439 // create new VTE with G3Pos corresponding to     432 // create new VTE with G3Pos corresponding to the
440 // envelope of divided volume                     433 // envelope of divided volume
441 {                                                 434 {
442     // verbose                                    435     // verbose
443     // G4cout << "  G3Division::CreateEnvelope    436     // G4cout << "  G3Division::CreateEnvelope " << "fIAaxis= " << fIAxis
444     //        << " hi= " << hi                    437     //        << " hi= " << hi
445     //        << " lo= " << lo                    438     //        << " lo= " << lo
446     //        << G4endl;                          439     //        << G4endl;
447                                                   440 
448     G4double *Rpar = new G4double[npar+2];        441     G4double *Rpar = new G4double[npar+2];
449     for (G4int i=0; i<npar; ++i){ Rpar[i] = pa    442     for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];}
450     G4double pos[3] = {0.,0.,0.};                 443     G4double pos[3] = {0.,0.,0.};
451                                                   444   
452     if ( shape == "BOX" ) {                       445     if ( shape == "BOX" ) {
453       Rpar[fIAxis-1] = (hi - lo)/2./cm;           446       Rpar[fIAxis-1] = (hi - lo)/2./cm;
454       pos [fIAxis-1] = (hi + lo)/2.;              447       pos [fIAxis-1] = (hi + lo)/2.;
455     }                                             448     }
456     else if ( shape == "TRD1" ) {                 449     else if ( shape == "TRD1" ) {
457       if ( fIAxis == 1 || fIAxis == 2  ) {        450       if ( fIAxis == 1 || fIAxis == 2  ) {
458         Exception("CreateEnvelope","TRD1-x,y")    451         Exception("CreateEnvelope","TRD1-x,y");
459       }                                           452       }
460       else if ( fIAxis == 3 ) {                   453       else if ( fIAxis == 3 ) {
461   // x = x1 + (c-z1)(x2 -x1)/(z2-z1)              454   // x = x1 + (c-z1)(x2 -x1)/(z2-z1)
462   G4double tn, x1, z1;                            455   G4double tn, x1, z1;
463         tn = (Rpar[1] - Rpar[0])/(2.* Rpar[3])    456         tn = (Rpar[1] - Rpar[0])/(2.* Rpar[3]); 
464         x1 = Rpar[0]; z1 = -Rpar[3];              457         x1 = Rpar[0]; z1 = -Rpar[3];
465         Rpar[0] = x1 + tn * (lo/cm - z1);         458         Rpar[0] = x1 + tn * (lo/cm - z1);
466         Rpar[1] = x1 + tn * (hi/cm - z1);         459         Rpar[1] = x1 + tn * (hi/cm - z1);
467         Rpar[3] = (hi - lo)/2./cm;                460         Rpar[3] = (hi - lo)/2./cm;
468         pos[2]  = (hi + lo)/2.;                   461         pos[2]  = (hi + lo)/2.;
469       }                                           462       }
470     }                                             463     }
471     else if ( shape == "TRD2" ) {                 464     else if ( shape == "TRD2" ) {
472       if ( fIAxis == 1 || fIAxis == 2) {          465       if ( fIAxis == 1 || fIAxis == 2) {
473         Exception("CreateEnvelope","TRD2-x,y")    466         Exception("CreateEnvelope","TRD2-x,y");
474       }                                           467       }
475       else if ( fIAxis == 3 ) {                   468       else if ( fIAxis == 3 ) {
476   // x = x1 + (c-z1)(x2 -x1)/(z2-z1)              469   // x = x1 + (c-z1)(x2 -x1)/(z2-z1)
477   // y = y1 + (c-z1)(y2 -y1)/(z2-z1)              470   // y = y1 + (c-z1)(y2 -y1)/(z2-z1)
478   G4double tn1, tn2, x1, y1, z1;                  471   G4double tn1, tn2, x1, y1, z1;
479         tn1 = (Rpar[1] - Rpar[0])/(2.* Rpar[4]    472         tn1 = (Rpar[1] - Rpar[0])/(2.* Rpar[4]); 
480         tn2 = (Rpar[3] - Rpar[2])/(2.* Rpar[4]    473         tn2 = (Rpar[3] - Rpar[2])/(2.* Rpar[4]); 
481         x1 = Rpar[0]; y1 = Rpar[2]; z1 = -Rpar    474         x1 = Rpar[0]; y1 = Rpar[2]; z1 = -Rpar[3];
482         Rpar[0] = x1 + tn1 * (lo/cm - z1);        475         Rpar[0] = x1 + tn1 * (lo/cm - z1);
483         Rpar[1] = x1 + tn1 * (hi/cm - z1);        476         Rpar[1] = x1 + tn1 * (hi/cm - z1);
484         Rpar[2] = y1 + tn2 * (lo/cm - z1);        477         Rpar[2] = y1 + tn2 * (lo/cm - z1);
485         Rpar[3] = y1 + tn2 * (hi/cm - z1);        478         Rpar[3] = y1 + tn2 * (hi/cm - z1);
486         Rpar[4] = (hi - lo)/2./cm;                479         Rpar[4] = (hi - lo)/2./cm;
487         pos[2]  = (hi + lo)/2.;                   480         pos[2]  = (hi + lo)/2.;
488       }                                           481       }
489     }                                             482     }
490     else if ( shape == "TRAP" ) {                 483     else if ( shape == "TRAP" ) {
491       Exception("CreateEnvelope","TRAP-x,y,z")    484       Exception("CreateEnvelope","TRAP-x,y,z");
492     }                                             485     }
493     else if ( shape == "TUBE" ) {                 486     else if ( shape == "TUBE" ) {
494       if ( fIAxis == 1 ) {                        487       if ( fIAxis == 1 ) {
495         Rpar[0] = lo/cm;                          488         Rpar[0] = lo/cm;
496         Rpar[1] = hi/cm;                          489         Rpar[1] = hi/cm;
497       }                                           490       }
498       else if ( fIAxis == 2 ) {                   491       else if ( fIAxis == 2 ) {
499         Rpar[3] = lo/deg;                         492         Rpar[3] = lo/deg;
500         Rpar[4] = hi/deg;                         493         Rpar[4] = hi/deg;
501         npar = npar + 2;                          494         npar = npar + 2;
502         shape = "TUBS";                           495         shape = "TUBS";
503       }                                           496       }
504       else if ( fIAxis == 3 ) {                   497       else if ( fIAxis == 3 ) {
505         Rpar[2] = (hi - lo)/2./cm;                498         Rpar[2] = (hi - lo)/2./cm;
506         pos [2] = (hi + lo)/2.;                   499         pos [2] = (hi + lo)/2.;
507       }                                           500       }
508     }                                             501     }
509     else if ( shape == "TUBS" ) {                 502     else if ( shape == "TUBS" ) {
510       if ( fIAxis == 1 ) {                        503       if ( fIAxis == 1 ) {
511         Rpar[0] = lo/cm;                          504         Rpar[0] = lo/cm;
512         Rpar[1] = hi/cm;                          505         Rpar[1] = hi/cm;
513       }                                           506       }
514       else if ( fIAxis == 2 ) {                   507       else if ( fIAxis == 2 ) {
515         Rpar[3] = lo/deg;                         508         Rpar[3] = lo/deg;
516         Rpar[4] = hi/deg;                         509         Rpar[4] = hi/deg;
517       }                                           510       }
518       else if ( fIAxis == 3 ) {                   511       else if ( fIAxis == 3 ) {
519         Rpar[2] = (hi - lo)/2./cm;                512         Rpar[2] = (hi - lo)/2./cm;
520         pos [2] = (hi + lo)/2.;                   513         pos [2] = (hi + lo)/2.;
521       }                                           514       }
522     }                                             515     }
523     else if ( shape == "CONE" ) {                 516     else if ( shape == "CONE" ) {
524       if ( fIAxis == 1) {                         517       if ( fIAxis == 1) {
525         Exception("CreateEnvelope","CONE-x,z")    518         Exception("CreateEnvelope","CONE-x,z");
526       }                                           519       }
527       else if ( fIAxis == 2 ) {                   520       else if ( fIAxis == 2 ) {
528         Rpar[5] = lo/deg;                         521         Rpar[5] = lo/deg;
529         Rpar[6] = hi/deg;                         522         Rpar[6] = hi/deg;
530         npar = npar + 2;                          523         npar = npar + 2;
531         shape = "CONS";                           524         shape = "CONS";
532       }                                           525       }
533       else if ( fIAxis == 3 ) {                   526       else if ( fIAxis == 3 ) {
534         G4double tn1, tn2, rmin, rmax, z1;        527         G4double tn1, tn2, rmin, rmax, z1;
535         tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]    528         tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); 
536         tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]    529         tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); 
537         rmin = Rpar[1]; rmax = Rpar[2]; z1 = -    530         rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0];
538         Rpar[1] = rmin + tn1 * (lo/cm - z1);      531         Rpar[1] = rmin + tn1 * (lo/cm - z1);
539         Rpar[3] = rmin + tn1 * (hi/cm - z1);      532         Rpar[3] = rmin + tn1 * (hi/cm - z1);
540         Rpar[2] = rmax + tn2 * (lo/cm - z1);      533         Rpar[2] = rmax + tn2 * (lo/cm - z1);
541         Rpar[4] = rmax + tn2 * (hi/cm - z1);      534         Rpar[4] = rmax + tn2 * (hi/cm - z1);
542         Rpar[0] = (hi - lo)/2./cm;                535         Rpar[0] = (hi - lo)/2./cm;
543         pos[2]  = (hi + lo)/2.;                   536         pos[2]  = (hi + lo)/2.;
544       }                                           537       }
545     }                                             538     }
546     else if ( shape == "CONS" ) {                 539     else if ( shape == "CONS" ) {
547       if ( fIAxis == 1 ) {                        540       if ( fIAxis == 1 ) {
548         Exception("CreateEnvelope","CONS-x");     541         Exception("CreateEnvelope","CONS-x");
549       }                                           542       }
550       else if ( fIAxis == 2 ) {                   543       else if ( fIAxis == 2 ) {
551         Rpar[5] = lo/deg;                         544         Rpar[5] = lo/deg;
552         Rpar[6] = hi/deg;                         545         Rpar[6] = hi/deg;
553       }                                           546       }
554       else if ( fIAxis == 3 ) {                   547       else if ( fIAxis == 3 ) {
555         G4double tn1, tn2, rmin, rmax, z1;        548         G4double tn1, tn2, rmin, rmax, z1;
556         tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]    549         tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); 
557         tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]    550         tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); 
558         rmin = Rpar[1]; rmax = Rpar[2]; z1 = -    551         rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0];
559         Rpar[1] = rmin + tn1 * (lo/cm - z1);      552         Rpar[1] = rmin + tn1 * (lo/cm - z1);
560         Rpar[3] = rmin + tn1 * (hi/cm - z1);      553         Rpar[3] = rmin + tn1 * (hi/cm - z1);
561         Rpar[2] = rmax + tn2 * (lo/cm - z1);      554         Rpar[2] = rmax + tn2 * (lo/cm - z1);
562         Rpar[4] = rmax + tn2 * (hi/cm - z1);      555         Rpar[4] = rmax + tn2 * (hi/cm - z1);
563         Rpar[0] = (hi - lo)/2./cm;                556         Rpar[0] = (hi - lo)/2./cm;
564         pos[2]  = (hi + lo)/2.;                   557         pos[2]  = (hi + lo)/2.;
565       }                                           558       }
566     }                                             559     }
567     else if ( shape == "SPHE" ) {                 560     else if ( shape == "SPHE" ) {
568       Exception("CreateEnvelope","SPHE-x,y,z")    561       Exception("CreateEnvelope","SPHE-x,y,z");                
569     }                                             562     }
570     else if ( shape == "PARA" ) {                 563     else if ( shape == "PARA" ) {
571       Exception("CreateEnvelope","PARA-x,y,z")    564       Exception("CreateEnvelope","PARA-x,y,z");
572     }                                             565     }
573     else if ( shape == "PGON" ) {                 566     else if ( shape == "PGON" ) {
574       if ( fIAxis == 2) {                         567       if ( fIAxis == 2) {
575   Rpar[0] = lo/deg;                               568   Rpar[0] = lo/deg;
576   Rpar[1] = hi/deg;                               569   Rpar[1] = hi/deg;
577   // rotm = ???                                   570   // rotm = ???
578       }                                           571       }
579       else {                                      572       else {
580         Exception("CreateEnvelope","PGON-x,z")    573         Exception("CreateEnvelope","PGON-x,z");
581       }                                           574       }
582     }                                             575     }
583     else if ( shape == "PCON" ) {                 576     else if ( shape == "PCON" ) {
584       if ( fIAxis == 2) {                         577       if ( fIAxis == 2) {
585   Rpar[0] = lo/deg;                               578   Rpar[0] = lo/deg;
586   Rpar[1] = hi/deg;                               579   Rpar[1] = hi/deg;
587   // rotm = ???                                   580   // rotm = ???
588       }                                           581       }
589       else {                                      582       else {
590         Exception("CreateEnvelope","PCON-x,z")    583         Exception("CreateEnvelope","PCON-x,z");
591       }                                           584       }
592     }                                             585     }
593     else {                                        586     else {
594        Exception("CreateEnvelope", "Unknown sh    587        Exception("CreateEnvelope", "Unknown shape" + shape);
595     }                                             588     }  
596                                                   589 
597     // create new VTE corresponding to envelop    590     // create new VTE corresponding to envelope
598     G4String envName = fVTE->GetName() + "_ENV    591     G4String envName = fVTE->GetName() + "_ENV"; 
599     G3VolTableEntry* envVTE                       592     G3VolTableEntry* envVTE 
600       = G4CreateVTE(envName, shape, fNmed, Rpa    593       = G4CreateVTE(envName, shape, fNmed, Rpar, npar);
601                                                   594 
602     // create a G3Pos object and add it to env    595     // create a G3Pos object and add it to envVTE
603     G4String motherName = fMVTE->GetMasterClon    596     G4String motherName = fMVTE->GetMasterClone()->GetName();
604     G4ThreeVector* offset = new G4ThreeVector(    597     G4ThreeVector* offset = new G4ThreeVector(pos[0],pos[1],pos[2]);    
605     G4String only = "ONLY";                       598     G4String only = "ONLY";
606     G3Pos* aG3Pos = new G3Pos(motherName, 1, o    599     G3Pos* aG3Pos = new G3Pos(motherName, 1, offset, 0, only);              
607     envVTE->AddG3Pos(aG3Pos);                     600     envVTE->AddG3Pos(aG3Pos);
608                                                   601 
609     delete [] Rpar;                               602     delete [] Rpar; 
610                                                   603 
611     return envVTE;                                604     return envVTE;
612 }                                                 605 }
613                                                   606 
614 void G3Division::CreateSolid(G4String shape, G    607 void G3Division::CreateSolid(G4String shape, G4double par[], G4int npar)
615 // create the solid corresponding to divided v    608 // create the solid corresponding to divided volume
616 // and set the fOffset for replica                609 // and set the fOffset for replica
617 {                                                 610 {
618     G4double *Rpar = new G4double[npar+2];        611     G4double *Rpar = new G4double[npar+2];
619     for (G4int i=0; i<npar; ++i){ Rpar[i] = pa    612     for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];}
620                                                   613 
621     // verbose                                    614     // verbose
622     // G4cout << "G3Division::CreateSolid volu    615     // G4cout << "G3Division::CreateSolid volume before: " 
623     //        << fVTE->GetName() << " " << sha    616     //        << fVTE->GetName() << " " << shape << G4endl;    
624     // G4cout << " npar,Rpar: " << npar;          617     // G4cout << " npar,Rpar: " << npar;
625     // for (G4int ii = 0; ii < npar; ++ii) G4c    618     // for (G4int ii = 0; ii < npar; ++ii) G4cout << " " << Rpar[ii];
626     // G4cout << G4endl;                          619     // G4cout << G4endl;
627                                                   620   
628     if ( shape == "BOX" ) {                       621     if ( shape == "BOX" ) {
629       if      ( fIAxis == 1 ) Rpar[0] = fWidth    622       if      ( fIAxis == 1 ) Rpar[0] = fWidth/2./cm;
630       else if ( fIAxis == 2 ) Rpar[1] = fWidth    623       else if ( fIAxis == 2 ) Rpar[1] = fWidth/2./cm; 
631       else if ( fIAxis == 3 ) Rpar[2] = fWidth    624       else if ( fIAxis == 3 ) Rpar[2] = fWidth/2./cm; 
632     }                                             625     }
633     else if ( shape == "TRD1" ) {                 626     else if ( shape == "TRD1" ) {
634       if ( fIAxis == 1 || fIAxis == 2 ) {         627       if ( fIAxis == 1 || fIAxis == 2 ) {
635         Exception("CreateSolid", "TRD1-x,y");     628         Exception("CreateSolid", "TRD1-x,y");
636       }                                           629       }
637       else if ( fIAxis == 3 ) {                   630       else if ( fIAxis == 3 ) {
638          Rpar[3] = fWidth/2./cm;                  631          Rpar[3] = fWidth/2./cm; 
639       }                                           632       }
640     }                                             633     }
641     else if ( shape == "TRD2" ) {                 634     else if ( shape == "TRD2" ) {
642       if ( fIAxis == 1 || fIAxis == 2 ) {         635       if ( fIAxis == 1 || fIAxis == 2 ) {
643         Exception("CreateSolid", "TRD2-x,y");     636         Exception("CreateSolid", "TRD2-x,y");
644       }                                           637       }
645       else if ( fIAxis == 3 ) {                   638       else if ( fIAxis == 3 ) {
646          Rpar[4] =  fWidth/2./cm;                 639          Rpar[4] =  fWidth/2./cm; 
647       }                                           640       }
648     }                                             641     }
649     else if ( shape == "TRAP" ) {                 642     else if ( shape == "TRAP" ) {
650       if ( fIAxis == 1 || fIAxis == 2) {          643       if ( fIAxis == 1 || fIAxis == 2) {
651         Exception("CreateSolid", "TRAP-x,y");     644         Exception("CreateSolid", "TRAP-x,y");
652       }                                           645       }
653       else if ( fIAxis == 3 ) {                   646       else if ( fIAxis == 3 ) {
654          Rpar[0] =  fWidth/2./cm;                 647          Rpar[0] =  fWidth/2./cm; 
655       }                                           648       }
656     }                                             649     }
657     else if ( shape == "TUBE" ) {                 650     else if ( shape == "TUBE" ) {
658       if ( fIAxis == 1 ) {                        651       if ( fIAxis == 1 ) {
659          Rpar[1] = Rpar[0] + fWidth/cm;           652          Rpar[1] = Rpar[0] + fWidth/cm;
660          fOffset = Rpar[0]*cm;                    653          fOffset = Rpar[0]*cm;
661       }                                           654       }
662       else if ( fIAxis == 2 ) {                   655       else if ( fIAxis == 2 ) {
663          Rpar[3] = 0.;                            656          Rpar[3] = 0.; 
664          Rpar[4] = fWidth/deg;                    657          Rpar[4] = fWidth/deg; 
665          shape = "TUBS";                          658          shape = "TUBS";
666          npar = npar + 2;                         659          npar = npar + 2;
667       }                                           660       }
668       else if ( fIAxis == 3 ) {                   661       else if ( fIAxis == 3 ) {
669          Rpar[2] = fWidth/2./cm;                  662          Rpar[2] = fWidth/2./cm; 
670       }                                           663       }
671     }                                             664     }
672     else if ( shape == "TUBS" ) {                 665     else if ( shape == "TUBS" ) {
673       if ( fIAxis == 1 ) {                        666       if ( fIAxis == 1 ) {
674         Rpar[1] = Rpar[0] + fWidth/cm;            667         Rpar[1] = Rpar[0] + fWidth/cm;
675         fOffset = Rpar[0]*cm;                     668         fOffset = Rpar[0]*cm;
676       }                                           669       }
677       else if ( fIAxis == 2 ) {                   670       else if ( fIAxis == 2 ) {
678          fOffset = Rpar[3]*deg;                   671          fOffset = Rpar[3]*deg; 
679          Rpar[3] = 0.;                            672          Rpar[3] = 0.;
680          Rpar[4] =  fWidth/deg;                   673          Rpar[4] =  fWidth/deg;
681       }                                           674       }
682       else if ( fIAxis == 3 ) {                   675       else if ( fIAxis == 3 ) {
683          Rpar[2] = fWidth/2./cm;                  676          Rpar[2] = fWidth/2./cm; 
684       }                                           677       }
685     }                                             678     }
686     else if ( shape == "CONE" ) {                 679     else if ( shape == "CONE" ) {
687       if ( fIAxis == 1 ) {                        680       if ( fIAxis == 1 ) {
688         Exception("CreateSolid", "CONE-x");       681         Exception("CreateSolid", "CONE-x"); 
689       }                                           682       }
690       else if ( fIAxis == 2 ) {                   683       else if ( fIAxis == 2 ) {
691          Rpar[5] = 0.;                            684          Rpar[5] = 0.;
692          Rpar[6] = fWidth/deg;                    685          Rpar[6] = fWidth/deg;
693          shape = "CONS";                          686          shape = "CONS";
694          npar = npar + 2;                         687          npar = npar + 2;
695       }                                           688       }
696       else if ( fIAxis == 3 ) {                   689       else if ( fIAxis == 3 ) {
697          Rpar[0] = fWidth/2./cm;                  690          Rpar[0] = fWidth/2./cm; 
698       }                                           691       }
699     }                                             692     }
700     else if ( shape == "CONS" ) {                 693     else if ( shape == "CONS" ) {
701       if ( fIAxis == 1 ) {                        694       if ( fIAxis == 1 ) {
702         Exception("CreateSolid", "CONS-x");       695         Exception("CreateSolid", "CONS-x"); 
703       }                                           696       }
704       else if ( fIAxis == 2 ) {                   697       else if ( fIAxis == 2 ) {
705          fOffset = Rpar[5]*deg;                   698          fOffset = Rpar[5]*deg;
706          Rpar[5] = 0.;                            699          Rpar[5] = 0.;
707          Rpar[6] = fWidth/deg;                    700          Rpar[6] = fWidth/deg;
708       }                                           701       }
709       else if ( fIAxis == 3 ) {                   702       else if ( fIAxis == 3 ) {
710          Rpar[0] = fWidth/2./cm;                  703          Rpar[0] = fWidth/2./cm; 
711       }                                           704       }
712     }                                             705     }
713     else if (shape == "PARA") {                   706     else if (shape == "PARA") {
714       if      ( fIAxis == 1 ) {                   707       if      ( fIAxis == 1 ) {
715          Rpar[0] = fWidth/2./cm;                  708          Rpar[0] = fWidth/2./cm;
716       }                                           709       }  
717       else if ( Rpar[4] == 0. && Rpar[5] == 0.    710       else if ( Rpar[4] == 0. && Rpar[5] == 0. ) {
718          // only special case for axis 2,3 is     711          // only special case for axis 2,3 is supported
719         if ( fIAxis == 2 ) {                      712         if ( fIAxis == 2 ) {
720           Rpar[1] = fWidth/2./cm;                 713           Rpar[1] = fWidth/2./cm;
721   }                                               714   }  
722     else if ( fIAxis == 3) {                      715     else if ( fIAxis == 3) {
723           Rpar[2] = fWidth/2./cm;                 716           Rpar[2] = fWidth/2./cm;
724   }                                               717   }
725       }                                           718       }   
726       else                                        719       else    
727          Exception("CreateSolid", shape);         720          Exception("CreateSolid", shape);
728     }                                             721     }  
729     else if (shape == "SPHE") {                   722     else if (shape == "SPHE") {
730       Exception("CreateSolid", shape);            723       Exception("CreateSolid", shape);
731     }                                             724     }
732     else if ( shape == "PGON" ) {                 725     else if ( shape == "PGON" ) {
733       if ( fIAxis == 2 ) {                        726       if ( fIAxis == 2 ) {
734          fOffset = Rpar[0]*deg;                   727          fOffset = Rpar[0]*deg;
735          Rpar[0] = 0.;                            728          Rpar[0] = 0.;
736          Rpar[1] = fWidth/deg;                    729          Rpar[1] = fWidth/deg;
737          Rpar[2] = 1.;                            730          Rpar[2] = 1.;
738       }                                           731       }
739       else                                        732       else
740        Exception("CreateSolid", shape);           733        Exception("CreateSolid", shape);
741     }                                             734     }
742     else if ( shape == "PCON" ) {                 735     else if ( shape == "PCON" ) {
743       if ( fIAxis == 2 ) {                        736       if ( fIAxis == 2 ) {
744          fOffset = Rpar[0]*deg;                   737          fOffset = Rpar[0]*deg;
745          Rpar[0] = 0.;                            738          Rpar[0] = 0.;
746          Rpar[1] = fWidth/deg;                    739          Rpar[1] = fWidth/deg;
747       }                                           740       }
748       else {                                      741       else {
749         Exception("CreateSolid", shape);          742         Exception("CreateSolid", shape);
750       }                                           743       } 
751     }                                             744     }
752     else {                                        745     else {
753        Exception("CreateSolid", "Unknown shape    746        Exception("CreateSolid", "Unknown shape" + shape);
754     }                                             747     }  
755                                                   748 
756     // create solid and set it to fVTE            749     // create solid and set it to fVTE
757     G4bool hasNegPars;                            750     G4bool hasNegPars;
758     G4bool deferred;                              751     G4bool deferred;   
759     G4bool okAxis[3];                             752     G4bool okAxis[3];
760     G4VSolid* solid                               753     G4VSolid* solid
761     = G3toG4MakeSolid(fVTE->GetName(), shape,     754     = G3toG4MakeSolid(fVTE->GetName(), shape, Rpar, npar, hasNegPars, deferred, okAxis);  
762                                                   755 
763     if (hasNegPars) {                             756     if (hasNegPars) {
764        G4String err_message = "CreateSolid VTE << 757        G4String name = fVTE->GetName();
765                             + " has negative p << 758        G4Exception("CreateSolid VTE " + name + " has negative parameters.");
766        G4Exception("G3Division::CreateSolid()" << 
767                    FatalException, err_message << 
768        return;                                 << 
769     }                                             759     }   
770                                                   760     
771     // update vte                                 761     // update vte
772     fVTE->SetSolid(solid);                        762     fVTE->SetSolid(solid);
773     fVTE->SetNRpar(npar, Rpar);                   763     fVTE->SetNRpar(npar, Rpar); 
774     fVTE->SetHasNegPars(hasNegPars);              764     fVTE->SetHasNegPars(hasNegPars);
775                                                   765 
776     // verbose                                    766     // verbose
777     // G4cout << "G3Division::CreateSolid volu    767     // G4cout << "G3Division::CreateSolid volume after: " 
778     //        << fVTE->GetName() << " " << sha    768     //        << fVTE->GetName() << " " << shape << G4endl;    
779     // G4cout << " npar,Rpar: " << npar;          769     // G4cout << " npar,Rpar: " << npar;
780     // for (G4int iii = 0; iii < npar; ++iii)     770     // for (G4int iii = 0; iii < npar; ++iii) G4cout << " " << Rpar[iii];
781     // G4cout << G4endl;                          771     // G4cout << G4endl;
782     delete [] Rpar;                            << 
783 }                                                 772 }
784                                                   773 
785                                                   774 
786 G3VolTableEntry* G3Division::Dvn()                775 G3VolTableEntry* G3Division::Dvn()
787 {                                                 776 {   
788   // no envelope need to be created               777   // no envelope need to be created 
789                                                   778 
790   // get parameters from mother                   779   // get parameters from mother
791   G4String shape = fMVTE->GetShape();             780   G4String shape = fMVTE->GetShape(); 
792   G4double* Rpar = fMVTE->GetRpar();              781   G4double* Rpar = fMVTE->GetRpar();
793   G4int     npar = fMVTE->GetNpar();              782   G4int     npar = fMVTE->GetNpar();
794                                                   783 
795   // set width for replica and create solid       784   // set width for replica and create solid
796   fWidth = (fHighRange - fLowRange)/fNofDivisi    785   fWidth = (fHighRange - fLowRange)/fNofDivisions;
797   CreateSolid(shape, Rpar, npar);                 786   CreateSolid(shape, Rpar, npar);
798                                                   787 
799   return 0;                                       788   return 0;       
800 }                                                 789 }
801                                                   790 
802 G3VolTableEntry* G3Division::Dvn2()               791 G3VolTableEntry* G3Division::Dvn2()
803 {                                                 792 {
804   // to be defined as const of this class         793   // to be defined as const of this class
805   G4double Rmin = 0.0001*cm;                      794   G4double Rmin = 0.0001*cm;
806                                                   795 
807   G4String shape = fMVTE->GetShape();             796   G4String shape = fMVTE->GetShape();
808   G4double* Rpar = fMVTE->GetRpar();              797   G4double* Rpar = fMVTE->GetRpar();
809   G4int     npar = fMVTE->GetNpar();              798   G4int     npar = fMVTE->GetNpar();
810                                                   799 
811   G4double c0 = fC0;                              800   G4double c0 = fC0;
812   if (fAxis == kPhi)  c0 = c0*deg;                801   if (fAxis == kPhi)  c0 = c0*deg;
813   else                c0 = c0*cm;                 802   else                c0 = c0*cm;
814                                                   803           
815   // create envelope (if needed)                  804   // create envelope (if needed)
816   G3VolTableEntry* envVTE = 0;                    805   G3VolTableEntry* envVTE = 0;
817   if( std::abs(c0 - fLowRange) > Rmin) {          806   if( std::abs(c0 - fLowRange) > Rmin) {
818     envVTE = CreateEnvelope(shape, fHighRange,    807     envVTE = CreateEnvelope(shape, fHighRange, c0, Rpar, npar);
819     Rpar = envVTE->GetRpar();                     808     Rpar = envVTE->GetRpar();
820     npar = envVTE->GetNpar();                     809     npar = envVTE->GetNpar();
821   }                                               810   }  
822                                                   811 
823   // set width for replica and create solid       812   // set width for replica and create solid
824   fWidth = (fHighRange - c0)/fNofDivisions;       813   fWidth = (fHighRange - c0)/fNofDivisions;
825   CreateSolid(shape, Rpar, npar);                 814   CreateSolid(shape, Rpar, npar);
826                                                   815 
827   return envVTE;                                  816   return envVTE;
828 }                                                 817 }
829                                                   818 
830 G3VolTableEntry* G3Division::Dvt()                819 G3VolTableEntry* G3Division::Dvt()
831 {                                                 820 {
832   // to be defined as const of this class         821   // to be defined as const of this class
833   G4double Rmin = 0.0001*cm;                      822   G4double Rmin = 0.0001*cm;
834                                                   823 
835   // get parameters from mother                   824   // get parameters from mother
836   G4String shape = fMVTE->GetShape();             825   G4String shape = fMVTE->GetShape();
837   G4double* Rpar = fMVTE->GetRpar();              826   G4double* Rpar = fMVTE->GetRpar();
838   G4int     npar = fMVTE->GetNpar();              827   G4int     npar = fMVTE->GetNpar();
839                                                   828 
840   // calculate the number of divisions            829   // calculate the number of divisions    
841   G4int ndvmx = fNofDivisions;                    830   G4int ndvmx = fNofDivisions;
842   G4double step = fStep;                          831   G4double step = fStep;
843                                                   832   
844   if (fAxis == kPhi)  step = step*deg;            833   if (fAxis == kPhi)  step = step*deg;
845   else                step = step*cm;             834   else                step = step*cm;
846                                                   835 
847   G4int ndiv = G4int((fHighRange - fLowRange +    836   G4int ndiv = G4int((fHighRange - fLowRange + Rmin)/step);
848   // to be added warning                          837   // to be added warning
849   if (ndvmx > 255) ndvmx = 255;                   838   if (ndvmx > 255) ndvmx = 255;
850   if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx    839   if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx;
851                                                   840 
852   // create envVTE (if needed)                    841   // create envVTE (if needed)
853   G3VolTableEntry* envVTE = 0;                    842   G3VolTableEntry* envVTE = 0;
854   G4double delta = std::abs((fHighRange - fLow    843   G4double delta = std::abs((fHighRange - fLowRange) - ndiv*step);
855   if (delta > Rmin) {                             844   if (delta > Rmin) {
856     envVTE                                        845     envVTE 
857        = CreateEnvelope(shape, fHighRange-delt    846        = CreateEnvelope(shape, fHighRange-delta/2., fLowRange+delta/2., 
858                         Rpar, npar);              847                         Rpar, npar);
859     Rpar = envVTE->GetRpar();                     848     Rpar = envVTE->GetRpar();
860     npar = envVTE->GetNpar();                     849     npar = envVTE->GetNpar();
861   }                                               850   }
862                                                   851 
863   // set width for replica and create solid       852   // set width for replica and create solid
864   fWidth = step;                                  853   fWidth = step;
865   fNofDivisions = ndiv;                           854   fNofDivisions = ndiv;
866   CreateSolid(shape, Rpar, npar);                 855   CreateSolid(shape, Rpar, npar);
867                                                   856 
868   return envVTE;                                  857   return envVTE;
869 }                                                 858 }
870                                                   859 
871 G3VolTableEntry* G3Division::Dvt2()               860 G3VolTableEntry* G3Division::Dvt2()
872 {                                                 861 {
873   // to be defined as const of this class         862   // to be defined as const of this class
874   G4double Rmin = 0.0001*cm;                      863   G4double Rmin = 0.0001*cm;
875                                                   864 
876   // get parameters from mother                   865   // get parameters from mother
877   G4String shape = fMVTE->GetShape();             866   G4String shape = fMVTE->GetShape();
878   G4double* Rpar = fMVTE->GetRpar();              867   G4double* Rpar = fMVTE->GetRpar();
879   G4int     npar = fMVTE->GetNpar();              868   G4int     npar = fMVTE->GetNpar();
880                                                   869 
881   // calculate the number of divisions            870   // calculate the number of divisions   
882   G4int ndvmx = fNofDivisions;                    871   G4int ndvmx = fNofDivisions;
883   G4double step = fStep;                          872   G4double step = fStep;
884   G4double c0 = fC0;                              873   G4double c0 = fC0;
885                                                   874 
886   if(fAxis == kPhi){                              875   if(fAxis == kPhi){
887     step = step*deg;                              876     step = step*deg;
888     c0 = c0*deg;                                  877     c0 = c0*deg;
889   }                                               878   } 
890   else {                                          879   else {
891     step = step*cm;                               880     step = step*cm;
892     c0 = c0*cm;                                   881     c0 = c0*cm;
893   }                                               882   }  
894                                                   883 
895   G4int ndiv = G4int((fHighRange - c0 + Rmin)/    884   G4int ndiv = G4int((fHighRange - c0 + Rmin)/step);
896   // to be added warning                          885   // to be added warning
897   if (ndvmx > 255) ndvmx = 255;                   886   if (ndvmx > 255) ndvmx = 255;
898   if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx    887   if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx;
899                                                   888 
900   // create envelope (if needed)                  889   // create envelope (if needed)
901   G3VolTableEntry* envVTE = 0;                    890   G3VolTableEntry* envVTE = 0;
902   G4double delta = std::abs((fHighRange - c0)     891   G4double delta = std::abs((fHighRange - c0) - ndiv*step);
903   if (std::abs(c0 - fLowRange) > Rmin) {          892   if (std::abs(c0 - fLowRange) > Rmin) {
904     envVTE                                        893     envVTE 
905       = CreateEnvelope(shape, fHighRange-delta    894       = CreateEnvelope(shape, fHighRange-delta/2., c0+delta/2., Rpar, npar);
906     Rpar = envVTE->GetRpar();                     895     Rpar = envVTE->GetRpar();
907     npar = envVTE->GetNpar();                     896     npar = envVTE->GetNpar();
908   }                                               897   }
909                                                   898 
910   // set with for replica and create solid        899   // set with for replica and create solid
911   fWidth = step;                                  900   fWidth = step;
912   fNofDivisions = ndiv;                           901   fNofDivisions = ndiv;
913   CreateSolid(shape, Rpar, npar);                 902   CreateSolid(shape, Rpar, npar);
914                                                   903 
915   return envVTE;                                  904   return envVTE;   
916 }                                                 905 }
917                                                   906