Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:11:48 

0001 #include "Fireworks/Geometry/interface/FWTGeoRecoGeometryESProducer.h"
0002 #include "Fireworks/Geometry/interface/FWTGeoRecoGeometry.h"
0003 #include "Fireworks/Geometry/interface/FWTGeoRecoGeometryRecord.h"
0004 
0005 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0006 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0007 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0008 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0009 #include "DataFormats/MuonDetId/interface/GEMDetId.h"
0010 #include "DataFormats/MuonDetId/interface/ME0DetId.h"
0011 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0012 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0013 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0014 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0015 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0016 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 
0019 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0020 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0021 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0022 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0023 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0024 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0025 #include "Geometry/CSCGeometry/interface/CSCChamber.h"
0026 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0027 #include "Geometry/DTGeometry/interface/DTChamber.h"
0028 #include "Geometry/DTGeometry/interface/DTLayer.h"
0029 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0030 #include "Geometry/GEMGeometry/interface/GEMEtaPartition.h"
0031 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0032 #include "Geometry/GEMGeometry/interface/ME0EtaPartition.h"
0033 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0034 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0035 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0036 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0037 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0038 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0039 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0040 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0041 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0042 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0043 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
0044 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
0045 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0046 
0047 #include "TGeoManager.h"
0048 #include "TGeoArb8.h"
0049 #include "TGeoMatrix.h"
0050 
0051 #include "TFile.h"
0052 #include "TTree.h"
0053 #include "TError.h"
0054 #include "TMath.h"
0055 
0056 #include "TEveVector.h"
0057 #include "TEveTrans.h"
0058 
0059 //std::map< const CaloCellGeometry::CCGFloat*, TGeoVolume*> caloShapeMap;
0060 namespace {
0061   typedef std::map<const float*, TGeoVolume*> CaloVolMap;
0062 }
0063 FWTGeoRecoGeometryESProducer::FWTGeoRecoGeometryESProducer(const edm::ParameterSet& pset) : m_dummyMedium(nullptr) {
0064   m_tracker = pset.getUntrackedParameter<bool>("Tracker", true);
0065   m_muon = pset.getUntrackedParameter<bool>("Muon", true);
0066   m_gem = pset.getUntrackedParameter<bool>("GEM", false);
0067   m_calo = pset.getUntrackedParameter<bool>("Calo", true);
0068 
0069   if (m_muon)
0070     m_gem = true;
0071 
0072   auto cc = setWhatProduced(this);
0073   if (m_tracker || m_muon || m_gem) {
0074     m_trackingGeomToken = cc.consumes();
0075   }
0076   if (m_tracker) {
0077     m_trackerTopologyToken = cc.consumes();
0078   }
0079   if (m_calo) {
0080     m_caloGeomToken = cc.consumes();
0081   }
0082 }
0083 
0084 FWTGeoRecoGeometryESProducer::~FWTGeoRecoGeometryESProducer(void) {}
0085 
0086 namespace {
0087 
0088   void AddLeafNode(TGeoVolume* mother, TGeoVolume* daughter, const char* name, TGeoMatrix* mtx) {
0089     int n = mother->GetNdaughters();
0090     mother->AddNode(daughter, 1, mtx);
0091     mother->GetNode(n)->SetName(name);
0092   }
0093 
0094   /** Create TGeo transformation of GeomDet */
0095   TGeoCombiTrans* createPlacement(const GeomDet* det) {
0096     // Position of the DetUnit's center
0097     float posx = det->surface().position().x();
0098     float posy = det->surface().position().y();
0099     float posz = det->surface().position().z();
0100 
0101     TGeoTranslation trans(posx, posy, posz);
0102 
0103     // Add the coeff of the rotation matrix
0104     // with a projection on the basis vectors
0105     TkRotation<float> detRot = det->surface().rotation();
0106 
0107     TGeoRotation rotation;
0108     const Double_t matrix[9] = {detRot.xx(),
0109                                 detRot.yx(),
0110                                 detRot.zx(),
0111                                 detRot.xy(),
0112                                 detRot.yy(),
0113                                 detRot.zy(),
0114                                 detRot.xz(),
0115                                 detRot.yz(),
0116                                 detRot.zz()};
0117     rotation.SetMatrix(matrix);
0118 
0119     return new TGeoCombiTrans(trans, rotation);
0120   }
0121 
0122 }  // namespace
0123 //______________________________________________________________________________
0124 
0125 TGeoVolume* FWTGeoRecoGeometryESProducer::GetDaughter(TGeoVolume* mother, const char* prefix, ERecoDet cidx, int id) {
0126   TGeoVolume* res = nullptr;
0127   if (mother->GetNdaughters()) {
0128     TGeoNode* n = mother->FindNode(Form("%s_%d_1", prefix, id));
0129     if (n)
0130       res = n->GetVolume();
0131   }
0132 
0133   if (!res) {
0134     res = new TGeoVolumeAssembly(Form("%s_%d", prefix, id));
0135     res->SetMedium(GetMedium(cidx));
0136     mother->AddNode(res, 1);
0137   }
0138 
0139   return res;
0140 }
0141 
0142 TGeoVolume* FWTGeoRecoGeometryESProducer::GetDaughter(TGeoVolume* mother, const char* prefix, ERecoDet cidx) {
0143   TGeoVolume* res = nullptr;
0144   if (mother->GetNdaughters()) {
0145     TGeoNode* n = mother->FindNode(Form("%s_1", prefix));
0146     if (n)
0147       res = n->GetVolume();
0148   }
0149 
0150   if (!res) {
0151     //      printf("GetDau... new holder %s for mother %s \n", mother->GetName(), prefix);
0152     res = new TGeoVolumeAssembly(prefix);
0153     res->SetMedium(GetMedium(cidx));
0154     mother->AddNode(res, 1);
0155   }
0156 
0157   return res;
0158 }
0159 
0160 TGeoVolume* FWTGeoRecoGeometryESProducer::GetTopHolder(const char* prefix, ERecoDet cidx) {
0161   //   printf("GetTopHolder res = %s \n", prefix);
0162   TGeoVolume* res = GetDaughter(gGeoManager->GetTopVolume(), prefix, cidx);
0163   return res;
0164 }
0165 
0166 namespace {
0167 
0168   enum GMCol {
0169     Green = 4,
0170     Blue0 = 13,
0171     Blue1 = 24,
0172     Blue2 = 6,
0173     Yellow0 = 3,
0174     Yellow1 = 16,
0175     Pink = 10,
0176     Red = 29,
0177     Orange0 = 79,
0178     Orange1 = 14,
0179     Magenta = 8,
0180     Gray = 12
0181   };
0182 
0183 }
0184 
0185 TGeoMedium* FWTGeoRecoGeometryESProducer::GetMedium(ERecoDet det) {
0186   std::map<ERecoDet, TGeoMedium*>::iterator it = m_recoMedium.find(det);
0187   if (it != m_recoMedium.end())
0188     return it->second;
0189 
0190   std::string name;
0191   int color;
0192 
0193   switch (det) {
0194     // TRACKER
0195     case kSiPixel:
0196       name = "SiPixel";
0197       color = GMCol::Green;
0198       break;
0199 
0200     case kSiStrip:
0201       name = "SiStrip";
0202       color = GMCol::Gray;
0203       break;
0204       // MUON
0205     case kMuonDT:
0206       name = "MuonDT";
0207       color = GMCol::Blue2;
0208       break;
0209 
0210     case kMuonRPC:
0211       name = "MuonRPC";
0212       color = GMCol::Red;
0213       break;
0214 
0215     case kMuonGEM:
0216       name = "MuonGEM";
0217       color = GMCol::Yellow1;
0218       break;
0219 
0220     case kMuonCSC:
0221       name = "MuonCSC";
0222       color = GMCol::Gray;
0223       break;
0224 
0225     case kMuonME0:
0226       name = "MuonME0";
0227       color = GMCol::Yellow0;
0228       break;
0229 
0230       // CALO
0231     case kECal:
0232       name = "ECal";
0233       color = GMCol::Blue2;
0234       break;
0235     case kHCal:
0236       name = "HCal";
0237       color = GMCol::Orange1;
0238       break;
0239     case kCaloTower:
0240       name = "CaloTower";
0241       color = GMCol::Green;
0242       break;
0243     case kHGCE:
0244       name = "HGCEE";
0245       color = GMCol::Blue2;
0246       break;
0247     case kHGCH:
0248       name = "HGCEH";
0249       color = GMCol::Blue1;
0250       break;
0251     default:
0252       printf("invalid medium id \n");
0253       return m_dummyMedium;
0254   }
0255 
0256   TGeoMaterial* mat = new TGeoMaterial(name.c_str(), 0, 0, 0);
0257   mat->SetZ(color);
0258   m_recoMedium[det] = new TGeoMedium(name.c_str(), 0, mat);
0259   mat->SetFillStyle(3000);  // tansparency 3000-3100
0260   mat->SetDensity(1);       // disable override of transparency in TGeoManager::DefaultColors()
0261 
0262   return m_recoMedium[det];
0263 }
0264 
0265 //______________________________________________________________________________
0266 
0267 std::unique_ptr<FWTGeoRecoGeometry> FWTGeoRecoGeometryESProducer::produce(const FWTGeoRecoGeometryRecord& record) {
0268   using namespace edm;
0269 
0270   auto fwTGeoRecoGeometry = std::make_unique<FWTGeoRecoGeometry>();
0271 
0272   if (m_calo) {
0273     m_caloGeom = &record.get(m_caloGeomToken);
0274   }
0275 
0276   TGeoManager* geom = new TGeoManager("cmsGeo", "CMS Detector");
0277   if (nullptr == gGeoIdentity) {
0278     gGeoIdentity = new TGeoIdentity("Identity");
0279   }
0280 
0281   fwTGeoRecoGeometry->manager(geom);
0282 
0283   // Default material is Vacuum
0284   TGeoMaterial* vacuum = new TGeoMaterial("Vacuum", 0, 0, 0);
0285   m_dummyMedium = new TGeoMedium("reco", 0, vacuum);
0286 
0287   TGeoVolume* top = geom->MakeBox("CMS", m_dummyMedium, 270., 270., 120.);
0288 
0289   if (nullptr == top) {
0290     return std::unique_ptr<FWTGeoRecoGeometry>();
0291   }
0292   geom->SetTopVolume(top);
0293   // ROOT chokes unless colors are assigned
0294   top->SetVisibility(kFALSE);
0295   top->SetLineColor(kBlue);
0296 
0297   if (m_tracker || m_muon || m_gem) {
0298     m_trackingGeom = &record.get(m_trackingGeomToken);
0299   }
0300 
0301   if (m_tracker) {
0302     DetId detId(DetId::Tracker, 0);
0303     m_trackerGeom = static_cast<const TrackerGeometry*>(m_trackingGeom->slaveGeometry(detId));
0304 
0305     m_trackerTopology = &record.get(m_trackerTopologyToken);
0306 
0307     addPixelBarrelGeometry();
0308     addPixelForwardGeometry();
0309 
0310     addTIBGeometry();
0311     addTIDGeometry();
0312     addTOBGeometry();
0313     addTECGeometry();
0314   }
0315 
0316   if (m_muon) {
0317     addDTGeometry();
0318     addCSCGeometry();
0319     addRPCGeometry();
0320     addME0Geometry();
0321   }
0322 
0323   if (m_gem) {
0324     addGEMGeometry();
0325   }
0326 
0327   if (m_calo) {
0328     addEcalCaloGeometry();
0329     addHcalCaloGeometryBarrel();
0330     addHcalCaloGeometryEndcap();
0331     addHcalCaloGeometryOuter();
0332     addHcalCaloGeometryForward();
0333     addCaloTowerGeometry();
0334   }
0335 
0336   geom->CloseGeometry();
0337 
0338   geom->DefaultColors();
0339   // printf("==== geo manager NNodes = %d \n", geom->GetNNodes());
0340   geom->CloseGeometry();
0341 
0342   return fwTGeoRecoGeometry;
0343 }
0344 
0345 /** Create TGeo shape for GeomDet */
0346 TGeoShape* FWTGeoRecoGeometryESProducer::createShape(const GeomDet* det) {
0347   TGeoShape* shape = nullptr;
0348 
0349   // Trapezoidal
0350   const Bounds* b = &((det->surface()).bounds());
0351   const TrapezoidalPlaneBounds* b2 = dynamic_cast<const TrapezoidalPlaneBounds*>(b);
0352   if (b2) {
0353     std::array<const float, 4> const& par = b2->parameters();
0354 
0355     // These parameters are half-lengths, as in CMSIM/GEANT3
0356     float hBottomEdge = par[0];
0357     float hTopEdge = par[1];
0358     float thickness = par[2];
0359     float apothem = par[3];
0360 
0361     std::stringstream s;
0362     s << "T_" << hBottomEdge << "_" << hTopEdge << "_" << thickness << "_" << apothem;
0363     std::string name = s.str();
0364 
0365     // Do not create identical shape,
0366     // if one already exists
0367     shape = m_nameToShape[name];
0368     if (nullptr == shape) {
0369       shape = new TGeoTrap(name.c_str(),
0370                            thickness,    //dz
0371                            0,            //theta
0372                            0,            //phi
0373                            apothem,      //dy1
0374                            hBottomEdge,  //dx1
0375                            hTopEdge,     //dx2
0376                            0,            //alpha1
0377                            apothem,      //dy2
0378                            hBottomEdge,  //dx3
0379                            hTopEdge,     //dx4
0380                            0);           //alpha2
0381 
0382       m_nameToShape[name] = shape;
0383     }
0384   }
0385   if (dynamic_cast<const RectangularPlaneBounds*>(b) != nullptr) {
0386     // Rectangular
0387     float length = det->surface().bounds().length();
0388     float width = det->surface().bounds().width();
0389     float thickness = det->surface().bounds().thickness();
0390 
0391     std::stringstream s;
0392     s << "R_" << width << "_" << length << "_" << thickness;
0393     std::string name = s.str();
0394 
0395     // Do not create identical shape,
0396     // if one already exists
0397     shape = m_nameToShape[name];
0398     if (nullptr == shape) {
0399       shape = new TGeoBBox(name.c_str(), width / 2., length / 2., thickness / 2.);  // dx, dy, dz
0400 
0401       m_nameToShape[name] = shape;
0402     }
0403   }
0404 
0405   return shape;
0406 }
0407 
0408 /** Create TGeo volume for GeomDet */
0409 TGeoVolume* FWTGeoRecoGeometryESProducer::createVolume(const std::string& name, const GeomDet* det, ERecoDet mid) {
0410   TGeoShape* solid = createShape(det);
0411 
0412   std::map<TGeoShape*, TGeoVolume*>::iterator vIt = m_shapeToVolume.find(solid);
0413   if (vIt != m_shapeToVolume.end())
0414     return vIt->second;
0415 
0416   TGeoVolume* volume = new TGeoVolume(name.c_str(), solid, GetMedium(mid));
0417 
0418   m_shapeToVolume[solid] = volume;
0419 
0420   return volume;
0421 }
0422 
0423 ///////////////////////////////////////////////////////////////////////////////////
0424 //////////////////////////////////////////////////////////////////////////////////
0425 ///////////////////////////////////////////////////////////////////////////////////
0426 //////////////////////////////////////////////////////////////////////////////////
0427 ///////////////////////////////////////////////////////////////////////////////////
0428 //////////////////////////////////////////////////////////////////////////////////
0429 ///////////////////////////////////////////////////////////////////////////////////
0430 //////////////////////////////////////////////////////////////////////////////////
0431 
0432 void FWTGeoRecoGeometryESProducer::addPixelBarrelGeometry() {
0433   TGeoVolume* tv = GetTopHolder("SiPixel", kSiPixel);
0434   TGeoVolume* assembly = GetDaughter(tv, "PXB", kSiPixel);
0435 
0436   for (auto it : m_trackerGeom->detsPXB()) {
0437     DetId detid = it->geographicalId();
0438     unsigned int layer = m_trackerTopology->pxbLayer(detid);
0439     unsigned int module = m_trackerTopology->pxbModule(detid);
0440     unsigned int ladder = m_trackerTopology->pxbLadder(detid);
0441 
0442     std::string name = Form("PXB Ly:%d, Md:%d Ld:%d ", layer, module, ladder);
0443     TGeoVolume* child = createVolume(name, it, kSiPixel);
0444 
0445     TGeoVolume* holder = GetDaughter(assembly, "Layer", kSiPixel, layer);
0446     holder = GetDaughter(holder, "Module", kSiPixel, module);
0447 
0448     AddLeafNode(holder, child, name.c_str(), createPlacement(it));
0449   }
0450 }
0451 
0452 void FWTGeoRecoGeometryESProducer::addPixelForwardGeometry() {
0453   TGeoVolume* tv = GetTopHolder("SiPixel", kSiPixel);
0454   TGeoVolume* assembly = GetDaughter(tv, "PXF", kSiPixel);
0455 
0456   for (auto it : m_trackerGeom->detsPXF()) {
0457     DetId detid = it->geographicalId();
0458     unsigned int disk = m_trackerTopology->pxfDisk(detid);
0459     unsigned int blade = m_trackerTopology->pxfBlade(detid);
0460     unsigned int panel = m_trackerTopology->pxfPanel(detid);
0461     unsigned int side = m_trackerTopology->side(detid);
0462 
0463     std::string name = Form("PXF D:%d, B:%d, P:%d, S:%d", disk, blade, panel, side);
0464     TGeoVolume* child = createVolume(name, it, kSiPixel);
0465 
0466     TGeoVolume* holder = GetDaughter(assembly, "Side", kSiPixel, side);
0467     holder = GetDaughter(holder, "Disk", kSiPixel, disk);
0468     holder = GetDaughter(holder, "Blade", kSiPixel, blade);
0469     holder = GetDaughter(holder, "Panel", kSiPixel, panel);
0470 
0471     AddLeafNode(holder, child, name.c_str(), createPlacement(it));
0472   }
0473 }
0474 
0475 void FWTGeoRecoGeometryESProducer::addTIBGeometry() {
0476   TGeoVolume* tv = GetTopHolder("SiStrip", kSiStrip);
0477   TGeoVolume* assembly = GetDaughter(tv, "TIB", kSiStrip);
0478 
0479   for (auto it : m_trackerGeom->detsTIB()) {
0480     DetId detid = it->geographicalId();
0481     unsigned int module = m_trackerTopology->tibModule(detid);
0482     unsigned int order = m_trackerTopology->tibOrder(detid);
0483     unsigned int side = m_trackerTopology->tibSide(detid);
0484 
0485     std::string name = m_trackerTopology->print(detid);
0486 
0487     TGeoVolume* child = createVolume(name, it, kSiStrip);
0488     TGeoVolume* holder = GetDaughter(assembly, "Module", kSiStrip, module);
0489     holder = GetDaughter(holder, "Order", kSiStrip, order);
0490     holder = GetDaughter(holder, "Side", kSiStrip, side);
0491     AddLeafNode(holder, child, name.c_str(), createPlacement(it));
0492   }
0493 }
0494 
0495 void FWTGeoRecoGeometryESProducer::addTIDGeometry() {
0496   TGeoVolume* tv = GetTopHolder("SiStrip", kSiStrip);
0497   TGeoVolume* assembly = GetDaughter(tv, "TID", kSiStrip);
0498 
0499   for (auto it : m_trackerGeom->detsTID()) {
0500     DetId detid = it->geographicalId();
0501     unsigned int side = m_trackerTopology->tidSide(detid);
0502     unsigned int wheel = m_trackerTopology->tidWheel(detid);
0503     unsigned int ring = m_trackerTopology->tidRing(detid);
0504 
0505     std::string name = m_trackerTopology->print(detid);
0506 
0507     TGeoVolume* child = createVolume(name, it, kSiStrip);
0508     TGeoVolume* holder = GetDaughter(assembly, "Side", kSiStrip, side);
0509     holder = GetDaughter(holder, "Wheel", kSiStrip, wheel);
0510     holder = GetDaughter(holder, "Ring", kSiStrip, ring);
0511     AddLeafNode(holder, child, name.c_str(), createPlacement(it));
0512   }
0513 }
0514 
0515 void FWTGeoRecoGeometryESProducer::addTOBGeometry() {
0516   TGeoVolume* tv = GetTopHolder("SiStrip", kSiStrip);
0517   TGeoVolume* assembly = GetDaughter(tv, "TOB", kSiStrip);
0518 
0519   for (auto it : m_trackerGeom->detsTOB()) {
0520     DetId detid = it->geographicalId();
0521     unsigned int rod = m_trackerTopology->tobRod(detid);
0522     unsigned int side = m_trackerTopology->tobSide(detid);
0523     unsigned int module = m_trackerTopology->tobModule(detid);
0524 
0525     std::string name = m_trackerTopology->print(detid);
0526 
0527     TGeoVolume* child = createVolume(name, it, kSiStrip);
0528     TGeoVolume* holder = GetDaughter(assembly, "Rod", kSiStrip, rod);
0529     holder = GetDaughter(holder, "Side", kSiStrip, side);
0530     holder = GetDaughter(holder, "Module", kSiStrip, module);
0531     AddLeafNode(holder, child, name.c_str(), createPlacement(it));
0532   }
0533 }
0534 
0535 void FWTGeoRecoGeometryESProducer::addTECGeometry() {
0536   TGeoVolume* tv = GetTopHolder("SiStrip", kSiStrip);
0537   TGeoVolume* assembly = GetDaughter(tv, "TEC", kSiStrip);
0538 
0539   for (auto it : m_trackerGeom->detsTEC()) {
0540     DetId detid = it->geographicalId();
0541     unsigned int order = m_trackerTopology->tecOrder(detid);
0542     unsigned int ring = m_trackerTopology->tecRing(detid);
0543     unsigned int module = m_trackerTopology->tecModule(detid);
0544 
0545     std::string name = m_trackerTopology->print(detid);
0546 
0547     TGeoVolume* child = createVolume(name, it, kSiStrip);
0548 
0549     TGeoVolume* holder = GetDaughter(assembly, "Order", kSiStrip, order);
0550     holder = GetDaughter(holder, "Ring", kSiStrip, ring);
0551     holder = GetDaughter(holder, "Module", kSiStrip, module);
0552     AddLeafNode(holder, child, name.c_str(), createPlacement(it));
0553   }
0554 }
0555 
0556 //==============================================================================
0557 //==============================================================================
0558 //=================================== MUON =====================================
0559 //==============================================================================
0560 
0561 void FWTGeoRecoGeometryESProducer::addDTGeometry() {
0562   TGeoVolume* tv = GetTopHolder("Muon", kMuonRPC);
0563   TGeoVolume* assemblyTop = GetDaughter(tv, "DT", kMuonDT);
0564 
0565   //
0566   // DT chambers geometry
0567   //
0568   {
0569     TGeoVolume* assembly = GetDaughter(assemblyTop, "DTChamber", kMuonDT);
0570     auto const& dtChamberGeom = m_trackingGeom->slaveGeometry(DTChamberId())->dets();
0571     for (auto it = dtChamberGeom.begin(), end = dtChamberGeom.end(); it != end; ++it) {
0572       if (auto chamber = dynamic_cast<const DTChamber*>(*it)) {
0573         DTChamberId detid = chamber->geographicalId();
0574         std::stringstream s;
0575         s << detid;
0576         std::string name = s.str();
0577 
0578         TGeoVolume* child = createVolume(name, chamber, kMuonDT);
0579         TGeoVolume* holder = GetDaughter(assembly, "Wheel", kMuonDT, detid.wheel());
0580         holder = GetDaughter(holder, "Station", kMuonDT, detid.station());
0581         holder = GetDaughter(holder, "Sector", kMuonDT, detid.sector());
0582 
0583         AddLeafNode(holder, child, name.c_str(), createPlacement(chamber));
0584       }
0585     }
0586   }
0587 
0588   // Fill in DT super layer parameters
0589   {
0590     TGeoVolume* assembly = GetDaughter(assemblyTop, "DTSuperLayer", kMuonDT);
0591     auto const& dtSuperLayerGeom = m_trackingGeom->slaveGeometry(DTSuperLayerId())->dets();
0592     for (auto it = dtSuperLayerGeom.begin(), end = dtSuperLayerGeom.end(); it != end; ++it) {
0593       if (auto* superlayer = dynamic_cast<const DTSuperLayer*>(*it)) {
0594         DTSuperLayerId detid(DetId(superlayer->geographicalId()));
0595         std::stringstream s;
0596         s << detid;
0597         std::string name = s.str();
0598 
0599         TGeoVolume* child = createVolume(name, superlayer, kMuonDT);
0600 
0601         TGeoVolume* holder = GetDaughter(assembly, "Wheel", kMuonDT, detid.wheel());
0602         holder = GetDaughter(holder, "Station", kMuonDT, detid.station());
0603         holder = GetDaughter(holder, "Sector", kMuonDT, detid.sector());
0604         holder = GetDaughter(holder, "SuperLayer", kMuonDT, detid.superlayer());
0605         AddLeafNode(holder, child, name.c_str(), createPlacement(superlayer));
0606       }
0607     }
0608   }
0609   // Fill in DT layer parameters
0610   {
0611     TGeoVolume* assembly = GetDaughter(assemblyTop, "DTLayer", kMuonDT);
0612     auto const& dtLayerGeom = m_trackingGeom->slaveGeometry(DTLayerId())->dets();
0613     for (auto it = dtLayerGeom.begin(), end = dtLayerGeom.end(); it != end; ++it) {
0614       if (auto layer = dynamic_cast<const DTLayer*>(*it)) {
0615         DTLayerId detid(DetId(layer->geographicalId()));
0616 
0617         std::stringstream s;
0618         s << detid;
0619         std::string name = s.str();
0620 
0621         TGeoVolume* child = createVolume(name, layer, kMuonDT);
0622 
0623         TGeoVolume* holder = GetDaughter(assembly, "Wheel", kMuonDT, detid.wheel());
0624         holder = GetDaughter(holder, "Station", kMuonDT, detid.station());
0625         holder = GetDaughter(holder, "Sector", kMuonDT, detid.sector());
0626         holder = GetDaughter(holder, "SuperLayer", kMuonDT, detid.superlayer());
0627         holder = GetDaughter(holder, "Layer", kMuonDT, detid.layer());
0628         AddLeafNode(holder, child, name.c_str(), createPlacement(layer));
0629       }
0630     }
0631   }
0632 }
0633 //______________________________________________________________________________
0634 
0635 void FWTGeoRecoGeometryESProducer::addCSCGeometry() {
0636   if (!m_trackingGeom->slaveGeometry(CSCDetId()))
0637     throw cms::Exception("FatalError") << "Cannnot find CSCGeometry\n";
0638 
0639   TGeoVolume* tv = GetTopHolder("Muon", kMuonRPC);
0640   TGeoVolume* assembly = GetDaughter(tv, "CSC", kMuonCSC);
0641 
0642   auto const& cscGeom = m_trackingGeom->slaveGeometry(CSCDetId())->dets();
0643   for (auto it = cscGeom.begin(), itEnd = cscGeom.end(); it != itEnd; ++it) {
0644     unsigned int rawid = (*it)->geographicalId();
0645     CSCDetId detId(rawid);
0646     std::stringstream s;
0647     s << "CSC" << detId;
0648     std::string name = s.str();
0649 
0650     TGeoVolume* child = nullptr;
0651 
0652     if (auto chamber = dynamic_cast<const CSCChamber*>(*it))
0653       child = createVolume(name, chamber, kMuonCSC);
0654     else if (auto* layer = dynamic_cast<const CSCLayer*>(*it))
0655       child = createVolume(name, layer, kMuonCSC);
0656 
0657     if (child) {
0658       TGeoVolume* holder = GetDaughter(assembly, "Endcap", kMuonCSC, detId.endcap());
0659       holder = GetDaughter(holder, "Station", kMuonCSC, detId.station());
0660       holder = GetDaughter(holder, "Ring", kMuonCSC, detId.ring());
0661       holder = GetDaughter(holder, "Chamber", kMuonCSC, detId.chamber());
0662 
0663       //   holder->AddNode(child, 1,  createPlacement( *it ));
0664       AddLeafNode(holder, child, name.c_str(), createPlacement(*it));
0665     }
0666   }
0667 }
0668 
0669 //______________________________________________________________________________
0670 
0671 void FWTGeoRecoGeometryESProducer::addGEMGeometry() {
0672   try {
0673     DetId detId(DetId::Muon, MuonSubdetId::GEM);
0674     const GEMGeometry* gemGeom = static_cast<const GEMGeometry*>(m_trackingGeom->slaveGeometry(detId));
0675 
0676     TGeoVolume* tv = GetTopHolder("Muon", kMuonRPC);
0677     TGeoVolume* assemblyTop = GetDaughter(tv, "GEM", kMuonGEM);
0678 
0679     {
0680       TGeoVolume* assembly = GetDaughter(assemblyTop, "GEMSuperChambers", kMuonGEM);
0681       for (auto it = gemGeom->superChambers().begin(), end = gemGeom->superChambers().end(); it != end; ++it) {
0682         const GEMSuperChamber* sc = (*it);
0683         if (sc) {
0684           GEMDetId detid = sc->geographicalId();
0685           std::stringstream s;
0686           s << detid;
0687           std::string name = s.str();
0688 
0689           TGeoVolume* child = createVolume(name, sc, kMuonGEM);
0690 
0691           TGeoVolume* holder = GetDaughter(assembly, "SuperChamber Region", kMuonGEM, detid.region());
0692           holder = GetDaughter(holder, "Ring", kMuonGEM, detid.ring());
0693           holder = GetDaughter(holder, "Station", kMuonGEM, detid.station());
0694           holder = GetDaughter(holder, "Chamber", kMuonGEM, detid.chamber());
0695 
0696           AddLeafNode(holder, child, name.c_str(), createPlacement(*it));
0697         }
0698       }
0699     }
0700 
0701     {
0702       TGeoVolume* assembly = GetDaughter(assemblyTop, "GEMetaPartitions", kMuonGEM);
0703       for (auto it = gemGeom->etaPartitions().begin(), end = gemGeom->etaPartitions().end(); it != end; ++it) {
0704         const GEMEtaPartition* roll = (*it);
0705         if (roll) {
0706           GEMDetId detid = roll->geographicalId();
0707           std::stringstream s;
0708           s << detid;
0709           std::string name = s.str();
0710 
0711           TGeoVolume* child = createVolume(name, roll, kMuonGEM);
0712 
0713           TGeoVolume* holder = GetDaughter(assembly, "ROLL Region", kMuonGEM, detid.region());
0714           holder = GetDaughter(holder, "Ring", kMuonGEM, detid.ring());
0715           holder = GetDaughter(holder, "Station", kMuonGEM, detid.station());
0716           holder = GetDaughter(holder, "Layer", kMuonGEM, detid.layer());
0717           holder = GetDaughter(holder, "Chamber", kMuonGEM, detid.chamber());
0718 
0719           AddLeafNode(holder, child, name.c_str(), createPlacement(*it));
0720         }
0721       }
0722     }
0723   } catch (cms::Exception& exception) {
0724     edm::LogInfo("FWRecoGeometry") << "failed to produce GEM geometry " << exception.what() << std::endl;
0725   }
0726 }
0727 
0728 //______________________________________________________________________________
0729 
0730 void FWTGeoRecoGeometryESProducer::addRPCGeometry() {
0731   TGeoVolume* tv = GetTopHolder("Muon", kMuonRPC);
0732   TGeoVolume* assembly = GetDaughter(tv, "RPC", kMuonRPC);
0733 
0734   DetId detId(DetId::Muon, MuonSubdetId::RPC);
0735   const RPCGeometry* rpcGeom = static_cast<const RPCGeometry*>(m_trackingGeom->slaveGeometry(detId));
0736   for (auto it = rpcGeom->rolls().begin(), end = rpcGeom->rolls().end(); it != end; ++it) {
0737     RPCRoll const* roll = (*it);
0738     if (roll) {
0739       RPCDetId detid = roll->geographicalId();
0740       std::stringstream s;
0741       s << detid;
0742       std::string name = s.str();
0743 
0744       TGeoVolume* child = createVolume(name, roll, kMuonRPC);
0745 
0746       TGeoVolume* holder = GetDaughter(assembly, "ROLL Region", kMuonRPC, detid.region());
0747       holder = GetDaughter(holder, "Ring", kMuonRPC, detid.ring());
0748       holder = GetDaughter(holder, "Station", kMuonRPC, detid.station());
0749       holder = GetDaughter(holder, "Sector", kMuonRPC, detid.sector());
0750       holder = GetDaughter(holder, "Layer", kMuonRPC, detid.layer());
0751       holder = GetDaughter(holder, "Subsector", kMuonRPC, detid.subsector());
0752 
0753       AddLeafNode(holder, child, name.c_str(), createPlacement(*it));
0754     }
0755   };
0756 }
0757 
0758 void FWTGeoRecoGeometryESProducer::addME0Geometry() {
0759   TGeoVolume* tv = GetTopHolder("Muon", kMuonME0);
0760   TGeoVolume* assembly = GetDaughter(tv, "ME0", kMuonME0);
0761 
0762   DetId detId(DetId::Muon, 5);
0763   try {
0764     const ME0Geometry* me0Geom = static_cast<const ME0Geometry*>(m_trackingGeom->slaveGeometry(detId));
0765 
0766     for (auto roll : me0Geom->etaPartitions()) {
0767       if (roll) {
0768         unsigned int rawid = roll->geographicalId().rawId();
0769         //std::cout << "AMT FWTTTTRecoGeometryES\n" << rawid ;
0770 
0771         ME0DetId detid(rawid);
0772         std::stringstream s;
0773         s << detid;
0774         std::string name = s.str();
0775         TGeoVolume* child = createVolume(name, roll, kMuonME0);
0776 
0777         TGeoVolume* holder = GetDaughter(assembly, "Region", kMuonME0, detid.region());
0778         holder = GetDaughter(holder, "Layer", kMuonME0, detid.layer());
0779         holder = GetDaughter(holder, "Chamber", kMuonME0, detid.chamber());
0780         AddLeafNode(holder, child, name.c_str(), createPlacement(roll));
0781       }
0782     }
0783   } catch (cms::Exception& exception) {
0784     edm::LogInfo("FWRecoGeometry") << "failed to produce ME0 geometry " << exception.what() << std::endl;
0785   }
0786 }
0787 
0788 //==============================================================================
0789 //================================= CALO =======================================
0790 //==============================================================================
0791 
0792 void FWTGeoRecoGeometryESProducer::addHcalCaloGeometryBarrel(void) {
0793   TGeoVolume* tv = GetTopHolder("HCal", kHCal);
0794   TGeoVolume* assembly = GetDaughter(tv, "HCalBarrel", kHCal);
0795 
0796   std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Hcal, HcalSubdetector::HcalBarrel);
0797 
0798   CaloVolMap caloShapeMapP;
0799   CaloVolMap caloShapeMapN;
0800   for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
0801     //HcalDetId detid = HcalDetId(it->rawId());
0802     HcalDetId detid(*it);
0803     const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
0804     const IdealObliquePrism* cell = dynamic_cast<const IdealObliquePrism*>(cellb);
0805 
0806     if (!cell) {
0807       printf("HB not oblique !!!\n");
0808       continue;
0809     }
0810 
0811     TGeoVolume* volume = nullptr;
0812     CaloVolMap& caloShapeMap = (cell->etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
0813     CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
0814     if (volIt == caloShapeMap.end()) {
0815       // printf("FIREWORKS NEW SHAPE BEGIN eta = %f etaPos = %f, phiPos %f >>>>>> \n", cell->eta(), cell->etaPos(), cell->phiPos());
0816       IdealObliquePrism::Pt3DVec lc(8);
0817       IdealObliquePrism::Pt3D ref;
0818       IdealObliquePrism::localCorners(lc, cell->param(), ref);
0819       HepGeom::Vector3D<float> lCenter;
0820       for (int c = 0; c < 8; ++c)
0821         lCenter += lc[c];
0822       lCenter *= 0.125;
0823 
0824       static const int arr[] = {1, 0, 3, 2, 5, 4, 7, 6};
0825       double points[16];
0826       for (int c = 0; c < 8; ++c) {
0827         if (cell->etaPos() > 0)
0828           points[c * 2 + 0] = -(lc[arr[c]].z() - lCenter.z());
0829         else
0830           points[c * 2 + 0] = (lc[arr[c]].z() - lCenter.z());
0831         points[c * 2 + 1] = (lc[arr[c]].y() - lCenter.y());
0832         // printf("AMT xy[%d] <=>[%d] = (%.4f, %.4f) \n", arr[c], c, points[c*2],  points[c*2+1]);
0833       }
0834 
0835       float dz = (lc[4].x() - lc[0].x()) * 0.5;
0836       TGeoShape* solid = new TGeoArb8(dz, &points[0]);
0837       volume = new TGeoVolume("hcal oblique prism", solid, GetMedium(kHCal));
0838       caloShapeMap[cell->param()] = volume;
0839     } else {
0840       volume = volIt->second;
0841     }
0842 
0843     HepGeom::Vector3D<float> gCenter;
0844     CaloCellGeometry::CornersVec const& gc = cell->getCorners();
0845     for (int c = 0; c < 8; ++c)
0846       gCenter += HepGeom::Vector3D<float>(gc[c].x(), gc[c].y(), gc[c].z());
0847     gCenter *= 0.125;
0848 
0849     TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
0850     TGeoRotation rot;
0851     rot.RotateY(90);
0852 
0853     TGeoRotation rotPhi;
0854     rotPhi.SetAngles(0, -cell->phiPos() * TMath::RadToDeg(), 0);
0855     rot.MultiplyBy(&rotPhi);
0856 
0857     TGeoVolume* holder = GetDaughter(assembly, "side", kHCal, detid.zside());
0858     holder = GetDaughter(holder, "ieta", kHCal, detid.ieta());
0859     std::stringstream nname;
0860     nname << detid;
0861     AddLeafNode(holder, volume, nname.str().c_str(), new TGeoCombiTrans(gtr, rot));
0862   }
0863 
0864   //  printf("HB map size P = %lu , N = %lu", caloShapeMapP.size(),caloShapeMapN.size() );
0865 }
0866 //______________________________________________________________________________
0867 
0868 void FWTGeoRecoGeometryESProducer::addHcalCaloGeometryEndcap(void) {
0869   CaloVolMap caloShapeMapP;
0870   CaloVolMap caloShapeMapN;
0871 
0872   TGeoVolume* tv = GetTopHolder("HCal", kHCal);
0873   TGeoVolume* assembly = GetDaughter(tv, "HCalEndcap", kHCal);
0874 
0875   std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Hcal, HcalSubdetector::HcalEndcap);
0876 
0877   for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
0878     HcalDetId detid = HcalDetId(it->rawId());
0879     const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
0880     const IdealObliquePrism* cell = dynamic_cast<const IdealObliquePrism*>(cellb);
0881 
0882     if (!cell) {
0883       printf("EC not oblique \n");
0884       continue;
0885     }
0886 
0887     TGeoVolume* volume = nullptr;
0888     CaloVolMap& caloShapeMap = (cell->etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
0889     CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
0890     if (volIt == caloShapeMap.end()) {
0891       IdealObliquePrism::Pt3DVec lc(8);
0892       IdealObliquePrism::Pt3D ref;
0893       IdealObliquePrism::localCorners(lc, cell->param(), ref);
0894       HepGeom::Vector3D<float> lCenter;
0895       for (int c = 0; c < 8; ++c)
0896         lCenter += lc[c];
0897       lCenter *= 0.125;
0898 
0899       //for( int c = 0; c < 8; ++c)
0900       //   printf("lc.push_back(TEveVector(%.4f, %.4f, %.4f));\n", lc[c].x(), lc[c].y(), lc[c].z() );
0901 
0902       static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
0903       static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
0904       const int* arr = (detid.ieta() > 0) ? &arrP[0] : &arrN[0];
0905 
0906       double points[16];
0907       for (int c = 0; c < 8; ++c) {
0908         points[c * 2 + 0] = lc[arr[c]].x() - lCenter.x();
0909         points[c * 2 + 1] = lc[arr[c]].y() - lCenter.y();
0910       }
0911 
0912       float dz = (lc[4].z() - lc[0].z()) * 0.5;
0913       TGeoShape* solid = new TGeoArb8(dz, &points[0]);
0914       volume = new TGeoVolume("ecal oblique prism", solid, GetMedium(kHCal));
0915       caloShapeMap[cell->param()] = volume;
0916     } else {
0917       volume = volIt->second;
0918     }
0919 
0920     HepGeom::Vector3D<float> gCenter;
0921     CaloCellGeometry::CornersVec const& gc = cell->getCorners();
0922     for (int c = 0; c < 8; ++c) {
0923       gCenter += HepGeom::Vector3D<float>(gc[c].x(), gc[c].y(), gc[c].z());
0924       //  printf("gc.push_back(TEveVector(%.4f, %.4f, %.4f));\n", gc[c].x(), gc[c].y(),gc[c].z() );
0925     }
0926     gCenter *= 0.125;
0927 
0928     TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
0929     TGeoRotation rot;
0930     rot.SetAngles(cell->phiPos() * TMath::RadToDeg(), 0, 0);
0931 
0932     TGeoVolume* holder = GetDaughter(assembly, "side", kHCal, detid.zside());
0933     holder = GetDaughter(holder, "ieta", kHCal, detid.ieta());
0934     std::stringstream nname;
0935     nname << detid;
0936     AddLeafNode(holder, volume, nname.str().c_str(), new TGeoCombiTrans(gtr, rot));
0937   }
0938 
0939   //   printf("HE map size P = %lu , N = %lu", caloShapeMapP.size(),caloShapeMapN.size() );
0940 }
0941 
0942 void FWTGeoRecoGeometryESProducer::addHcalCaloGeometryOuter() {
0943   CaloVolMap caloShapeMapP;
0944   CaloVolMap caloShapeMapN;
0945 
0946   TGeoVolume* tv = GetTopHolder("HCal", kHCal);
0947   TGeoVolume* assembly = GetDaughter(tv, "HCalOuter", kHCal);
0948 
0949   std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Hcal, HcalSubdetector::HcalOuter);
0950 
0951   for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
0952     HcalDetId detid = HcalDetId(it->rawId());
0953     const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
0954     const IdealObliquePrism* cell = dynamic_cast<const IdealObliquePrism*>(cellb);
0955 
0956     if (!cell) {
0957       printf("EC not oblique \n");
0958       continue;
0959     }
0960 
0961     TGeoVolume* volume = nullptr;
0962     CaloVolMap& caloShapeMap = (cell->etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
0963     CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
0964     if (volIt == caloShapeMap.end()) {
0965       IdealObliquePrism::Pt3DVec lc(8);
0966       IdealObliquePrism::Pt3D ref;
0967       IdealObliquePrism::localCorners(lc, cell->param(), ref);
0968       HepGeom::Vector3D<float> lCenter;
0969       for (int c = 0; c < 8; ++c)
0970         lCenter += lc[c];
0971       lCenter *= 0.125;
0972       static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
0973       static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
0974       const int* arr = (detid.ieta() > 0) ? &arrP[0] : &arrN[0];
0975 
0976       double points[16];
0977       for (int c = 0; c < 8; ++c) {
0978         points[c * 2 + 0] = lc[arr[c]].x() - lCenter.x();
0979         points[c * 2 + 1] = lc[arr[c]].y() - lCenter.y();
0980       }
0981 
0982       float dz = (lc[4].z() - lc[0].z()) * 0.5;
0983       TGeoShape* solid = new TGeoArb8(dz, &points[0]);
0984       volume = new TGeoVolume("ecal oblique prism", solid, GetMedium(kHCal));
0985       caloShapeMap[cell->param()] = volume;
0986     } else {
0987       volume = volIt->second;
0988     }
0989     HepGeom::Vector3D<float> gCenter;
0990     CaloCellGeometry::CornersVec const& gc = cell->getCorners();
0991     for (int c = 0; c < 8; ++c) {
0992       gCenter += HepGeom::Vector3D<float>(gc[c].x(), gc[c].y(), gc[c].z());
0993     }
0994     gCenter *= 0.125;
0995 
0996     TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
0997     TGeoRotation rot;
0998     rot.SetAngles(cell->phiPos() * TMath::RadToDeg(), 0, 0);
0999 
1000     TGeoVolume* holder = GetDaughter(assembly, "side", kHCal, detid.zside());
1001     holder = GetDaughter(holder, "ieta", kHCal, detid.ieta());
1002     std::stringstream nname;
1003     nname << detid;
1004     AddLeafNode(holder, volume, nname.str().c_str(), new TGeoCombiTrans(gtr, rot));
1005   }
1006 }
1007 
1008 void FWTGeoRecoGeometryESProducer::addHcalCaloGeometryForward() {
1009   CaloVolMap caloShapeMapP;
1010   CaloVolMap caloShapeMapN;
1011 
1012   TGeoVolume* tv = GetTopHolder("HCal", kHCal);
1013   TGeoVolume* assembly = GetDaughter(tv, "HCalForward", kHCal);
1014 
1015   std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Hcal, HcalSubdetector::HcalForward);
1016 
1017   for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
1018     HcalDetId detid = HcalDetId(it->rawId());
1019     const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
1020     const IdealZPrism* cell = dynamic_cast<const IdealZPrism*>(cellb);
1021 
1022     if (!cell) {
1023       printf("EC not Z prism \n");
1024       continue;
1025     }
1026 
1027     TGeoVolume* volume = nullptr;
1028     CaloVolMap& caloShapeMap = (cell->etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
1029     CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
1030     if (volIt == caloShapeMap.end()) {
1031       IdealZPrism::Pt3DVec lc(8);
1032       IdealZPrism::Pt3D ref;
1033       IdealZPrism::localCorners(lc, cell->param(), ref);
1034       HepGeom::Vector3D<float> lCenter;
1035       for (int c = 0; c < 8; ++c)
1036         lCenter += lc[c];
1037       lCenter *= 0.125;
1038       static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
1039       static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
1040       const int* arr = (detid.ieta() > 0) ? &arrP[0] : &arrN[0];
1041 
1042       double points[16];
1043       for (int c = 0; c < 8; ++c) {
1044         points[c * 2 + 0] = lc[arr[c]].x() - lCenter.x();
1045         points[c * 2 + 1] = lc[arr[c]].y() - lCenter.y();
1046       }
1047 
1048       float dz = (lc[4].z() - lc[0].z()) * 0.5;
1049       TGeoShape* solid = new TGeoArb8(dz, &points[0]);
1050       volume = new TGeoVolume("ecal oblique prism", solid, GetMedium(kHCal));
1051       caloShapeMap[cell->param()] = volume;
1052     } else {
1053       volume = volIt->second;
1054     }
1055     HepGeom::Vector3D<float> gCenter;
1056     CaloCellGeometry::CornersVec const& gc = cell->getCorners();
1057     for (int c = 0; c < 8; ++c) {
1058       gCenter += HepGeom::Vector3D<float>(gc[c].x(), gc[c].y(), gc[c].z());
1059     }
1060     gCenter *= 0.125;
1061 
1062     TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
1063     TGeoRotation rot;
1064     rot.SetAngles(cell->phiPos() * TMath::RadToDeg(), 0, 0);
1065 
1066     TGeoVolume* holder = GetDaughter(assembly, "side", kHCal, detid.zside());
1067     holder = GetDaughter(holder, "ieta", kHCal, detid.ieta());
1068     std::stringstream nname;
1069     nname << detid;
1070     AddLeafNode(holder, volume, nname.str().c_str(), new TGeoCombiTrans(gtr, rot));
1071   }
1072 }
1073 
1074 void FWTGeoRecoGeometryESProducer::addCaloTowerGeometry() {
1075   CaloVolMap caloShapeMapP;
1076   CaloVolMap caloShapeMapN;
1077 
1078   TGeoVolume* tv = GetTopHolder("CaloTower", kCaloTower);
1079   TGeoVolume* assembly = GetDaughter(tv, "CaloTower", kCaloTower);
1080 
1081   std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Calo, CaloTowerDetId::SubdetId);
1082   for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
1083     CaloTowerDetId detid = CaloTowerDetId(it->rawId());
1084     const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
1085     const IdealObliquePrism* cell = dynamic_cast<const IdealObliquePrism*>(cellb);
1086     if (!cell) {
1087       printf("EC not oblique \n");
1088       continue;
1089     }
1090     TGeoVolume* volume = nullptr;
1091     CaloVolMap& caloShapeMap = (cell->etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
1092     CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
1093     if (volIt == caloShapeMap.end()) {
1094       IdealObliquePrism::Pt3DVec lc(8);
1095       IdealObliquePrism::Pt3D ref;
1096       IdealObliquePrism::localCorners(lc, cell->param(), ref);
1097       HepGeom::Vector3D<float> lCenter;
1098       for (int c = 0; c < 8; ++c)
1099         lCenter += lc[c];
1100       lCenter *= 0.125;
1101 
1102       static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
1103       static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
1104       const int* arr = (detid.ieta() > 0) ? &arrP[0] : &arrN[0];
1105 
1106       double points[16];
1107       for (int c = 0; c < 8; ++c) {
1108         points[c * 2 + 0] = lc[arr[c]].x() - lCenter.x();
1109         points[c * 2 + 1] = lc[arr[c]].y() - lCenter.y();
1110       }
1111 
1112       float dz = (lc[4].z() - lc[0].z()) * 0.5;
1113       TGeoShape* solid = new TGeoArb8(dz, &points[0]);
1114       volume = new TGeoVolume("ecal oblique prism", solid, GetMedium(kCaloTower));
1115       caloShapeMap[cell->param()] = volume;
1116     } else {
1117       volume = volIt->second;
1118     }
1119 
1120     HepGeom::Vector3D<float> gCenter;
1121     CaloCellGeometry::CornersVec const& gc = cell->getCorners();
1122     for (int c = 0; c < 8; ++c) {
1123       gCenter += HepGeom::Vector3D<float>(gc[c].x(), gc[c].y(), gc[c].z());
1124     }
1125     gCenter *= 0.125;
1126 
1127     TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
1128     TGeoRotation rot;
1129     rot.SetAngles(cell->phiPos() * TMath::RadToDeg(), 0, 0);
1130 
1131     TGeoVolume* holder = GetDaughter(assembly, "side", kCaloTower, detid.zside());
1132     holder = GetDaughter(holder, "ieta", kCaloTower, detid.ieta());
1133     std::stringstream nname;
1134     nname << detid;
1135     AddLeafNode(holder, volume, nname.str().c_str(), new TGeoCombiTrans(gtr, rot));
1136   }
1137 }
1138 
1139 //______________________________________________________________________________
1140 
1141 TGeoHMatrix* getEcalTrans(CaloCellGeometry::CornersVec const& gc) {
1142   TEveVector gCenter;
1143   for (int i = 0; i < 8; ++i)
1144     gCenter += TEveVector(gc[i].x(), gc[i].y(), gc[i].z());
1145   gCenter *= 0.125;
1146 
1147   TEveVector tgCenter;  // top center 4 corners
1148   for (int i = 4; i < 8; ++i)
1149     tgCenter += TEveVector(gc[i].x(), gc[i].y(), gc[i].z());
1150   tgCenter *= 0.25;
1151 
1152   TEveVector axis = tgCenter - gCenter;
1153   axis.Normalize();
1154 
1155   TEveTrans tr;
1156   TVector3 v1t;
1157   tr.GetBaseVec(1, v1t);
1158 
1159   TEveVector v1(v1t.x(), v1t.y(), v1t.z());
1160   double dot13 = axis.Dot(v1);
1161   TEveVector gd = axis;
1162   gd *= dot13;
1163   v1 -= gd;
1164   v1.Normalize();
1165   TEveVector v2;
1166   TMath::Cross(v1.Arr(), axis.Arr(), v2.Arr());
1167   TMath::Cross(axis.Arr(), v1.Arr(), v2.Arr());
1168   v2.Normalize();
1169 
1170   tr.SetBaseVec(1, v1.fX, v1.fY, v1.fZ);
1171   tr.SetBaseVec(2, v2.fX, v2.fY, v2.fZ);
1172   tr.SetBaseVec(3, axis.fX, axis.fY, axis.fZ);
1173   tr.Move3PF(gCenter.fX, gCenter.fY, gCenter.fZ);
1174 
1175   TGeoHMatrix* out = new TGeoHMatrix();
1176   tr.SetGeoHMatrix(*out);
1177   return out;
1178 }
1179 
1180 TGeoShape* makeEcalShape(const TruncatedPyramid* cell) {
1181   // printf("BEGIN EE SHAPE --------------------------------\n");
1182   // std::cout << detid << std::endl;
1183   const HepGeom::Transform3D idtr;
1184   TruncatedPyramid::Pt3DVec co(8);
1185   TruncatedPyramid::Pt3D ref;
1186   TruncatedPyramid::localCorners(co, cell->param(), ref);
1187   //for( int c = 0; c < 8; ++c)
1188   //   printf("lc.push_back(TEveVector(%.4f, %.4f, %.4f));\n", co[c].x(),co[c].y(),co[c].z() );
1189 
1190   double points[16];
1191   for (int c = 0; c < 8; ++c) {
1192     points[c * 2] = co[c].x();
1193     points[c * 2 + 1] = co[c].y();
1194   }
1195   TGeoShape* solid = new TGeoArb8(cell->param()[0], points);
1196   return solid;
1197 }
1198 
1199 //______________________________________________________________________________
1200 
1201 void FWTGeoRecoGeometryESProducer::addEcalCaloGeometry(void) {
1202   TGeoVolume* tv = GetTopHolder("ECal", kECal);
1203   CaloVolMap caloShapeMap;
1204 
1205   {
1206     TGeoVolume* assembly = GetDaughter(tv, "ECalBarrel", kECal);
1207 
1208     std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Ecal, EcalSubdetector::EcalBarrel);
1209     for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
1210       EBDetId detid(*it);
1211       const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
1212       const TruncatedPyramid* cell = dynamic_cast<const TruncatedPyramid*>(cellb);
1213       if (!cell) {
1214         printf("ecalBarrel cell not a TruncatedPyramid !!\n");
1215         return;
1216       }
1217 
1218       TGeoVolume* volume = nullptr;
1219       CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
1220       if (volIt == caloShapeMap.end()) {
1221         volume = new TGeoVolume("EE TruncatedPyramid", makeEcalShape(cell), GetMedium(kECal));
1222         caloShapeMap[cell->param()] = volume;
1223       } else {
1224         volume = volIt->second;
1225       }
1226       TGeoHMatrix* mtx = getEcalTrans(cell->getCorners());
1227       TGeoVolume* holder = GetDaughter(assembly, "side", kECal, detid.zside());
1228       holder = GetDaughter(holder, "ieta", kECal, detid.ieta());
1229       std::stringstream nname;
1230       nname << detid;
1231       AddLeafNode(holder, volume, nname.str().c_str(), mtx);
1232     }
1233   }
1234 
1235   {
1236     TGeoVolume* assembly = GetDaughter(tv, "ECalEndcap", kECal);
1237 
1238     std::vector<DetId> vid = m_caloGeom->getValidDetIds(DetId::Ecal, EcalSubdetector::EcalEndcap);
1239     for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
1240       EEDetId detid(*it);
1241       const CaloCellGeometry* cellb = (m_caloGeom->getGeometry(*it)).get();
1242       const TruncatedPyramid* cell = dynamic_cast<const TruncatedPyramid*>(cellb);
1243       if (!cell) {
1244         printf("ecalEndcap cell not a TruncatedPyramid !!\n");
1245         continue;
1246       }
1247 
1248       TGeoVolume* volume = nullptr;
1249       CaloVolMap::iterator volIt = caloShapeMap.find(cell->param());
1250       if (volIt == caloShapeMap.end()) {
1251         volume = new TGeoVolume("EE TruncatedPyramid", makeEcalShape(cell), GetMedium(kECal));
1252         caloShapeMap[cell->param()] = volume;
1253       } else {
1254         volume = volIt->second;
1255       }
1256       TGeoHMatrix* mtx = getEcalTrans(cell->getCorners());
1257       TGeoVolume* holder = GetDaughter(assembly, "side", kECal, detid.zside());
1258       holder = GetDaughter(holder, "ix", kECal, detid.ix());
1259       std::stringstream nname;
1260       nname << detid;
1261       AddLeafNode(holder, volume, nname.str().c_str(), mtx);
1262     }
1263   }
1264 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!