G3Division.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // by I.Hrivnacova, V.Berejnoi 13.10.99
00030 
00031 #include <assert.h>
00032 
00033 #include "G3Division.hh"
00034 #include "G3VolTableEntry.hh"
00035 #include "G3toG4MakeSolid.hh"
00036 #include "G4Para.hh"
00037 #include "G3Pos.hh"
00038 #include "G4SystemOfUnits.hh"
00039 #include "G4LogicalVolume.hh"
00040 #include "G4VPhysicalVolume.hh"
00041 #include "G4PVPlacement.hh"
00042 #include "G4PVReplica.hh"
00043 #ifndef G3G4_NO_REFLECTION
00044 #include "G4ReflectionFactory.hh"
00045 #endif
00046 
00047 G3VolTableEntry* G4CreateVTE(G4String vname, G4String shape, G4int nmed,
00048                                G4double Rpar[], G4int npar);
00049 
00050 G3Division::G3Division(G3DivType type, G3VolTableEntry* vte, 
00051                 G3VolTableEntry* mvte, G4int nofDivisions, 
00052                 G4int iaxis, G4int nmed, G4double c0, G4double step)
00053   : fType(type),
00054     fVTE(vte),
00055     fMVTE(mvte),
00056     fNofDivisions(nofDivisions),
00057     fIAxis(iaxis),
00058     fNmed(nmed),
00059     fC0(c0),
00060     fStep(step),
00061     fLowRange(0.),
00062     fHighRange(0.),
00063     fWidth(0.),
00064     fOffset(0.),
00065     fAxis(kXAxis)
00066 {
00067   fVTE->SetHasNegPars(true);
00068 }
00069 
00070 G3Division::G3Division(G3VolTableEntry* vte, G3VolTableEntry* mvte,
00071                        const G3Division& division)
00072   : fVTE(vte),
00073     fMVTE(mvte)
00074 {
00075   // only "input" parameters are copied from division
00076   fType = division.fType;
00077   fNofDivisions = division.fNofDivisions;
00078   fIAxis = division.fIAxis;
00079   fNmed = division.fNmed;
00080   fC0 = division.fC0;
00081   fStep = division.fStep;
00082 
00083   // other parameters are set as in standard constructor
00084   fLowRange = 0.;
00085   fHighRange = 0.;
00086   fWidth = 0.;
00087   fOffset = 0.;
00088   fAxis = kXAxis;
00089   fVTE->SetHasNegPars(true);
00090 }
00091 
00092 G3Division::~G3Division()
00093 {}
00094 
00095 // public methods
00096 
00097 void G3Division::UpdateVTE()
00098 {
00099   if (fVTE->HasNegPars() && !(fMVTE->HasNegPars())) {
00100 
00101     // set nmed from mother 
00102     if (fNmed == 0) fNmed = fMVTE->GetNmed();
00103     fVTE->SetNmed(fNmed);
00104    
00105     SetRangeAndAxis();
00106     
00107     // create envelope (if necessary)
00108     // and solid    
00109     G3VolTableEntry* envVTE = 0;
00110     if      (fType == kDvn)  envVTE = Dvn(); 
00111     else if (fType == kDvn2) envVTE = Dvn2(); 
00112     else if (fType == kDvt)  envVTE = Dvt(); 
00113     else if (fType == kDvt2) envVTE = Dvt2();
00114     
00115     if (envVTE) {
00116       // reset mother <-> daughter
00117       fMVTE->ReplaceDaughter(fVTE, envVTE);
00118       fVTE->ReplaceMother(fMVTE, envVTE);
00119       envVTE->AddDaughter(fVTE);
00120       envVTE->AddMother(fMVTE);
00121 
00122       // replace mother with envelope
00123       fMVTE = envVTE;
00124     }
00125   }  
00126 }
00127 
00128 void G3Division::CreatePVReplica()
00129 {
00130   G4String name = fVTE->GetName();
00131   G4LogicalVolume* lv =  fVTE->GetLV();
00132   G4LogicalVolume* mlv = fMVTE->GetLV();
00133   
00134   G4String shape = fMVTE->GetShape();
00135   if (shape == "PARA") {
00136     // The para volume cannot be replicated using G4PVReplica.
00137     // (Replicating a volume along a cartesian axis means "slicing" it
00138     // with slices -perpendicular- to that axis.)
00139     
00140     // position the replicated elements    
00141     for (G4int i=0; i<fNofDivisions; i++) {
00142        G4ThreeVector position = G4ThreeVector(); 
00143        position[fIAxis-1] = fLowRange + fWidth/2. + i*fWidth;
00144        if (position.y()!=0.) 
00145          position.setX(position.y()*((G4Para*)lv->GetSolid())->GetTanAlpha());
00146 
00147        #ifndef G3G4_NO_REFLECTION
00148        G4ReflectionFactory::Instance()
00149          ->Place(G4Translate3D(position), name, lv, mlv, 0, i);
00150 
00151        #else  
00152        new G4PVPlacement(0, position, lv, name, mlv, 0, i);
00153 
00154        #endif
00155     }
00156     
00157     // G4PVReplica cannot be created
00158     return;   
00159   }     
00160   
00161   #ifdef G3G4DEBUG
00162     G4cout << "Create G4PVReplica name " << name << " logical volume name " 
00163            << lv->GetName() << " mother logical volme name "
00164            << mlv->GetName() << " axis " << fAxis << " ndivisions " 
00165            << fNofDivisions << " width " << fWidth << " Offset "
00166            << fOffset << G4endl;
00167   #endif
00168 
00169   #ifndef G3G4_NO_REFLECTION
00170   G4ReflectionFactory::Instance()
00171     ->Replicate(name, lv, mlv, fAxis, fNofDivisions, fWidth, fOffset);
00172 
00173   #else    
00174   new G4PVReplica(name, lv, mlv, fAxis, fNofDivisions, fWidth, fOffset);
00175 
00176   #endif
00177 }
00178 
00179 // private methods
00180 
00181 void G3Division::Exception(G4String where, G4String what)
00182 {
00183   G4String err_message = "G3Division::" + where + " for "
00184                        + what + " is not implemented";
00185   G4Exception("G3Division::Exception()", "G3toG40004",
00186               FatalException, err_message);
00187   return;
00188 }  
00189 
00190 void G3Division::SetRangeAndAxis()
00191 // set fHighRange, fLowRange, fAxis
00192 {
00193     G4String shape = fMVTE->GetShape();
00194     G4double *Rpar = fMVTE->GetRpar();
00195     
00196     switch (fIAxis) {
00197       case 1: fAxis = kXAxis;
00198               break;
00199       case 2: fAxis = kYAxis;
00200               break;
00201       case 3: fAxis = kZAxis;
00202               break;
00203       default: G4Exception("G3Division::SetRangeAndAxis()", "G3toG40005",
00204                             FatalException, "Wrong iaxis defenition!");
00205     }
00206 
00207     if ( shape == "BOX" ) {
00208       fHighRange = Rpar[fIAxis-1]*cm;
00209       fLowRange = -fHighRange;
00210     }
00211     else if ( shape == "TRD1" ) {
00212       if (fIAxis == 1){
00213         fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm);
00214       }
00215       else if( fIAxis == 2) {
00216        fHighRange = Rpar[2]*cm;
00217       }
00218       else if( fIAxis == 3) {
00219        fHighRange = Rpar[3]*cm;
00220       }
00221       fLowRange = - fHighRange;
00222     }
00223     else if ( shape == "TRD2" ) {
00224       if (fIAxis == 1){
00225         fHighRange = std::max(Rpar[0]*cm, Rpar[1]*cm);
00226       }
00227       else if( fIAxis == 2) {
00228         fHighRange = std::max(Rpar[2]*cm, Rpar[3]*cm);
00229       }
00230       else if( fIAxis == 3) {
00231        fHighRange = Rpar[4]*cm;
00232       }
00233     }
00234     else if ( shape == "TRAP" ) {
00235       if ( fIAxis == 3 ) fHighRange = Rpar[0]*cm;
00236       else               fHighRange = 0.;
00237       fLowRange = -fHighRange;
00238     }
00239     else if ( shape == "TUBE" ) {
00240       if (fIAxis == 1){
00241         fHighRange = Rpar[1]*cm;
00242         fLowRange = Rpar[0]*cm;
00243         fAxis = kRho;
00244       }
00245       else if( fIAxis == 2) {
00246         fHighRange = 360.*deg;
00247         fLowRange = 0.;
00248         fAxis = kPhi;
00249       }
00250       else if( fIAxis == 3) {
00251        fHighRange = Rpar[2]*cm;
00252        fLowRange = -fHighRange;
00253       }
00254     }
00255     else if ( shape == "TUBS" ) {
00256       if (fIAxis == 1){
00257         fHighRange = Rpar[1]*cm;
00258         fLowRange = Rpar[0]*cm;
00259         fAxis = kRho;
00260       }
00261       else if( fIAxis == 2) {
00262 
00263        fLowRange = Rpar[3]*deg;
00264        fHighRange = Rpar[4]*deg - fLowRange;
00265        if ( Rpar[4]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg;
00266        fHighRange = fHighRange + fLowRange;
00267        fAxis = kPhi;
00268       }
00269       else if( fIAxis == 3) {
00270        fHighRange = Rpar[2]*cm;
00271        fLowRange = -fHighRange;
00272       }
00273     }
00274     else if ( shape == "CONE" ) {
00275       if (fIAxis == 1){
00276         fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm);
00277         fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm);
00278         fAxis = kRho;
00279       }
00280       else if( fIAxis == 2) {
00281 
00282        fLowRange = 0.;
00283        fHighRange = 360.*deg;
00284        fAxis = kPhi;
00285       }
00286       else if( fIAxis == 3) {
00287        fHighRange = Rpar[0]*cm;
00288        fLowRange = -fHighRange;
00289       }
00290     }
00291     else if ( shape == "CONS" ) {
00292       if (fIAxis == 1){
00293         fHighRange = std::max(Rpar[2]*cm,Rpar[4]*cm);
00294         fLowRange = std::max(Rpar[1]*cm,Rpar[3]*cm);
00295         fAxis = kRho;
00296       }
00297       else if( fIAxis == 2) {
00298 
00299        fLowRange = Rpar[5]*deg;
00300        fHighRange = Rpar[6]*deg - fLowRange;
00301        if ( Rpar[6]*deg <= fLowRange )fHighRange = fHighRange + 360.*deg;
00302        fHighRange = fHighRange + fLowRange;
00303        fAxis = kPhi;
00304       }
00305       else if( fIAxis == 3) {
00306        fHighRange = Rpar[2]*cm;
00307        fLowRange = -fHighRange;
00308       }
00309     }
00310     else if ( shape == "SPHE" ) {
00311       if (fIAxis == 1){
00312         fHighRange = Rpar[1]*cm;
00313         fLowRange = Rpar[0]*cm;
00314         fAxis = kRho;
00315       }
00316       else if( fIAxis == 2) {
00317        fLowRange = std::min(Rpar[2]*deg,Rpar[3]*deg);
00318        fHighRange = std::max(Rpar[2]*deg,Rpar[3]*deg);
00319        fAxis = kPhi;
00320       }
00321       else if( fIAxis == 3) {
00322        fLowRange = std::min(Rpar[4]*deg,Rpar[5]*deg);
00323        fHighRange = std::max(Rpar[4]*deg,Rpar[5]*deg);
00324        fAxis = kPhi; // ?????? 
00325       }
00326     }
00327     else if ( shape == "PARA" ) {
00328       fHighRange = Rpar[fIAxis-1]*cm;
00329       fLowRange = -fHighRange;
00330     }
00331     else if ( shape == "PGON" ) {
00332         G4int i;
00333         G4int nz = G4int(Rpar[3]);
00334 
00335         G4double pPhi1 = Rpar[0]*deg;
00336         G4double dPhi  = Rpar[1]*deg;
00337     
00338         G4double *DzArray = new G4double[nz];
00339         G4double *Rmax    = new G4double[nz];
00340         G4double *Rmin    = new G4double[nz];
00341         G4double rangehi[3], rangelo[3];
00342         rangehi[0] = -kInfinity  ;
00343         rangelo[0] =  kInfinity ;
00344         rangehi[2] = -kInfinity ;
00345         rangelo[2] =  kInfinity ;
00346 
00347         for(i=0; i<nz; i++) 
00348         {
00349             G4int i4=3*i+4;
00350             G4int i5=i4+1;
00351             G4int i6=i4+2;
00352             
00353             DzArray[i] = Rpar[i4]*cm;
00354             Rmin[i] = Rpar[i5]*cm;
00355             Rmax[i] = Rpar[i6]*cm;
00356             rangelo[0] = std::min(rangelo[0], Rmin[i]);
00357             rangehi[0] = std::max(rangehi[0], Rmax[i]);
00358             rangelo[2] = std::min(rangelo[2], DzArray[i]);
00359             rangehi[2] = std::max(rangehi[2], DzArray[i]);
00360         }
00361         for (i=0;i<nz;i++){
00362             assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]);
00363         }
00364         rangehi[1] = pPhi1 + dPhi;
00365         rangelo[1] = pPhi1;
00366         fHighRange = rangehi[fIAxis-1];
00367         fLowRange = rangelo[fIAxis-1];
00368         if      (fIAxis == 1)fAxis = kRho;
00369         else if (fIAxis == 2)fAxis = kPhi;
00370         else if (fIAxis == 3)fAxis = kZAxis;
00371 
00372         delete [] DzArray;
00373         delete [] Rmin;
00374         delete [] Rmax;
00375 
00376     }
00377     else if ( shape == "PCON" ) {
00378 
00379         G4int i;
00380         G4double pPhi1 = Rpar[0]*deg;
00381         G4double dPhi  = Rpar[1]*deg;    
00382         G4int nz = G4int(Rpar[2]);
00383     
00384         G4double *DzArray = new G4double[nz];
00385         G4double *Rmax    = new G4double[nz];
00386         G4double *Rmin    = new G4double[nz];
00387         G4double rangehi[3],rangelo[3];
00388 
00389         rangehi[0] = -kInfinity  ;
00390         rangelo[0] =  kInfinity ;
00391         rangehi[2] = -kInfinity ;
00392         rangelo[2] =  kInfinity ;
00393         
00394         for(i=0; i<nz; i++){
00395             G4int i4=3*i+3;
00396             G4int i5=i4+1;
00397             G4int i6=i4+2;
00398             
00399             DzArray[i] = Rpar[i4]*cm;
00400             Rmin[i] = Rpar[i5]*cm;
00401             Rmax[i] = Rpar[i6]*cm;
00402             rangelo[0] = std::min(rangelo[0], Rmin[i]);
00403             rangehi[0] = std::max(rangehi[0], Rmax[i]);
00404             rangelo[2] = std::min(rangelo[2], DzArray[i]);
00405             rangehi[2] = std::max(rangehi[2], DzArray[i]);
00406         }
00407         for (i=0;i<nz;i++){
00408             assert(Rmin[i]>=0 && Rmax[i]>=Rmin[i]);
00409         }
00410         rangehi[1] = pPhi1 + dPhi;
00411         rangelo[1] = pPhi1;
00412         fHighRange = rangehi[fIAxis-1];
00413         fLowRange = rangelo[fIAxis-1];
00414         if      (fIAxis == 1)fAxis = kRho;
00415         else if (fIAxis == 2)fAxis = kPhi;
00416         else if (fIAxis == 3)fAxis = kZAxis;
00417 
00418 
00419         delete [] DzArray;
00420         delete [] Rmin;
00421         delete [] Rmax;
00422     }
00423     else if ( shape == "ELTU" ||  shape == "HYPE" || shape == "GTRA" ||
00424          shape == "CTUB") {
00425        Exception("SetRangeAndAxis", shape);
00426     }
00427     else {
00428        Exception("SetRangeAndAxis", "Unknown shape" + shape);
00429     }  
00430 
00431     // verbose
00432     #ifdef G3G4DEBUG
00433       G4cout << "Shape " << shape << " SetRangeAndAxis: " 
00434              << fLowRange << " " << fHighRange << " " << fAxis << G4endl;
00435     #endif
00436 }
00437 
00438 G3VolTableEntry* G3Division::CreateEnvelope(G4String shape, G4double hi, 
00439                                G4double lo, G4double par[], G4int npar)
00440 // create new VTE with G3Pos corresponding to the
00441 // envelope of divided volume
00442 {
00443     // verbose
00444     // G4cout << "  G3Division::CreateEnvelope " << "fIAaxis= " << fIAxis
00445     //        << " hi= " << hi
00446     //        << " lo= " << lo
00447     //        << G4endl;
00448 
00449     G4double *Rpar = new G4double[npar+2];
00450     for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];}
00451     G4double pos[3] = {0.,0.,0.};
00452   
00453     if ( shape == "BOX" ) {
00454       Rpar[fIAxis-1] = (hi - lo)/2./cm;
00455       pos [fIAxis-1] = (hi + lo)/2.;
00456     }
00457     else if ( shape == "TRD1" ) {
00458       if ( fIAxis == 1 || fIAxis == 2  ) {
00459         Exception("CreateEnvelope","TRD1-x,y");
00460       }
00461       else if ( fIAxis == 3 ) {
00462         // x = x1 + (c-z1)(x2 -x1)/(z2-z1)
00463         G4double tn, x1, z1;
00464         tn = (Rpar[1] - Rpar[0])/(2.* Rpar[3]); 
00465         x1 = Rpar[0]; z1 = -Rpar[3];
00466         Rpar[0] = x1 + tn * (lo/cm - z1);
00467         Rpar[1] = x1 + tn * (hi/cm - z1);
00468         Rpar[3] = (hi - lo)/2./cm;
00469         pos[2]  = (hi + lo)/2.;
00470       }
00471     }
00472     else if ( shape == "TRD2" ) {
00473       if ( fIAxis == 1 || fIAxis == 2) {
00474         Exception("CreateEnvelope","TRD2-x,y");
00475       }
00476       else if ( fIAxis == 3 ) {
00477         // x = x1 + (c-z1)(x2 -x1)/(z2-z1)
00478         // y = y1 + (c-z1)(y2 -y1)/(z2-z1)
00479         G4double tn1, tn2, x1, y1, z1;
00480         tn1 = (Rpar[1] - Rpar[0])/(2.* Rpar[4]); 
00481         tn2 = (Rpar[3] - Rpar[2])/(2.* Rpar[4]); 
00482         x1 = Rpar[0]; y1 = Rpar[2]; z1 = -Rpar[3];
00483         Rpar[0] = x1 + tn1 * (lo/cm - z1);
00484         Rpar[1] = x1 + tn1 * (hi/cm - z1);
00485         Rpar[2] = y1 + tn2 * (lo/cm - z1);
00486         Rpar[3] = y1 + tn2 * (hi/cm - z1);
00487         Rpar[4] = (hi - lo)/2./cm;
00488         pos[2]  = (hi + lo)/2.;
00489       }
00490     }
00491     else if ( shape == "TRAP" ) {
00492       Exception("CreateEnvelope","TRAP-x,y,z");
00493     }
00494     else if ( shape == "TUBE" ) {
00495       if ( fIAxis == 1 ) {
00496         Rpar[0] = lo/cm;
00497         Rpar[1] = hi/cm;
00498       }
00499       else if ( fIAxis == 2 ) {
00500         Rpar[3] = lo/deg;
00501         Rpar[4] = hi/deg;
00502         npar = npar + 2;
00503         shape = "TUBS";
00504       }
00505       else if ( fIAxis == 3 ) {
00506         Rpar[2] = (hi - lo)/2./cm;
00507         pos [2] = (hi + lo)/2.;
00508       }
00509     }
00510     else if ( shape == "TUBS" ) {
00511       if ( fIAxis == 1 ) {
00512         Rpar[0] = lo/cm;
00513         Rpar[1] = hi/cm;
00514       }
00515       else if ( fIAxis == 2 ) {
00516         Rpar[3] = lo/deg;
00517         Rpar[4] = hi/deg;
00518       }
00519       else if ( fIAxis == 3 ) {
00520         Rpar[2] = (hi - lo)/2./cm;
00521         pos [2] = (hi + lo)/2.;
00522       }
00523     }
00524     else if ( shape == "CONE" ) {
00525       if ( fIAxis == 1) {
00526         Exception("CreateEnvelope","CONE-x,z");
00527       }
00528       else if ( fIAxis == 2 ) {
00529         Rpar[5] = lo/deg;
00530         Rpar[6] = hi/deg;
00531         npar = npar + 2;
00532         shape = "CONS";
00533       }
00534       else if ( fIAxis == 3 ) {
00535         G4double tn1, tn2, rmin, rmax, z1;
00536         tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); 
00537         tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); 
00538         rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0];
00539         Rpar[1] = rmin + tn1 * (lo/cm - z1);
00540         Rpar[3] = rmin + tn1 * (hi/cm - z1);
00541         Rpar[2] = rmax + tn2 * (lo/cm - z1);
00542         Rpar[4] = rmax + tn2 * (hi/cm - z1);
00543         Rpar[0] = (hi - lo)/2./cm;
00544         pos[2]  = (hi + lo)/2.;
00545       }
00546     }
00547     else if ( shape == "CONS" ) {
00548       if ( fIAxis == 1 ) {
00549         Exception("CreateEnvelope","CONS-x");
00550       }
00551       else if ( fIAxis == 2 ) {
00552         Rpar[5] = lo/deg;
00553         Rpar[6] = hi/deg;
00554       }
00555       else if ( fIAxis == 3 ) {
00556         G4double tn1, tn2, rmin, rmax, z1;
00557         tn1 = (Rpar[3] - Rpar[1])/(2.* Rpar[0]); 
00558         tn2 = (Rpar[4] - Rpar[2])/(2.* Rpar[0]); 
00559         rmin = Rpar[1]; rmax = Rpar[2]; z1 = -Rpar[0];
00560         Rpar[1] = rmin + tn1 * (lo/cm - z1);
00561         Rpar[3] = rmin + tn1 * (hi/cm - z1);
00562         Rpar[2] = rmax + tn2 * (lo/cm - z1);
00563         Rpar[4] = rmax + tn2 * (hi/cm - z1);
00564         Rpar[0] = (hi - lo)/2./cm;
00565         pos[2]  = (hi + lo)/2.;
00566       }
00567     }
00568     else if ( shape == "SPHE" ) {
00569       Exception("CreateEnvelope","SPHE-x,y,z");                
00570     }
00571     else if ( shape == "PARA" ) {
00572       Exception("CreateEnvelope","PARA-x,y,z");
00573     }
00574     else if ( shape == "PGON" ) {
00575       if ( fIAxis == 2) {
00576         Rpar[0] = lo/deg;
00577         Rpar[1] = hi/deg;
00578         // rotm = ???
00579       }
00580       else {
00581         Exception("CreateEnvelope","PGON-x,z");
00582       }
00583     }
00584     else if ( shape == "PCON" ) {
00585       if ( fIAxis == 2) {
00586         Rpar[0] = lo/deg;
00587         Rpar[1] = hi/deg;
00588         // rotm = ???
00589       }
00590       else {
00591         Exception("CreateEnvelope","PCON-x,z");
00592       }
00593     }
00594     else {
00595        Exception("CreateEnvelope", "Unknown shape" + shape);
00596     }  
00597 
00598     // create new VTE corresponding to envelope
00599     G4String envName = fVTE->GetName() + "_ENV"; 
00600     G3VolTableEntry* envVTE 
00601       = G4CreateVTE(envName, shape, fNmed, Rpar, npar);
00602 
00603     // create a G3Pos object and add it to envVTE
00604     G4String motherName = fMVTE->GetMasterClone()->GetName();
00605     G4ThreeVector* offset = new G4ThreeVector(pos[0],pos[1],pos[2]);    
00606     G4String only = "ONLY";
00607     G3Pos* aG3Pos = new G3Pos(motherName, 1, offset, 0, only);              
00608     envVTE->AddG3Pos(aG3Pos);
00609 
00610     delete [] Rpar; 
00611 
00612     return envVTE;
00613 }
00614 
00615 void G3Division::CreateSolid(G4String shape, G4double par[], G4int npar)
00616 // create the solid corresponding to divided volume
00617 // and set the fOffset for replica
00618 {
00619     G4double *Rpar = new G4double[npar+2];
00620     for (G4int i=0; i<npar; ++i){ Rpar[i] = par[i];}
00621 
00622     // verbose
00623     // G4cout << "G3Division::CreateSolid volume before: " 
00624     //        << fVTE->GetName() << " " << shape << G4endl;    
00625     // G4cout << " npar,Rpar: " << npar;
00626     // for (G4int ii = 0; ii < npar; ++ii) G4cout << " " << Rpar[ii];
00627     // G4cout << G4endl;
00628   
00629     if ( shape == "BOX" ) {
00630       if      ( fIAxis == 1 ) Rpar[0] = fWidth/2./cm;
00631       else if ( fIAxis == 2 ) Rpar[1] = fWidth/2./cm; 
00632       else if ( fIAxis == 3 ) Rpar[2] = fWidth/2./cm; 
00633     }
00634     else if ( shape == "TRD1" ) {
00635       if ( fIAxis == 1 || fIAxis == 2 ) {
00636         Exception("CreateSolid", "TRD1-x,y");
00637       }
00638       else if ( fIAxis == 3 ) {
00639          Rpar[3] = fWidth/2./cm; 
00640       }
00641     }
00642     else if ( shape == "TRD2" ) {
00643       if ( fIAxis == 1 || fIAxis == 2 ) {
00644         Exception("CreateSolid", "TRD2-x,y");
00645       }
00646       else if ( fIAxis == 3 ) {
00647          Rpar[4] =  fWidth/2./cm; 
00648       }
00649     }
00650     else if ( shape == "TRAP" ) {
00651       if ( fIAxis == 1 || fIAxis == 2) {
00652         Exception("CreateSolid", "TRAP-x,y");
00653       }
00654       else if ( fIAxis == 3 ) {
00655          Rpar[0] =  fWidth/2./cm; 
00656       }
00657     }
00658     else if ( shape == "TUBE" ) {
00659       if ( fIAxis == 1 ) {
00660          Rpar[1] = Rpar[0] + fWidth/cm;
00661          fOffset = Rpar[0]*cm;
00662       }
00663       else if ( fIAxis == 2 ) {
00664          Rpar[3] = 0.; 
00665          Rpar[4] = fWidth/deg; 
00666          shape = "TUBS";
00667          npar = npar + 2;
00668       }
00669       else if ( fIAxis == 3 ) {
00670          Rpar[2] = fWidth/2./cm; 
00671       }
00672     }
00673     else if ( shape == "TUBS" ) {
00674       if ( fIAxis == 1 ) {
00675         Rpar[1] = Rpar[0] + fWidth/cm;
00676         fOffset = Rpar[0]*cm;
00677       }
00678       else if ( fIAxis == 2 ) {
00679          fOffset = Rpar[3]*deg; 
00680          Rpar[3] = 0.;
00681          Rpar[4] =  fWidth/deg;
00682       }
00683       else if ( fIAxis == 3 ) {
00684          Rpar[2] = fWidth/2./cm; 
00685       }
00686     }
00687     else if ( shape == "CONE" ) {
00688       if ( fIAxis == 1 ) {
00689         Exception("CreateSolid", "CONE-x"); 
00690       }
00691       else if ( fIAxis == 2 ) {
00692          Rpar[5] = 0.;
00693          Rpar[6] = fWidth/deg;
00694          shape = "CONS";
00695          npar = npar + 2;
00696       }
00697       else if ( fIAxis == 3 ) {
00698          Rpar[0] = fWidth/2./cm; 
00699       }
00700     }
00701     else if ( shape == "CONS" ) {
00702       if ( fIAxis == 1 ) {
00703         Exception("CreateSolid", "CONS-x"); 
00704       }
00705       else if ( fIAxis == 2 ) {
00706          fOffset = Rpar[5]*deg;
00707          Rpar[5] = 0.;
00708          Rpar[6] = fWidth/deg;
00709       }
00710       else if ( fIAxis == 3 ) {
00711          Rpar[0] = fWidth/2./cm; 
00712       }
00713     }
00714     else if (shape == "PARA") {
00715       if      ( fIAxis == 1 ) {
00716          Rpar[0] = fWidth/2./cm;
00717       }  
00718       else if ( Rpar[4] == 0. && Rpar[5] == 0. ) {
00719          // only special case for axis 2,3 is supported
00720         if ( fIAxis == 2 ) {
00721           Rpar[1] = fWidth/2./cm;
00722         }  
00723         else if ( fIAxis == 3) {
00724           Rpar[2] = fWidth/2./cm;
00725         }
00726       }   
00727       else    
00728          Exception("CreateSolid", shape);
00729     }    
00730     else if (shape == "SPHE") {
00731       Exception("CreateSolid", shape);
00732     }
00733     else if ( shape == "PGON" ) {
00734       if ( fIAxis == 2 ) {
00735          fOffset = Rpar[0]*deg;
00736          Rpar[0] = 0.;
00737          Rpar[1] = fWidth/deg;
00738          Rpar[2] = 1.;
00739       }
00740       else
00741        Exception("CreateSolid", shape);
00742     }
00743     else if ( shape == "PCON" ) {
00744       if ( fIAxis == 2 ) {
00745          fOffset = Rpar[0]*deg;
00746          Rpar[0] = 0.;
00747          Rpar[1] = fWidth/deg;
00748       }
00749       else {
00750         Exception("CreateSolid", shape);
00751       } 
00752     }
00753     else {
00754        Exception("CreateSolid", "Unknown shape" + shape);
00755     }  
00756 
00757     // create solid and set it to fVTE
00758     G4bool hasNegPars;
00759     G4bool deferred;   
00760     G4bool okAxis[3];
00761     G4VSolid* solid
00762     = G3toG4MakeSolid(fVTE->GetName(), shape, Rpar, npar, hasNegPars, deferred, okAxis);  
00763 
00764     if (hasNegPars) {
00765        G4String err_message = "CreateSolid VTE " + fVTE->GetName()
00766                             + " has negative parameters.";
00767        G4Exception("G3Division::CreateSolid()", "G3toG40006",
00768                    FatalException, err_message);
00769        return;
00770     }   
00771     
00772     // update vte
00773     fVTE->SetSolid(solid);
00774     fVTE->SetNRpar(npar, Rpar); 
00775     fVTE->SetHasNegPars(hasNegPars);
00776 
00777     // verbose
00778     // G4cout << "G3Division::CreateSolid volume after: " 
00779     //        << fVTE->GetName() << " " << shape << G4endl;    
00780     // G4cout << " npar,Rpar: " << npar;
00781     // for (G4int iii = 0; iii < npar; ++iii) G4cout << " " << Rpar[iii];
00782     // G4cout << G4endl;
00783     delete [] Rpar;
00784 }
00785 
00786 
00787 G3VolTableEntry* G3Division::Dvn()
00788 {   
00789   // no envelope need to be created 
00790 
00791   // get parameters from mother
00792   G4String shape = fMVTE->GetShape(); 
00793   G4double* Rpar = fMVTE->GetRpar();
00794   G4int     npar = fMVTE->GetNpar();
00795 
00796   // set width for replica and create solid
00797   fWidth = (fHighRange - fLowRange)/fNofDivisions;
00798   CreateSolid(shape, Rpar, npar);
00799 
00800   return 0;                       
00801 }
00802 
00803 G3VolTableEntry* G3Division::Dvn2()
00804 {
00805   // to be defined as const of this class
00806   G4double Rmin = 0.0001*cm;
00807 
00808   G4String shape = fMVTE->GetShape();
00809   G4double* Rpar = fMVTE->GetRpar();
00810   G4int     npar = fMVTE->GetNpar();
00811 
00812   G4double c0 = fC0;
00813   if (fAxis == kPhi)  c0 = c0*deg;
00814   else                c0 = c0*cm;
00815           
00816   // create envelope (if needed)
00817   G3VolTableEntry* envVTE = 0;
00818   if( std::abs(c0 - fLowRange) > Rmin) {
00819     envVTE = CreateEnvelope(shape, fHighRange, c0, Rpar, npar);
00820     Rpar = envVTE->GetRpar();
00821     npar = envVTE->GetNpar();
00822   }  
00823 
00824   // set width for replica and create solid
00825   fWidth = (fHighRange - c0)/fNofDivisions;
00826   CreateSolid(shape, Rpar, npar);
00827 
00828   return envVTE;
00829 }
00830 
00831 G3VolTableEntry* G3Division::Dvt()
00832 {
00833   // to be defined as const of this class
00834   G4double Rmin = 0.0001*cm;
00835 
00836   // get parameters from mother
00837   G4String shape = fMVTE->GetShape();
00838   G4double* Rpar = fMVTE->GetRpar();
00839   G4int     npar = fMVTE->GetNpar();
00840 
00841   // calculate the number of divisions    
00842   G4int ndvmx = fNofDivisions;
00843   G4double step = fStep;
00844   
00845   if (fAxis == kPhi)  step = step*deg;
00846   else                step = step*cm;
00847 
00848   G4int ndiv = G4int((fHighRange - fLowRange + Rmin)/step);
00849   // to be added warning
00850   if (ndvmx > 255) ndvmx = 255;
00851   if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx;
00852 
00853   // create envVTE (if needed)
00854   G3VolTableEntry* envVTE = 0;
00855   G4double delta = std::abs((fHighRange - fLowRange) - ndiv*step);
00856   if (delta > Rmin) {
00857     envVTE 
00858        = CreateEnvelope(shape, fHighRange-delta/2., fLowRange+delta/2., 
00859                         Rpar, npar);
00860     Rpar = envVTE->GetRpar();
00861     npar = envVTE->GetNpar();
00862   }
00863 
00864   // set width for replica and create solid
00865   fWidth = step;
00866   fNofDivisions = ndiv;
00867   CreateSolid(shape, Rpar, npar);
00868 
00869   return envVTE;
00870 }
00871 
00872 G3VolTableEntry* G3Division::Dvt2()
00873 {
00874   // to be defined as const of this class
00875   G4double Rmin = 0.0001*cm;
00876 
00877   // get parameters from mother
00878   G4String shape = fMVTE->GetShape();
00879   G4double* Rpar = fMVTE->GetRpar();
00880   G4int     npar = fMVTE->GetNpar();
00881 
00882   // calculate the number of divisions   
00883   G4int ndvmx = fNofDivisions;
00884   G4double step = fStep;
00885   G4double c0 = fC0;
00886 
00887   if(fAxis == kPhi){
00888     step = step*deg;
00889     c0 = c0*deg;
00890   } 
00891   else {
00892     step = step*cm;
00893     c0 = c0*cm;
00894   }  
00895 
00896   G4int ndiv = G4int((fHighRange - c0 + Rmin)/step);
00897   // to be added warning
00898   if (ndvmx > 255) ndvmx = 255;
00899   if (ndiv > ndvmx && ndvmx > 0 ) ndiv = ndvmx;
00900 
00901   // create envelope (if needed)
00902   G3VolTableEntry* envVTE = 0;
00903   G4double delta = std::abs((fHighRange - c0) - ndiv*step);
00904   if (std::abs(c0 - fLowRange) > Rmin) {
00905     envVTE 
00906       = CreateEnvelope(shape, fHighRange-delta/2., c0+delta/2., Rpar, npar);
00907     Rpar = envVTE->GetRpar();
00908     npar = envVTE->GetNpar();
00909   }
00910 
00911   // set with for replica and create solid
00912   fWidth = step;
00913   fNofDivisions = ndiv;
00914   CreateSolid(shape, Rpar, npar);
00915 
00916   return envVTE;         
00917 }

Generated on Mon May 27 17:47:35 2013 for Geant4 by  doxygen 1.4.7