47 #include "TGeoManager.h"
49 #include "TGeoMatrix.h"
56 #include "TEveVector.h"
57 #include "TEveTrans.h"
61 typedef std::map<const float*, TGeoVolume*> CaloVolMap;
64 m_tracker =
pset.getUntrackedParameter<
bool>(
"Tracker",
true);
65 m_muon =
pset.getUntrackedParameter<
bool>(
"Muon",
true);
66 m_calo =
pset.getUntrackedParameter<
bool>(
"Calo",
true);
84 void AddLeafNode(TGeoVolume* mother, TGeoVolume* daughter,
const char*
name, TGeoMatrix* mtx) {
85 int n = mother->GetNdaughters();
86 mother->AddNode(daughter, 1, mtx);
87 mother->GetNode(
n)->SetName(
name);
97 TGeoTranslation trans(posx, posy, posz);
104 const Double_t
matrix[9] = {detRot.
xx(),
115 return new TGeoCombiTrans(trans,
rotation);
122 TGeoVolume*
res =
nullptr;
123 if (mother->GetNdaughters()) {
124 TGeoNode*
n = mother->FindNode(Form(
"%s_%d_1",
prefix,
id));
126 res =
n->GetVolume();
130 res =
new TGeoVolumeAssembly(Form(
"%s_%d",
prefix,
id));
132 mother->AddNode(
res, 1);
139 TGeoVolume*
res =
nullptr;
140 if (mother->GetNdaughters()) {
141 TGeoNode*
n = mother->FindNode(Form(
"%s_1",
prefix));
143 res =
n->GetVolume();
150 mother->AddNode(
res, 1);
182 std::map<ERecoDet, TGeoMedium*>::iterator it =
m_recoMedium.find(det);
193 color = GMCol::Green;
203 color = GMCol::Blue2;
213 color = GMCol::Yellow1;
223 color = GMCol::Yellow0;
229 color = GMCol::Blue2;
233 color = GMCol::Orange1;
237 color = GMCol::Green;
241 color = GMCol::Blue2;
245 color = GMCol::Blue1;
248 printf(
"invalid medium id \n");
252 TGeoMaterial* mat =
new TGeoMaterial(
name.c_str(), 0, 0, 0);
255 mat->SetFillStyle(3000);
266 auto fwTGeoRecoGeometry = std::make_unique<FWTGeoRecoGeometry>();
272 TGeoManager*
geom =
new TGeoManager(
"cmsGeo",
"CMS Detector");
273 if (
nullptr == gGeoIdentity) {
274 gGeoIdentity =
new TGeoIdentity(
"Identity");
277 fwTGeoRecoGeometry->manager(
geom);
280 TGeoMaterial* vacuum =
new TGeoMaterial(
"Vacuum", 0, 0, 0);
285 if (
nullptr == top) {
286 return std::unique_ptr<FWTGeoRecoGeometry>();
288 geom->SetTopVolume(top);
290 top->SetVisibility(kFALSE);
291 top->SetLineColor(kBlue);
329 geom->CloseGeometry();
331 geom->DefaultColors();
333 geom->CloseGeometry();
335 return fwTGeoRecoGeometry;
340 TGeoShape* shape =
nullptr;
346 std::array<const float, 4>
const& par =
b2->parameters();
349 float hBottomEdge = par[0];
350 float hTopEdge = par[1];
352 float apothem = par[3];
355 s <<
"T_" << hBottomEdge <<
"_" << hTopEdge <<
"_" <<
thickness <<
"_" << apothem;
361 if (
nullptr == shape) {
362 shape =
new TGeoTrap(
name.c_str(),
378 if (dynamic_cast<const RectangularPlaneBounds*>(
b) !=
nullptr) {
391 if (
nullptr == shape) {
405 std::map<TGeoShape*, TGeoVolume*>::iterator vIt =
m_shapeToVolume.find(solid);
409 TGeoVolume* volume =
new TGeoVolume(
name.c_str(), solid,
GetMedium(mid));
430 DetId detid = it->geographicalId();
450 DetId detid = it->geographicalId();
456 std::string name = Form(
"PXF D:%d, B:%d, P:%d, S:%d", disk, blade, panel, side);
473 DetId detid = it->geographicalId();
493 DetId detid = it->geographicalId();
513 DetId detid = it->geographicalId();
533 DetId detid = it->geographicalId();
564 for (
auto it = dtChamberGeom.begin(),
end = dtChamberGeom.end(); it !=
end; ++it) {
565 if (
auto chamber = dynamic_cast<const DTChamber*>(*it)) {
585 for (
auto it = dtSuperLayerGeom.begin(),
end = dtSuperLayerGeom.end(); it !=
end; ++it) {
586 if (
auto* superlayer = dynamic_cast<const DTSuperLayer*>(*it)) {
606 for (
auto it = dtLayerGeom.begin(),
end = dtLayerGeom.end(); it !=
end; ++it) {
607 if (
auto layer = dynamic_cast<const DTLayer*>(*it)) {
630 throw cms::Exception(
"FatalError") <<
"Cannnot find CSCGeometry\n";
636 for (
auto it = cscGeom.begin(), itEnd = cscGeom.end(); it != itEnd; ++it) {
637 unsigned int rawid = (*it)->geographicalId();
643 TGeoVolume*
child =
nullptr;
645 if (
auto chamber = dynamic_cast<const CSCChamber*>(*it))
647 else if (
auto* layer = dynamic_cast<const CSCLayer*>(*it))
684 TGeoVolume* holder =
GetDaughter(assembly,
"SuperChamber Region",
kMuonGEM, detid.region());
717 edm::LogInfo(
"FWRecoGeometry") <<
"failed to produce GEM geometry " <<
exception.what() << std::endl;
729 for (
auto it = rpcGeom->
rolls().begin(),
end = rpcGeom->
rolls().end(); it !=
end; ++it) {
761 unsigned int rawid = roll->geographicalId().rawId();
777 edm::LogInfo(
"FWRecoGeometry") <<
"failed to produce ME0 geometry " <<
exception.what() << std::endl;
791 CaloVolMap caloShapeMapP;
792 CaloVolMap caloShapeMapN;
793 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
800 printf(
"HB not oblique !!!\n");
804 TGeoVolume* volume =
nullptr;
805 CaloVolMap& caloShapeMap = (cell->
etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
806 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
807 if (volIt == caloShapeMap.end()) {
812 HepGeom::Vector3D<float> lCenter;
813 for (
int c = 0;
c < 8; ++
c)
817 static const int arr[] = {1, 0, 3, 2, 5, 4, 7, 6};
819 for (
int c = 0;
c < 8; ++
c) {
821 points[
c * 2 + 0] = -(lc[arr[
c]].
z() - lCenter.z());
823 points[
c * 2 + 0] = (lc[arr[
c]].z() - lCenter.z());
824 points[
c * 2 + 1] = (lc[arr[
c]].y() - lCenter.y());
828 float dz = (lc[4].x() - lc[0].x()) * 0.5;
829 TGeoShape* solid =
new TGeoArb8(
dz, &
points[0]);
830 volume =
new TGeoVolume(
"hcal oblique prism", solid,
GetMedium(
kHCal));
831 caloShapeMap[cell->
param()] = volume;
833 volume = volIt->second;
836 HepGeom::Vector3D<float> gCenter;
838 for (
int c = 0;
c < 8; ++
c)
839 gCenter += HepGeom::Vector3D<float>(gc[
c].
x(), gc[
c].
y(), gc[
c].
z());
842 TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
847 rotPhi.SetAngles(0, -cell->
phiPos() * TMath::RadToDeg(), 0);
848 rot.MultiplyBy(&rotPhi);
852 std::stringstream nname;
854 AddLeafNode(holder, volume, nname.str().c_str(),
new TGeoCombiTrans(gtr,
rot));
862 CaloVolMap caloShapeMapP;
863 CaloVolMap caloShapeMapN;
870 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
876 printf(
"EC not oblique \n");
880 TGeoVolume* volume =
nullptr;
881 CaloVolMap& caloShapeMap = (cell->
etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
882 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
883 if (volIt == caloShapeMap.end()) {
887 HepGeom::Vector3D<float> lCenter;
888 for (
int c = 0;
c < 8; ++
c)
895 static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
896 static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
897 const int* arr = (detid.
ieta() > 0) ? &arrP[0] : &arrN[0];
900 for (
int c = 0;
c < 8; ++
c) {
901 points[
c * 2 + 0] = lc[arr[
c]].x() - lCenter.x();
902 points[
c * 2 + 1] = lc[arr[
c]].y() - lCenter.y();
905 float dz = (lc[4].z() - lc[0].z()) * 0.5;
906 TGeoShape* solid =
new TGeoArb8(
dz, &
points[0]);
907 volume =
new TGeoVolume(
"ecal oblique prism", solid,
GetMedium(
kHCal));
908 caloShapeMap[cell->
param()] = volume;
910 volume = volIt->second;
913 HepGeom::Vector3D<float> gCenter;
915 for (
int c = 0;
c < 8; ++
c) {
916 gCenter += HepGeom::Vector3D<float>(gc[
c].
x(), gc[
c].
y(), gc[
c].
z());
921 TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
923 rot.SetAngles(cell->
phiPos() * TMath::RadToDeg(), 0, 0);
927 std::stringstream nname;
929 AddLeafNode(holder, volume, nname.str().c_str(),
new TGeoCombiTrans(gtr,
rot));
936 CaloVolMap caloShapeMapP;
937 CaloVolMap caloShapeMapN;
944 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
950 printf(
"EC not oblique \n");
954 TGeoVolume* volume =
nullptr;
955 CaloVolMap& caloShapeMap = (cell->
etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
956 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
957 if (volIt == caloShapeMap.end()) {
961 HepGeom::Vector3D<float> lCenter;
962 for (
int c = 0;
c < 8; ++
c)
965 static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
966 static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
967 const int* arr = (detid.
ieta() > 0) ? &arrP[0] : &arrN[0];
970 for (
int c = 0;
c < 8; ++
c) {
971 points[
c * 2 + 0] = lc[arr[
c]].x() - lCenter.x();
972 points[
c * 2 + 1] = lc[arr[
c]].y() - lCenter.y();
975 float dz = (lc[4].z() - lc[0].z()) * 0.5;
976 TGeoShape* solid =
new TGeoArb8(
dz, &
points[0]);
977 volume =
new TGeoVolume(
"ecal oblique prism", solid,
GetMedium(
kHCal));
978 caloShapeMap[cell->
param()] = volume;
980 volume = volIt->second;
982 HepGeom::Vector3D<float> gCenter;
984 for (
int c = 0;
c < 8; ++
c) {
985 gCenter += HepGeom::Vector3D<float>(gc[
c].
x(), gc[
c].
y(), gc[
c].
z());
989 TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
991 rot.SetAngles(cell->
phiPos() * TMath::RadToDeg(), 0, 0);
995 std::stringstream nname;
997 AddLeafNode(holder, volume, nname.str().c_str(),
new TGeoCombiTrans(gtr,
rot));
1002 CaloVolMap caloShapeMapP;
1003 CaloVolMap caloShapeMapN;
1010 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
1013 const IdealZPrism* cell = dynamic_cast<const IdealZPrism*>(cellb);
1016 printf(
"EC not Z prism \n");
1020 TGeoVolume* volume =
nullptr;
1021 CaloVolMap& caloShapeMap = (cell->
etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
1022 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
1023 if (volIt == caloShapeMap.end()) {
1027 HepGeom::Vector3D<float> lCenter;
1028 for (
int c = 0;
c < 8; ++
c)
1031 static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
1032 static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
1033 const int* arr = (detid.
ieta() > 0) ? &arrP[0] : &arrN[0];
1036 for (
int c = 0;
c < 8; ++
c) {
1037 points[
c * 2 + 0] = lc[arr[
c]].x() - lCenter.x();
1038 points[
c * 2 + 1] = lc[arr[
c]].y() - lCenter.y();
1041 float dz = (lc[4].z() - lc[0].z()) * 0.5;
1042 TGeoShape* solid =
new TGeoArb8(
dz, &
points[0]);
1043 volume =
new TGeoVolume(
"ecal oblique prism", solid,
GetMedium(
kHCal));
1044 caloShapeMap[cell->
param()] = volume;
1046 volume = volIt->second;
1048 HepGeom::Vector3D<float> gCenter;
1050 for (
int c = 0;
c < 8; ++
c) {
1051 gCenter += HepGeom::Vector3D<float>(gc[
c].
x(), gc[
c].
y(), gc[
c].
z());
1055 TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
1057 rot.SetAngles(cell->
phiPos() * TMath::RadToDeg(), 0, 0);
1061 std::stringstream nname;
1063 AddLeafNode(holder, volume, nname.str().c_str(),
new TGeoCombiTrans(gtr,
rot));
1068 CaloVolMap caloShapeMapP;
1069 CaloVolMap caloShapeMapN;
1075 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
1080 printf(
"EC not oblique \n");
1083 TGeoVolume* volume =
nullptr;
1084 CaloVolMap& caloShapeMap = (cell->
etaPos() > 0) ? caloShapeMapP : caloShapeMapN;
1085 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
1086 if (volIt == caloShapeMap.end()) {
1090 HepGeom::Vector3D<float> lCenter;
1091 for (
int c = 0;
c < 8; ++
c)
1095 static const int arrP[] = {3, 2, 1, 0, 7, 6, 5, 4};
1096 static const int arrN[] = {7, 6, 5, 4, 3, 2, 1, 0};
1097 const int* arr = (detid.
ieta() > 0) ? &arrP[0] : &arrN[0];
1100 for (
int c = 0;
c < 8; ++
c) {
1101 points[
c * 2 + 0] = lc[arr[
c]].x() - lCenter.x();
1102 points[
c * 2 + 1] = lc[arr[
c]].y() - lCenter.y();
1105 float dz = (lc[4].z() - lc[0].z()) * 0.5;
1106 TGeoShape* solid =
new TGeoArb8(
dz, &
points[0]);
1108 caloShapeMap[cell->
param()] = volume;
1110 volume = volIt->second;
1113 HepGeom::Vector3D<float> gCenter;
1115 for (
int c = 0;
c < 8; ++
c) {
1116 gCenter += HepGeom::Vector3D<float>(gc[
c].
x(), gc[
c].
y(), gc[
c].
z());
1120 TGeoTranslation gtr(gCenter.x(), gCenter.y(), gCenter.z());
1122 rot.SetAngles(cell->
phiPos() * TMath::RadToDeg(), 0, 0);
1126 std::stringstream nname;
1128 AddLeafNode(holder, volume, nname.str().c_str(),
new TGeoCombiTrans(gtr,
rot));
1136 for (
int i = 0;
i < 8; ++
i)
1137 gCenter += TEveVector(gc[
i].
x(), gc[
i].
y(), gc[
i].
z());
1140 TEveVector tgCenter;
1141 for (
int i = 4;
i < 8; ++
i)
1142 tgCenter += TEveVector(gc[
i].
x(), gc[
i].
y(), gc[
i].
z());
1145 TEveVector axis = tgCenter - gCenter;
1150 tr.GetBaseVec(1, v1t);
1152 TEveVector v1(v1t.x(), v1t.y(), v1t.z());
1153 double dot13 = axis.Dot(v1);
1154 TEveVector gd = axis;
1159 TMath::Cross(v1.Arr(), axis.Arr(), v2.Arr());
1160 TMath::Cross(axis.Arr(), v1.Arr(), v2.Arr());
1163 tr.SetBaseVec(1, v1.fX, v1.fY, v1.fZ);
1164 tr.SetBaseVec(2, v2.fX, v2.fY, v2.fZ);
1165 tr.SetBaseVec(3, axis.fX, axis.fY, axis.fZ);
1166 tr.Move3PF(gCenter.fX, gCenter.fY, gCenter.fZ);
1168 TGeoHMatrix*
out =
new TGeoHMatrix();
1169 tr.SetGeoHMatrix(*
out);
1176 const HepGeom::Transform3D idtr;
1184 for (
int c = 0;
c < 8; ++
c) {
1188 TGeoShape* solid =
new TGeoArb8(cell->
param()[0],
points);
1196 CaloVolMap caloShapeMap;
1202 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
1205 const TruncatedPyramid* cell = dynamic_cast<const TruncatedPyramid*>(cellb);
1207 printf(
"ecalBarrel cell not a TruncatedPyramid !!\n");
1211 TGeoVolume* volume =
nullptr;
1212 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
1213 if (volIt == caloShapeMap.end()) {
1215 caloShapeMap[cell->
param()] = volume;
1217 volume = volIt->second;
1222 std::stringstream nname;
1224 AddLeafNode(holder, volume, nname.str().c_str(), mtx);
1232 for (std::vector<DetId>::const_iterator it =
vid.begin(),
end =
vid.end(); it !=
end; ++it) {
1235 const TruncatedPyramid* cell = dynamic_cast<const TruncatedPyramid*>(cellb);
1237 printf(
"ecalEndcap cell not a TruncatedPyramid !!\n");
1241 TGeoVolume* volume =
nullptr;
1242 CaloVolMap::iterator volIt = caloShapeMap.find(cell->
param());
1243 if (volIt == caloShapeMap.end()) {
1245 caloShapeMap[cell->
param()] = volume;
1247 volume = volIt->second;
1252 std::stringstream nname;
1254 AddLeafNode(holder, volume, nname.str().c_str(), mtx);