source: trunk/source/g3tog4/src/G3Division.cc @ 1202

Last change on this file since 1202 was 965, checked in by garnier, 15 years ago

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