CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
DisplayGeom.cc
Go to the documentation of this file.
1 
5 
12 
14 
17 
19 
21 
23 
25 
26 #include <boost/algorithm/string/case_conv.hpp>
27 
28 #include "TEveManager.h"
29 #include "TEveTrack.h"
30 #include "TEveTrackPropagator.h"
31 
32 #include "TGeoManager.h"
33 #include "TGeoMatrix.h"
34 #include "TEveGeoNode.h"
35 #include "TEveTrans.h"
36 #include "TEveScene.h"
37 #include "TEveViewer.h"
38 #include "TVector.h"
39 #include "TGLScenePad.h"
40 #include "TGLRnrCtx.h"
41 #include "TEvePointSet.h"
42 #include "TRandom.h"
43 #include "TEveUtil.h"
44 
45 #include "TEveQuadSet.h"
46 #include "TEveStraightLineSet.h"
47 #include "TEveRGBAPalette.h"
48 #include "TSystem.h"
49 #include "TStyle.h"
50 // class decleration
51 //
52 
54 public:
55  explicit DisplayGeom(const edm::ParameterSet&);
56  ~DisplayGeom() override;
57 
58 protected:
59  TEveGeoTopNode* make_node(const TString& path, Int_t vis_level, Bool_t global_cs);
60 
61 private:
62  void beginJob() override;
63  void analyze(const edm::Event&, const edm::EventSetup&) override;
64 
65  void endJob() override;
66 
68 
69  TEveElement* m_geomList;
70 
71  int m_level;
72 
73  bool m_MF;
75  std::vector<double> m_MF_plane_d0;
76  std::vector<double> m_MF_plane_d1;
77  std::vector<double> m_MF_plane_d2;
82 
84 
87 
88  void remakeGeometry(const DisplayGeomRecord& dgRec);
89 };
90 
92 
94  : m_eve(),
95  m_geomList(nullptr),
96  m_MF_component(0),
97  m_geomWatcher(this, &DisplayGeom::remakeGeometry),
98  m_displayGeomToken(esConsumes()) {
99  m_level = iConfig.getUntrackedParameter<int>("level", 2); // Geometry level to visualize
100 
101  m_MF = iConfig.getUntrackedParameter<int>("MF", false); // Show the MF geometry, instead of detector geometry
102 
103  std::string component = iConfig.getUntrackedParameter<std::string>("MF_component", "NONE");
104  boost::algorithm::to_upper(component);
105 
106  if (component == "NONE") {
107  m_MF_component = -1;
108  } else if (component == "ABSBZ") {
109  m_MF_component = 1;
110  } else if (component == "ABSBR") {
111  m_MF_component = 2;
112  } else if (component == "ABSBPHI") {
113  m_MF_component = 3;
114  } else if (component == "BR") {
115  m_MF_component = 4;
116  } else if (component == "BPHI") {
117  m_MF_component = 5;
118  } else { // Anything else -> |B|
119  m_MF_component = 0;
120  }
121 
122  if (m_MF_component != -1) {
124  }
125 
126  if (m_MF_component == 0) {
127  m_MF_plane_d0 = iConfig.getUntrackedParameter<std::vector<double> >("MF_plane_d0", std::vector<double>(3, 0.0));
128  m_MF_plane_d1 = iConfig.getParameter<std::vector<double> >("MF_plane_d1");
129  m_MF_plane_d2 = iConfig.getParameter<std::vector<double> >("MF_plane_d2");
130 
131  m_MF_plane_N1 = iConfig.getUntrackedParameter<UInt_t>("MF_plane_N", 100);
132  m_MF_plane_N2 = iConfig.getUntrackedParameter<UInt_t>("MF_plane_N2", m_MF_plane_N1);
133 
134  m_MF_plane_draw_dir = iConfig.getUntrackedParameter<int>("MF_plane_draw_dir", true);
135  m_MF_isPickable = iConfig.getUntrackedParameter<bool>("MF_pickable", true);
136  }
137 }
138 
140 
141 //==============================================================================
142 // Protected helpers
143 //==============================================================================
144 
145 TEveGeoTopNode* DisplayGeom::make_node(const TString& path, Int_t vis_level, Bool_t global_cs) {
146  if (!gGeoManager->cd(path)) {
147  Warning("make_node", "Path '%s' not found.", path.Data());
148  return nullptr;
149  }
150 
151  TEveGeoTopNode* tn = new TEveGeoTopNode(gGeoManager, gGeoManager->GetCurrentNode());
152  tn->SetVisLevel(vis_level);
153  if (global_cs) {
154  tn->RefMainTrans().SetFrom(*gGeoManager->GetCurrentMatrix());
155  }
156  m_geomList->AddElement(tn);
157  gEve->AddToListTree(tn, true);
158  return tn;
159 }
160 
161 //==============================================================================
162 // member functions
163 //==============================================================================
164 
166  if (m_eve) {
167  // Remake geometry if it has changed.
168  m_geomWatcher.check(iSetup);
169 
170  if (m_MF_component != -1) {
171  MagneticField const& field = iSetup.getData(m_magFieldToken);
172 
173  gStyle->SetPalette(1, nullptr);
174 
175  int minval = 0;
176  int maxval = 4000;
177  if (m_MF_component == 1) { //AbsBZ
178  minval = 0, maxval = 4000;
179  } else if (m_MF_component == 2) { //AbsBR
180  minval = 0, maxval = 4000;
181  } else if (m_MF_component == 3) { //AbsBphi
182  minval = 0, maxval = 1000;
183  } else if (m_MF_component == 4) { //BR
184  minval = -4000, maxval = 4000;
185  } else if (m_MF_component == 5) { //Bphi
186  minval = -1000, maxval = 1000;
187  }
188 
189  TEveRGBAPalette* pal = new TEveRGBAPalette(minval, maxval);
190 
191  TEveStraightLineSet* ls = nullptr;
192  if (m_MF_plane_draw_dir) {
193  ls = new TEveStraightLineSet("MF_line_direction");
194  ls->SetPickable(false);
195  ls->SetLineColor(kGreen);
196  ls->SetMarkerColor(kGreen);
197  ls->SetMarkerStyle(1);
198  }
199 
200  TEveQuadSet* q = new TEveQuadSet("MF_quad_values");
201  q->Reset(TEveQuadSet::kQT_RectangleXY, kFALSE, 32);
202  q->SetOwnIds(kTRUE);
203  q->SetAlwaysSecSelect(true);
204  q->SetPickable(m_MF_isPickable);
205  q->SetPalette(pal);
206 
207  TEveVectorD v0(m_MF_plane_d0[0], m_MF_plane_d0[1], m_MF_plane_d0[2]);
208  TEveVectorD v01(m_MF_plane_d1[0], m_MF_plane_d1[1], m_MF_plane_d1[2]);
209  TEveVectorD v02(m_MF_plane_d2[0], m_MF_plane_d2[1], m_MF_plane_d2[2]);
210 
211  TEveVectorD b01 = (v01 - v0);
212  TEveVectorD b02 = (v02 - v0);
213  TEveVectorD b03 = b01.Cross(b02);
214 
215  TEveTrans trans;
216  trans.SetBaseVec(1, b01.fX, b01.fY, b01.fZ);
217  trans.SetBaseVec(2, b02.fX, b02.fY, b02.fZ);
218  trans.SetBaseVec(3, b03.fX, b03.fY, b03.fZ);
219  trans.SetPos(v0.Arr());
220  trans.OrtoNorm3();
221  q->SetTransMatrix(trans.Array());
222 
223  double w_step = b01.Mag() / m_MF_plane_N1;
224  double h_step = b02.Mag() / m_MF_plane_N2;
225 
226  q->SetDefWidth(w_step);
227  q->SetDefHeight(h_step);
228  TEveVectorD d1;
229  trans.GetBaseVec(1).GetXYZ(d1);
230  d1 *= w_step;
231  TEveVectorD d2;
232  trans.GetBaseVec(2).GetXYZ(d2);
233  d2 *= h_step;
234 
235  //d1.Print();
236  d2.Dump();
237  double line_step_size = TMath::Min(w_step, h_step);
238 
239  for (int i = 0; i < m_MF_plane_N1; i++) {
240  for (int j = 0; j < m_MF_plane_N2; j++) {
241  TEveVectorD p = d1 * Double_t(i) + d2 * Double_t(j) + v0;
242  GlobalPoint pos(p.fX, p.fY, p.fZ);
243  GlobalVector b = field.inTesla(pos) * 1000.; // in mT
244  float value = 0.;
245  if (m_MF_component == 0) { //BMOD
246  value = b.mag();
247  } else if (m_MF_component == 1) { //BZ
248  value = fabs(b.z());
249  } else if (m_MF_component == 2) { //ABSBR
250  value = fabs(b.x() * cos(pos.phi()) + b.y() * sin(pos.phi()));
251  } else if (m_MF_component == 3) { //ABSBPHI
252  value = fabs(-b.x() * sin(pos.phi()) + b.y() * cos(pos.phi()));
253  } else if (m_MF_component == 2) { //BR
254  value = b.x() * cos(pos.phi()) + b.y() * sin(pos.phi());
255  } else if (m_MF_component == 5) { //BPHI
256  value = -b.x() * sin(pos.phi()) + b.y() * cos(pos.phi());
257  }
258 
259  q->AddQuad(w_step * i, h_step * j);
260  q->QuadValue(value);
261  if (m_MF_isPickable)
262  q->QuadId(new TNamed(Form("Mag (%f, %f, %f) val = %f", b.x(), b.y(), b.z(), b.mag()), "Dong!"));
263 
264  if (ls) {
265  if (b.mag() > 1e-6) {
266  b.unit();
267  b *= line_step_size;
268  ls->AddLine(p.fX, p.fY, p.fZ, p.fX + b.x(), p.fY + b.y(), p.fZ + b.z());
269  } else {
270  ls->AddLine(p.fX, p.fY, p.fZ, p.fX + b.x(), p.fY + b.y(), p.fZ + b.z());
271  }
272 
273  ls->AddMarker(ls->GetLinePlex().Size() - 1, 0);
274  ls->AddMarker(i * m_MF_plane_N1 + j, 0);
275  }
276  }
277  }
278 
279  TEveScene* eps = gEve->SpawnNewScene("FillStyleScene");
280  gEve->GetDefaultViewer()->AddScene(eps);
281  eps->GetGLScene()->SetStyle(TGLRnrCtx::kFill);
282  eps->AddElement(q);
283  if (ls)
284  m_eve->AddElement(ls);
285  } else {
286  // // Add a test obj
287  // if (!gRandom)
288  // gRandom = new TRandom(0);
289  // TRandom& r= *gRandom;
290 
291  // Float_t s = 100;
292 
293  // TEvePointSet* ps = new TEvePointSet();
294  // ps->SetOwnIds(kTRUE);
295  // for(Int_t i = 0; i< 100; i++)
296  // {
297  // ps->SetNextPoint(r.Uniform(-s,s), r.Uniform(-s,s), r.Uniform(-s,s));
298  // ps->SetPointId(new TNamed(Form("Point %d", i), ""));
299  // }
300 
301  // ps->SetMarkerColor(TMath::Nint(r.Uniform(2, 9)));
302  // ps->SetMarkerSize(r.Uniform(1, 2));
303  // ps->SetMarkerStyle(4);
304  // m_eve->AddElement(ps);
305  }
306  }
307 }
308 
309 // ------------ method called once each job just before starting event loop ------------
311  if (m_eve) {
312  m_geomList = new TEveElementList("Display Geom");
314  // m_eve->getManager()->GetGlobalScene()->GetGLScene()->SetStyle(TGLRnrCtx::kWireFrame);
315  }
316 }
317 
318 // ------------ method called once each job just after ending the event loop ------------
320 
321 //------------------------------------------------------------------------------
323  m_geomList->DestroyElements();
324 
325  TEveGeoManagerHolder _tgeo(const_cast<TGeoManager*>(&dgRec.get(m_displayGeomToken)));
326 
327  // To have a full one, all detectors in one top-node:
328  // make_node("/cms:World_1/cms:CMSE_1", 4, kTRUE);
329 
330  if (m_MF) {
331  make_node("/cms:World_1", m_level, kTRUE);
332  } else {
333  make_node("/cms:World_1/cms:CMSE_1/tracker:Tracker_1", m_level, kTRUE);
334  make_node("/cms:World_1/cms:CMSE_1/caloBase:CALO_1", m_level, kTRUE);
335  make_node("/cms:World_1/cms:CMSE_1/muonBase:MUON_1", m_level, kTRUE);
336  }
337 }
T getUntrackedParameter(std::string const &, T const &) const
int m_MF_plane_N2
Definition: DisplayGeom.cc:79
void remakeGeometry(const DisplayGeomRecord &dgRec)
Definition: DisplayGeom.cc:322
WorkSpace int float eps
void beginJob() override
Definition: DisplayGeom.cc:310
DisplayGeom(const edm::ParameterSet &)
Definition: DisplayGeom.cc:93
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
def ls
Definition: eostools.py:349
T y() const
Definition: PV3DBase.h:60
int m_MF_component
Definition: DisplayGeom.cc:74
int m_MF_plane_draw_dir
Definition: DisplayGeom.cc:80
bool m_MF_isPickable
Definition: DisplayGeom.cc:81
T Min(T a, T b)
Definition: MathUtil.h:39
std::vector< double > m_MF_plane_d1
Definition: DisplayGeom.cc:76
edm::Service< EveService > m_eve
Definition: DisplayGeom.cc:67
bool getData(T &iHolder) const
Definition: EventSetup.h:128
TEveGeoTopNode * make_node(const TString &path, Int_t vis_level, Bool_t global_cs)
Definition: DisplayGeom.cc:145
int iEvent
Definition: GenABIO.cc:224
T mag() const
Definition: PV3DBase.h:64
std::vector< double > m_MF_plane_d2
Definition: DisplayGeom.cc:77
ProductT const & get(ESGetToken< ProductT, DepRecordT > const &iToken) const
T z() const
Definition: PV3DBase.h:61
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > m_magFieldToken
Definition: DisplayGeom.cc:85
Vector3DBase unit() const
Definition: Vector3DBase.h:54
~DisplayGeom() override
Definition: DisplayGeom.cc:139
void endJob() override
Definition: DisplayGeom.cc:319
std::vector< double > m_MF_plane_d0
Definition: DisplayGeom.cc:75
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
double b
Definition: hdecay.h:118
bool check(const edm::EventSetup &iSetup)
Definition: ESWatcher.h:57
int m_MF_plane_N1
Definition: DisplayGeom.cc:78
void AddElement(TEveElement *el)
Definition: EveService.cc:261
edm::ESWatcher< DisplayGeomRecord > m_geomWatcher
Definition: DisplayGeom.cc:83
static constexpr float d1
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: DisplayGeom.cc:165
T x() const
Definition: PV3DBase.h:59
ESGetTokenH3DDVariant esConsumes(std::string const &Reccord, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
TEveElement * m_geomList
Definition: DisplayGeom.cc:69
void AddGlobalElement(TEveElement *el)
Definition: EveService.cc:263
const edm::ESGetToken< TGeoManager, DisplayGeomRecord > m_displayGeomToken
Definition: DisplayGeom.cc:86