CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_6/src/Fireworks/Geometry/src/TGeoMgrFromDdd.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:     Geometry
00004 // Class  :     TGeoMgrFromDdd
00005 // 
00006 // Implementation:
00007 //     [Notes on implementation]
00008 //
00009 // Original Author:  
00010 //         Created:  Fri Jul  2 16:11:42 CEST 2010
00011 // $Id: TGeoMgrFromDdd.cc,v 1.13 2012/05/30 14:02:07 yana Exp $
00012 //
00013 
00014 #include "Fireworks/Geometry/interface/TGeoMgrFromDdd.h"
00015 
00016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00017 #include "FWCore/ServiceRegistry/interface/ActivityRegistry.h"
00018 #include "FWCore/Framework/interface/EventSetup.h"
00019 #include "FWCore/Framework/interface/ESTransientHandle.h"
00020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00021 
00022 #include "DetectorDescription/Core/interface/DDCompactView.h"
00023 #include "DetectorDescription/Core/interface/DDSolid.h"
00024 #include "DetectorDescription/Core/interface/DDMaterial.h"
00025 
00026 #include "Geometry/Records/interface/IdealGeometryRecord.h"
00027 
00028 #include "Fireworks/Geometry/interface/DisplayGeomRecord.h"
00029 
00030 #include "TGeoManager.h"
00031 #include "TGeoMatrix.h"
00032 #include "TGeoCompositeShape.h"
00033 #include "TGeoPcon.h"
00034 #include "TGeoPgon.h"
00035 #include "TGeoCone.h"
00036 #include "TGeoBoolNode.h"
00037 #include "TGeoTube.h"
00038 #include "TGeoArb8.h"
00039 #include "TGeoTrd2.h"
00040 #include "TGeoTorus.h"
00041 
00042 #include "Math/GenVector/RotationX.h"
00043 #include "Math/GenVector/RotationZ.h"
00044 
00045 #include "CLHEP/Units/GlobalSystemOfUnits.h"
00046 #include <math.h>
00047 
00048 TGeoMgrFromDdd::TGeoMgrFromDdd( const edm::ParameterSet& pset )
00049   : m_level( pset.getUntrackedParameter<int> ( "level", 10 )),
00050     m_verbose( pset.getUntrackedParameter<bool>( "verbose", false ))
00051 {
00052    // The following line is needed to tell the framework what data is
00053    // being produced.
00054    setWhatProduced( this );
00055 }
00056 
00057 TGeoMgrFromDdd::~TGeoMgrFromDdd( void )
00058 {}
00059 
00060 //==============================================================================
00061 // Local helpers
00062 //==============================================================================
00063 
00064 namespace
00065 {
00066    TGeoCombiTrans* createPlacement(const DDRotationMatrix& iRot,
00067                                    const DDTranslation&    iTrans)
00068    {
00069       //  std::cout << "in createPlacement" << std::endl;
00070       double elements[9];
00071       iRot.GetComponents(elements);
00072       TGeoRotation r;
00073       r.SetMatrix(elements);
00074 
00075       TGeoTranslation t(iTrans.x()/cm,
00076                         iTrans.y()/cm,
00077                         iTrans.z()/cm);
00078 
00079       return new TGeoCombiTrans(t,r);
00080    }
00081 }
00082 
00083 
00084 //==============================================================================
00085 // public member functions
00086 //==============================================================================
00087 
00088 TGeoMgrFromDdd::ReturnType
00089 TGeoMgrFromDdd::produce(const DisplayGeomRecord& iRecord)
00090 {
00091    using namespace edm;
00092 
00093    ESTransientHandle<DDCompactView> viewH;
00094    iRecord.getRecord<IdealGeometryRecord>().get(viewH);
00095 
00096    if ( ! viewH.isValid()) {
00097       return boost::shared_ptr<TGeoManager>();
00098    }
00099 
00100    TGeoManager *geo_mgr = new TGeoManager("cmsGeo","CMS Detector");
00101    // NOTE: the default constructor does not create the identity matrix
00102    if (gGeoIdentity == 0)
00103    {
00104       gGeoIdentity = new TGeoIdentity("Identity");
00105    }
00106 
00107    std::cout << "about to initialize the DDCompactView walker" << std::endl;
00108    DDCompactView::walker_type             walker(viewH->graph());
00109    DDCompactView::walker_type::value_type info = walker.current();
00110 
00111    // The top most item is actually the volume holding both the
00112    // geometry AND the magnetic field volumes!
00113    walker.firstChild();
00114    if( ! walker.firstChild()) {
00115       return boost::shared_ptr<TGeoManager>();
00116    }
00117 
00118    TGeoVolume *top = createVolume(info.first.name().fullname(),
00119                                   info.first.solid(),
00120                                   info.first.material());
00121    if (top == 0) {
00122       return boost::shared_ptr<TGeoManager>();
00123    }
00124 
00125    geo_mgr->SetTopVolume(top);
00126    // ROOT chokes unless colors are assigned
00127    top->SetVisibility(kFALSE);
00128    top->SetLineColor(kBlue);
00129 
00130    std::vector<TGeoVolume*> parentStack;
00131    parentStack.push_back(top);
00132 
00133    do
00134    {
00135       DDCompactView::walker_type::value_type info = walker.current();
00136 
00137       if (m_verbose)
00138       {
00139          for(unsigned int i=0; i<parentStack.size();++i) {
00140             std::cout <<" ";
00141          }
00142          std::cout << info.first.name()<<" "<<info.second->copyno_<<" "
00143                    << DDSolidShapesName::name(info.first.solid().shape())<<std::endl;
00144       }
00145 
00146       bool childAlreadyExists = (0 != nameToVolume_[info.first.name().fullname()]);
00147       TGeoVolume *child = createVolume(info.first.name().fullname(),
00148                                        info.first.solid(),
00149                                        info.first.material());
00150       if (0!=child && info.second != 0)
00151       {
00152          parentStack.back()->AddNode(child,
00153                                  info.second->copyno_,
00154                                  createPlacement(info.second->rotation(),
00155                                                  info.second->translation()));
00156          child->SetLineColor(kBlue);
00157       }
00158       else
00159       {
00160         if ( info.second == 0 ) {
00161           break;
00162         }
00163       }
00164       if (0 == child || childAlreadyExists || m_level == int(parentStack.size()))
00165       {
00166          if (0!=child)
00167          {
00168             child->SetLineColor(kRed);
00169          }
00170          //stop descending
00171          if ( ! walker.nextSibling())
00172          {
00173             while (walker.parent())
00174             {
00175                parentStack.pop_back();
00176                if (walker.nextSibling()) {
00177                   break;
00178                }
00179             }
00180          }
00181       }
00182       else
00183       {
00184          if (walker.firstChild())
00185          {
00186             parentStack.push_back(child);
00187          }
00188          else
00189          {          
00190             if ( ! walker.nextSibling())
00191             {
00192                while (walker.parent())
00193                {
00194                   parentStack.pop_back();
00195                   if (walker.nextSibling()) {
00196                      break;
00197                   }
00198                }
00199             }
00200          }
00201       }
00202    } while ( ! parentStack.empty());
00203 
00204    geo_mgr->CloseGeometry();
00205 
00206    geo_mgr->DefaultColors();
00207 
00208    nameToShape_.clear();
00209    nameToVolume_.clear();
00210    nameToMaterial_.clear();
00211    nameToMedium_.clear();
00212 
00213    return boost::shared_ptr<TGeoManager>(geo_mgr);
00214 }
00215 
00216 
00217 //==============================================================================
00218 // private member functions
00219 //==============================================================================
00220 
00221 TGeoShape* 
00222 TGeoMgrFromDdd::createShape(const std::string& iName,
00223                             const DDSolid&     iSolid)
00224 {
00225    LogDebug("TGeoMgrFromDdd::createShape") << "with name: " << iName << " and solid: " << iSolid;
00226 
00227    DDBase<DDName, DDI::Solid*>::def_type defined( iSolid.isDefined());
00228    if( !defined.first ) throw cms::Exception("TGeoMgrFromDdd::createShape * solid " + iName + " is not declared * " );
00229    if( !defined.second ) throw cms::Exception("TGeoMgrFromDdd::createShape * solid " + defined.first->name() + " is not defined *" );
00230    
00231    TGeoShape* rSolid= nameToShape_[iName];
00232    if (rSolid == 0)
00233    {
00234       const std::vector<double>& params = iSolid.parameters();
00235       switch(iSolid.shape())
00236       {
00237          case ddbox:
00238             rSolid = new TGeoBBox(
00239                                   iName.c_str(),
00240                                   params[0]/cm,
00241                                   params[1]/cm,
00242                                   params[2]/cm);
00243             break;
00244          case ddcons:
00245             rSolid = new TGeoConeSeg(
00246                                      iName.c_str(),
00247                                      params[0]/cm,
00248                                      params[1]/cm,
00249                                      params[2]/cm,
00250                                      params[3]/cm,
00251                                      params[4]/cm,
00252                                      params[5]/deg,
00253                                      params[6]/deg+params[5]/deg
00254                                      );
00255             break;
00256          case ddtubs:
00257             //Order in params is  zhalf,rIn,rOut,startPhi,deltaPhi
00258             rSolid= new TGeoTubeSeg(
00259                                     iName.c_str(),
00260                                     params[1]/cm,
00261                                     params[2]/cm,
00262                                     params[0]/cm,
00263                                     params[3]/deg,
00264                                     params[3]/deg + params[4]/deg);
00265             break;
00266          case ddtrap:
00267             rSolid =new TGeoTrap(
00268                                  iName.c_str(),
00269                                  params[0]/cm,  //dz
00270                                  params[1]/deg, //theta
00271                                  params[2]/deg, //phi
00272                                  params[3]/cm,  //dy1
00273                                  params[4]/cm,  //dx1
00274                                  params[5]/cm,  //dx2
00275                                  params[6]/deg, //alpha1
00276                                  params[7]/cm,  //dy2
00277                                  params[8]/cm,  //dx3
00278                                  params[9]/cm,  //dx4
00279                                  params[10]/deg);//alpha2
00280             break;
00281          case ddpolycone_rrz:    
00282             rSolid = new TGeoPcon(
00283                                   iName.c_str(),
00284                                   params[0]/deg,
00285                                   params[1]/deg,
00286                                   (params.size()-2)/3) ;
00287             {
00288                std::vector<double> temp(params.size()+1);
00289                temp.reserve(params.size()+1);
00290                temp[0]=params[0]/deg;
00291                temp[1]=params[1]/deg;
00292                temp[2]=(params.size()-2)/3;
00293                std::copy(params.begin()+2,params.end(),temp.begin()+3);
00294                for(std::vector<double>::iterator it=temp.begin()+3;
00295                    it != temp.end();
00296                    ++it) {
00297                   *it /=cm;
00298                }               
00299                rSolid->SetDimensions(&(*(temp.begin())));
00300             }
00301             break;
00302          case ddpolyhedra_rrz:
00303             rSolid = new TGeoPgon(
00304                                   iName.c_str(),
00305                                   params[1]/deg,
00306                                   params[2]/deg,
00307                                   static_cast<int>(params[0]),
00308                                   (params.size()-3)/3);
00309             {
00310                std::vector<double> temp(params.size()+1);
00311                temp[0]=params[1]/deg;
00312                temp[1]=params[2]/deg;
00313                temp[2]=params[0];
00314                temp[3]=(params.size()-3)/3;
00315                std::copy(params.begin()+3,params.end(),temp.begin()+4);
00316                for(std::vector<double>::iterator it=temp.begin()+4;
00317                    it != temp.end();
00318                    ++it) {
00319                   *it /=cm;
00320                }
00321                rSolid->SetDimensions(&(*(temp.begin())));
00322             }
00323             break;
00324          case ddpseudotrap:
00325          {
00326             //implementation taken from SimG4Core/Geometry/src/DDG4SolidConverter.cc
00327             static DDRotationMatrix s_rot( ROOT::Math::RotationX( 90.*deg ));
00328             DDPseudoTrap pt( iSolid );
00329 
00330             double r = pt.radius();
00331             bool atMinusZ = pt.atMinusZ();
00332             double x = 0;
00333             double h = 0;
00334             bool intersec = false; // union or intersection solid
00335 
00336             if( atMinusZ )
00337             {
00338                x = pt.x1(); // tubs radius
00339             }
00340             else
00341             {
00342                x = pt.x2(); // tubs radius
00343             }
00344             double halfOpeningAngle = asin( x / abs( r ))/deg;
00345             double displacement = 0;
00346             double startPhi = 0;
00347             /* calculate the displacement of the tubs w.r.t. to the trap,
00348                determine the opening angle of the tubs */
00349             double delta = sqrt( r * r - x * x );
00350 
00351             if( r < 0 && abs( r ) >= x )
00352             {
00353               intersec = true; // intersection solid
00354               h = pt.y1() < pt.y2() ? pt.y2() : pt.y1(); // tubs half height
00355               h += h/20.; // enlarge a bit - for subtraction solid
00356               if( atMinusZ )
00357               {
00358                 displacement = - pt.halfZ() - delta; 
00359                 startPhi = 90. - halfOpeningAngle;
00360               }
00361               else
00362               {
00363                 displacement =   pt.halfZ() + delta;
00364                 startPhi = -90.- halfOpeningAngle;
00365               }
00366             }
00367             else if( r > 0 && abs( r ) >= x )
00368             {
00369               if( atMinusZ )
00370               {
00371                 displacement = - pt.halfZ() + delta;
00372                 startPhi = 270.- halfOpeningAngle;
00373                 h = pt.y1();
00374               }
00375               else
00376               {
00377                 displacement =   pt.halfZ() - delta; 
00378                 startPhi = 90. - halfOpeningAngle;
00379                 h = pt.y2();
00380               }    
00381             }
00382             else
00383             {
00384               throw cms::Exception( "Check parameters of the PseudoTrap! name=" + pt.name().name());   
00385             }
00386 
00387             std::auto_ptr<TGeoShape> trap( new TGeoTrd2( pt.name().name().c_str(),
00388                                                          pt.x1()/cm,
00389                                                          pt.x2()/cm,
00390                                                          pt.y1()/cm,
00391                                                          pt.y2()/cm,
00392                                                          pt.halfZ()/cm ));
00393               
00394             std::auto_ptr<TGeoShape> tubs( new TGeoTubeSeg( pt.name().name().c_str(),
00395                                                             0.,
00396                                                             abs(r)/cm, // radius cannot be negative!!!
00397                                                             h/cm,
00398                                                             startPhi,
00399                                                             startPhi + halfOpeningAngle * 2. ));
00400             if( intersec )
00401             {
00402               TGeoSubtraction* sub = new TGeoSubtraction( trap.release(),
00403                                                           tubs.release(),
00404                                                           0,
00405                                                           createPlacement( s_rot,
00406                                                                            DDTranslation( 0.,
00407                                                                                           0.,
00408                                                                                           displacement )));
00409               rSolid = new TGeoCompositeShape( iName.c_str(),
00410                                                sub );
00411             }
00412             else
00413             {
00414               std::auto_ptr<TGeoShape> box( new TGeoBBox( 1.1*x/cm, 1.1*h/cm, sqrt(r*r-x*x)/cm ));
00415               
00416               TGeoSubtraction* sub = new TGeoSubtraction( tubs.release(),
00417                                                           box.release(),
00418                                                           0,
00419                                                           createPlacement( s_rot,
00420                                                                            DDTranslation( 0.,
00421                                                                                           0.,
00422                                                                                           0. )));
00423               
00424               std::auto_ptr<TGeoShape> tubicCap( new TGeoCompositeShape( iName.c_str(), sub ));
00425                                                 
00426               TGeoUnion* boolS = new TGeoUnion( trap.release(),
00427                                                 tubicCap.release(),
00428                                                 0,
00429                                                 createPlacement( s_rot,
00430                                                                  DDTranslation( 0.,
00431                                                                                 0.,
00432                                                                                 displacement )));
00433               
00434               rSolid = new TGeoCompositeShape( iName.c_str(),
00435                                                boolS );
00436             }
00437             
00438             break;
00439          }
00440          case ddtorus:
00441          {
00442             DDTorus solid( iSolid );
00443             rSolid = new TGeoTorus( iName.c_str(),
00444                                     solid.rTorus()/cm,
00445                                     solid.rMin()/cm,
00446                                     solid.rMax()/cm,
00447                                     solid.startPhi()/deg,
00448                                     solid.deltaPhi()/deg);
00449             break;
00450          }      
00451          case ddsubtraction:
00452          {
00453             DDBooleanSolid boolSolid(iSolid);
00454             if(!boolSolid) {
00455                throw cms::Exception("GeomConvert") <<"conversion to DDBooleanSolid failed";
00456             }
00457             
00458             std::auto_ptr<TGeoShape> left( createShape(boolSolid.solidA().name().fullname(),
00459                                                        boolSolid.solidA()) );
00460             std::auto_ptr<TGeoShape> right( createShape(boolSolid.solidB().name().fullname(),
00461                                                         boolSolid.solidB()));
00462             if( 0 != left.get() &&
00463                 0 != right.get() ) {
00464                TGeoSubtraction* sub = new TGeoSubtraction(left.release(),right.release(),
00465                                                           0,
00466                                                           createPlacement(
00467                                                                           *(boolSolid.rotation().matrix()),
00468                                                                           boolSolid.translation()));
00469                rSolid = new TGeoCompositeShape(iName.c_str(),
00470                                                sub);
00471             }
00472             break;
00473          }
00474          case ddtrunctubs:
00475          {
00476             DDTruncTubs tt( iSolid );
00477             if( !tt )
00478             {
00479                throw cms::Exception( "GeomConvert" ) << "conversion to DDTruncTubs failed";
00480             }
00481             double rIn( tt.rIn());
00482             double rOut( tt.rOut());
00483             double zHalf( tt.zHalf());
00484             double startPhi( tt.startPhi());
00485             double deltaPhi( tt.deltaPhi());
00486             double cutAtStart( tt.cutAtStart());
00487             double cutAtDelta( tt.cutAtDelta());
00488             bool cutInside( bool( tt.cutInside()));
00489             std::string name = tt.name().name();
00490             
00491             // check the parameters
00492             if( rIn <= 0 || rOut <=0 || cutAtStart <=0 || cutAtDelta <= 0 )
00493             {
00494               std::string s = "TruncTubs " + std::string( tt.name().fullname()) + ": 0 <= rIn,cutAtStart,rOut,cutAtDelta,rOut violated!";
00495               throw cms::Exception( s );
00496             }
00497             if( rIn >= rOut )
00498             {
00499               std::string s = "TruncTubs " + std::string( tt.name().fullname()) + ": rIn<rOut violated!";
00500               throw cms::Exception(s);
00501             }
00502             if( startPhi != 0. )
00503             {
00504               std::string s= "TruncTubs " + std::string( tt.name().fullname()) + ": startPhi != 0 not supported!";
00505               throw cms::Exception( s );
00506             }
00507             
00508             startPhi = 0.;
00509             double r( cutAtStart );
00510             double R( cutAtDelta );
00511             
00512             // Note: startPhi is always 0.0
00513             std::auto_ptr<TGeoShape> tubs( new TGeoTubeSeg( name.c_str(), rIn/cm, rOut/cm, zHalf/cm, startPhi, deltaPhi/deg ));
00514 
00515             double boxX( rOut ), boxY( rOut ); // exaggerate dimensions - does not matter, it's subtracted!
00516             
00517             // width of the box > width of the tubs
00518             double boxZ( 1.1 * zHalf );
00519             
00520             // angle of the box w.r.t. tubs
00521             double cath = r - R * cos( deltaPhi );
00522             double hypo = sqrt( r * r + R * R - 2. * r * R * cos( deltaPhi ));
00523             double cos_alpha = cath / hypo;
00524             double alpha = -acos( cos_alpha );
00525             
00526             // rotationmatrix of box w.r.t. tubs
00527             TGeoRotation rot;
00528             rot.RotateX( 90 );
00529             rot.RotateZ( alpha/deg );
00530             
00531             // center point of the box
00532             double xBox;
00533             if( !cutInside )
00534             {
00535               xBox = r + boxX / sin( fabs( alpha ));
00536             }
00537             else
00538             {
00539               xBox = - ( boxX / sin( fabs( alpha )) - r );
00540             }
00541             std::auto_ptr<TGeoShape> box( new TGeoBBox( name.c_str(), boxX/cm, boxZ/cm, boxY/cm ));
00542 
00543             TGeoTranslation trans( xBox/cm, 0., 0.);
00544 
00545             TGeoSubtraction* sub = new TGeoSubtraction( tubs.release(),
00546                                                         box.release(),
00547                                                         0, new TGeoCombiTrans( trans, rot ));
00548 
00549             rSolid = new TGeoCompositeShape( iName.c_str(),
00550                                              sub );
00551             break;   
00552          }
00553          case ddunion:
00554          {
00555             DDBooleanSolid boolSolid(iSolid);
00556             if(!boolSolid) {
00557                throw cms::Exception("GeomConvert") <<"conversion to DDBooleanSolid failed";
00558             }
00559             
00560             std::auto_ptr<TGeoShape> left( createShape(boolSolid.solidA().name().fullname(),
00561                                                        boolSolid.solidA()) );
00562             std::auto_ptr<TGeoShape> right( createShape(boolSolid.solidB().name().fullname(),
00563                                                         boolSolid.solidB()));
00564             //DEBUGGING
00565             //break;
00566             if( 0 != left.get() &&
00567                 0 != right.get() ) {
00568                TGeoUnion* boolS = new TGeoUnion(left.release(),right.release(),
00569                                                 0,
00570                                                 createPlacement(
00571                                                                 *(boolSolid.rotation().matrix()),
00572                                                                 boolSolid.translation()));
00573                rSolid = new TGeoCompositeShape(iName.c_str(),
00574                                                boolS);
00575             }
00576             break;
00577          }
00578          case ddintersection:
00579          {
00580             DDBooleanSolid boolSolid(iSolid);
00581             if(!boolSolid) {
00582                throw cms::Exception("GeomConvert") <<"conversion to DDBooleanSolid failed";
00583             }
00584             
00585             std::auto_ptr<TGeoShape> left( createShape(boolSolid.solidA().name().fullname(),
00586                                                        boolSolid.solidA()) );
00587             std::auto_ptr<TGeoShape> right( createShape(boolSolid.solidB().name().fullname(),
00588                                                         boolSolid.solidB()));
00589             if( 0 != left.get() &&
00590                 0 != right.get() ) {
00591                TGeoIntersection* boolS = new TGeoIntersection(left.release(),
00592                                                               right.release(),
00593                                                               0,
00594                                                               createPlacement(
00595                                                                               *(boolSolid.rotation().matrix()),
00596                                                                               boolSolid.translation()));
00597                rSolid = new TGeoCompositeShape(iName.c_str(),
00598                                                boolS);
00599             }
00600             break;
00601          }
00602          default:
00603             break;
00604       }
00605       nameToShape_[iName]=rSolid;
00606    }
00607    if (rSolid == 0)
00608    {
00609       std::cerr <<"COULD NOT MAKE "<<iName<<" of a shape "<<iSolid<<std::endl;
00610    }
00611 
00612    LogDebug("TGeoMgrFromDdd::createShape") << "solid " << iName << " has been created.";
00613 
00614    return rSolid;
00615 }
00616 
00617 TGeoVolume* 
00618 TGeoMgrFromDdd::createVolume(const std::string& iName,
00619                        const DDSolid& iSolid,
00620                        const DDMaterial& iMaterial)
00621 {
00622    TGeoVolume* v=nameToVolume_[iName];
00623    if (v == 0)
00624    {
00625       TGeoShape* solid     = createShape(iSolid.name().fullname(),
00626                                          iSolid);
00627       std::string mat_name = iMaterial.name().fullname();
00628       TGeoMedium *geo_med  = nameToMedium_[mat_name];
00629       if (geo_med == 0)
00630       {
00631          TGeoMaterial *geo_mat = createMaterial(iMaterial);
00632          geo_med = new TGeoMedium(mat_name.c_str(), 0, geo_mat);
00633          nameToMedium_[mat_name] = geo_med;
00634       }
00635       if (solid)
00636       {
00637          v = new TGeoVolume(iName.c_str(),
00638                             solid,
00639                             geo_med);
00640       }
00641       nameToVolume_[iName] = v;
00642    }
00643    return v;
00644 }
00645 
00646 TGeoMaterial*
00647 TGeoMgrFromDdd::createMaterial(const DDMaterial& iMaterial)
00648 {
00649    std::string   mat_name = iMaterial.name().fullname();
00650    TGeoMaterial *mat      = nameToMaterial_[mat_name];
00651 
00652    if (mat == 0)
00653    {
00654       if (iMaterial.noOfConstituents() > 0)
00655       {
00656          TGeoMixture *mix = new TGeoMixture(mat_name.c_str(),
00657                                             iMaterial.noOfConstituents(),
00658                                             iMaterial.density()*cm3/g);
00659          for (int i = 0; i < iMaterial.noOfConstituents(); ++i)
00660          {
00661             mix->AddElement(createMaterial(iMaterial.constituent(i).first),
00662                             iMaterial.constituent(i).second);
00663          }
00664          mat = mix;
00665       }
00666       else
00667       {
00668          mat = new TGeoMaterial(mat_name.c_str(),
00669                                 iMaterial.a()*mole/g, iMaterial.z(),
00670                                 iMaterial.density()*cm3/g);
00671       }
00672       nameToMaterial_[mat_name] = mat;
00673    }
00674 
00675    return mat;
00676 }
00677 
00678 //
00679 // const member functions
00680 //
00681 
00682 //
00683 // static member functions
00684 //