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   static void fillDescriptions(edm::ConfigurationDescriptions&);
0059
0060 protected:
0061   TEveGeoTopNode* make_node(const TString& path, Int_t vis_level, Bool_t global_cs);
0062
0063 private:
0064   void beginJob() override;
0065   void analyze(const edm::Event&, const edm::EventSetup&) override;
0066
0067   void endJob() override;
0068
0069   edm::Service<EveService> m_eve;
0070
0071   TEveElement* m_geomList;
0072
0073   int m_level;
0074
0075   typedef std::vector<std::string> vstring;
0076   vstring m_nodes;
0077
0078   int m_MF_component;
0079   std::vector<double> m_MF_plane_d0;
0080   std::vector<double> m_MF_plane_d1;
0081   std::vector<double> m_MF_plane_d2;
0082   int m_MF_plane_N1;
0083   int m_MF_plane_N2;
0084   int m_MF_plane_draw_dir;
0085   bool m_MF_isPickable;
0086
0087   edm::ESWatcher<DisplayGeomRecord> m_geomWatcher;
0088
0089   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> m_magFieldToken;
0090   const edm::ESGetToken<TGeoManager, DisplayGeomRecord> m_displayGeomToken;
0091
0092   void remakeGeometry(const DisplayGeomRecord& dgRec);
0093 };
0094
0095 DEFINE_FWK_MODULE(DisplayGeom);
0096
0097 DisplayGeom::DisplayGeom(const edm::ParameterSet& iConfig)
0098     : m_eve(),
0099       m_geomList(nullptr),
0100       m_level(iConfig.getUntrackedParameter<int>("level")),
0101       m_nodes(iConfig.getUntrackedParameter<vstring>("nodes")),
0102       m_MF_component(0),
0103       m_geomWatcher(this, &DisplayGeom::remakeGeometry),
0104       m_displayGeomToken(esConsumes()) {
0105   std::string component = iConfig.getUntrackedParameter<std::string>("MF_component");
0106   boost::algorithm::to_upper(component);
0107
0108   if (component == "NONE") {
0109     m_MF_component = -1;
0110   } else if (component == "ABSBZ") {
0111     m_MF_component = 1;
0112   } else if (component == "ABSBR") {
0113     m_MF_component = 2;
0114   } else if (component == "ABSBPHI") {
0115     m_MF_component = 3;
0116   } else if (component == "BR") {
0117     m_MF_component = 4;
0118   } else if (component == "BPHI") {
0119     m_MF_component = 5;
0120   } else {  // Anything else -> |B|
0121     m_MF_component = 0;
0122   }
0123
0124   if (m_MF_component != -1) {
0125     m_magFieldToken = esConsumes();
0126
0127     m_MF_plane_d0 = iConfig.getUntrackedParameter<std::vector<double> >("MF_plane_d0");
0128     m_MF_plane_d1 = iConfig.getUntrackedParameter<std::vector<double> >("MF_plane_d1");
0129     m_MF_plane_d2 = iConfig.getUntrackedParameter<std::vector<double> >("MF_plane_d2");
0130
0131     m_MF_plane_N1 = iConfig.getUntrackedParameter<int>("MF_plane_N");
0132     m_MF_plane_N2 = iConfig.getUntrackedParameter<int>("MF_plane_N2");
0133     if (m_MF_plane_N2 < 0)
0134       m_MF_plane_N2 = m_MF_plane_N1;
0135
0136     m_MF_plane_draw_dir = iConfig.getUntrackedParameter<int>("MF_plane_draw_dir");
0137     m_MF_isPickable = iConfig.getUntrackedParameter<bool>("MF_pickable");
0138   }
0139 }
0140
0141 DisplayGeom::~DisplayGeom() {}
0142
0143 //==============================================================================
0144 // Protected helpers
0145 //==============================================================================
0146
0147 TEveGeoTopNode* DisplayGeom::make_node(const TString& path, Int_t vis_level, Bool_t global_cs) {
0148   if (!gGeoManager->cd(path)) {
0149     Warning("make_node", "Path '%s' not found.", path.Data());
0150     return nullptr;
0151   }
0152
0153   TEveGeoTopNode* tn = new TEveGeoTopNode(gGeoManager, gGeoManager->GetCurrentNode());
0154   tn->SetVisLevel(vis_level);
0155   if (global_cs) {
0156     tn->RefMainTrans().SetFrom(*gGeoManager->GetCurrentMatrix());
0157   }
0158   m_geomList->AddElement(tn);
0159   gEve->AddToListTree(tn, true);
0160   return tn;
0161 }
0162
0163 //==============================================================================
0164 // member functions
0165 //==============================================================================
0166
0167 void DisplayGeom::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0168   if (m_eve) {
0169     // Remake geometry if it has changed.
0170     m_geomWatcher.check(iSetup);
0171
0172     if (m_MF_component != -1) {
0173       MagneticField const& field = iSetup.getData(m_magFieldToken);
0174
0175       gStyle->SetPalette(1, nullptr);
0176
0177       int minval = 0;
0178       int maxval = 4000;
0179       if (m_MF_component == 1) {  //AbsBZ
0180         minval = 0, maxval = 4000;
0181       } else if (m_MF_component == 2) {  //AbsBR
0182         minval = 0, maxval = 4000;
0183       } else if (m_MF_component == 3) {  //AbsBphi
0184         minval = 0, maxval = 1000;
0185       } else if (m_MF_component == 4) {  //BR
0186         minval = -4000, maxval = 4000;
0187       } else if (m_MF_component == 5) {  //Bphi
0188         minval = -200, maxval = 200;
0189       }
0190
0191       TEveRGBAPalette* pal = new TEveRGBAPalette(minval, maxval);
0192
0193       TEveStraightLineSet* ls = nullptr;
0194       if (m_MF_plane_draw_dir) {
0195         ls = new TEveStraightLineSet("MF_line_direction");
0196         ls->SetPickable(false);
0197         ls->SetLineColor(kGreen);
0198         ls->SetMarkerColor(kGreen);
0199         ls->SetMarkerStyle(1);
0200       }
0201
0202       TEveQuadSet* q = new TEveQuadSet("MF_quad_values");
0203       q->Reset(TEveQuadSet::kQT_RectangleXY, kFALSE, 32);
0204       q->SetOwnIds(kTRUE);
0205       q->SetAlwaysSecSelect(true);
0206       q->SetPickable(m_MF_isPickable);
0207       q->SetPalette(pal);
0208
0209       TEveVectorD v0(m_MF_plane_d0[0], m_MF_plane_d0[1], m_MF_plane_d0[2]);
0210       TEveVectorD v01(m_MF_plane_d1[0], m_MF_plane_d1[1], m_MF_plane_d1[2]);
0211       TEveVectorD v02(m_MF_plane_d2[0], m_MF_plane_d2[1], m_MF_plane_d2[2]);
0212
0213       TEveVectorD b01 = (v01 - v0);
0214       TEveVectorD b02 = (v02 - v0);
0215       TEveVectorD b03 = b01.Cross(b02);
0216
0217       TEveTrans trans;
0218       trans.SetBaseVec(1, b01.fX, b01.fY, b01.fZ);
0219       trans.SetBaseVec(2, b02.fX, b02.fY, b02.fZ);
0220       trans.SetBaseVec(3, b03.fX, b03.fY, b03.fZ);
0221       trans.SetPos(v0.Arr());
0222       trans.OrtoNorm3();
0223       q->SetTransMatrix(trans.Array());
0224
0225       double w_step = b01.Mag() / m_MF_plane_N1;
0226       double h_step = b02.Mag() / m_MF_plane_N2;
0227
0228       q->SetDefWidth(w_step);
0229       q->SetDefHeight(h_step);
0230       TEveVectorD d1;
0231       trans.GetBaseVec(1).GetXYZ(d1);
0232       d1 *= w_step;
0233       TEveVectorD d2;
0234       trans.GetBaseVec(2).GetXYZ(d2);
0235       d2 *= h_step;
0236
0237       //d1.Print();
0238       d2.Dump();
0239       double line_step_size = TMath::Min(w_step, h_step);
0240
0241       for (int i = 0; i < m_MF_plane_N1; i++) {
0242         for (int j = 0; j < m_MF_plane_N2; j++) {
0243           TEveVectorD p = d1 * Double_t(i) + d2 * Double_t(j) + v0;
0244           GlobalPoint pos(p.fX, p.fY, p.fZ);
0245           GlobalVector b = field.inTesla(pos) * 1000.;  // in mT
0246           float value = 0.;
0247           if (m_MF_component == 0) {  //BMOD
0248             value = b.mag();
0249           } else if (m_MF_component == 1) {  //BZ
0250             value = fabs(b.z());
0251           } else if (m_MF_component == 2) {  //ABSBR
0252             value = fabs(b.x() * cos(pos.phi()) + b.y() * sin(pos.phi()));
0253           } else if (m_MF_component == 3) {  //ABSBPHI
0254             value = fabs(-b.x() * sin(pos.phi()) + b.y() * cos(pos.phi()));
0255           } else if (m_MF_component == 2) {  //BR
0256             value = b.x() * cos(pos.phi()) + b.y() * sin(pos.phi());
0257           } else if (m_MF_component == 5) {  //BPHI
0258             value = -b.x() * sin(pos.phi()) + b.y() * cos(pos.phi());
0259           }
0260
0261           q->AddQuad(w_step * i, h_step * j);
0262           q->QuadValue(value);
0263           if (m_MF_isPickable)
0264             q->QuadId(new TNamed(Form("Mag (%f, %f, %f) val = %f", b.x(), b.y(), b.z(), b.mag()), "Dong!"));
0265
0266           if (ls) {
0267             if (b.mag() > 1e-6) {
0268               b.unit();
0269               b *= line_step_size;
0270               ls->AddLine(p.fX, p.fY, p.fZ, p.fX + b.x(), p.fY + b.y(), p.fZ + b.z());
0271             } else {
0272               ls->AddLine(p.fX, p.fY, p.fZ, p.fX + b.x(), p.fY + b.y(), p.fZ + b.z());
0273             }
0274
0275             ls->AddMarker(ls->GetLinePlex().Size() - 1, 0);
0276             ls->AddMarker(i * m_MF_plane_N1 + j, 0);
0277           }
0278         }
0279       }
0280
0281       TEveScene* eps = gEve->SpawnNewScene("MF Map");
0282       gEve->GetDefaultViewer()->AddScene(eps);
0283       eps->GetGLScene()->SetStyle(TGLRnrCtx::kFill);
0284       eps->AddElement(q);
0285       if (ls)
0286         m_eve->AddElement(ls);
0287     } else {
0288       //     // Add a test obj
0289       //     if (!gRandom)
0290       //       gRandom = new TRandom(0);
0291       //     TRandom& r= *gRandom;
0292
0293       //     Float_t s = 100;
0294
0295       //     TEvePointSet* ps = new TEvePointSet();
0296       //     ps->SetOwnIds(kTRUE);
0297       //     for(Int_t i = 0; i< 100; i++)
0298       //       {
0299       //         ps->SetNextPoint(r.Uniform(-s,s), r.Uniform(-s,s), r.Uniform(-s,s));
0300       //         ps->SetPointId(new TNamed(Form("Point %d", i), ""));
0301       //       }
0302
0303       //     ps->SetMarkerColor(TMath::Nint(r.Uniform(2, 9)));
0304       //     ps->SetMarkerSize(r.Uniform(1, 2));
0305       //     ps->SetMarkerStyle(4);
0306       //     m_eve->AddElement(ps);
0307     }
0308   }
0309   m_eve->getManager()->FullRedraw3D(true, true);
0310 }
0311
0312 // ------------ method called once each job just before starting event loop  ------------
0313 void DisplayGeom::beginJob() {
0314   if (m_eve) {
0315     m_geomList = new TEveElementList("Display Geom");
0316     m_eve->AddGlobalElement(m_geomList);
0317     //      m_eve->getManager()->GetGlobalScene()->GetGLScene()->SetStyle(TGLRnrCtx::kWireFrame);
0318   }
0319 }
0320
0321 // ------------ method called once each job just after ending the event loop  ------------
0322 void DisplayGeom::endJob() {}
0323
0324 //------------------------------------------------------------------------------
0325 void DisplayGeom::remakeGeometry(const DisplayGeomRecord& dgRec) {
0326   m_geomList->DestroyElements();
0327
0328   TEveGeoManagerHolder _tgeo(const_cast<TGeoManager*>(&dgRec.get(m_displayGeomToken)));
0329
0330   for (std::string& aNode : m_nodes) {
0331     make_node(aNode, m_level, kTRUE);
0332   }
0333 }
0334
0335 void DisplayGeom::fillDescriptions(edm::ConfigurationDescriptions& conf) {
0336   edm::ParameterSetDescription desc;
0337   desc.addUntracked<vstring>("nodes", vstring{"tracker:Tracker_1", "muonBase:MUON_1", "caloBase:CALO_1"})
0338       ->setComment("List of nodes to visualize");
0339   ;
0340   desc.addUntracked<int>("level", 4)->setComment("Levels into the geometry hierarchy visualized at startup");
0341   desc.addUntracked<std::string>("MF_component", "None")
0342       ->setComment("Component of the MF to show: 'None', 'B', 'AbsBZ', 'AbsBR', 'AbsBphi', 'BR', 'Bphi'");
0343   desc.addUntracked<std::vector<double> >("MF_plane_d0", std::vector<double>{0., -900., -2400.})
0344       ->setComment("1st corner of MF map");
0345   desc.addUntracked<std::vector<double> >("MF_plane_d1", std::vector<double>{0., -900., 2400.})
0346       ->setComment("2nd corner of MF map");
0347   desc.addUntracked<std::vector<double> >("MF_plane_d2", std::vector<double>{0., 900., -2400.})
0348       ->setComment("3rd corner of MF map");
0349   desc.addUntracked<int>("MF_plane_N", 200)->setComment("Number of bins for the MF map in 1st coord");
0350   desc.addUntracked<int>("MF_plane_N2", -1)->setComment("Number of bins for the MF map in 2nd coord");
0351   desc.addUntracked<int>("MF_plane_draw_dir", true)->setComment("Draw MF direction arrows (slow)");
0352   desc.addUntracked<bool>("MF_pickable", false)->setComment("MF values are pickable (slow)");
0353
0354   conf.add("DisplayGeom", desc);
0355 }