ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/FastOpenGlDisplayer/src/OpenGLDisplayer.cc
Revision: 1.10
Committed: Fri Mar 21 19:40:46 2008 UTC (17 years, 1 month ago) by querten
Content type: text/plain
Branch: MAIN
Changes since 1.9: +3 -21 lines
Log Message:
Include Geometry of Calorimter

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 querten 1.3 // $Id: OpenGLDisplayer.cc,v 1.2 2008/03/15 17:02:20 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    
58    
59 querten 1.1 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
60     #include "SimDataFormats/TrackingHit/interface/PSimHit.h"
61     #include "SimDataFormats/CrossingFrame/interface/CrossingFrame.h"
62     #include "SimDataFormats/CrossingFrame/interface/MixCollection.h"
63     #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
64     #include "SimDataFormats/Track/interface/SimTrack.h"
65     #include "SimDataFormats/Vertex/interface/SimVertex.h"
66    
67     #include "DataFormats/Provenance/interface/BranchDescription.h"
68     #include "DataFormats/Provenance/interface/Provenance.h"
69     #include "DataFormats/Candidate/interface/Candidate.h"
70     #include "DataFormats/HepMCCandidate/interface/GenParticleCandidate.h"
71     #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
72     #include "DataFormats/TrackReco/interface/Track.h"
73    
74     #include "DataFormats/Math/interface/Point3D.h"
75 querten 1.7 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
76     #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
77 querten 1.1
78     #include "Visualisation/OpenGLDisplayer/interface/SimEvent.h"
79 querten 1.3 #include "Visualisation/OpenGLDisplayer/interface/Geometry.h"
80 querten 1.1
81     using namespace edm;
82     using namespace std;
83    
84     //
85     // class decleration
86     //
87    
88     class OpenGLDisplayer : public edm::EDAnalyzer {
89     public:
90     explicit OpenGLDisplayer(const edm::ParameterSet&);
91     ~OpenGLDisplayer();
92    
93    
94     private:
95 querten 1.3 virtual void beginJob(const edm::EventSetup& iSetup);
96     virtual void analyze (const edm::Event&, const edm::EventSetup&);
97     virtual void endJob ();
98 querten 1.1
99    
100     MySimEvents* MyEvents;
101 querten 1.3 Geometry* Geom;
102 querten 1.1
103    
104     std::vector<std::string> SimHitSubdetectors;
105    
106     std::string OutputFile;
107    
108    
109     // ----------member data ---------------------------
110     };
111    
112     //
113     // constructors and destructor
114     //
115     OpenGLDisplayer::OpenGLDisplayer(const edm::ParameterSet& iConfig)
116     {
117     SimHitSubdetectors = iConfig.getParameter<std::vector<std::string> >("SimHitSubdetectors");
118    
119     OutputFile = iConfig.getParameter<std::string >("OutputFile");
120     }
121    
122    
123     OpenGLDisplayer::~OpenGLDisplayer()
124     {
125     }
126    
127     // ------------ method called once each job just before starting event loop ------------
128     void
129 querten 1.3 OpenGLDisplayer::beginJob(const edm::EventSetup& iSetup)
130 querten 1.1 {
131 querten 1.9 DetId Detid;
132     int SubDet;
133 querten 1.1
134 querten 1.3 MyEvents = new MySimEvents();
135     Geom = new Geometry();
136    
137     edm::ESHandle<TrackerGeometry> tkGeom;
138     iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
139 querten 1.9 vector<GeomDet*> TrackerDets = tkGeom->dets();
140 querten 1.3
141 querten 1.9 for(unsigned int i=0;i<TrackerDets.size();i++){
142     Detid = TrackerDets[i]->geographicalId();
143     // SubDet = Detid.subdetId();
144 querten 1.4
145 querten 1.9 GeomDet* DetUnit = TrackerDets[i];
146 querten 1.3 if(!DetUnit)continue;
147 querten 1.7 const BoundPlane plane = DetUnit->surface();
148     const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
149     const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
150 querten 1.3
151 querten 1.9 float width = 0;
152     float length = 0;
153     float thickness = 0;
154     float TrapezoidalParam = 0;
155 querten 1.7
156     if(trapezoidalBounds)
157     {
158     std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
159 querten 1.8 width = parameters[0]*2;
160     length = parameters[3]*2;
161 querten 1.7 thickness = (*trapezoidalBounds).thickness();
162     TrapezoidalParam = parameters[1]/parameters[0];
163     }else if(rectangularBounds){
164     width = DetUnit->surface().bounds().width();
165     length = DetUnit->surface().bounds().length();
166     thickness = DetUnit->surface().bounds().thickness();
167     TrapezoidalParam = 1;
168     }
169 querten 1.4
170     Surface::GlobalPoint WidthVector = plane.toGlobal( LocalPoint(width/2, 0, 0) );
171     Surface::GlobalPoint LengthVector = plane.toGlobal( LocalPoint(0, length/2, 0) );
172     Surface::GlobalPoint ThickVector = plane.toGlobal( LocalPoint(0, 0, thickness/2) );
173    
174 querten 1.7 GlobalVector Pos = GlobalVector(DetUnit->position().basicVector());
175 querten 1.4
176 querten 1.7 Geom->Add_TrackerDet(Detid.rawId(), TrapezoidalParam,
177     Pos.x(), Pos.y(), Pos.z(),
178     WidthVector.x() -Pos.x(), WidthVector.y() -Pos.y(), WidthVector.z() -Pos.z(),
179     LengthVector.x()-Pos.x(), LengthVector.y()-Pos.y(), LengthVector.z()-Pos.z(),
180     ThickVector.x() -Pos.x(), ThickVector.y() -Pos.y(), ThickVector.z() -Pos.z());
181    
182 querten 1.3
183     }
184 querten 1.9
185    
186     // ### CALO GEOMETRY ###
187    
188     edm::ESHandle<CaloGeometry> CaloGeom;
189     iSetup.get<IdealGeometryRecord>().get( CaloGeom );
190     const vector<DetId> CaloDets = CaloGeom->getValidDetIds();
191    
192     for(unsigned int i=0;i<CaloDets.size();i++)
193     {
194     Detid = CaloDets[i];
195     SubDet = Detid.subdetId();
196    
197 querten 1.10 if(Detid.det()!=DetId::Ecal) continue;
198 querten 1.9
199     //EcalBarrel=1, EcalEndcap=2, EcalPreshower=3, EcalTriggerTower=4, EcalLaserPnDiode=5
200 querten 1.10 if(SubDet!=EcalBarrel)continue;
201 querten 1.9
202     const CaloCellGeometry* CellGeom = CaloGeom->getGeometry(Detid);
203     GlobalPoint CellPos = CellGeom->getPosition();
204    
205     const CaloCellGeometry::CornersVec CellCorners = CellGeom->getCorners();
206 querten 1.10
207 querten 1.9 printf("Det = %i Subdet = %i --> %6.2f %6.2f %6.2f - %i\n", Detid.det(), SubDet, CellPos.x(), CellPos.y(), CellPos.z(), CellCorners.size());
208     }
209    
210    
211 querten 1.3 Geom->Save("Tracker.geom");
212 querten 1.1
213    
214     }
215    
216     // ------------ method called once each job just after ending the event loop ------------
217     void
218     OpenGLDisplayer::endJob() {
219    
220     MyEvents->Save((char*) OutputFile.c_str());
221     MyEvents->Load((char*) OutputFile.c_str());
222     }
223    
224    
225    
226     //
227     // member functions
228     //
229    
230     // ------------ method called to for each event ------------
231     void
232     OpenGLDisplayer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
233     {
234    
235     edm::Handle<std::vector< SimTrack > > h_SimTracks;
236     iEvent.getByLabel("g4SimHits", h_SimTracks);
237     std::vector< SimTrack > SimTrackColl = *h_SimTracks.product();
238    
239     edm::Handle<std::vector< SimVertex > > h_Vertex;
240     iEvent.getByLabel("g4SimHits", h_Vertex);
241     std::vector< SimVertex > VertexColl = *h_Vertex.product();
242    
243    
244    
245     // access the tracker
246     edm::ESHandle<TrackerGeometry> theTrackerGeometry;
247     iSetup.get<TrackerDigiGeometryRecord>().get(theTrackerGeometry);
248     const TrackerGeometry& theTracker(*theTrackerGeometry);
249    
250     // the DetUnits
251     TrackingGeometry::DetContainer theDetUnits = theTracker.dets();
252    
253    
254    
255    
256     // Save the data
257    
258     MySimEvent* MyEvent = new MySimEvent;
259    
260     for ( unsigned int a = 0; a < SimTrackColl.size(); ++a ) {
261     MySimTrack MysimTrack;
262     SimTrack simTrack =SimTrackColl[a];
263    
264     MysimTrack.track_id =simTrack.trackId();
265     MysimTrack.Type =simTrack.type();
266     MysimTrack.parent_vertex=simTrack.vertIndex();
267     MysimTrack.Px =simTrack.momentum().x();
268     MysimTrack.Py =simTrack.momentum().y();
269     MysimTrack.Pz =simTrack.momentum().z();
270     MysimTrack.E =simTrack.momentum().e();
271     MysimTrack.charge =simTrack.charge();
272    
273     MyEvent->MySimTrackCollection.push_back(MysimTrack);
274     }
275    
276    
277     for (unsigned int b = 0; b < VertexColl.size(); ++b ) {
278     MySimVertex MyVertex;
279     SimVertex Vertex = VertexColl[b];
280    
281     MyVertex.parentTrack_id =Vertex.parentIndex ();
282     MyVertex.x =Vertex.position().x();
283     MyVertex.y =Vertex.position().y();
284     MyVertex.z =Vertex.position().z();
285    
286     MyEvent->MySimVertexCollection.push_back(MyVertex);
287     }
288    
289     for(unsigned int i=0; i<SimHitSubdetectors.size(); i++)
290     {
291     edm::Handle<std::vector< PSimHit > > h_Hits;
292     iEvent.getByLabel("g4SimHits",SimHitSubdetectors[i].c_str(), h_Hits);
293     std::vector< PSimHit > Hits = *h_Hits.product();
294    
295     for(unsigned int h=0; h<Hits.size(); h++)
296     {
297     DetId theDetUnitId(Hits[h].detUnitId());
298     const GeomDet * theDet = theTracker.idToDet(theDetUnitId);
299    
300     MyPSimHit Hit;
301     Hit.x = theDet->surface().toGlobal(Hits[h].localPosition()).x();
302     Hit.y = theDet->surface().toGlobal(Hits[h].localPosition()).y();
303     Hit.z = theDet->surface().toGlobal(Hits[h].localPosition()).z();
304     Hit.ProcessType = Hits[h].processType();
305     Hit.dEdX = Hits[h].energyLoss();
306    
307     MyEvent->MyPSimHitCollection.push_back(Hit);
308    
309     }
310     }
311    
312     edm::Handle<std::vector< reco::Track > > h_Tracks;
313     iEvent.getByLabel("ctfWithMaterialTracks", h_Tracks);
314     std::vector< reco::Track > TrackColl = *h_Tracks.product();
315    
316     for ( unsigned int t = 0; t < TrackColl.size(); ++t ) {
317     MyRecoTrack MyrecoTrack;
318     reco::Track recoTrack =TrackColl[t];
319    
320 querten 1.2 // MyrecoTrack.Hits = new MyRecoHit[MyrecoTrack.N];
321 querten 1.1
322 querten 1.2 for(unsigned int h=0;h<recoTrack.recHitsSize();h++){
323 querten 1.1 TrackingRecHitRef h_it = recoTrack.recHit(h);
324     if(!h_it->isValid() )continue;
325     DetId detId = h_it->geographicalId();
326     const GeomDet * theDet = theTracker.idToDet(detId);
327     LocalPoint localPos = h_it->localPosition();
328 querten 1.2
329     MyRecoHit hit;
330     hit.x = theDet->surface().toGlobal(localPos).x();
331     hit.y = theDet->surface().toGlobal(localPos).y();
332     hit.z = theDet->surface().toGlobal(localPos).z();
333     hit.DetId = detId.rawId();
334     hit.Charge = -1;
335     MyrecoTrack.Hits.push_back(hit);
336    
337    
338     // printf("%8.2f %8.2f %8.2f\n",theDet->surface().toGlobal(localPos).x(),theDet->surface().toGlobal(localPos).y(),theDet->surface().toGlobal(localPos).z());
339 querten 1.1 }
340     MyEvent->MyRecoTrackCollection.push_back(MyrecoTrack);
341     }
342    
343     MyEvents->Events.push_back(MyEvent);
344    
345     }
346    
347    
348     //define this as a plug-in
349     DEFINE_FWK_MODULE(OpenGLDisplayer);
350    
351