CMS 3D CMS Logo

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