Geometry/plugins/DisplayGeom.cc

0001
0002 #include "FWCore/Framework/interface/Frameworkfwd.h"
0003 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008 #include "FWCore/Framework/interface/Event.h"
0009 #include "FWCore/Framework/interface/EventSetup.h"
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/ESWatcher.h"
0012
0013 #include "FWCore/ServiceRegistry/interface/Service.h"
0014
0015 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0016 #include "MagneticField/Engine/interface/MagneticField.h"
0017
0018 #include "Fireworks/Eve/interface/EveService.h"
0019
0020 #include "Fireworks/Geometry/interface/DisplayGeomRecord.h"
0021
0022 #include "DataFormats/TrackReco/interface/Track.h"
0023
0024 #include "Fireworks/Tracks/interface/TrackUtils.h"
0025
0026 #include <boost/algorithm/string/case_conv.hpp>
0027
0028 #include "TEveManager.h"
0029 #include "TEveTrack.h"
0030 #include "TEveTrackPropagator.h"
0031
0032 #include "TGeoManager.h"
0033 #include "TGeoMatrix.h"
0034 #include "TEveGeoNode.h"
0035 #include "TEveTrans.h"
0036 #include "TEveScene.h"
0037 #include "TEveViewer.h"
0038 #include "TVector.h"
0039 #include "TGLScenePad.h"
0040 #include "TGLRnrCtx.h"
0041 #include "TEvePointSet.h"
0042 #include "TRandom.h"
0043 #include "TEveUtil.h"
0044
0045 #include "TEveQuadSet.h"
0046 #include "TEveStraightLineSet.h"
0047 #include "TEveRGBAPalette.h"
0048 #include "TSystem.h"
0049 #include "TStyle.h"
0050 // class decleration
0051 //
0052
0053 class DisplayGeom : public edm::one::EDAnalyzer<> {
0054 public:
0055   explicit DisplayGeom(const edm::ParameterSet&);
0056   ~DisplayGeom() override;
0057
0058 protected:
0059   TEveGeoTopNode* make_node(const TString& path, Int_t vis_level, Bool_t global_cs);
0060
0061 private:
0062   void beginJob() override;
0063   void analyze(const edm::Event&, const edm::EventSetup&) override;
0064
0065   void endJob() override;
0066
0067   edm::Service<EveService> m_eve;
0068
0069   TEveElement* m_geomList;
0070
0071   int m_level;
0072
0073   bool m_MF;
0074   int m_MF_component;
0075   std::vector<double> m_MF_plane_d0;
0076   std::vector<double> m_MF_plane_d1;
0077   std::vector<double> m_MF_plane_d2;
0078   int m_MF_plane_N1;
0079   int m_MF_plane_N2;
0080   int m_MF_plane_draw_dir;
0081   bool m_MF_isPickable;
0082
0083   edm::ESWatcher<DisplayGeomRecord> m_geomWatcher;
0084
0085   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_magFieldToken;
0086   const edm::ESGetToken<TGeoManager, DisplayGeomRecord> m_displayGeomToken;
0087
0088   void remakeGeometry(const DisplayGeomRecord& dgRec);
0089 };
0090
0091 DEFINE_FWK_MODULE(DisplayGeom);
0092
0093 DisplayGeom::DisplayGeom(const edm::ParameterSet& iConfig)
0094     : m_eve(),
0095       m_geomList(nullptr),
0096       m_MF_component(0),
0097       m_geomWatcher(this, &DisplayGeom::remakeGeometry),
0098       m_displayGeomToken(esConsumes()) {
0099   m_level = iConfig.getUntrackedParameter<int>("level", 2);  // Geometry level to visualize
0100
0101   m_MF = iConfig.getUntrackedParameter<int>("MF", false);  // Show the MF geometry, instead of detector geometry
0102
0103   std::string component = iConfig.getUntrackedParameter<std::string>("MF_component", "NONE");
0104   boost::algorithm::to_upper(component);
0105
0106   if (component == "NONE") {
0107     m_MF_component = -1;
0108   } else if (component == "ABSBZ") {
0109     m_MF_component = 1;
0110   } else if (component == "ABSBR") {
0111     m_MF_component = 2;
0112   } else if (component == "ABSBPHI") {
0113     m_MF_component = 3;
0114   } else if (component == "BR") {
0115     m_MF_component = 4;
0116   } else if (component == "BPHI") {
0117     m_MF_component = 5;
0118   } else {  // Anything else -> |B|
0119     m_MF_component = 0;
0120   }
0121
0122   if (m_MF_component != -1) {
0123     m_magFieldToken = esConsumes();
0124   }
0125
0126   if (m_MF_component == 0) {
0127     m_MF_plane_d0 = iConfig.getUntrackedParameter<std::vector<double> >("MF_plane_d0", std::vector<double>(3, 0.0));
0128     m_MF_plane_d1 = iConfig.getParameter<std::vector<double> >("MF_plane_d1");
0129     m_MF_plane_d2 = iConfig.getParameter<std::vector<double> >("MF_plane_d2");
0130
0131     m_MF_plane_N1 = iConfig.getUntrackedParameter<UInt_t>("MF_plane_N", 100);
0132     m_MF_plane_N2 = iConfig.getUntrackedParameter<UInt_t>("MF_plane_N2", m_MF_plane_N1);
0133
0134     m_MF_plane_draw_dir = iConfig.getUntrackedParameter<int>("MF_plane_draw_dir", true);
0135     m_MF_isPickable = iConfig.getUntrackedParameter<bool>("MF_pickable", true);
0136   }
0137 }
0138
0139 DisplayGeom::~DisplayGeom() {}
0140
0141 //==============================================================================
0142 // Protected helpers
0143 //==============================================================================
0144
0145 TEveGeoTopNode* DisplayGeom::make_node(const TString& path, Int_t vis_level, Bool_t global_cs) {
0146   if (!gGeoManager->cd(path)) {
0147     Warning("make_node", "Path '%s' not found.", path.Data());
0148     return nullptr;
0149   }
0150
0151   TEveGeoTopNode* tn = new TEveGeoTopNode(gGeoManager, gGeoManager->GetCurrentNode());
0152   tn->SetVisLevel(vis_level);
0153   if (global_cs) {
0154     tn->RefMainTrans().SetFrom(*gGeoManager->GetCurrentMatrix());
0155   }
0156   m_geomList->AddElement(tn);
0157   gEve->AddToListTree(tn, true);
0158   return tn;
0159 }
0160
0161 //==============================================================================
0162 // member functions
0163 //==============================================================================
0164
0165 void DisplayGeom::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0166   if (m_eve) {
0167     // Remake geometry if it has changed.
0168     m_geomWatcher.check(iSetup);
0169
0170     if (m_MF_component != -1) {
0171       MagneticField const& field = iSetup.getData(m_magFieldToken);
0172
0173       gStyle->SetPalette(1, nullptr);
0174
0175       int minval = 0;
0176       int maxval = 4000;
0177       if (m_MF_component == 1) {  //AbsBZ
0178         minval = 0, maxval = 4000;
0179       } else if (m_MF_component == 2) {  //AbsBR
0180         minval = 0, maxval = 4000;
0181       } else if (m_MF_component == 3) {  //AbsBphi
0182         minval = 0, maxval = 1000;
0183       } else if (m_MF_component == 4) {  //BR
0184         minval = -4000, maxval = 4000;
0185       } else if (m_MF_component == 5) {  //Bphi
0186         minval = -1000, maxval = 1000;
0187       }
0188
0189       TEveRGBAPalette* pal = new TEveRGBAPalette(minval, maxval);
0190
0191       TEveStraightLineSet* ls = nullptr;
0192       if (m_MF_plane_draw_dir) {
0193         ls = new TEveStraightLineSet("MF_line_direction");
0194         ls->SetPickable(false);
0195         ls->SetLineColor(kGreen);
0196         ls->SetMarkerColor(kGreen);
0197         ls->SetMarkerStyle(1);
0198       }
0199
0200       TEveQuadSet* q = new TEveQuadSet("MF_quad_values");
0201       q->Reset(TEveQuadSet::kQT_RectangleXY, kFALSE, 32);
0202       q->SetOwnIds(kTRUE);
0203       q->SetAlwaysSecSelect(true);
0204       q->SetPickable(m_MF_isPickable);
0205       q->SetPalette(pal);
0206
0207       TEveVectorD v0(m_MF_plane_d0[0], m_MF_plane_d0[1], m_MF_plane_d0[2]);
0208       TEveVectorD v01(m_MF_plane_d1[0], m_MF_plane_d1[1], m_MF_plane_d1[2]);
0209       TEveVectorD v02(m_MF_plane_d2[0], m_MF_plane_d2[1], m_MF_plane_d2[2]);
0210
0211       TEveVectorD b01 = (v01 - v0);
0212       TEveVectorD b02 = (v02 - v0);
0213       TEveVectorD b03 = b01.Cross(b02);
0214
0215       TEveTrans trans;
0216       trans.SetBaseVec(1, b01.fX, b01.fY, b01.fZ);
0217       trans.SetBaseVec(2, b02.fX, b02.fY, b02.fZ);
0218       trans.SetBaseVec(3, b03.fX, b03.fY, b03.fZ);
0219       trans.SetPos(v0.Arr());
0220       trans.OrtoNorm3();
0221       q->SetTransMatrix(trans.Array());
0222
0223       double w_step = b01.Mag() / m_MF_plane_N1;
0224       double h_step = b02.Mag() / m_MF_plane_N2;
0225
0226       q->SetDefWidth(w_step);
0227       q->SetDefHeight(h_step);
0228       TEveVectorD d1;
0229       trans.GetBaseVec(1).GetXYZ(d1);
0230       d1 *= w_step;
0231       TEveVectorD d2;
0232       trans.GetBaseVec(2).GetXYZ(d2);
0233       d2 *= h_step;
0234
0235       //d1.Print();
0236       d2.Dump();
0237       double line_step_size = TMath::Min(w_step, h_step);
0238
0239       for (int i = 0; i < m_MF_plane_N1; i++) {
0240         for (int j = 0; j < m_MF_plane_N2; j++) {
0241           TEveVectorD p = d1 * Double_t(i) + d2 * Double_t(j) + v0;
0242           GlobalPoint pos(p.fX, p.fY, p.fZ);
0243           GlobalVector b = field.inTesla(pos) * 1000.;  // in mT
0244           float value = 0.;
0245           if (m_MF_component == 0) {  //BMOD
0246             value = b.mag();
0247           } else if (m_MF_component == 1) {  //BZ
0248             value = fabs(b.z());
0249           } else if (m_MF_component == 2) {  //ABSBR
0250             value = fabs(b.x() * cos(pos.phi()) + b.y() * sin(pos.phi()));
0251           } else if (m_MF_component == 3) {  //ABSBPHI
0252             value = fabs(-b.x() * sin(pos.phi()) + b.y() * cos(pos.phi()));
0253           } else if (m_MF_component == 2) {  //BR
0254             value = b.x() * cos(pos.phi()) + b.y() * sin(pos.phi());
0255           } else if (m_MF_component == 5) {  //BPHI
0256             value = -b.x() * sin(pos.phi()) + b.y() * cos(pos.phi());
0257           }
0258
0259           q->AddQuad(w_step * i, h_step * j);
0260           q->QuadValue(value);
0261           if (m_MF_isPickable)
0262             q->QuadId(new TNamed(Form("Mag (%f, %f, %f) val = %f", b.x(), b.y(), b.z(), b.mag()), "Dong!"));
0263
0264           if (ls) {
0265             if (b.mag() > 1e-6) {
0266               b.unit();
0267               b *= line_step_size;
0268               ls->AddLine(p.fX, p.fY, p.fZ, p.fX + b.x(), p.fY + b.y(), p.fZ + b.z());
0269             } else {
0270               ls->AddLine(p.fX, p.fY, p.fZ, p.fX + b.x(), p.fY + b.y(), p.fZ + b.z());
0271             }
0272
0273             ls->AddMarker(ls->GetLinePlex().Size() - 1, 0);
0274             ls->AddMarker(i * m_MF_plane_N1 + j, 0);
0275           }
0276         }
0277       }
0278
0279       TEveScene* eps = gEve->SpawnNewScene("FillStyleScene");
0280       gEve->GetDefaultViewer()->AddScene(eps);
0281       eps->GetGLScene()->SetStyle(TGLRnrCtx::kFill);
0282       eps->AddElement(q);
0283       if (ls)
0284         m_eve->AddElement(ls);
0285     } else {
0286       //     // Add a test obj
0287       //     if (!gRandom)
0288       //       gRandom = new TRandom(0);
0289       //     TRandom& r= *gRandom;
0290
0291       //     Float_t s = 100;
0292
0293       //     TEvePointSet* ps = new TEvePointSet();
0294       //     ps->SetOwnIds(kTRUE);
0295       //     for(Int_t i = 0; i< 100; i++)
0296       //       {
0297       //         ps->SetNextPoint(r.Uniform(-s,s), r.Uniform(-s,s), r.Uniform(-s,s));
0298       //         ps->SetPointId(new TNamed(Form("Point %d", i), ""));
0299       //       }
0300
0301       //     ps->SetMarkerColor(TMath::Nint(r.Uniform(2, 9)));
0302       //     ps->SetMarkerSize(r.Uniform(1, 2));
0303       //     ps->SetMarkerStyle(4);
0304       //     m_eve->AddElement(ps);
0305     }
0306   }
0307 }
0308
0309 // ------------ method called once each job just before starting event loop  ------------
0310 void DisplayGeom::beginJob() {
0311   if (m_eve) {
0312     m_geomList = new TEveElementList("Display Geom");
0313     m_eve->AddGlobalElement(m_geomList);
0314     //      m_eve->getManager()->GetGlobalScene()->GetGLScene()->SetStyle(TGLRnrCtx::kWireFrame);
0315   }
0316 }
0317
0318 // ------------ method called once each job just after ending the event loop  ------------
0319 void DisplayGeom::endJob() {}
0320
0321 //------------------------------------------------------------------------------
0322 void DisplayGeom::remakeGeometry(const DisplayGeomRecord& dgRec) {
0323   m_geomList->DestroyElements();
0324
0325   TEveGeoManagerHolder _tgeo(const_cast<TGeoManager*>(&dgRec.get(m_displayGeomToken)));
0326
0327   // To have a full one, all detectors in one top-node:
0328   // make_node("/cms:World_1/cms:CMSE_1", 4, kTRUE);
0329
0330   if (m_MF) {
0331     make_node("/cms:World_1", m_level, kTRUE);
0332   } else {
0333     make_node("/cms:World_1/cms:CMSE_1/tracker:Tracker_1", m_level, kTRUE);
0334     make_node("/cms:World_1/cms:CMSE_1/caloBase:CALO_1", m_level, kTRUE);
0335     make_node("/cms:World_1/cms:CMSE_1/muonBase:MUON_1", m_level, kTRUE);
0336   }
0337 }