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

# User Rev Content
1 querten 1.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 roberfro 1.14 // $Id: OpenGLDisplayer.cc,v 1.13 2008/03/25 13:12:31 querten Exp $
17 querten 1.1 //
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 querten 1.9 #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 querten 1.12 #include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
58     #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
59     #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
60 querten 1.9
61 querten 1.1 #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 querten 1.7 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
78     #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
79 querten 1.1
80     #include "Visualisation/OpenGLDisplayer/interface/SimEvent.h"
81 querten 1.3 #include "Visualisation/OpenGLDisplayer/interface/Geometry.h"
82 querten 1.1
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 querten 1.3 virtual void beginJob(const edm::EventSetup& iSetup);
98     virtual void analyze (const edm::Event&, const edm::EventSetup&);
99     virtual void endJob ();
100 querten 1.1
101    
102     MySimEvents* MyEvents;
103 querten 1.11 Geometry* Geom_Tracker;
104     Geometry* Geom_ECAL;
105     Geometry* Geom_HCAL;
106 querten 1.1
107    
108     std::vector<std::string> SimHitSubdetectors;
109    
110     std::string OutputFile;
111 roberfro 1.14 std::string trackProducer;
112 querten 1.1
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 roberfro 1.14
126     trackProducer = iConfig.getParameter<std::string >("TrackProducer");
127 querten 1.1 }
128    
129    
130     OpenGLDisplayer::~OpenGLDisplayer()
131     {
132     }
133    
134     // ------------ method called once each job just before starting event loop ------------
135     void
136 querten 1.3 OpenGLDisplayer::beginJob(const edm::EventSetup& iSetup)
137 querten 1.1 {
138 querten 1.9 DetId Detid;
139     int SubDet;
140 querten 1.1
141 querten 1.11 MyEvents = new MySimEvents();
142     Geom_Tracker = new Geometry();
143     Geom_ECAL = new Geometry();
144     Geom_HCAL = new Geometry();
145    
146 querten 1.3
147     edm::ESHandle<TrackerGeometry> tkGeom;
148     iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
149 querten 1.9 vector<GeomDet*> TrackerDets = tkGeom->dets();
150 querten 1.3
151 querten 1.9 for(unsigned int i=0;i<TrackerDets.size();i++){
152     Detid = TrackerDets[i]->geographicalId();
153     // SubDet = Detid.subdetId();
154 querten 1.4
155 querten 1.9 GeomDet* DetUnit = TrackerDets[i];
156 querten 1.3 if(!DetUnit)continue;
157 querten 1.7 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 querten 1.3
161 querten 1.9 float width = 0;
162     float length = 0;
163     float thickness = 0;
164     float TrapezoidalParam = 0;
165 querten 1.7
166     if(trapezoidalBounds)
167     {
168     std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
169 querten 1.8 width = parameters[0]*2;
170     length = parameters[3]*2;
171 querten 1.7 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 querten 1.4
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 querten 1.7 GlobalVector Pos = GlobalVector(DetUnit->position().basicVector());
185 querten 1.4
186 querten 1.11 Geom_Tracker->Add_TrackerDet(Detid.rawId(), TrapezoidalParam,
187 querten 1.7 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 querten 1.3
193     }
194 querten 1.9
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 querten 1.11 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 querten 1.9
226     const CaloCellGeometry* CellGeom = CaloGeom->getGeometry(Detid);
227     GlobalPoint CellPos = CellGeom->getPosition();
228 querten 1.11 const CaloCellGeometry::CornersVec CellCorners = CellGeom->getCorners();
229 querten 1.9
230 querten 1.11 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 querten 1.9 }
241    
242 querten 1.11 Geom_Tracker->Save("Tracker.geom");
243     Geom_ECAL->Save("ECAL.geom");
244     Geom_HCAL->Save("HCAL.geom");
245 querten 1.9
246 querten 1.1
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 querten 1.12 // MyEvents->Load((char*) OutputFile.c_str());
256 querten 1.1 }
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 roberfro 1.14 // access primary tracks
347 querten 1.1 edm::Handle<std::vector< reco::Track > > h_Tracks;
348 roberfro 1.14 iEvent.getByLabel(trackProducer, h_Tracks);
349 querten 1.1 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 querten 1.2 // MyrecoTrack.Hits = new MyRecoHit[MyrecoTrack.N];
356 querten 1.1
357 querten 1.2 for(unsigned int h=0;h<recoTrack.recHitsSize();h++){
358 querten 1.1 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 querten 1.2
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 querten 1.1 }
375     MyEvent->MyRecoTrackCollection.push_back(MyrecoTrack);
376     }
377    
378 querten 1.12
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 querten 1.1 MyEvents->Events.push_back(MyEvent);
474    
475     }
476    
477    
478     //define this as a plug-in
479     DEFINE_FWK_MODULE(OpenGLDisplayer);
480    
481