00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
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
00053
00054 setWhatProduced( this );
00055 }
00056
00057 TGeoMgrFromDdd::~TGeoMgrFromDdd( void )
00058 {}
00059
00060
00061
00062
00063
00064 namespace
00065 {
00066 TGeoCombiTrans* createPlacement(const DDRotationMatrix& iRot,
00067 const DDTranslation& iTrans)
00068 {
00069
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
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
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
00112
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
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
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
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
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,
00270 params[1]/deg,
00271 params[2]/deg,
00272 params[3]/cm,
00273 params[4]/cm,
00274 params[5]/cm,
00275 params[6]/deg,
00276 params[7]/cm,
00277 params[8]/cm,
00278 params[9]/cm,
00279 params[10]/deg);
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
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;
00335
00336 if( atMinusZ )
00337 {
00338 x = pt.x1();
00339 }
00340 else
00341 {
00342 x = pt.x2();
00343 }
00344 double halfOpeningAngle = asin( x / abs( r ))/deg;
00345 double displacement = 0;
00346 double startPhi = 0;
00347
00348
00349 double delta = sqrt( r * r - x * x );
00350
00351 if( r < 0 && abs( r ) >= x )
00352 {
00353 intersec = true;
00354 h = pt.y1() < pt.y2() ? pt.y2() : pt.y1();
00355 h += h/20.;
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,
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
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
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 );
00516
00517
00518 double boxZ( 1.1 * zHalf );
00519
00520
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
00527 TGeoRotation rot;
00528 rot.RotateX( 90 );
00529 rot.RotateZ( alpha/deg );
00530
00531
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
00565
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
00680
00681
00682
00683
00684
1.7.3