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
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
00052
00053 setWhatProduced( this );
00054 }
00055
00056 TGeoMgrFromDdd::~TGeoMgrFromDdd( void )
00057 {}
00058
00059
00060
00061
00062
00063 namespace
00064 {
00065 TGeoCombiTrans* createPlacement(const DDRotationMatrix& iRot,
00066 const DDTranslation& iTrans)
00067 {
00068
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
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
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
00111
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
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
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
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
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
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,
00264 params[1]/deg,
00265 params[2]/deg,
00266 params[3]/cm,
00267 params[4]/cm,
00268 params[5]/cm,
00269 params[6]/deg,
00270 params[7]/cm,
00271 params[8]/cm,
00272 params[9]/cm,
00273 params[10]/deg);
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
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;
00329
00330 if( atMinusZ )
00331 {
00332 x = pt.x1();
00333 }
00334 else
00335 {
00336 x = pt.x2();
00337 }
00338 double halfOpeningAngle = asin( x / abs( r ))/deg;
00339 double displacement = 0;
00340 double startPhi = 0;
00341
00342
00343 double delta = sqrt( r * r - x * x );
00344
00345 if( r < 0 && abs( r ) >= x )
00346 {
00347 intersec = true;
00348 h = pt.y1() < pt.y2() ? pt.y2() : pt.y1();
00349 h += h/20.;
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,
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
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
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 );
00510
00511
00512 double boxZ( 1.1 * zHalf );
00513
00514
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
00521 TGeoRotation rot;
00522 rot.RotateX( 90 );
00523 rot.RotateZ( alpha/deg );
00524
00525
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
00559
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
00671
00672
00673
00674
00675