ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FastOpenGlDisplayer/src/OpenGLDisplayer.cc
Revision: 1.18
Committed: Fri Apr 18 11:31:21 2008 UTC (17 years ago) by querten
Content type: text/plain
Branch: MAIN
Changes since 1.17: +68 -4 lines
Log Message:
Add DT Geometry In

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: OpenGLDisplayer
4 // Class: OpenGLDisplayer
5 //
6 /**\class OpenGLDisplayer OpenGLDisplayer.cc Visualisation/OpenGLDisplayer/src/OpenGLDisplayer.cc
7
8 Description: <one line class summary>
9
10 Implementation:
11 <Notes on implementation>
12 */
13 //
14 // Original Author: Loic QUERTENMONT
15 // Created: Fri Oct 26 07:22:12 CEST 2007
16 // $Id: OpenGLDisplayer.cc,v 1.17 2008/04/16 08:21:47 querten Exp $
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23
24 // user include files
25 #include "FWCore/Framework/interface/Frameworkfwd.h"
26 #include "FWCore/Framework/interface/EDAnalyzer.h"
27
28 #include "FWCore/Framework/interface/Event.h"
29 #include "FWCore/Framework/interface/MakerMacros.h"
30
31 #include "FWCore/ParameterSet/interface/ParameterSet.h"
32 #include "FWCore/ServiceRegistry/interface/Service.h"
33 #include "FWCore/Framework/interface/ESHandle.h"
34
35 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
36 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
37 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
38 #include "Geometry/TrackerNumberingBuilder/interface/GeometricDet.h"
39 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
40 #include "Geometry/CommonTopologies/interface/StripTopology.h"
41 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetType.h"
42 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
43 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
44 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
45 #include "DataFormats/GeometrySurface/interface/BoundSurface.h"
46 #include "DataFormats/DetId/interface/DetId.h"
47
48 #include "Geometry/CaloEventSetup/interface/CaloTopologyRecord.h"
49 #include "Geometry/Records/interface/IdealGeometryRecord.h"
50 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
51 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
52 #include "Geometry/EcalBarrelAlgo/interface/EcalBarrelGeometry.h"
53 #include "Geometry/EcalEndcapAlgo/interface/EcalEndcapGeometry.h"
54 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
55 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
56
57 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
58 #include "Geometry/Records/interface/MuonGeometryRecord.h"
59 #include "Geometry/DTGeometry/interface/DTGeometry.h"
60
61
62 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
63 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
64 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
65
66 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
67 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
68 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
69 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
70 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
71 #include "SimDataFormats/Track/interface/SimTrack.h"
72 #include "SimDataFormats/Vertex/interface/SimVertex.h"
73
74 #include "DataFormats/Provenance/interface/BranchDescription.h"
75 #include "DataFormats/Provenance/interface/Provenance.h"
76 #include "DataFormats/Candidate/interface/Candidate.h"
77 #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
78 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
79 #include "DataFormats/TrackReco/interface/Track.h"
80
81 #include "DataFormats/Math/interface/Point3D.h"
82 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
83 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
84
85 #include "Visualisation/OpenGLDisplayer/interface/SimEvent.h"
86 #include "Visualisation/OpenGLDisplayer/interface/Geometry.h"
87
88 using namespace edm;
89 using namespace std;
90
91 //
92 // class decleration
93 //
94
95 class OpenGLDisplayer : public edm::EDAnalyzer {
96 public:
97 explicit OpenGLDisplayer(const edm::ParameterSet&);
98 ~OpenGLDisplayer();
99
100
101 private:
102 virtual void beginJob(const edm::EventSetup& iSetup);
103 virtual void analyze (const edm::Event&, const edm::EventSetup&);
104 virtual void endJob ();
105
106
107 MySimEvents* MyEvents;
108 Geometry* Geom_Tracker;
109 Geometry* Geom_ECAL;
110 Geometry* Geom_HCAL;
111 Geometry* Geom_Muon;
112
113
114 std::string OutputFile;
115
116 std::vector<InputTag> SimTrackProducers;
117 std::vector<InputTag> SimVertexProducers;
118 std::vector<InputTag> SimHitProducers;
119
120 std::vector<InputTag> TrackProducers;
121 std::vector<InputTag> EcalRecHitProducers;
122 std::vector<InputTag> HcalHBHERecHitProducers;
123 std::vector<InputTag> HcalHORecHitProducers;
124 std::vector<InputTag> HcalHFRecHitProducers;
125
126
127
128
129 // ----------member data ---------------------------
130 };
131
132 //
133 // constructors and destructor
134 //
135 OpenGLDisplayer::OpenGLDisplayer(const edm::ParameterSet& iConfig)
136 {
137 OutputFile = iConfig.getParameter<std::string >("OutputFile");
138
139 SimTrackProducers = iConfig.getParameter<std::vector<InputTag> >("SimTrackProducers");
140 SimVertexProducers = iConfig.getParameter<std::vector<InputTag> >("SimVertexProducers");
141 SimHitProducers = iConfig.getParameter<std::vector<InputTag> >("SimHitProducers");
142
143 TrackProducers = iConfig.getParameter<std::vector<InputTag> >("TrackProducers");
144 EcalRecHitProducers = iConfig.getParameter<std::vector<InputTag> >("EcalRecHitProducers");
145 HcalHBHERecHitProducers = iConfig.getParameter<std::vector<InputTag> >("HcalHBHERecHitProducers");
146 HcalHORecHitProducers = iConfig.getParameter<std::vector<InputTag> >("HcalHORecHitProducers");
147 HcalHFRecHitProducers = iConfig.getParameter<std::vector<InputTag> >("HcalHFRecHitProducers");
148
149 }
150
151
152 OpenGLDisplayer::~OpenGLDisplayer()
153 {
154 }
155
156 // ------------ method called once each job just before starting event loop ------------
157 void
158 OpenGLDisplayer::beginJob(const edm::EventSetup& iSetup)
159 {
160 DetId Detid;
161 int SubDet;
162
163 MyEvents = new MySimEvents();
164 Geom_Tracker = new Geometry();
165 Geom_ECAL = new Geometry();
166 Geom_HCAL = new Geometry();
167 Geom_Muon = new Geometry();
168
169
170 // ### TRACKER GEOMETRY ###
171
172 edm::ESHandle<TrackerGeometry> tkGeom;
173 iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
174 vector<GeomDet*> TrackerDets = tkGeom->dets();
175
176 for(unsigned int i=0;i<TrackerDets.size();i++){
177 Detid = TrackerDets[i]->geographicalId();
178 // SubDet = Detid.subdetId();
179
180 GeomDet* DetUnit = TrackerDets[i];
181 if(!DetUnit)continue;
182 const BoundPlane plane = DetUnit->surface();
183 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
184 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
185
186 float width = 0;
187 float length = 0;
188 float thickness = 0;
189 float TrapezoidalParam = 0;
190
191 if(trapezoidalBounds)
192 {
193 std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
194 width = parameters[0]*2;
195 length = parameters[3]*2;
196 thickness = (*trapezoidalBounds).thickness();
197 TrapezoidalParam = parameters[1]/parameters[0];
198 }else if(rectangularBounds){
199 width = DetUnit->surface().bounds().width();
200 length = DetUnit->surface().bounds().length();
201 thickness = DetUnit->surface().bounds().thickness();
202 TrapezoidalParam = 1;
203 }
204
205 Surface::GlobalPoint WidthVector = plane.toGlobal( LocalPoint(width/2, 0, 0) );
206 Surface::GlobalPoint LengthVector = plane.toGlobal( LocalPoint(0, length/2, 0) );
207 Surface::GlobalPoint ThickVector = plane.toGlobal( LocalPoint(0, 0, thickness/2) );
208
209 GlobalVector Pos = GlobalVector(DetUnit->position().basicVector());
210
211 Geom_Tracker->Add_TrackerDet(Detid.rawId(), TrapezoidalParam,
212 Pos.x(), Pos.y(), Pos.z(),
213 WidthVector.x() -Pos.x(), WidthVector.y() -Pos.y(), WidthVector.z() -Pos.z(),
214 LengthVector.x()-Pos.x(), LengthVector.y()-Pos.y(), LengthVector.z()-Pos.z(),
215 ThickVector.x() -Pos.x(), ThickVector.y() -Pos.y(), ThickVector.z() -Pos.z());
216
217
218 }
219
220
221 // ### CALO GEOMETRY ###
222
223 edm::ESHandle<CaloGeometry> CaloGeom;
224 iSetup.get<IdealGeometryRecord>().get( CaloGeom );
225 const vector<DetId> CaloDets = CaloGeom->getValidDetIds();
226
227 for(unsigned int i=0;i<CaloDets.size();i++)
228 {
229 Detid = CaloDets[i];
230 SubDet = Detid.subdetId();
231
232 Geometry* Geom_temp = NULL;
233 if(Detid.det()==DetId::Ecal){
234 Geom_temp = Geom_ECAL;
235
236 if(SubDet<1 || SubDet>3){
237 printf("Don't Save EcalTriggerTower or EcalLaserPnDiode\n");
238 continue;
239 }
240 }else if(Detid.det()==DetId::Hcal){
241 Geom_temp = Geom_HCAL;
242
243 if(SubDet<1 || SubDet>4){
244 printf("Don't Save HcalEmpty, HcalTriggerTower, HcalOther\n");
245 continue;
246 }
247 }else{
248 continue;
249 }
250
251 const CaloCellGeometry* CellGeom = CaloGeom->getGeometry(Detid);
252 GlobalPoint CellPos = CellGeom->getPosition();
253 const CaloCellGeometry::CornersVec CellCorners = CellGeom->getCorners();
254
255 Geom_temp->Add_CaloDet(Detid.rawId(),
256 CellPos.x() , CellPos.y() , CellPos.z(),
257 CellCorners[0].x() , CellCorners[0].y() , CellCorners[0].z() ,
258 CellCorners[1].x() , CellCorners[1].y() , CellCorners[1].z() ,
259 CellCorners[2].x() , CellCorners[2].y() , CellCorners[2].z() ,
260 CellCorners[3].x() , CellCorners[3].y() , CellCorners[3].z() ,
261 CellCorners[4].x() , CellCorners[4].y() , CellCorners[4].z() ,
262 CellCorners[5].x() , CellCorners[5].y() , CellCorners[5].z() ,
263 CellCorners[6].x() , CellCorners[6].y() , CellCorners[6].z() ,
264 CellCorners[7].x() , CellCorners[7].y() , CellCorners[7].z() );
265 }
266
267 // ### MUON GEOMETRY ###
268
269 edm::ESHandle<DTGeometry> DtGeom;
270 // iSetup.get<IdealGeometryRecord>().get( DtGeom );
271 iSetup.get<MuonGeometryRecord>().get( DtGeom );
272 const vector<GeomDet*> DtDets = DtGeom->dets();
273
274
275 for(unsigned int i=0;i<DtDets.size();i++)
276 {
277 Detid = DetId(DtDets[i]->geographicalId());
278 SubDet = Detid.subdetId();
279
280 GeomDet* DetUnit = DtDets[i];
281 if(!DetUnit)continue;
282 const BoundPlane plane = DetUnit->surface();
283 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
284 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
285
286 float width = 0;
287 float length = 0;
288 float thickness = 0;
289 float TrapezoidalParam = 0;
290
291 if(trapezoidalBounds)
292 {
293 std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
294 width = parameters[0]*2;
295 length = parameters[3]*2;
296 thickness = (*trapezoidalBounds).thickness();
297 TrapezoidalParam = parameters[1]/parameters[0];
298 }else if(rectangularBounds){
299 width = DetUnit->surface().bounds().width();
300 length = DetUnit->surface().bounds().length();
301 thickness = DetUnit->surface().bounds().thickness();
302 TrapezoidalParam = 1;
303 }
304
305 Surface::GlobalPoint WidthVector = plane.toGlobal( LocalPoint(width/2, 0, 0) );
306 Surface::GlobalPoint LengthVector = plane.toGlobal( LocalPoint(0, length/2, 0) );
307 Surface::GlobalPoint ThickVector = plane.toGlobal( LocalPoint(0, 0, thickness/2) );
308
309 GlobalVector Pos = GlobalVector(DetUnit->position().basicVector());
310
311 printf("Det = %i\n",Detid.det());
312 Geom_Muon->Add_TrackerDet(Detid.rawId(), TrapezoidalParam,
313 Pos.x(), Pos.y(), Pos.z(),
314 WidthVector.x() -Pos.x(), WidthVector.y() -Pos.y(), WidthVector.z() -Pos.z(),
315 LengthVector.x()-Pos.x(), LengthVector.y()-Pos.y(), LengthVector.z()-Pos.z(),
316 ThickVector.x() -Pos.x(), ThickVector.y() -Pos.y(), ThickVector.z() -Pos.z());
317 }
318
319
320 // ### Save .geom ###
321
322 Geom_Tracker->Save("Tracker.geom");
323 Geom_ECAL ->Save("Ecal.geom");
324 Geom_HCAL ->Save("Hcal.geom");
325 Geom_Muon ->Save("Muon.geom");
326
327 }
328
329 // ------------ method called once each job just after ending the event loop ------------
330 void
331 OpenGLDisplayer::endJob() {
332
333 MyEvents->Save((char*) OutputFile.c_str());
334 // MyEvents->Load((char*) OutputFile.c_str());
335 }
336
337
338
339 //
340 // member functions
341 //
342
343 // ------------ method called to for each event ------------
344 void
345 OpenGLDisplayer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
346 {
347 MySimEvent* MyEvent = new MySimEvent;
348
349
350
351 // access the tracker
352 edm::ESHandle<TrackerGeometry> theTrackerGeometry;
353 iSetup.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
354 const TrackerGeometry& theTracker(*theTrackerGeometry);
355
356 // the DetUnits
357 TrackingGeometry::DetContainer theDetUnits = theTracker.dets();
358
359
360 // Save Data
361 for(unsigned int i=0;i<SimTrackProducers.size();i++){
362 edm::Handle<std::vector< SimTrack > > h_SimTracks;
363 iEvent.getByLabel(SimTrackProducers[i], h_SimTracks);
364 std::vector< SimTrack > SimTrackColl = *h_SimTracks.product();
365
366 for ( unsigned int a = 0; a < SimTrackColl.size(); ++a ) {
367 MySimTrack MysimTrack;
368 SimTrack simTrack =SimTrackColl[a];
369
370 MysimTrack.track_id =simTrack.trackId();
371 MysimTrack.Type =simTrack.type();
372 MysimTrack.parent_vertex=simTrack.vertIndex();
373 MysimTrack.Px =simTrack.momentum().x();
374 MysimTrack.Py =simTrack.momentum().y();
375 MysimTrack.Pz =simTrack.momentum().z();
376 MysimTrack.E =simTrack.momentum().e();
377 MysimTrack.charge =simTrack.charge();
378
379 MyEvent->MySimTrackCollection.push_back(MysimTrack);
380 }
381 }
382
383
384 for(unsigned int i=0;i<SimVertexProducers.size();i++){
385 edm::Handle<std::vector< SimVertex > > h_Vertex;
386 iEvent.getByLabel(SimVertexProducers[i], h_Vertex);
387 std::vector< SimVertex > VertexColl = *h_Vertex.product();
388
389 for (unsigned int b = 0; b < VertexColl.size(); ++b ) {
390 MySimVertex MyVertex;
391 SimVertex Vertex = VertexColl[b];
392
393 MyVertex.parentTrack_id =Vertex.parentIndex ();
394 MyVertex.x =Vertex.position().x();
395 MyVertex.y =Vertex.position().y();
396 MyVertex.z =Vertex.position().z();
397
398 MyEvent->MySimVertexCollection.push_back(MyVertex);
399 }
400 }
401
402 for(unsigned int i=0;i<SimHitProducers.size();i++){
403 edm::Handle<std::vector< PSimHit > > h_Hits;
404 iEvent.getByLabel(SimHitProducers[i], h_Hits);
405 std::vector< PSimHit > Hits = *h_Hits.product();
406
407 for(unsigned int h=0; h<Hits.size(); h++)
408 {
409 DetId theDetUnitId(Hits[h].detUnitId());
410 const GeomDet * theDet = theTracker.idToDet(theDetUnitId);
411
412 MyPSimHit Hit;
413 Hit.x = theDet->surface().toGlobal(Hits[h].localPosition()).x();
414 Hit.y = theDet->surface().toGlobal(Hits[h].localPosition()).y();
415 Hit.z = theDet->surface().toGlobal(Hits[h].localPosition()).z();
416 Hit.ProcessType = Hits[h].processType();
417 Hit.dEdX = Hits[h].energyLoss();
418
419 MyEvent->MyPSimHitCollection.push_back(Hit);
420 }
421 }
422
423
424 for(unsigned int i=0;i<TrackProducers.size();i++){
425 edm::Handle<std::vector< reco::Track > > h_Tracks;
426 iEvent.getByLabel(TrackProducers[i], h_Tracks);
427 std::vector< reco::Track > TrackColl = *h_Tracks.product();
428
429 for ( unsigned int t = 0; t < TrackColl.size(); ++t ) {
430 MyRecoTrack MyrecoTrack;
431 reco::Track recoTrack =TrackColl[t];
432
433 for(unsigned int h=0;h<recoTrack.recHitsSize();h++){
434 TrackingRecHitRef h_it = recoTrack.recHit(h);
435 if(!h_it->isValid() )continue;
436 DetId detId = h_it->geographicalId();
437 const GeomDet * theDet = theTracker.idToDet(detId);
438 LocalPoint localPos = h_it->localPosition();
439
440 MyRecoHit hit;
441 hit.x = theDet->surface().toGlobal(localPos).x();
442 hit.y = theDet->surface().toGlobal(localPos).y();
443 hit.z = theDet->surface().toGlobal(localPos).z();
444 hit.DetId = detId.rawId();
445 hit.Charge = -1;
446 MyrecoTrack.Hits.push_back(hit);
447 }
448 MyEvent->MyRecoTrackCollection.push_back(MyrecoTrack);
449 }
450 }
451
452
453 for(unsigned int i=0;i<EcalRecHitProducers.size();i++){
454 edm::Handle<EcalRecHitCollection > h_Ecal_RecHits;
455 iEvent.getByLabel(EcalRecHitProducers[i], h_Ecal_RecHits);
456 EcalRecHitCollection Ecal_RecHits = *h_Ecal_RecHits.product();
457
458 for(unsigned int eh=0;eh<Ecal_RecHits.size();eh++){
459 MyCaloHit temp_EcalHit;
460 temp_EcalHit.E = Ecal_RecHits[eh].energy();
461 temp_EcalHit.t = Ecal_RecHits[eh].time();
462 temp_EcalHit.DetId = (Ecal_RecHits[eh].detid()).rawId();
463
464 MyEvent->MyEcalCaloHitCollection.push_back(temp_EcalHit);
465 }
466 }
467
468
469 for(unsigned int i=0;i<HcalHBHERecHitProducers.size();i++){
470 edm::Handle<HBHERecHitCollection > h_HcalHBHE_RecHits;
471 iEvent.getByLabel(HcalHBHERecHitProducers[i], h_HcalHBHE_RecHits);
472 HBHERecHitCollection HcalHBHE_RecHits = *h_HcalHBHE_RecHits.product();
473
474 for(unsigned int hh=0;hh<HcalHBHE_RecHits.size();hh++){
475 MyCaloHit temp_HcalHBHEHit;
476 temp_HcalHBHEHit.E = HcalHBHE_RecHits[hh].energy();
477 temp_HcalHBHEHit.t = HcalHBHE_RecHits[hh].time();
478 temp_HcalHBHEHit.DetId = (HcalHBHE_RecHits[hh].detid()).rawId();
479
480 MyEvent->MyHcalCaloHitCollection.push_back(temp_HcalHBHEHit);
481 }
482 }
483
484
485 for(unsigned int i=0;i<HcalHORecHitProducers.size();i++){
486 edm::Handle<HORecHitCollection > h_HcalHO_RecHits;
487 iEvent.getByLabel(HcalHORecHitProducers[i], h_HcalHO_RecHits);
488 HORecHitCollection HcalHO_RecHits = *h_HcalHO_RecHits.product();
489
490 for(unsigned int hh=0;hh<HcalHO_RecHits.size();hh++){
491 MyCaloHit temp_HcalHOHit;
492 temp_HcalHOHit.E = HcalHO_RecHits[hh].energy();
493 temp_HcalHOHit.t = HcalHO_RecHits[hh].time();
494 temp_HcalHOHit.DetId = (HcalHO_RecHits[hh].detid()).rawId();
495
496 MyEvent->MyHcalCaloHitCollection.push_back(temp_HcalHOHit);
497 }
498 }
499
500
501 for(unsigned int i=0;i<HcalHFRecHitProducers.size();i++){
502 edm::Handle<HFRecHitCollection > h_HcalHF_RecHits;
503 iEvent.getByLabel(HcalHFRecHitProducers[i], h_HcalHF_RecHits);
504 HFRecHitCollection HcalHF_RecHits = *h_HcalHF_RecHits.product();
505
506 for(unsigned int hh=0;hh<HcalHF_RecHits.size();hh++){
507 MyCaloHit temp_HcalHFHit;
508 temp_HcalHFHit.E = HcalHF_RecHits[hh].energy();
509 temp_HcalHFHit.t = HcalHF_RecHits[hh].time();
510 temp_HcalHFHit.DetId = (HcalHF_RecHits[hh].detid()).rawId();
511
512 MyEvent->MyHcalCaloHitCollection.push_back(temp_HcalHFHit);
513 }
514 }
515
516 MyEvents->Events.push_back(MyEvent);
517 }
518
519
520 //define this as a plug-in
521 DEFINE_FWK_MODULE(OpenGLDisplayer);
522
523