ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FastOpenGlDisplayer/src/OpenGLDisplayer.cc
Revision: 1.14
Committed: Thu Apr 3 12:18:29 2008 UTC (17 years, 1 month ago) by roberfro
Content type: text/plain
Branch: MAIN
Changes since 1.13: +6 -2 lines
Log Message:
Possibility to merge track and add nuclear interaction sec tracks

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.13 2008/03/25 13:12:31 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 "DataFormats/CaloRecHit/interface/CaloRecHit.h"
58 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
59 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
60
61 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
62 #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
63 #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
64 #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
65 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
66 #include "SimDataFormats/Track/interface/SimTrack.h"
67 #include "SimDataFormats/Vertex/interface/SimVertex.h"
68
69 #include "DataFormats/Provenance/interface/BranchDescription.h"
70 #include "DataFormats/Provenance/interface/Provenance.h"
71 #include "DataFormats/Candidate/interface/Candidate.h"
72 #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
73 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
74 #include "DataFormats/TrackReco/interface/Track.h"
75
76 #include "DataFormats/Math/interface/Point3D.h"
77 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
78 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
79
80 #include "Visualisation/OpenGLDisplayer/interface/SimEvent.h"
81 #include "Visualisation/OpenGLDisplayer/interface/Geometry.h"
82
83 using namespace edm;
84 using namespace std;
85
86 //
87 // class decleration
88 //
89
90 class OpenGLDisplayer : public edm::EDAnalyzer {
91 public:
92 explicit OpenGLDisplayer(const edm::ParameterSet&);
93 ~OpenGLDisplayer();
94
95
96 private:
97 virtual void beginJob(const edm::EventSetup& iSetup);
98 virtual void analyze (const edm::Event&, const edm::EventSetup&);
99 virtual void endJob ();
100
101
102 MySimEvents* MyEvents;
103 Geometry* Geom_Tracker;
104 Geometry* Geom_ECAL;
105 Geometry* Geom_HCAL;
106
107
108 std::vector<std::string> SimHitSubdetectors;
109
110 std::string OutputFile;
111 std::string trackProducer;
112
113
114 // ----------member data ---------------------------
115 };
116
117 //
118 // constructors and destructor
119 //
120 OpenGLDisplayer::OpenGLDisplayer(const edm::ParameterSet& iConfig)
121 {
122 SimHitSubdetectors = iConfig.getParameter<std::vector<std::string> >("SimHitSubdetectors");
123
124 OutputFile = iConfig.getParameter<std::string >("OutputFile");
125
126 trackProducer = iConfig.getParameter<std::string >("TrackProducer");
127 }
128
129
130 OpenGLDisplayer::~OpenGLDisplayer()
131 {
132 }
133
134 // ------------ method called once each job just before starting event loop ------------
135 void
136 OpenGLDisplayer::beginJob(const edm::EventSetup& iSetup)
137 {
138 DetId Detid;
139 int SubDet;
140
141 MyEvents = new MySimEvents();
142 Geom_Tracker = new Geometry();
143 Geom_ECAL = new Geometry();
144 Geom_HCAL = new Geometry();
145
146
147 edm::ESHandle<TrackerGeometry> tkGeom;
148 iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
149 vector<GeomDet*> TrackerDets = tkGeom->dets();
150
151 for(unsigned int i=0;i<TrackerDets.size();i++){
152 Detid = TrackerDets[i]->geographicalId();
153 // SubDet = Detid.subdetId();
154
155 GeomDet* DetUnit = TrackerDets[i];
156 if(!DetUnit)continue;
157 const BoundPlane plane = DetUnit->surface();
158 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
159 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
160
161 float width = 0;
162 float length = 0;
163 float thickness = 0;
164 float TrapezoidalParam = 0;
165
166 if(trapezoidalBounds)
167 {
168 std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
169 width = parameters[0]*2;
170 length = parameters[3]*2;
171 thickness = (*trapezoidalBounds).thickness();
172 TrapezoidalParam = parameters[1]/parameters[0];
173 }else if(rectangularBounds){
174 width = DetUnit->surface().bounds().width();
175 length = DetUnit->surface().bounds().length();
176 thickness = DetUnit->surface().bounds().thickness();
177 TrapezoidalParam = 1;
178 }
179
180 Surface::GlobalPoint WidthVector = plane.toGlobal( LocalPoint(width/2, 0, 0) );
181 Surface::GlobalPoint LengthVector = plane.toGlobal( LocalPoint(0, length/2, 0) );
182 Surface::GlobalPoint ThickVector = plane.toGlobal( LocalPoint(0, 0, thickness/2) );
183
184 GlobalVector Pos = GlobalVector(DetUnit->position().basicVector());
185
186 Geom_Tracker->Add_TrackerDet(Detid.rawId(), TrapezoidalParam,
187 Pos.x(), Pos.y(), Pos.z(),
188 WidthVector.x() -Pos.x(), WidthVector.y() -Pos.y(), WidthVector.z() -Pos.z(),
189 LengthVector.x()-Pos.x(), LengthVector.y()-Pos.y(), LengthVector.z()-Pos.z(),
190 ThickVector.x() -Pos.x(), ThickVector.y() -Pos.y(), ThickVector.z() -Pos.z());
191
192
193 }
194
195
196 // ### CALO GEOMETRY ###
197
198 edm::ESHandle<CaloGeometry> CaloGeom;
199 iSetup.get<IdealGeometryRecord>().get( CaloGeom );
200 const vector<DetId> CaloDets = CaloGeom->getValidDetIds();
201
202 for(unsigned int i=0;i<CaloDets.size();i++)
203 {
204 Detid = CaloDets[i];
205 SubDet = Detid.subdetId();
206
207 Geometry* Geom_temp = NULL;
208 if(Detid.det()==DetId::Ecal){
209 Geom_temp = Geom_ECAL;
210
211 if(SubDet<1 || SubDet>3){
212 printf("Don't Save EcalTriggerTower or EcalLaserPnDiode\n");
213 continue;
214 }
215 }else if(Detid.det()==DetId::Hcal){
216 Geom_temp = Geom_HCAL;
217
218 if(SubDet<1 || SubDet>4){
219 printf("Don't Save HcalEmpty, HcalTriggerTower, HcalOther\n");
220 continue;
221 }
222 }else{
223 continue;
224 }
225
226 const CaloCellGeometry* CellGeom = CaloGeom->getGeometry(Detid);
227 GlobalPoint CellPos = CellGeom->getPosition();
228 const CaloCellGeometry::CornersVec CellCorners = CellGeom->getCorners();
229
230 Geom_temp->Add_CaloDet(Detid.rawId(),
231 CellPos.x() , CellPos.y() , CellPos.z(),
232 CellCorners[0].x() , CellCorners[0].y() , CellCorners[0].z() ,
233 CellCorners[1].x() , CellCorners[1].y() , CellCorners[1].z() ,
234 CellCorners[2].x() , CellCorners[2].y() , CellCorners[2].z() ,
235 CellCorners[3].x() , CellCorners[3].y() , CellCorners[3].z() ,
236 CellCorners[4].x() , CellCorners[4].y() , CellCorners[4].z() ,
237 CellCorners[5].x() , CellCorners[5].y() , CellCorners[5].z() ,
238 CellCorners[6].x() , CellCorners[6].y() , CellCorners[6].z() ,
239 CellCorners[7].x() , CellCorners[7].y() , CellCorners[7].z() );
240 }
241
242 Geom_Tracker->Save("Tracker.geom");
243 Geom_ECAL->Save("ECAL.geom");
244 Geom_HCAL->Save("HCAL.geom");
245
246
247
248 }
249
250 // ------------ method called once each job just after ending the event loop ------------
251 void
252 OpenGLDisplayer::endJob() {
253
254 MyEvents->Save((char*) OutputFile.c_str());
255 // MyEvents->Load((char*) OutputFile.c_str());
256 }
257
258
259
260 //
261 // member functions
262 //
263
264 // ------------ method called to for each event ------------
265 void
266 OpenGLDisplayer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
267 {
268
269 edm::Handle<std::vector< SimTrack > > h_SimTracks;
270 iEvent.getByLabel("g4SimHits", h_SimTracks);
271 std::vector< SimTrack > SimTrackColl = *h_SimTracks.product();
272
273 edm::Handle<std::vector< SimVertex > > h_Vertex;
274 iEvent.getByLabel("g4SimHits", h_Vertex);
275 std::vector< SimVertex > VertexColl = *h_Vertex.product();
276
277
278
279 // access the tracker
280 edm::ESHandle<TrackerGeometry> theTrackerGeometry;
281 iSetup.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
282 const TrackerGeometry& theTracker(*theTrackerGeometry);
283
284 // the DetUnits
285 TrackingGeometry::DetContainer theDetUnits = theTracker.dets();
286
287
288
289
290 // Save the data
291
292 MySimEvent* MyEvent = new MySimEvent;
293
294 for ( unsigned int a = 0; a < SimTrackColl.size(); ++a ) {
295 MySimTrack MysimTrack;
296 SimTrack simTrack =SimTrackColl[a];
297
298 MysimTrack.track_id =simTrack.trackId();
299 MysimTrack.Type =simTrack.type();
300 MysimTrack.parent_vertex=simTrack.vertIndex();
301 MysimTrack.Px =simTrack.momentum().x();
302 MysimTrack.Py =simTrack.momentum().y();
303 MysimTrack.Pz =simTrack.momentum().z();
304 MysimTrack.E =simTrack.momentum().e();
305 MysimTrack.charge =simTrack.charge();
306
307 MyEvent->MySimTrackCollection.push_back(MysimTrack);
308 }
309
310
311 for (unsigned int b = 0; b < VertexColl.size(); ++b ) {
312 MySimVertex MyVertex;
313 SimVertex Vertex = VertexColl[b];
314
315 MyVertex.parentTrack_id =Vertex.parentIndex ();
316 MyVertex.x =Vertex.position().x();
317 MyVertex.y =Vertex.position().y();
318 MyVertex.z =Vertex.position().z();
319
320 MyEvent->MySimVertexCollection.push_back(MyVertex);
321 }
322
323 for(unsigned int i=0; i<SimHitSubdetectors.size(); i++)
324 {
325 edm::Handle<std::vector< PSimHit > > h_Hits;
326 iEvent.getByLabel("g4SimHits",SimHitSubdetectors[i].c_str(), h_Hits);
327 std::vector< PSimHit > Hits = *h_Hits.product();
328
329 for(unsigned int h=0; h<Hits.size(); h++)
330 {
331 DetId theDetUnitId(Hits[h].detUnitId());
332 const GeomDet * theDet = theTracker.idToDet(theDetUnitId);
333
334 MyPSimHit Hit;
335 Hit.x = theDet->surface().toGlobal(Hits[h].localPosition()).x();
336 Hit.y = theDet->surface().toGlobal(Hits[h].localPosition()).y();
337 Hit.z = theDet->surface().toGlobal(Hits[h].localPosition()).z();
338 Hit.ProcessType = Hits[h].processType();
339 Hit.dEdX = Hits[h].energyLoss();
340
341 MyEvent->MyPSimHitCollection.push_back(Hit);
342
343 }
344 }
345
346 // access primary tracks
347 edm::Handle<std::vector< reco::Track > > h_Tracks;
348 iEvent.getByLabel(trackProducer, h_Tracks);
349 std::vector< reco::Track > TrackColl = *h_Tracks.product();
350
351 for ( unsigned int t = 0; t < TrackColl.size(); ++t ) {
352 MyRecoTrack MyrecoTrack;
353 reco::Track recoTrack =TrackColl[t];
354
355 // MyrecoTrack.Hits = new MyRecoHit[MyrecoTrack.N];
356
357 for(unsigned int h=0;h<recoTrack.recHitsSize();h++){
358 TrackingRecHitRef h_it = recoTrack.recHit(h);
359 if(!h_it->isValid() )continue;
360 DetId detId = h_it->geographicalId();
361 const GeomDet * theDet = theTracker.idToDet(detId);
362 LocalPoint localPos = h_it->localPosition();
363
364 MyRecoHit hit;
365 hit.x = theDet->surface().toGlobal(localPos).x();
366 hit.y = theDet->surface().toGlobal(localPos).y();
367 hit.z = theDet->surface().toGlobal(localPos).z();
368 hit.DetId = detId.rawId();
369 hit.Charge = -1;
370 MyrecoTrack.Hits.push_back(hit);
371
372
373 // printf("%8.2f %8.2f %8.2f\n",theDet->surface().toGlobal(localPos).x(),theDet->surface().toGlobal(localPos).y(),theDet->surface().toGlobal(localPos).z());
374 }
375 MyEvent->MyRecoTrackCollection.push_back(MyrecoTrack);
376 }
377
378
379
380 edm::Handle<EcalRecHitCollection > h_EcalEB_RecHits;
381 iEvent.getByLabel("ecalRecHit","EcalRecHitsEB", h_EcalEB_RecHits);
382 EcalRecHitCollection EcalEB_RecHits = *h_EcalEB_RecHits.product();
383
384 for(unsigned int eh=0;eh<EcalEB_RecHits.size();eh++){
385 MyCaloHit temp_EcalEBHit;
386 temp_EcalEBHit.E = EcalEB_RecHits[eh].energy();
387 temp_EcalEBHit.t = EcalEB_RecHits[eh].time();
388 temp_EcalEBHit.DetId = (EcalEB_RecHits[eh].detid()).rawId();
389
390 MyEvent->MyEcalCaloHitCollection.push_back(temp_EcalEBHit);
391 }
392
393
394 edm::Handle<EcalRecHitCollection > h_EcalEE_RecHits;
395 iEvent.getByLabel("ecalRecHit","EcalRecHitsEE", h_EcalEE_RecHits);
396 EcalRecHitCollection EcalEE_RecHits = *h_EcalEE_RecHits.product();
397
398 for(unsigned int eh=0;eh<EcalEE_RecHits.size();eh++){
399 MyCaloHit temp_EcalEEHit;
400 temp_EcalEEHit.E = EcalEE_RecHits[eh].energy();
401 temp_EcalEEHit.t = EcalEE_RecHits[eh].time();
402 temp_EcalEEHit.DetId = (EcalEE_RecHits[eh].detid()).rawId();
403
404 MyEvent->MyEcalCaloHitCollection.push_back(temp_EcalEEHit);
405 }
406
407
408 edm::Handle<EcalRecHitCollection > h_EcalES_RecHits;
409 iEvent.getByLabel("ecalPreshowerRecHit","EcalRecHitsES", h_EcalES_RecHits);
410 EcalRecHitCollection EcalES_RecHits = *h_EcalES_RecHits.product();
411
412 for(unsigned int eh=0;eh<EcalES_RecHits.size();eh++){
413 MyCaloHit temp_EcalESHit;
414 temp_EcalESHit.E = EcalES_RecHits[eh].energy();
415 temp_EcalESHit.t = EcalES_RecHits[eh].time();
416 temp_EcalESHit.DetId = (EcalES_RecHits[eh].detid()).rawId();
417
418 MyEvent->MyEcalCaloHitCollection.push_back(temp_EcalESHit);
419 }
420
421
422
423 edm::Handle<HBHERecHitCollection > h_HcalHBHE_RecHits;
424 iEvent.getByLabel("hbhereco", h_HcalHBHE_RecHits);
425 HBHERecHitCollection HcalHBHE_RecHits = *h_HcalHBHE_RecHits.product();
426
427 for(unsigned int hh=0;hh<HcalHBHE_RecHits.size();hh++){
428 MyCaloHit temp_HcalHBHEHit;
429 temp_HcalHBHEHit.E = HcalHBHE_RecHits[hh].energy();
430 temp_HcalHBHEHit.t = HcalHBHE_RecHits[hh].time();
431 temp_HcalHBHEHit.DetId = (HcalHBHE_RecHits[hh].detid()).rawId();
432
433 MyEvent->MyHcalCaloHitCollection.push_back(temp_HcalHBHEHit);
434 }
435
436
437 edm::Handle<HORecHitCollection > h_HcalHO_RecHits;
438 iEvent.getByLabel("horeco", h_HcalHO_RecHits);
439 HORecHitCollection HcalHO_RecHits = *h_HcalHO_RecHits.product();
440
441 for(unsigned int hh=0;hh<HcalHO_RecHits.size();hh++){
442 MyCaloHit temp_HcalHOHit;
443 temp_HcalHOHit.E = HcalHO_RecHits[hh].energy();
444 temp_HcalHOHit.t = HcalHO_RecHits[hh].time();
445 temp_HcalHOHit.DetId = (HcalHO_RecHits[hh].detid()).rawId();
446
447 MyEvent->MyHcalCaloHitCollection.push_back(temp_HcalHOHit);
448 }
449
450 edm::Handle<HFRecHitCollection > h_HcalHF_RecHits;
451 iEvent.getByLabel("hfreco", h_HcalHF_RecHits);
452 HFRecHitCollection HcalHF_RecHits = *h_HcalHF_RecHits.product();
453
454 for(unsigned int hh=0;hh<HcalHF_RecHits.size();hh++){
455 MyCaloHit temp_HcalHFHit;
456 temp_HcalHFHit.E = HcalHF_RecHits[hh].energy();
457 temp_HcalHFHit.t = HcalHF_RecHits[hh].time();
458 temp_HcalHFHit.DetId = (HcalHF_RecHits[hh].detid()).rawId();
459
460 MyEvent->MyHcalCaloHitCollection.push_back(temp_HcalHFHit);
461 }
462
463
464
465
466
467
468
469
470
471
472
473 MyEvents->Events.push_back(MyEvent);
474
475 }
476
477
478 //define this as a plug-in
479 DEFINE_FWK_MODULE(OpenGLDisplayer);
480
481