CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiPixelErrorEstimation.cc
Go to the documentation of this file.
1 
2 // Package: SiPixelErrorEstimation
3 // Class: SiPixelErrorEstimation
4 // Original Author: Gavril Giurgiu (JHU)
5 // Created: Fri May 4 17:48:24 CDT 2007
6 
7 #include <iostream>
8 #include "CalibTracker/SiPixelErrorEstimation/interface/SiPixelErrorEstimation.h"
9 
10 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
11 #include "DataFormats/GeometryVector/interface/LocalVector.h"
12 
13 #include "DataFormats/TrackReco/interface/Track.h"
14 #include "DataFormats/DetId/interface/DetId.h"
15 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
16 
17 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
18 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
19 #include "Geometry/Records/interface/IdealGeometryRecord.h"
20 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
21 #include "DataFormats/TrackReco/interface/TrackFwd.h"
22 
23 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
24 
25 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
26 
27 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
28 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
29 
30 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
31 #include "FWCore/ParameterSet/interface/FileInPath.h"
32 
33 using namespace std;
34 using namespace edm;
35 
36 SiPixelErrorEstimation::SiPixelErrorEstimation(const edm::ParameterSet& ps):tfile_(0), ttree_all_hits_(0),
37  ttree_track_hits_(0), ttree_track_hits_strip_(0)
38 {
39  //Read config file
40  outputFile_ = ps.getUntrackedParameter<string>( "outputFile", "SiPixelErrorEstimation_Ntuple.root" );
41 
42  // Replace "ctfWithMaterialTracks" with "generalTracks"
43  //src_ = ps.getUntrackedParameter<std::string>( "src", "ctfWithMaterialTracks" );
44  src_ = ps.getUntrackedParameter<std::string>( "src", "generalTracks" );
45 
46  checkType_ = ps.getParameter<bool>( "checkType" );
47  genType_ = ps.getParameter<int>( "genType" );
48  include_trk_hits_ = ps.getParameter<bool>( "include_trk_hits" );
49 
50  tTrajectory = consumes<std::vector<Trajectory>> (src_);
51  tPixRecHitCollection = consumes<SiPixelRecHitCollection>(edm::InputTag( "siPixelRecHits"));
52  tSimTrackContainer = consumes <edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
53  tTrackCollection = consumes<reco::TrackCollection>(src_);
54 
55 }
56 
57 SiPixelErrorEstimation::~SiPixelErrorEstimation()
58 {}
59 
60 void SiPixelErrorEstimation::beginJob()
61 {
62  int bufsize = 64000;
63 
64 
65  tfile_ = new TFile ( outputFile_.c_str() , "RECREATE");
66 
67 
68  ttree_track_hits_strip_ = new TTree("TrackHitNtupleStrip", "TrackHitNtupleStrip");
69 
70  ttree_track_hits_strip_->Branch("strip_rechitx", &strip_rechitx, "strip_rechitx/F" , bufsize);
71  ttree_track_hits_strip_->Branch("strip_rechity", &strip_rechity, "strip_rechity/F" , bufsize);
72  ttree_track_hits_strip_->Branch("strip_rechitz", &strip_rechitz, "strip_rechitz/F" , bufsize);
73 
74  ttree_track_hits_strip_->Branch("strip_rechiterrx", &strip_rechiterrx, "strip_rechiterrx/F" , bufsize);
75  ttree_track_hits_strip_->Branch("strip_rechiterry", &strip_rechiterry, "strip_rechiterry/F" , bufsize);
76 
77  ttree_track_hits_strip_->Branch("strip_rechitresx", &strip_rechitresx, "strip_rechitresx/F" , bufsize);
78 
79  ttree_track_hits_strip_->Branch("strip_rechitresx2", &strip_rechitresx2, "strip_rechitresx2/F" , bufsize);
80 
81 
82  ttree_track_hits_strip_->Branch("strip_rechitresy", &strip_rechitresy, "strip_rechitresy/F" , bufsize);
83 
84  ttree_track_hits_strip_->Branch("strip_rechitpullx", &strip_rechitpullx, "strip_rechitpullx/F", bufsize);
85  ttree_track_hits_strip_->Branch("strip_rechitpully", &strip_rechitpully, "strip_rechitpully/F", bufsize);
86 
87  ttree_track_hits_strip_->Branch("strip_is_stereo", &strip_is_stereo, "strip_is_stereo/I", bufsize);
88  ttree_track_hits_strip_->Branch("strip_hit_type" , &strip_hit_type , "strip_hit_type/I" , bufsize);
89  ttree_track_hits_strip_->Branch("detector_type" , &detector_type , "detector_type/I" , bufsize);
90 
91  ttree_track_hits_strip_->Branch("strip_trk_pt" , &strip_trk_pt , "strip_trk_pt/F" , bufsize);
92  ttree_track_hits_strip_->Branch("strip_cotalpha" , &strip_cotalpha , "strip_cotalpha/F" , bufsize);
93  ttree_track_hits_strip_->Branch("strip_cotbeta" , &strip_cotbeta , "strip_cotbeta/F" , bufsize);
94  ttree_track_hits_strip_->Branch("strip_locbx" , &strip_locbx , "strip_locbx/F" , bufsize);
95  ttree_track_hits_strip_->Branch("strip_locby" , &strip_locby , "strip_locby/F" , bufsize);
96  ttree_track_hits_strip_->Branch("strip_locbz" , &strip_locbz , "strip_locbz/F" , bufsize);
97  ttree_track_hits_strip_->Branch("strip_charge" , &strip_charge , "strip_charge/F" , bufsize);
98  ttree_track_hits_strip_->Branch("strip_size" , &strip_size , "strip_size/I" , bufsize);
99 
100 
101  ttree_track_hits_strip_->Branch("strip_edge" , &strip_edge , "strip_edge/I" , bufsize);
102  ttree_track_hits_strip_->Branch("strip_nsimhit", &strip_nsimhit, "strip_nsimhit/I" , bufsize);
103  ttree_track_hits_strip_->Branch("strip_pidhit" , &strip_pidhit , "strip_pidhit/I" , bufsize);
104  ttree_track_hits_strip_->Branch("strip_simproc", &strip_simproc, "strip_simproc/I" , bufsize);
105 
106 
107 ttree_track_hits_strip_->Branch("strip_subdet_id" , &strip_subdet_id , "strip_subdet_id/I" , bufsize);
108 
109 ttree_track_hits_strip_->Branch("strip_tib_layer" , &strip_tib_layer , "strip_tib_layer/I" , bufsize);
110 ttree_track_hits_strip_->Branch("strip_tib_module" , &strip_tib_module , "strip_tib_module/I" , bufsize);
111 ttree_track_hits_strip_->Branch("strip_tib_order" , &strip_tib_order , "strip_tib_order/I" , bufsize);
112 ttree_track_hits_strip_->Branch("strip_tib_side" , &strip_tib_side , "strip_tib_side/I" , bufsize);
113 ttree_track_hits_strip_->Branch("strip_tib_is_double_side" , &strip_tib_is_double_side , "strip_tib_is_double_side/I" , bufsize);
114 ttree_track_hits_strip_->Branch("strip_tib_is_z_plus_side" , &strip_tib_is_z_plus_side , "strip_tib_is_z_plus_side/I" , bufsize);
115 ttree_track_hits_strip_->Branch("strip_tib_is_z_minus_side" , &strip_tib_is_z_minus_side , "strip_tib_is_z_minus_side/I" , bufsize);
116 ttree_track_hits_strip_->Branch("strip_tib_layer_number" , &strip_tib_layer_number , "strip_tib_layer_number/I" , bufsize);
117 ttree_track_hits_strip_->Branch("strip_tib_string_number" , &strip_tib_string_number , "strip_tib_string_number/I" , bufsize);
118 ttree_track_hits_strip_->Branch("strip_tib_module_number" , &strip_tib_module_number ,"strip_tib_module_number/I" , bufsize);
119 ttree_track_hits_strip_->Branch("strip_tib_is_internal_string", &strip_tib_is_internal_string,"strip_tib_is_internal_string/I", bufsize);
120 ttree_track_hits_strip_->Branch("strip_tib_is_external_string", &strip_tib_is_external_string,"strip_tib_is_external_string/I", bufsize);
121 ttree_track_hits_strip_->Branch("strip_tib_is_rphi" , &strip_tib_is_rphi , "strip_tib_is_rphi/I" , bufsize);
122 ttree_track_hits_strip_->Branch("strip_tib_is_stereo" , &strip_tib_is_stereo , "strip_tib_is_stereo/I" , bufsize);
123 ttree_track_hits_strip_->Branch("strip_tob_layer" , &strip_tob_layer , "strip_tob_layer/I" , bufsize);
124 ttree_track_hits_strip_->Branch("strip_tob_module" , &strip_tob_module , "strip_tob_module/I" , bufsize);
125 ttree_track_hits_strip_->Branch("strip_tob_side" , &strip_tob_side , "strip_tob_side/I" , bufsize);
126 ttree_track_hits_strip_->Branch("strip_tob_is_double_side" , &strip_tob_is_double_side , "strip_tob_is_double_side/I" , bufsize);
127 ttree_track_hits_strip_->Branch("strip_tob_is_z_plus_side" , &strip_tob_is_z_plus_side , "strip_tob_is_z_plus_side/I" , bufsize);
128 ttree_track_hits_strip_->Branch("strip_tob_is_z_minus_side" , &strip_tob_is_z_minus_side , "strip_tob_is_z_minus_side/I" , bufsize);
129 ttree_track_hits_strip_->Branch("strip_tob_layer_number" , &strip_tob_layer_number , "strip_tob_layer_number/I" , bufsize);
130 ttree_track_hits_strip_->Branch("strip_tob_rod_number" , &strip_tob_rod_number , "strip_tob_rod_number/I" , bufsize);
131 ttree_track_hits_strip_->Branch("strip_tob_module_number" , &strip_tob_module_number , "strip_tob_module_number/I" , bufsize);
132 
133 
134 ttree_track_hits_strip_->Branch("strip_prob", &strip_prob, "strip_prob/F" , bufsize);
135 ttree_track_hits_strip_->Branch("strip_qbin", &strip_qbin, "strip_qbin/I", bufsize);
136 
137 ttree_track_hits_strip_->Branch("strip_nprm", &strip_nprm, "strip_nprm/I", bufsize);
138 
139 ttree_track_hits_strip_->Branch("strip_pidhit1" , &strip_pidhit1 , "strip_pidhit1/I" , bufsize);
140 ttree_track_hits_strip_->Branch("strip_simproc1", &strip_simproc1, "strip_simproc1/I" , bufsize);
141 
142 ttree_track_hits_strip_->Branch("strip_pidhit2" , &strip_pidhit2 , "strip_pidhit2/I" , bufsize);
143 ttree_track_hits_strip_->Branch("strip_simproc2", &strip_simproc2, "strip_simproc2/I" , bufsize);
144 
145 ttree_track_hits_strip_->Branch("strip_pidhit3" , &strip_pidhit3 , "strip_pidhit3/I" , bufsize);
146 ttree_track_hits_strip_->Branch("strip_simproc3", &strip_simproc3, "strip_simproc3/I" , bufsize);
147 
148 ttree_track_hits_strip_->Branch("strip_pidhit4" , &strip_pidhit4 , "strip_pidhit4/I" , bufsize);
149 ttree_track_hits_strip_->Branch("strip_simproc4", &strip_simproc4, "strip_simproc4/I" , bufsize);
150 
151 ttree_track_hits_strip_->Branch("strip_pidhit5" , &strip_pidhit5 , "strip_pidhit5/I" , bufsize);
152 ttree_track_hits_strip_->Branch("strip_simproc5", &strip_simproc5, "strip_simproc5/I" , bufsize);
153 
154 ttree_track_hits_strip_->Branch("strip_split", &strip_split, "strip_split/I" , bufsize);
155 
156 ttree_track_hits_strip_->Branch("strip_clst_err_x", &strip_clst_err_x, "strip_clst_err_x/F" , bufsize);
157 ttree_track_hits_strip_->Branch("strip_clst_err_y", &strip_clst_err_y, "strip_clst_err_y/F" , bufsize);
158 
159  if ( include_trk_hits_ )
160  {
161  //tfile_ = new TFile ("SiPixelErrorEstimation_Ntuple.root" , "RECREATE");
162  //const char* tmp_name = outputFile_.c_str();
163 
164 
165  ttree_track_hits_ = new TTree("TrackHitNtuple", "TrackHitNtuple");
166 
167  ttree_track_hits_->Branch("evt", &evt, "evt/I", bufsize);
168  ttree_track_hits_->Branch("run", &run, "run/I", bufsize);
169 
170  ttree_track_hits_->Branch("subdetId", &subdetId, "subdetId/I", bufsize);
171 
172  ttree_track_hits_->Branch("layer" , &layer , "layer/I" , bufsize);
173  ttree_track_hits_->Branch("ladder", &ladder, "ladder/I", bufsize);
174  ttree_track_hits_->Branch("mod" , &mod , "mod/I" , bufsize);
175 
176  ttree_track_hits_->Branch("side" , &side , "side/I" , bufsize);
177  ttree_track_hits_->Branch("disk" , &disk , "disk/I" , bufsize);
178  ttree_track_hits_->Branch("blade" , &blade , "blade/I" , bufsize);
179  ttree_track_hits_->Branch("panel" , &panel , "panel/I" , bufsize);
180  ttree_track_hits_->Branch("plaq" , &plaq , "plaq/I" , bufsize);
181 
182  ttree_track_hits_->Branch("half" , &half , "half/I" , bufsize);
183  ttree_track_hits_->Branch("flipped", &flipped, "flipped/I", bufsize);
184 
185  ttree_track_hits_->Branch("rechitx", &rechitx, "rechitx/F" , bufsize);
186  ttree_track_hits_->Branch("rechity", &rechity, "rechity/F" , bufsize);
187  ttree_track_hits_->Branch("rechitz", &rechitz, "rechitz/F" , bufsize);
188 
189  ttree_track_hits_->Branch("rechiterrx", &rechiterrx, "rechiterrx/F" , bufsize);
190  ttree_track_hits_->Branch("rechiterry", &rechiterry, "rechiterry/F" , bufsize);
191 
192  ttree_track_hits_->Branch("rechitresx", &rechitresx, "rechitresx/F" , bufsize);
193  ttree_track_hits_->Branch("rechitresy", &rechitresy, "rechitresy/F" , bufsize);
194 
195  ttree_track_hits_->Branch("rechitpullx", &rechitpullx, "rechitpullx/F", bufsize);
196  ttree_track_hits_->Branch("rechitpully", &rechitpully, "rechitpully/F", bufsize);
197 
198 
199  ttree_track_hits_->Branch("npix" , &npix , "npix/I" , bufsize);
200  ttree_track_hits_->Branch("nxpix" , &nxpix , "nxpix/I" , bufsize);
201  ttree_track_hits_->Branch("nypix" , &nypix , "nypix/I" , bufsize);
202  ttree_track_hits_->Branch("charge", &charge, "charge/F", bufsize);
203 
204  ttree_track_hits_->Branch("edgex", &edgex, "edgex/I", bufsize);
205  ttree_track_hits_->Branch("edgey", &edgey, "edgey/I", bufsize);
206 
207  ttree_track_hits_->Branch("bigx", &bigx, "bigx/I", bufsize);
208  ttree_track_hits_->Branch("bigy", &bigy, "bigy/I", bufsize);
209 
210  ttree_track_hits_->Branch("alpha", &alpha, "alpha/F", bufsize);
211  ttree_track_hits_->Branch("beta" , &beta , "beta/F" , bufsize);
212 
213  ttree_track_hits_->Branch("trk_alpha", &trk_alpha, "trk_alpha/F", bufsize);
214  ttree_track_hits_->Branch("trk_beta" , &trk_beta , "trk_beta/F" , bufsize);
215 
216  ttree_track_hits_->Branch("phi", &phi, "phi/F", bufsize);
217  ttree_track_hits_->Branch("eta", &eta, "eta/F", bufsize);
218 
219  ttree_track_hits_->Branch("simhitx", &simhitx, "simhitx/F", bufsize);
220  ttree_track_hits_->Branch("simhity", &simhity, "simhity/F", bufsize);
221 
222  ttree_track_hits_->Branch("nsimhit", &nsimhit, "nsimhit/I", bufsize);
223  ttree_track_hits_->Branch("pidhit" , &pidhit , "pidhit/I" , bufsize);
224  ttree_track_hits_->Branch("simproc", &simproc, "simproc/I", bufsize);
225 
226  ttree_track_hits_->Branch("pixel_split", &pixel_split, "pixel_split/I", bufsize);
227 
228  ttree_track_hits_->Branch("pixel_clst_err_x", &pixel_clst_err_x, "pixel_clst_err_x/F" , bufsize);
229  ttree_track_hits_->Branch("pixel_clst_err_y", &pixel_clst_err_y, "pixel_clst_err_y/F" , bufsize);
230 
231  } // if ( include_trk_hits_ )
232 
233  // ----------------------------------------------------------------------
234 
235  ttree_all_hits_ = new TTree("AllHitNtuple", "AllHitNtuple");
236 
237  ttree_all_hits_->Branch("evt", &evt, "evt/I", bufsize);
238  ttree_all_hits_->Branch("run", &run, "run/I", bufsize);
239 
240  ttree_all_hits_->Branch("subdetid", &all_subdetid, "subdetid/I", bufsize);
241 
242  ttree_all_hits_->Branch("layer" , &all_layer , "layer/I" , bufsize);
243  ttree_all_hits_->Branch("ladder", &all_ladder, "ladder/I", bufsize);
244  ttree_all_hits_->Branch("mod" , &all_mod , "mod/I" , bufsize);
245 
246  ttree_all_hits_->Branch("side" , &all_side , "side/I" , bufsize);
247  ttree_all_hits_->Branch("disk" , &all_disk , "disk/I" , bufsize);
248  ttree_all_hits_->Branch("blade" , &all_blade , "blade/I" , bufsize);
249  ttree_all_hits_->Branch("panel" , &all_panel , "panel/I" , bufsize);
250  ttree_all_hits_->Branch("plaq" , &all_plaq , "plaq/I" , bufsize);
251 
252  ttree_all_hits_->Branch("half" , &all_half , "half/I" , bufsize);
253  ttree_all_hits_->Branch("flipped", &all_flipped, "flipped/I", bufsize);
254 
255  ttree_all_hits_->Branch("cols", &all_cols, "cols/I", bufsize);
256  ttree_all_hits_->Branch("rows", &all_rows, "rows/I", bufsize);
257 
258  ttree_all_hits_->Branch("rechitx" , &all_rechitx , "rechitx/F" , bufsize);
259  ttree_all_hits_->Branch("rechity" , &all_rechity , "rechity/F" , bufsize);
260  ttree_all_hits_->Branch("rechitz" , &all_rechitz , "rechitz/F" , bufsize);
261 
262  ttree_all_hits_->Branch("rechiterrx" , &all_rechiterrx , "rechiterrx/F" , bufsize);
263  ttree_all_hits_->Branch("rechiterry" , &all_rechiterry , "rechiterry/F" , bufsize);
264 
265  ttree_all_hits_->Branch("rechitresx" , &all_rechitresx , "rechitresx/F" , bufsize);
266  ttree_all_hits_->Branch("rechitresy" , &all_rechitresy , "rechitresy/F" , bufsize);
267 
268  ttree_all_hits_->Branch("rechitpullx", &all_rechitpullx, "rechitpullx/F", bufsize);
269  ttree_all_hits_->Branch("rechitpully", &all_rechitpully, "rechitpully/F", bufsize);
270 
271  ttree_all_hits_->Branch("npix" , &all_npix , "npix/I" , bufsize);
272  ttree_all_hits_->Branch("nxpix" , &all_nxpix , "nxpix/I" , bufsize);
273  ttree_all_hits_->Branch("nypix" , &all_nypix , "nypix/I" , bufsize);
274 
275  ttree_all_hits_->Branch("edgex", &all_edgex, "edgex/I", bufsize);
276  ttree_all_hits_->Branch("edgey", &all_edgey, "edgey/I", bufsize);
277 
278  ttree_all_hits_->Branch("bigx", &all_bigx, "bigx/I", bufsize);
279  ttree_all_hits_->Branch("bigy", &all_bigy, "bigy/I", bufsize);
280 
281  ttree_all_hits_->Branch("alpha", &all_alpha, "alpha/F", bufsize);
282  ttree_all_hits_->Branch("beta" , &all_beta , "beta/F" , bufsize);
283 
284  ttree_all_hits_->Branch("simhitx", &all_simhitx, "simhitx/F", bufsize);
285  ttree_all_hits_->Branch("simhity", &all_simhity, "simhity/F", bufsize);
286 
287  ttree_all_hits_->Branch("nsimhit", &all_nsimhit, "nsimhit/I", bufsize);
288  ttree_all_hits_->Branch("pidhit" , &all_pidhit , "pidhit/I" , bufsize);
289  ttree_all_hits_->Branch("simproc", &all_simproc, "simproc/I", bufsize);
290 
291  ttree_all_hits_->Branch("vtxr", &all_vtxr, "vtxr/F", bufsize);
292  ttree_all_hits_->Branch("vtxz", &all_vtxz, "vtxz/F", bufsize);
293 
294  ttree_all_hits_->Branch("simpx", &all_simpx, "simpx/F", bufsize);
295  ttree_all_hits_->Branch("simpy", &all_simpy, "simpy/F", bufsize);
296  ttree_all_hits_->Branch("simpz", &all_simpz, "simpz/F", bufsize);
297 
298  ttree_all_hits_->Branch("eloss", &all_eloss, "eloss/F", bufsize);
299 
300  ttree_all_hits_->Branch("simphi", &all_simphi, "simphi/F", bufsize);
301  ttree_all_hits_->Branch("simtheta", &all_simtheta, "simtheta/F", bufsize);
302 
303  ttree_all_hits_->Branch("trkid", &all_trkid, "trkid/I", bufsize);
304 
305  ttree_all_hits_->Branch("x1", &all_x1, "x1/F", bufsize);
306  ttree_all_hits_->Branch("x2", &all_x2, "x2/F", bufsize);
307  ttree_all_hits_->Branch("y1", &all_x1, "y1/F", bufsize);
308  ttree_all_hits_->Branch("y2", &all_x2, "y2/F", bufsize);
309  ttree_all_hits_->Branch("z1", &all_x1, "z1/F", bufsize);
310  ttree_all_hits_->Branch("z2", &all_x2, "z2/F", bufsize);
311 
312  ttree_all_hits_->Branch("row1", &all_row1, "row1/F", bufsize);
313  ttree_all_hits_->Branch("row2", &all_row2, "row2/F", bufsize);
314  ttree_all_hits_->Branch("col1", &all_col1, "col1/F", bufsize);
315  ttree_all_hits_->Branch("col2", &all_col2, "col2/F", bufsize);
316 
317  ttree_all_hits_->Branch("gx1", &all_gx1, "gx1/F", bufsize);
318  ttree_all_hits_->Branch("gx2", &all_gx2, "gx2/F", bufsize);
319  ttree_all_hits_->Branch("gy1", &all_gx1, "gy1/F", bufsize);
320  ttree_all_hits_->Branch("gy2", &all_gx2, "gy2/F", bufsize);
321  ttree_all_hits_->Branch("gz1", &all_gx1, "gz1/F", bufsize);
322  ttree_all_hits_->Branch("gz2", &all_gx2, "gz2/F", bufsize);
323 
324  ttree_all_hits_->Branch("simtrketa", &all_simtrketa, "simtrketa/F", bufsize);
325  ttree_all_hits_->Branch("simtrkphi", &all_simtrkphi, "simtrkphi/F", bufsize);
326 
327  ttree_all_hits_->Branch("clust_row", &all_clust_row, "clust_row/F", bufsize);
328  ttree_all_hits_->Branch("clust_col", &all_clust_col, "clust_col/F", bufsize);
329 
330  ttree_all_hits_->Branch("clust_x", &all_clust_x, "clust_x/F", bufsize);
331  ttree_all_hits_->Branch("clust_y", &all_clust_y, "clust_y/F", bufsize);
332 
333  ttree_all_hits_->Branch("clust_q", &all_clust_q, "clust_q/F", bufsize);
334 
335  ttree_all_hits_->Branch("clust_maxpixcol", &all_clust_maxpixcol, "clust_maxpixcol/I", bufsize);
336  ttree_all_hits_->Branch("clust_maxpixrow", &all_clust_maxpixrow, "clust_maxpixrow/I", bufsize);
337  ttree_all_hits_->Branch("clust_minpixcol", &all_clust_minpixcol, "clust_minpixcol/I", bufsize);
338  ttree_all_hits_->Branch("clust_minpixrow", &all_clust_minpixrow, "clust_minpixrow/I", bufsize);
339 
340  ttree_all_hits_->Branch("clust_geoid", &all_clust_geoid, "clust_geoid/I", bufsize);
341 
342  ttree_all_hits_->Branch("clust_alpha", &all_clust_alpha, "clust_alpha/F", bufsize);
343  ttree_all_hits_->Branch("clust_beta" , &all_clust_beta , "clust_beta/F" , bufsize);
344 
345  ttree_all_hits_->Branch("rowpix", all_pixrow, "row[npix]/F", bufsize);
346  ttree_all_hits_->Branch("colpix", all_pixcol, "col[npix]/F", bufsize);
347  ttree_all_hits_->Branch("adc", all_pixadc, "adc[npix]/F", bufsize);
348  ttree_all_hits_->Branch("xpix", all_pixx, "x[npix]/F", bufsize);
349  ttree_all_hits_->Branch("ypix", all_pixy, "y[npix]/F", bufsize);
350  ttree_all_hits_->Branch("gxpix", all_pixgx, "gx[npix]/F", bufsize);
351  ttree_all_hits_->Branch("gypix", all_pixgy, "gy[npix]/F", bufsize);
352  ttree_all_hits_->Branch("gzpix", all_pixgz, "gz[npix]/F", bufsize);
353 
354  ttree_all_hits_->Branch("hit_probx", &all_hit_probx, "hit_probx/F" , bufsize);
355  ttree_all_hits_->Branch("hit_proby", &all_hit_proby, "hit_proby/F" , bufsize);
356 
357  ttree_all_hits_->Branch("all_pixel_split", &all_pixel_split, "all_pixel_split/I" , bufsize);
358 
359  ttree_all_hits_->Branch("all_pixel_clst_err_x", &all_pixel_clst_err_x, "all_pixel_clst_err_x/F" , bufsize);
360  ttree_all_hits_->Branch("all_pixel_clst_err_y", &all_pixel_clst_err_y, "all_pixel_clst_err_y/F" , bufsize);
361 
362 }
363 
364 void SiPixelErrorEstimation::endJob()
365 {
366  tfile_->Write();
367  tfile_->Close();
368 }
369 
370 void
371 SiPixelErrorEstimation::analyze(const edm::Event& e, const edm::EventSetup& es)
372 {
373  //Retrieve tracker topology from geometry
374  edm::ESHandle<TrackerTopology> tTopoHandle;
375  es.get<IdealGeometryRecord>().get(tTopoHandle);
376  const TrackerTopology* const tTopo = tTopoHandle.product();
377 
378  using namespace edm;
379 
380  run = e.id().run();
381  evt = e.id().event();
382 
383  if ( evt%1000 == 0 )
384  cout << "evt = " << evt << endl;
385 
386  float math_pi = 3.14159265;
387  float radtodeg = 180.0 / math_pi;
388 
389  LocalPoint position;
390  LocalError error;
391  float mindist = 999999.9;
392 
393  std::vector<PSimHit> matched;
394  TrackerHitAssociator associate(e);
395 
396  edm::ESHandle<TrackerGeometry> pDD;
397  es.get<TrackerDigiGeometryRecord> ().get (pDD);
398  const TrackerGeometry* tracker = &(* pDD);
399 
400 
401  //cout << "...1..." << endl;
402 
403 
404  edm::ESHandle<MagneticField> magneticField;
405  es.get<IdealMagneticFieldRecord>().get( magneticField );
406  //const MagneticField* magField_ = magFieldHandle.product();
407 
408  edm::FileInPath FileInPath_;
409 
410 
411 
412  // Strip hits ==============================================================================================================
413 
414 
415  edm::Handle<vector<Trajectory> > trajCollectionHandle;
416 
417  e.getByToken( tTrajectory, trajCollectionHandle);
418  //e.getByLabel( "generalTracks", trajCollectionHandle);
419 
420  for ( vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end(); ++it )
421  {
422 
423  vector<TrajectoryMeasurement> tmColl = it->measurements();
424  for ( vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj!=tmColl.end(); ++itTraj )
425  {
426 
427  if ( !itTraj->updatedState().isValid() )
428  continue;
429 
430  strip_rechitx = -9999999.9;
431  strip_rechity = -9999999.9;
432  strip_rechitz = -9999999.9;
433  strip_rechiterrx = -9999999.9;
434  strip_rechiterry = -9999999.9;
435  strip_rechitresx = -9999999.9;
436  strip_rechitresx2 = -9999999.9;
437 
438 
439  strip_rechitresy = -9999999.9;
440  strip_rechitpullx = -9999999.9;
441  strip_rechitpully = -9999999.9;
442  strip_is_stereo = -9999999 ;
443  strip_hit_type = -9999999 ;
444  detector_type = -9999999 ;
445 
446  strip_trk_pt = -9999999.9;
447  strip_cotalpha = -9999999.9;
448  strip_cotbeta = -9999999.9;
449  strip_locbx = -9999999.9;
450  strip_locby = -9999999.9;
451  strip_locbz = -9999999.9;
452  strip_charge = -9999999.9;
453  strip_size = -9999999 ;
454 
455  strip_edge = -9999999 ;
456  strip_nsimhit = -9999999 ;
457  strip_pidhit = -9999999 ;
458  strip_simproc = -9999999 ;
459 
460 
461  strip_subdet_id = -9999999;
462 
463  strip_tib_layer = -9999999;
464  strip_tib_module = -9999999;
465  strip_tib_order = -9999999;
466  strip_tib_side = -9999999;
467  strip_tib_is_double_side = -9999999;
468  strip_tib_is_z_plus_side = -9999999;
469  strip_tib_is_z_minus_side = -9999999;
470  strip_tib_layer_number = -9999999;
471  strip_tib_string_number = -9999999;
472  strip_tib_module_number = -9999999;
473  strip_tib_is_internal_string= -9999999;
474  strip_tib_is_external_string= -9999999;
475  strip_tib_is_rphi = -9999999;
476  strip_tib_is_stereo = -9999999;
477 
478  strip_tob_layer = -9999999;
479  strip_tob_module = -9999999;
480  strip_tob_side = -9999999;
481  strip_tob_is_double_side = -9999999;
482  strip_tob_is_z_plus_side = -9999999;
483  strip_tob_is_z_minus_side = -9999999;
484  strip_tob_layer_number = -9999999;
485  strip_tob_rod_number = -9999999;
486  strip_tob_module_number = -9999999;
487 
488  strip_tob_is_rphi = -9999999;
489  strip_tob_is_stereo = -9999999;
490 
491  strip_prob = -9999999.9;
492  strip_qbin = -9999999;
493 
494  strip_nprm = -9999999;
495 
496  strip_pidhit1 = -9999999 ;
497  strip_simproc1 = -9999999 ;
498 
499  strip_pidhit2 = -9999999 ;
500  strip_simproc2 = -9999999 ;
501 
502  strip_pidhit3 = -9999999 ;
503  strip_simproc3 = -9999999 ;
504 
505  strip_pidhit4 = -9999999 ;
506  strip_simproc4 = -9999999 ;
507 
508  strip_pidhit5 = -9999999 ;
509  strip_simproc5 = -9999999 ;
510 
511  strip_split = -9999999;
512  strip_clst_err_x = -9999999.9;
513  strip_clst_err_y = -9999999.9;
514 
515  const TransientTrackingRecHit::ConstRecHitPointer trans_trk_rec_hit_point = itTraj->recHit();
516 
517  if ( trans_trk_rec_hit_point == NULL )
518  continue;
519 
520  const TrackingRecHit *trk_rec_hit = (*trans_trk_rec_hit_point).hit();
521 
522  if ( trk_rec_hit == NULL )
523  continue;
524 
525  DetId detid = (trk_rec_hit)->geographicalId();
526 
527  strip_subdet_id = (int)detid.subdetId();
528 
529  if ( (int)detid.subdetId() != (int)(StripSubdetector::TIB) && (int)detid.subdetId() != (int)(StripSubdetector::TOB) )
530  continue;
531 
532  const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>( (*trans_trk_rec_hit_point).hit() );
533  const SiStripRecHit2D * hit2d = dynamic_cast<const SiStripRecHit2D *>( (*trans_trk_rec_hit_point).hit() );
534  const SiStripRecHit1D * hit1d = dynamic_cast<const SiStripRecHit1D *>( (*trans_trk_rec_hit_point).hit() );
535 
536  if ( !matchedhit && !hit2d && !hit1d )
537  continue;
538 
539  position = (trk_rec_hit)->localPosition();
540  error = (trk_rec_hit)->localPositionError();
541 
542  strip_rechitx = position.x();
543  strip_rechity = position.y();
544  strip_rechitz = position.z();
545  strip_rechiterrx = sqrt( error.xx() );
546  strip_rechiterry = sqrt( error.yy() );
547 
548 
549  //cout << "strip_rechitx = " << strip_rechitx << endl;
550  //cout << "strip_rechity = " << strip_rechity << endl;
551  //cout << "strip_rechitz = " << strip_rechitz << endl;
552 
553  //cout << "strip_rechiterrx = " << strip_rechiterrx << endl;
554  //cout << "strip_rechiterry = " << strip_rechiterry << endl;
555 
556  TrajectoryStateOnSurface tsos = itTraj->updatedState();
557 
558  strip_trk_pt = tsos.globalMomentum().perp();
559 
560  LocalTrajectoryParameters ltp = tsos.localParameters();
561 
562  LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
563 
564  float locx = localDir.x();
565  float locy = localDir.y();
566  float locz = localDir.z();
567 
568  //alpha_ = atan2( locz, locx );
569  //beta_ = atan2( locz, locy );
570 
571  strip_cotalpha = locx/locz;
572  strip_cotbeta = locy/locz;
573 
574 
575  StripSubdetector StripSubdet = (StripSubdetector)detid;
576 
577  if ( StripSubdet.stereo() == 0 )
578  strip_is_stereo = 0;
579  else
580  strip_is_stereo = 1;
581 
582 
583  SiStripDetId si_strip_det_id = SiStripDetId( detid );
584 
585  //const StripGeomDetUnit* strip_geom_det_unit = dynamic_cast<const StripGeomDetUnit*> ( tracker->idToDet( detid ) );
586  const StripGeomDetUnit* strip_geom_det_unit = (const StripGeomDetUnit*)tracker->idToDetUnit( detid );
587 
588  if ( strip_geom_det_unit != NULL )
589  {
590  LocalVector lbfield
591  = (strip_geom_det_unit->surface()).toLocal( (*magneticField).inTesla( strip_geom_det_unit->surface().position() ) );
592 
593  strip_locbx = lbfield.x();
594  strip_locby = lbfield.y();
595  strip_locbz = lbfield.z();
596  }
597 
598 
599  //enum ModuleGeometry {UNKNOWNGEOMETRY, IB1, IB2, OB1, OB2, W1A, W2A, W3A, W1B, W2B, W3B, W4, W5, W6, W7};
600 
601  if ( si_strip_det_id.moduleGeometry() == 1 )
602  {
603  detector_type = 1;
604  //cout << "si_strip_det_id.moduleGeometry() = IB1" << endl;
605  //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
606  }
607  else if ( si_strip_det_id.moduleGeometry() == 2 )
608  {
609  detector_type = 2;
610  //cout << "si_strip_det_id.moduleGeometry() = IB2" << endl;
611  //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
612  }
613  else if ( si_strip_det_id.moduleGeometry() == 3 )
614  {
615  detector_type = 3;
616  //cout << "si_strip_det_id.moduleGeometry() = OB1" << endl;
617  //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
618  }
619  else if ( si_strip_det_id.moduleGeometry() == 4 )
620  {
621  detector_type = 4;
622  //cout << "si_strip_det_id.moduleGeometry() = OB2" << endl;
623  //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
624  }
625 
626 
627  // Store ntuple variables
628 
629  if ( (int)detid.subdetId() == int(StripSubdetector::TIB) )
630  {
631 
632 
633 
634  strip_tib_layer = (int)tTopo->tibLayer( detid );
635  strip_tib_module = (int)tTopo->tibModule( detid );
636  strip_tib_order = (int)tTopo->tibOrder( detid );
637  strip_tib_side = (int)tTopo->tibSide( detid );
638  strip_tib_is_double_side = (int)tTopo->tibIsDoubleSide( detid );
639  strip_tib_is_z_plus_side = (int)tTopo->tibIsZPlusSide( detid );
640  strip_tib_is_z_minus_side = (int)tTopo->tibIsZMinusSide( detid );
641  strip_tib_layer_number = (int)tTopo->tibLayer( detid );
642  strip_tib_string_number = (int)tTopo->tibString( detid ) ;
643  strip_tib_module_number = (int)tTopo->tibModule( detid );
644  strip_tib_is_internal_string = (int)tTopo->tibIsInternalString( detid );
645  strip_tib_is_external_string = (int)tTopo->tibIsExternalString( detid );
646  strip_tib_is_rphi = (int)tTopo->tibIsRPhi( detid );
647  strip_tib_is_stereo = (int)tTopo->tibIsStereo( detid );
648  }
649 
650 
651  if ( (int)detid.subdetId() == int(StripSubdetector::TOB) )
652  {
653 
654 
655 
656  strip_tob_layer = (int)tTopo->tobLayer( detid );
657  strip_tob_module = (int)tTopo->tobModule( detid );
658 
659  strip_tob_side = (int)tTopo->tobSide( detid );
660  strip_tob_is_double_side = (int)tTopo->tobIsDoubleSide( detid );
661  strip_tob_is_z_plus_side = (int)tTopo->tobIsZPlusSide( detid );
662  strip_tob_is_z_minus_side = (int)tTopo->tobIsZMinusSide( detid );
663  strip_tob_layer_number = (int)tTopo->tobLayer( detid );
664  strip_tob_rod_number = (int)tTopo->tobRod( detid );
665  strip_tob_module_number = (int)tTopo->tobModule( detid );
666 
667 
668  strip_tob_is_rphi = (int)tTopo->tobIsRPhi( detid );
669  strip_tob_is_stereo = (int)tTopo->tobIsStereo( detid );
670 
671  }
672 
673 
674  if ( matchedhit )
675  {
676  //cout << endl << endl << endl;
677  //cout << "Found a SiStripMatchedRecHit2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
678  //cout << endl << endl << endl;
679  strip_hit_type = 0;
680 
681  } // if ( matchedhit )
682 
683 
684  if ( hit1d )
685  {
686  strip_hit_type = 1;
687 
688  const SiStripRecHit1D::ClusterRef & cluster = hit1d->cluster();
689 
690  if ( cluster->getSplitClusterError() > 0.0 )
691  strip_split = 1;
692  else
693  strip_split = 0;
694 
695  strip_clst_err_x = cluster->getSplitClusterError();
696  //strip_clst_err_y = ...
697 
698  // Get cluster total charge
699  const std::vector<uint8_t>& stripCharges = cluster->amplitudes();
700  uint16_t charge = 0;
701  for (unsigned int i = 0; i < stripCharges.size(); ++i)
702  {
703  charge += stripCharges.at(i);
704  }
705 
706  strip_charge = (float)charge;
707  strip_size = (int)( (cluster->amplitudes()).size() );
708 
709 
710  // Association of the rechit to the simhit
711  float mindist = 999999.9;
712  matched.clear();
713 
714  matched = associate.associateHit(*hit1d);
715 
716  strip_nsimhit = (int)matched.size();
717 
718  if ( !matched.empty())
719  {
720  PSimHit closest;
721 
722  // Get the closest simhit
723 
724  int strip_nprimaries = 0;
725  int current_index = 0;
726 
727  for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
728  {
729  ++current_index;
730 
731  if ( (*m).processType() == 2 )
732  ++strip_nprimaries;
733 
734  if ( current_index == 1 )
735  {
736  strip_pidhit1 = (*m).particleType();
737  strip_simproc1 = (*m).processType();
738  }
739  else if ( current_index == 2 )
740  {
741  strip_pidhit2 = (*m).particleType();
742  strip_simproc2 = (*m).processType();
743  }
744  else if ( current_index == 3 )
745  {
746  strip_pidhit3 = (*m).particleType();
747  strip_simproc3 = (*m).processType();
748  }
749  else if ( current_index == 4 )
750  {
751  strip_pidhit4 = (*m).particleType();
752  strip_simproc4 = (*m).processType();
753  }
754  else if ( current_index == 5 )
755  {
756  strip_pidhit5 = (*m).particleType();
757  strip_simproc5 = (*m).processType();
758  }
759 
760 
761  float dist = abs( (hit1d)->localPosition().x() - (*m).localPosition().x() );
762 
763  if ( dist<mindist )
764  {
765  mindist = dist;
766  closest = (*m);
767  }
768 
769  } // for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
770 
771  strip_nprm = strip_nprimaries;
772 
773  strip_rechitresx = strip_rechitx - closest.localPosition().x();
774  strip_rechitresy = strip_rechity - closest.localPosition().y();
775 
776  strip_rechitresx2 = strip_rechitx - matched[0].localPosition().x();
777 
778  strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
779  strip_rechitpully = strip_rechitresy / strip_rechiterry;
780 
781  strip_pidhit = closest.particleType();
782  strip_simproc = closest.processType();
783 
784 
785  } // if( !matched.empty())
786 
787 
788  //strip_prob = (hit1d)->getTemplProb();
789  //strip_qbin = (hit1d)->getTemplQbin();
790 
791  //cout << endl;
792  //cout << "SiPixelErrorEstimation 1d hit: " << endl;
793  //cout << "prob 1d = " << strip_prob << endl;
794  //cout << "qbin 1d = " << strip_qbin << endl;
795  //cout << endl;
796 
797 
798  // Is the cluster on edge ?
799  /*
800  SiStripDetInfoFileReader* reader;
801 
802  FileInPath_("CalibTracker/SiStripCommon/data/SiStripDetInfo.dat");
803 
804  reader = new SiStripDetInfoFileReader( FileInPath_.fullPath() );
805 
806  uint16_t firstStrip = cluster->firstStrip();
807  uint16_t lastStrip = firstStrip + (cluster->amplitudes()).size() -1;
808  unsigned short Nstrips;
809  Nstrips = reader->getNumberOfApvsAndStripLength(id1).first*128;
810 
811  if ( firstStrip == 0 || lastStrip == (Nstrips-1) )
812  strip_edge = 1;
813  else
814  strip_edge = 0;
815  */
816 
817  } // if ( hit1d )
818 
819 
820  if ( hit2d )
821  {
822  strip_hit_type = 2;
823 
824  const SiStripRecHit1D::ClusterRef & cluster = hit2d->cluster();
825 
826  //if ( cluster->getSplitClusterError() > 0.0 )
827  //strip_split = 1;
828  //else
829  //strip_split = 0;
830 
831  //strip_clst_err_x = cluster->getSplitClusterError();
832  //strip_clst_err_y = ...
833 
834  // Get cluster total charge
835  const std::vector<uint8_t>& stripCharges = cluster->amplitudes();
836  uint16_t charge = 0;
837  for (unsigned int i = 0; i < stripCharges.size(); ++i)
838  {
839  charge += stripCharges.at(i);
840  }
841 
842  strip_charge = (float)charge;
843  strip_size = (int)( (cluster->amplitudes()).size() );
844 
845  // Association of the rechit to the simhit
846  float mindist = 999999.9;
847  matched.clear();
848 
849  matched = associate.associateHit(*hit2d);
850 
851  strip_nsimhit = (int)matched.size();
852 
853  if ( !matched.empty())
854  {
855  PSimHit closest;
856 
857  // Get the closest simhit
858 
859  for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
860  {
861  float dist = abs( (hit2d)->localPosition().x() - (*m).localPosition().x() );
862 
863  if ( dist<mindist )
864  {
865  mindist = dist;
866  closest = (*m);
867  }
868  }
869 
870  strip_rechitresx = strip_rechitx - closest.localPosition().x();
871  strip_rechitresy = strip_rechity - closest.localPosition().y();
872 
873  strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
874  strip_rechitpully = strip_rechitresy / strip_rechiterry;
875 
876  strip_pidhit = closest.particleType();
877  strip_simproc = closest.processType();
878 
879 
880  } // if( !matched.empty())
881 
882 
883  //strip_prob = (hit2d)->getTemplProb();
884  //strip_qbin = (hit2d)->getTemplQbin();
885 
886  //cout << endl;
887  //cout << "SiPixelErrorEstimation 2d hit: " << endl;
888  //cout << "prob 2d = " << strip_prob << endl;
889  //cout << "qbin 2d = " << strip_qbin << endl;
890  //cout << endl;
891 
892 
893  } // if ( hit2d )
894 
895 
896  ttree_track_hits_strip_->Fill();
897 
898 
899  } // for ( vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj!=tmColl.end(); ++itTraj )
900 
901  } // for( vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end(); ++it)
902 
903 
904 
905 
906 
907 
908 
909 
910  //cout << "...2..." << endl;
911 
912 
913 
914 
915 
916 
917  // --------------------------------------- all hits -----------------------------------------------------------
918  edm::Handle<SiPixelRecHitCollection> recHitColl;
919  e.getByToken(tPixRecHitCollection, recHitColl);
920 
921  Handle<edm::SimTrackContainer> simtracks;
922  e.getByToken(tSimTrackContainer, simtracks);
923 
924  //-----Iterate over detunits
925  for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin(); it != pDD->dets().end(); it++)
926  {
927  DetId detId = ((*it)->geographicalId());
928 
929  SiPixelRecHitCollection::const_iterator dsmatch = recHitColl->find(detId);
930  if (dsmatch == recHitColl->end()) continue;
931 
932  SiPixelRecHitCollection::DetSet pixelrechitRange = *dsmatch;
933  SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
934  SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
935  SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
936  std::vector<PSimHit> matched;
937 
938  //----Loop over rechits for this detId
939  for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter)
940  {
941  matched.clear();
942  matched = associate.associateHit(*pixeliter);
943  // only consider rechits that have associated simhit
944  // otherwise cannot determine residiual
945  if ( matched.empty() )
946  {
947  cout << "SiPixelErrorEstimation::analyze: rechits without associated simhit !!!!!!!"
948  << endl;
949  continue;
950  }
951 
952  all_subdetid = -9999;
953 
954  all_layer = -9999;
955  all_ladder = -9999;
956  all_mod = -9999;
957 
958  all_side = -9999;
959  all_disk = -9999;
960  all_blade = -9999;
961  all_panel = -9999;
962  all_plaq = -9999;
963 
964  all_half = -9999;
965  all_flipped = -9999;
966 
967  all_cols = -9999;
968  all_rows = -9999;
969 
970  all_rechitx = -9999;
971  all_rechity = -9999;
972  all_rechitz = -9999;
973 
974  all_simhitx = -9999;
975  all_simhity = -9999;
976 
977  all_rechiterrx = -9999;
978  all_rechiterry = -9999;
979 
980  all_rechitresx = -9999;
981  all_rechitresy = -9999;
982 
983  all_rechitpullx = -9999;
984  all_rechitpully = -9999;
985 
986  all_npix = -9999;
987  all_nxpix = -9999;
988  all_nypix = -9999;
989 
990  all_edgex = -9999;
991  all_edgey = -9999;
992 
993  all_bigx = -9999;
994  all_bigy = -9999;
995 
996  all_alpha = -9999;
997  all_beta = -9999;
998 
999  all_simphi = -9999;
1000  all_simtheta = -9999;
1001 
1002  all_simhitx = -9999;
1003  all_simhity = -9999;
1004 
1005  all_nsimhit = -9999;
1006  all_pidhit = -9999;
1007  all_simproc = -9999;
1008 
1009  all_vtxr = -9999;
1010  all_vtxz = -9999;
1011 
1012  all_simpx = -9999;
1013  all_simpy = -9999;
1014  all_simpz = -9999;
1015 
1016  all_eloss = -9999;
1017 
1018  all_trkid = -9999;
1019 
1020  all_x1 = -9999;
1021  all_x2 = -9999;
1022  all_y1 = -9999;
1023  all_y2 = -9999;
1024  all_z1 = -9999;
1025  all_z2 = -9999;
1026 
1027  all_row1 = -9999;
1028  all_row2 = -9999;
1029  all_col1 = -9999;
1030  all_col2 = -9999;
1031 
1032  all_gx1 = -9999;
1033  all_gx2 = -9999;
1034  all_gy1 = -9999;
1035  all_gy2 = -9999;
1036  all_gz1 = -9999;
1037  all_gz2 = -9999;
1038 
1039  all_simtrketa = -9999;
1040  all_simtrkphi = -9999;
1041 
1042  all_clust_row = -9999;
1043  all_clust_col = -9999;
1044 
1045  all_clust_x = -9999;
1046  all_clust_y = -9999;
1047 
1048  all_clust_q = -9999;
1049 
1050  all_clust_maxpixcol = -9999;
1051  all_clust_maxpixrow = -9999;
1052  all_clust_minpixcol = -9999;
1053  all_clust_minpixrow = -9999;
1054 
1055  all_clust_geoid = -9999;
1056 
1057  all_clust_alpha = -9999;
1058  all_clust_beta = -9999;
1059 
1060  /*
1061  for (int i=0; i<all_npix; ++i)
1062  {
1063  all_pixrow[i] = -9999;
1064  all_pixcol[i] = -9999;
1065  all_pixadc[i] = -9999;
1066  all_pixx[i] = -9999;
1067  all_pixy[i] = -9999;
1068  all_pixgx[i] = -9999;
1069  all_pixgy[i] = -9999;
1070  all_pixgz[i] = -9999;
1071  }
1072  */
1073 
1074  all_hit_probx = -9999;
1075  all_hit_proby = -9999;
1076  all_hit_cprob0 = -9999;
1077  all_hit_cprob1 = -9999;
1078  all_hit_cprob2 = -9999;
1079 
1080  all_pixel_split = -9999;
1081  all_pixel_clst_err_x = -9999.9;
1082  all_pixel_clst_err_y = -9999.9;
1083 
1084  all_nsimhit = (int)matched.size();
1085 
1086  all_subdetid = (int)detId.subdetId();
1087  // only consider rechits in pixel barrel and pixel forward
1088  if ( !(all_subdetid==1 || all_subdetid==2) )
1089  {
1090  cout << "SiPixelErrorEstimation::analyze: Not in a pixel detector !!!!!" << endl;
1091  continue;
1092  }
1093 
1094  const PixelGeomDetUnit* theGeomDet
1095  = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(detId) );
1096 
1097  const PixelTopology* topol = &(theGeomDet->specificTopology());
1098 
1099  if ( pixeliter->cluster()->getSplitClusterErrorX() > 0.0 &&
1100  pixeliter->cluster()->getSplitClusterErrorY() > 0.0 )
1101  {
1102  all_pixel_split = 1;
1103  }
1104  else
1105  {
1106  all_pixel_split = 0;
1107  }
1108 
1109  all_pixel_clst_err_x = pixeliter->cluster()->getSplitClusterErrorX();
1110  all_pixel_clst_err_y = pixeliter->cluster()->getSplitClusterErrorY();
1111 
1112 
1113  const int maxPixelCol = pixeliter->cluster()->maxPixelCol();
1114  const int maxPixelRow = pixeliter->cluster()->maxPixelRow();
1115  const int minPixelCol = pixeliter->cluster()->minPixelCol();
1116  const int minPixelRow = pixeliter->cluster()->minPixelRow();
1117 
1118  //all_hit_probx = (float)pixeliter->probabilityX();
1119  //all_hit_proby = (float)pixeliter->probabilityY();
1120  all_hit_cprob0 = (float)pixeliter->clusterProbability(0);
1121  all_hit_cprob1 = (float)pixeliter->clusterProbability(1);
1122  all_hit_cprob2 = (float)pixeliter->clusterProbability(2);
1123 
1124  // check whether the cluster is at the module edge
1125  if ( topol->isItEdgePixelInX( minPixelRow ) ||
1126  topol->isItEdgePixelInX( maxPixelRow ) )
1127  all_edgex = 1;
1128  else
1129  all_edgex = 0;
1130 
1131  if ( topol->isItEdgePixelInY( minPixelCol ) ||
1132  topol->isItEdgePixelInY( maxPixelCol ) )
1133  all_edgey = 1;
1134  else
1135  all_edgey = 0;
1136 
1137  // check whether this rechit contains big pixels
1138  if ( topol->containsBigPixelInX(minPixelRow, maxPixelRow) )
1139  all_bigx = 1;
1140  else
1141  all_bigx = 0;
1142 
1143  if ( topol->containsBigPixelInY(minPixelCol, maxPixelCol) )
1144  all_bigy = 1;
1145  else
1146  all_bigy = 0;
1147 
1148  if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel )
1149  {
1150 
1151  all_layer = tTopo->pxbLayer(detId);
1152  all_ladder = tTopo->pxbLadder(detId);
1153  all_mod = tTopo->pxbModule(detId);
1154 
1155  int tmp_nrows = theGeomDet->specificTopology().nrows();
1156  if ( tmp_nrows == 80 )
1157  all_half = 1;
1158  else if ( tmp_nrows == 160 )
1159  all_half = 0;
1160  else
1161  cout << "-------------------------------------------------- Wrong module size !!!" << endl;
1162 
1163  float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
1164  float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
1165 
1166  if ( tmp2<tmp1 )
1167  all_flipped = 1;
1168  else
1169  all_flipped = 0;
1170  }
1171  else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
1172  {
1173 
1174  all_side = tTopo->pxfSide(detId);
1175  all_disk = tTopo->pxfDisk(detId);
1176  all_blade = tTopo->pxfBlade(detId);
1177  all_panel = tTopo->pxfPanel(detId);
1178  all_plaq = tTopo->pxfModule(detId); // also known as plaquette
1179 
1180  } // else if ( detId.subdetId()==PixelSubdetector::PixelEndcap )
1181  else std::cout << "We are not in the pixel detector" << (int)detId.subdetId() << endl;
1182 
1183  all_cols = theGeomDet->specificTopology().ncolumns();
1184  all_rows = theGeomDet->specificTopology().nrows();
1185 
1186  LocalPoint lp = pixeliter->localPosition();
1187  // gavril: change this name
1188  all_rechitx = lp.x();
1189  all_rechity = lp.y();
1190  all_rechitz = lp.z();
1191 
1192  LocalError le = pixeliter->localPositionError();
1193  all_rechiterrx = sqrt( le.xx() );
1194  all_rechiterry = sqrt( le.yy() );
1195 
1196  bool found_hit_from_generated_particle = false;
1197 
1198  //---Loop over sim hits, fill closest
1199  float closest_dist = 99999.9;
1200  std::vector<PSimHit>::const_iterator closest_simhit = matched.begin();
1201 
1202  for (std::vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); m++)
1203  {
1204  if ( checkType_ )
1205  {
1206  int pid = (*m).particleType();
1207  if ( abs(pid) != genType_ )
1208  continue;
1209  }
1210 
1211  float simhitx = 0.5 * ( (*m).entryPoint().x() + (*m).exitPoint().x() );
1212  float simhity = 0.5 * ( (*m).entryPoint().y() + (*m).exitPoint().y() );
1213 
1214  float x_res = simhitx - rechitx;
1215  float y_res = simhity - rechity;
1216 
1217  float dist = sqrt( x_res*x_res + y_res*y_res );
1218 
1219  if ( dist < closest_dist )
1220  {
1221  closest_dist = dist;
1222  closest_simhit = m;
1223  found_hit_from_generated_particle = true;
1224  }
1225  } // end sim hit loop
1226 
1227  // If this recHit does not have any simHit with the same particleType as the particles generated
1228  // ignore it as most probably comes from delta rays.
1229  if ( checkType_ && !found_hit_from_generated_particle )
1230  continue;
1231 
1232  all_x1 = (*closest_simhit).entryPoint().x(); // width (row index, in col direction)
1233  all_y1 = (*closest_simhit).entryPoint().y(); // length (col index, in row direction)
1234  all_z1 = (*closest_simhit).entryPoint().z();
1235  all_x2 = (*closest_simhit).exitPoint().x();
1236  all_y2 = (*closest_simhit).exitPoint().y();
1237  all_z2 = (*closest_simhit).exitPoint().z();
1238  GlobalPoint GP1 =
1239  theGeomDet->surface().toGlobal( Local3DPoint( (*closest_simhit).entryPoint().x(),
1240  (*closest_simhit).entryPoint().y(),
1241  (*closest_simhit).entryPoint().z() ) );
1242  GlobalPoint GP2 =
1243  theGeomDet->surface().toGlobal (Local3DPoint( (*closest_simhit).exitPoint().x(),
1244  (*closest_simhit).exitPoint().y(),
1245  (*closest_simhit).exitPoint().z() ) );
1246  all_gx1 = GP1.x();
1247  all_gx2 = GP2.x();
1248  all_gy1 = GP1.y();
1249  all_gy2 = GP2.y();
1250  all_gz1 = GP1.z();
1251  all_gz2 = GP2.z();
1252 
1253  MeasurementPoint mp1 =
1254  topol->measurementPosition( LocalPoint( (*closest_simhit).entryPoint().x(),
1255  (*closest_simhit).entryPoint().y(),
1256  (*closest_simhit).entryPoint().z() ) );
1257  MeasurementPoint mp2 =
1258  topol->measurementPosition( LocalPoint( (*closest_simhit).exitPoint().x(),
1259  (*closest_simhit).exitPoint().y(),
1260  (*closest_simhit).exitPoint().z() ) );
1261  all_row1 = mp1.x();
1262  all_col1 = mp1.y();
1263  all_row2 = mp2.x();
1264  all_col2 = mp2.y();
1265 
1266  all_simhitx = 0.5*(all_x1+all_x2);
1267  all_simhity = 0.5*(all_y1+all_y2);
1268 
1269  all_rechitresx = all_rechitx - all_simhitx;
1270  all_rechitresy = all_rechity - all_simhity;
1271 
1272  all_rechitpullx = all_rechitresx / all_rechiterrx;
1273  all_rechitpully = all_rechitresy / all_rechiterry;
1274 
1275  SiPixelRecHit::ClusterRef const& clust = pixeliter->cluster();
1276 
1277  all_npix = clust->size();
1278  all_nxpix = clust->sizeX();
1279  all_nypix = clust->sizeY();
1280 
1281  all_clust_row = clust->x();
1282  all_clust_col = clust->y();
1283 
1284  LocalPoint lp2 = topol->localPosition( MeasurementPoint( all_clust_row, all_clust_col ) );
1285  all_clust_x = lp2.x();
1286  all_clust_y = lp2.y();
1287 
1288  all_clust_q = clust->charge();
1289 
1290  all_clust_maxpixcol = clust->maxPixelCol();
1291  all_clust_maxpixrow = clust->maxPixelRow();
1292  all_clust_minpixcol = clust->minPixelCol();
1293  all_clust_minpixrow = clust->minPixelRow();
1294 
1295  all_clust_geoid = 0; // never set!
1296 
1297  all_simpx = (*closest_simhit).momentumAtEntry().x();
1298  all_simpy = (*closest_simhit).momentumAtEntry().y();
1299  all_simpz = (*closest_simhit).momentumAtEntry().z();
1300  all_eloss = (*closest_simhit).energyLoss();
1301  all_simphi = (*closest_simhit).phiAtEntry();
1302  all_simtheta = (*closest_simhit).thetaAtEntry();
1303  all_pidhit = (*closest_simhit).particleType();
1304  all_trkid = (*closest_simhit).trackId();
1305 
1306  //--- Fill alpha and beta -- more useful for exploring the residuals...
1307  all_beta = atan2(all_simpz, all_simpy);
1308  all_alpha = atan2(all_simpz, all_simpx);
1309 
1310  all_simproc = (int)closest_simhit->processType();
1311 
1312  const edm::SimTrackContainer& trks = *(simtracks.product());
1313  SimTrackContainer::const_iterator trksiter;
1314  for (trksiter = trks.begin(); trksiter != trks.end(); trksiter++)
1315  if ( (int)trksiter->trackId() == all_trkid )
1316  {
1317  all_simtrketa = trksiter->momentum().eta();
1318  all_simtrkphi = trksiter->momentum().phi();
1319  }
1320 
1321  all_vtxz = theGeomDet->surface().position().z();
1322  all_vtxr = theGeomDet->surface().position().perp();
1323 
1324  //computeAnglesFromDetPosition(clust,
1325  // theGeomDet,
1326  // all_clust_alpha, all_clust_beta )
1327 
1328  const std::vector<SiPixelCluster::Pixel>& pixvector = clust->pixels();
1329  for ( int i=0; i<(int)pixvector.size(); ++i)
1330  {
1331  SiPixelCluster::Pixel holdpix = pixvector[i];
1332  all_pixrow[i] = holdpix.x;
1333  all_pixcol[i] = holdpix.y;
1334  all_pixadc[i] = holdpix.adc;
1335  LocalPoint lp = topol->localPosition(MeasurementPoint(holdpix.x, holdpix.y));
1336  all_pixx[i] = lp.x();
1337  all_pixy[i]= lp.y();
1338  GlobalPoint GP = theGeomDet->surface().toGlobal(Local3DPoint(lp.x(),lp.y(),lp.z()));
1339  all_pixgx[i] = GP.x();
1340  all_pixgy[i]= GP.y();
1341  all_pixgz[i]= GP.z();
1342  }
1343 
1344  ttree_all_hits_->Fill();
1345 
1346  } // for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter)
1347 
1348  } // for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin(); it != pDD->dets().end(); it++)
1349 
1350  // ------------------------------------------------ all hits ---------------------------------------------------------------
1351 
1352 
1353  //cout << "...3..." << endl;
1354 
1355 
1356  // ------------------------------------------------ track hits only --------------------------------------------------------
1357 
1358 
1359 
1360  if ( include_trk_hits_ )
1361  {
1362  // Get tracks
1363  edm::Handle<reco::TrackCollection> trackCollection;
1364  e.getByToken(tTrackCollection, trackCollection);
1365  const reco::TrackCollection *tracks = trackCollection.product();
1366  reco::TrackCollection::const_iterator tciter;
1367 
1368  if ( tracks->size() > 0 )
1369  {
1370  // Loop on tracks
1371  for ( tciter=tracks->begin(); tciter!=tracks->end(); ++tciter)
1372  {
1373  // First loop on hits: find matched hits
1374  for ( trackingRecHit_iterator it = tciter->recHitsBegin(); it != tciter->recHitsEnd(); ++it)
1375  {
1376  const TrackingRecHit &trk_rec_hit = **it;
1377  // Is it a matched hit?
1378  const SiPixelRecHit* matchedhit = dynamic_cast<const SiPixelRecHit*>(&trk_rec_hit);
1379 
1380  if ( matchedhit )
1381  {
1382  rechitx = -9999.9;
1383  rechity = -9999.9;
1384  rechitz = -9999.9;
1385  rechiterrx = -9999.9;
1386  rechiterry = -9999.9;
1387  rechitresx = -9999.9;
1388  rechitresy = -9999.9;
1389  rechitpullx = -9999.9;
1390  rechitpully = -9999.9;
1391 
1392  npix = -9999;
1393  nxpix = -9999;
1394  nypix = -9999;
1395  charge = -9999.9;
1396 
1397  edgex = -9999;
1398  edgey = -9999;
1399 
1400  bigx = -9999;
1401  bigy = -9999;
1402 
1403  alpha = -9999.9;
1404  beta = -9999.9;
1405 
1406  phi = -9999.9;
1407  eta = -9999.9;
1408 
1409  subdetId = -9999;
1410 
1411  layer = -9999;
1412  ladder = -9999;
1413  mod = -9999;
1414  side = -9999;
1415  disk = -9999;
1416  blade = -9999;
1417  panel = -9999;
1418  plaq = -9999;
1419 
1420  half = -9999;
1421  flipped = -9999;
1422 
1423  nsimhit = -9999;
1424  pidhit = -9999;
1425  simproc = -9999;
1426 
1427  simhitx = -9999.9;
1428  simhity = -9999.9;
1429 
1430  hit_probx = -9999.9;
1431  hit_proby = -9999.9;
1432  hit_cprob0 = -9999.9;
1433  hit_cprob1 = -9999.9;
1434  hit_cprob2 = -9999.9;
1435 
1436  pixel_split = -9999;
1437 
1438  pixel_clst_err_x = -9999.9;
1439  pixel_clst_err_y = -9999.9;
1440 
1441  position = (*it)->localPosition();
1442  error = (*it)->localPositionError();
1443 
1444  rechitx = position.x();
1445  rechity = position.y();
1446  rechitz = position.z();
1447  rechiterrx = sqrt(error.xx());
1448  rechiterry = sqrt(error.yy());
1449 
1450  npix = matchedhit->cluster()->size();
1451  nxpix = matchedhit->cluster()->sizeX();
1452  nypix = matchedhit->cluster()->sizeY();
1453  charge = matchedhit->cluster()->charge();
1454 
1455  if ( matchedhit->cluster()->getSplitClusterErrorX() > 0.0 &&
1456  matchedhit->cluster()->getSplitClusterErrorY() > 0.0 )
1457  pixel_split = 1;
1458  else
1459  pixel_split = 0;
1460 
1461  pixel_clst_err_x = matchedhit->cluster()->getSplitClusterErrorX();
1462  pixel_clst_err_y = matchedhit->cluster()->getSplitClusterErrorY();
1463 
1464  //hit_probx = (float)matchedhit->probabilityX();
1465  //hit_proby = (float)matchedhit->probabilityY();
1466  hit_cprob0 = (float)matchedhit->clusterProbability(0);
1467  hit_cprob1 = (float)matchedhit->clusterProbability(1);
1468  hit_cprob2 = (float)matchedhit->clusterProbability(2);
1469 
1470 
1471  //Association of the rechit to the simhit
1472  matched.clear();
1473  matched = associate.associateHit(*matchedhit);
1474 
1475  nsimhit = (int)matched.size();
1476 
1477  if ( !matched.empty() )
1478  {
1479  mindist = 999999.9;
1480  float distx, disty, dist;
1481  bool found_hit_from_generated_particle = false;
1482 
1483  int n_assoc_muon = 0;
1484 
1485  vector<PSimHit>::const_iterator closestit = matched.begin();
1486  for (vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
1487  {
1488  if ( checkType_ )
1489  { // only consider associated simhits with the generated pid (muons)
1490  int pid = (*m).particleType();
1491  if ( abs(pid) != genType_ )
1492  continue;
1493  }
1494 
1495  float simhitx = 0.5 * ( (*m).entryPoint().x() + (*m).exitPoint().x() );
1496  float simhity = 0.5 * ( (*m).entryPoint().y() + (*m).exitPoint().y() );
1497 
1498  distx = fabs(rechitx - simhitx);
1499  disty = fabs(rechity - simhity);
1500  dist = sqrt( distx*distx + disty*disty );
1501 
1502  if ( dist < mindist )
1503  {
1504  n_assoc_muon++;
1505 
1506  mindist = dist;
1507  closestit = m;
1508  found_hit_from_generated_particle = true;
1509  }
1510  } // for (vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); m++)
1511 
1512  // This recHit does not have any simHit with the same particleType as the particles generated
1513  // Ignore it as most probably come from delta rays.
1514  if ( checkType_ && !found_hit_from_generated_particle )
1515  continue;
1516 
1517  //if ( n_assoc_muon > 1 )
1518  //{
1519  // // cout << " ----- This is not good: n_assoc_muon = " << n_assoc_muon << endl;
1520  // // cout << "evt = " << evt << endl;
1521  //}
1522 
1523  DetId detId = (*it)->geographicalId();
1524 
1525  const PixelGeomDetUnit* theGeomDet =
1526  dynamic_cast<const PixelGeomDetUnit*> ((*tracker).idToDet(detId) );
1527 
1528  const PixelTopology* theTopol = &(theGeomDet->specificTopology());
1529 
1530  pidhit = (*closestit).particleType();
1531  simproc = (int)(*closestit).processType();
1532 
1533  simhitx = 0.5*( (*closestit).entryPoint().x() + (*closestit).exitPoint().x() );
1534  simhity = 0.5*( (*closestit).entryPoint().y() + (*closestit).exitPoint().y() );
1535 
1536  rechitresx = rechitx - simhitx;
1537  rechitresy = rechity - simhity;
1538  rechitpullx = ( rechitx - simhitx ) / sqrt(error.xx());
1539  rechitpully = ( rechity - simhity ) / sqrt(error.yy());
1540 
1541  float simhitpx = (*closestit).momentumAtEntry().x();
1542  float simhitpy = (*closestit).momentumAtEntry().y();
1543  float simhitpz = (*closestit).momentumAtEntry().z();
1544 
1545  beta = atan2(simhitpz, simhitpy) * radtodeg;
1546  alpha = atan2(simhitpz, simhitpx) * radtodeg;
1547 
1548  //beta = fabs(atan2(simhitpz, simhitpy)) * radtodeg;
1549  //alpha = fabs(atan2(simhitpz, simhitpx)) * radtodeg;
1550 
1551  // calculate alpha and beta exactly as in PixelCPEBase.cc
1552  float locx = simhitpx;
1553  float locy = simhitpy;
1554  float locz = simhitpz;
1555 
1556  bool isFlipped = false;
1557  float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
1558  float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
1559  if ( tmp2<tmp1 )
1560  isFlipped = true;
1561  else
1562  isFlipped = false;
1563 
1564  trk_alpha = acos(locx/sqrt(locx*locx+locz*locz)) * radtodeg;
1565  if ( isFlipped ) // &&& check for FPIX !!!
1566  trk_alpha = 180.0 - trk_alpha ;
1567 
1568  trk_beta = acos(locy/sqrt(locy*locy+locz*locz)) * radtodeg;
1569 
1570 
1571  phi = tciter->momentum().phi() / math_pi*180.0;
1572  eta = tciter->momentum().eta();
1573 
1574  const int maxPixelCol = (*matchedhit).cluster()->maxPixelCol();
1575  const int maxPixelRow = (*matchedhit).cluster()->maxPixelRow();
1576  const int minPixelCol = (*matchedhit).cluster()->minPixelCol();
1577  const int minPixelRow = (*matchedhit).cluster()->minPixelRow();
1578 
1579  // check whether the cluster is at the module edge
1580  if ( theTopol->isItEdgePixelInX( minPixelRow ) ||
1581  theTopol->isItEdgePixelInX( maxPixelRow ) )
1582  edgex = 1;
1583  else
1584  edgex = 0;
1585 
1586  if ( theTopol->isItEdgePixelInY( minPixelCol ) ||
1587  theTopol->isItEdgePixelInY( maxPixelCol ) )
1588  edgey = 1;
1589  else
1590  edgey = 0;
1591 
1592  // check whether this rechit contains big pixels
1593  if ( theTopol->containsBigPixelInX(minPixelRow, maxPixelRow) )
1594  bigx = 1;
1595  else
1596  bigx = 0;
1597 
1598  if ( theTopol->containsBigPixelInY(minPixelCol, maxPixelCol) )
1599  bigy = 1;
1600  else
1601  bigy = 0;
1602 
1603  subdetId = (int)detId.subdetId();
1604 
1605  if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel )
1606  {
1607 
1608  int tmp_nrows = theGeomDet->specificTopology().nrows();
1609  if ( tmp_nrows == 80 )
1610  half = 1;
1611  else if ( tmp_nrows == 160 )
1612  half = 0;
1613  else
1614  cout << "-------------------------------------------------- Wrong module size !!!" << endl;
1615 
1616  float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
1617  float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
1618 
1619  if ( tmp2<tmp1 )
1620  flipped = 1;
1621  else
1622  flipped = 0;
1623 
1624 
1625  layer = tTopo->pxbLayer(detId); // Layer: 1,2,3.
1626  ladder = tTopo->pxbLadder(detId); // Ladder: 1-20, 32, 44.
1627  mod = tTopo->pxbModule(detId); // Mod: 1-8.
1628  }
1629  else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
1630  {
1631 
1632  side = tTopo->pxfSide(detId);
1633  disk = tTopo->pxfDisk(detId);
1634  blade = tTopo->pxfBlade(detId);
1635  panel = tTopo->pxfPanel(detId);
1636  plaq = tTopo->pxfModule(detId); // also known as plaquette
1637 
1638  float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
1639  float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
1640 
1641  if ( tmp2<tmp1 )
1642  flipped = 1;
1643  else
1644  flipped = 0;
1645 
1646  } // else if ( detId.subdetId()==PixelSubdetector::PixelEndcap )
1647  //else std::// cout << "We are not in the pixel detector. detId.subdetId() = " << (int)detId.subdetId() << endl;
1648 
1649  ttree_track_hits_->Fill();
1650 
1651  } // if ( !matched.empty() )
1652  else
1653  cout << "---------------- RecHit with no associated SimHit !!! -------------------------- " << endl;
1654 
1655  } // if ( matchedhit )
1656 
1657  } // end of loop on hits
1658 
1659  } //end of loop on track
1660 
1661  } // tracks > 0.
1662 
1663  } // if ( include_trk_hits_ )
1664 
1665 
1666 
1667  // ----------------------------------------------- track hits only -----------------------------------------------------------
1668 
1669 }
1670 
1671 void SiPixelErrorEstimation::
1672 computeAnglesFromDetPosition(const SiPixelCluster & cl,
1673  const GeomDetUnit & det,
1674  float& alpha, float& beta )
1675 {
1676  //--- This is a new det unit, so cache it
1677  const PixelGeomDetUnit* theDet = dynamic_cast<const PixelGeomDetUnit*>( &det );
1678  if (! theDet)
1679  {
1680  cout << "---------------------------------------------- Not a pixel detector !!!!!!!!!!!!!!" << endl;
1681  assert(0);
1682  }
1683 
1684  const PixelTopology* theTopol = &(theDet->specificTopology());
1685 
1686  // get cluster center of gravity (of charge)
1687  float xcenter = cl.x();
1688  float ycenter = cl.y();
1689 
1690  // get the cluster position in local coordinates (cm)
1691  LocalPoint lp = theTopol->localPosition( MeasurementPoint(xcenter, ycenter) );
1692  //float lp_mod = sqrt( lp.x()*lp.x() + lp.y()*lp.y() + lp.z()*lp.z() );
1693 
1694  // get the cluster position in global coordinates (cm)
1695  GlobalPoint gp = theDet->surface().toGlobal( lp );
1696  float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() );
1697 
1698  // normalize
1699  float gpx = gp.x()/gp_mod;
1700  float gpy = gp.y()/gp_mod;
1701  float gpz = gp.z()/gp_mod;
1702 
1703  // make a global vector out of the global point; this vector will point from the
1704  // origin of the detector to the cluster
1705  GlobalVector gv(gpx, gpy, gpz);
1706 
1707  // make local unit vector along local X axis
1708  const Local3DVector lvx(1.0, 0.0, 0.0);
1709 
1710  // get the unit X vector in global coordinates/
1711  GlobalVector gvx = theDet->surface().toGlobal( lvx );
1712 
1713  // make local unit vector along local Y axis
1714  const Local3DVector lvy(0.0, 1.0, 0.0);
1715 
1716  // get the unit Y vector in global coordinates
1717  GlobalVector gvy = theDet->surface().toGlobal( lvy );
1718 
1719  // make local unit vector along local Z axis
1720  const Local3DVector lvz(0.0, 0.0, 1.0);
1721 
1722  // get the unit Z vector in global coordinates
1723  GlobalVector gvz = theDet->surface().toGlobal( lvz );
1724 
1725  // calculate the components of gv (the unit vector pointing to the cluster)
1726  // in the local coordinate system given by the basis {gvx, gvy, gvz}
1727  // note that both gv and the basis {gvx, gvy, gvz} are given in global coordinates
1728  float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z();
1729  float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z();
1730  float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z();
1731 
1732  // calculate angles
1733  alpha = atan2( gv_dot_gvz, gv_dot_gvx );
1734  beta = atan2( gv_dot_gvz, gv_dot_gvy );
1735 
1736  // calculate cotalpha and cotbeta
1737  // cotalpha_ = 1.0/tan(alpha_);
1738  // cotbeta_ = 1.0/tan(beta_ );
1739  // or like this
1740  //cotalpha_ = gv_dot_gvx / gv_dot_gvz;
1741  //cotbeta_ = gv_dot_gvy / gv_dot_gvz;
1742 }
1743 
1744 //define this as a plug-in
1745 DEFINE_FWK_MODULE(SiPixelErrorEstimation);
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:39
SiStripRecHit1D::cluster
ClusterRef cluster() const
Definition: SiStripRecHit1D.h:28
Surface::toGlobal
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:41
edm::ParameterSet::getUntrackedParameter
T getUntrackedParameter(std::string const &, T const &) const
i
int i
Definition: DBlmapReader.cc:9
LocalError::xx
float xx() const
Definition: LocalError.h:24
edmNew::DetSetVector::const_iterator
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
Definition: DetSetVectorNew.h:119
TrackerTopology::tobIsDoubleSide
bool tobIsDoubleSide(const DetId &id) const
Definition: TrackerTopology.h:222
SiPixelRecHit::clusterProbability
float clusterProbability(unsigned int flags=0) const
Definition: SiPixelRecHit.cc:10
PV2DBase::y
T y() const
Definition: PV2DBase.h:46
TrackerTopology::tibIsDoubleSide
bool tibIsDoubleSide(const DetId &id) const
Definition: TrackerTopology.h:224
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:72
TrackerTopology::tibLayer
unsigned int tibLayer(const DetId &id) const
Definition: TrackerTopology.h:121
TrackerTopology::tibString
unsigned int tibString(const DetId &id) const
Definition: TrackerTopology.h:323
Topology::localPosition
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
TrackerTopology::tobIsStereo
bool tobIsStereo(const DetId &id) const
Definition: TrackerTopology.h:240
PixelTopology::ncolumns
virtual int ncolumns() const =0
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:446
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:17
TrackerTopology::pxfDisk
unsigned int pxfDisk(const DetId &id) const
Definition: TrackerTopology.h:359
SiPixelErrorEstimation::tTrackCollection
edm::EDGetTokenT< reco::TrackCollection > tTrackCollection
Definition: SiPixelErrorEstimation.h:84
TrackerTopology::pxbLadder
unsigned int pxbLadder(const DetId &id) const
Definition: TrackerTopology.h:127
reco::TrackCollection
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:13
PV3DBase::y
T y() const
Definition: PV3DBase.h:63
PixelTopology::nrows
virtual int nrows() const =0
NULL
#define NULL
Definition: scimark2.h:8
TrackerTopology::pxbModule
unsigned int pxbModule(const DetId &id) const
Definition: TrackerTopology.h:132
toLocal
LocalVector toLocal(const reco::Track::Vector &v, const Surface &s)
Definition: ConversionProducer.h:207
TrackerTopology::tobIsRPhi
bool tobIsRPhi(const DetId &id) const
Definition: TrackerTopology.h:302
GeomDet::surface
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:40
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
TrackerTopology::tibIsZPlusSide
bool tibIsZPlusSide(const DetId &id) const
Definition: TrackerTopology.h:230
TrackerTopology::tibSide
unsigned int tibSide(const DetId &id) const
Definition: TrackerTopology.h:161
MeasurementPoint
Measurement2DPoint MeasurementPoint
Measurement points are two-dimensional by default.
Definition: MeasurementPoint.h:12
PixelTopology::isItEdgePixelInX
virtual bool isItEdgePixelInX(int ixbin) const =0
TrackerTopology::tibIsExternalString
bool tibIsExternalString(const DetId &id) const
Definition: TrackerTopology.h:339
PV3DBase::mag
T mag() const
Definition: PV3DBase.h:67
TrackerTopology::tibIsRPhi
bool tibIsRPhi(const DetId &id) const
Definition: TrackerTopology.h:304
SiPixelErrorEstimation::tSimTrackContainer
edm::EDGetTokenT< edm::SimTrackContainer > tSimTrackContainer
Definition: SiPixelErrorEstimation.h:83
LocalError::yy
float yy() const
Definition: LocalError.h:26
PixelTopology::containsBigPixelInY
virtual bool containsBigPixelInY(int iymin, int iymax) const =0
TrackerTopology::tibIsZMinusSide
bool tibIsZMinusSide(const DetId &id) const
Definition: TrackerTopology.h:231
PSimHit::localPosition
Local3DPoint localPosition() const
Definition: PSimHit.h:44
TrackingRecHit::ConstRecHitPointer
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
Definition: TrackingRecHit.h:31
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
TrackerTopology::tobIsZPlusSide
bool tobIsZPlusSide(const DetId &id) const
Definition: TrackerTopology.h:227
TrackerTopology::tobSide
unsigned int tobSide(const DetId &id) const
Definition: TrackerTopology.h:153
TrackerTopology::tobIsZMinusSide
bool tobIsZMinusSide(const DetId &id) const
Definition: TrackerTopology.h:228
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
PixelTopology::containsBigPixelInX
virtual bool containsBigPixelInX(int ixmin, int ixmax) const =0
SiPixelErrorEstimation::SiPixelErrorEstimation
SiPixelErrorEstimation(const edm::ParameterSet &)
Definition: SiPixelErrorEstimation.cc:36
patCandidatesForDimuonsSequences_cff.tracker
tuple tracker
Definition: patCandidatesForDimuonsSequences_cff.py:69
SiStripRecHit2D::cluster
ClusterRef cluster() const
Definition: SiStripRecHit2D.h:29
Topology::measurementPosition
virtual MeasurementPoint measurementPosition(const LocalPoint &) const =0
LocalTrajectoryParameters::momentum
LocalVector momentum() const
Momentum vector in the local frame.
Definition: LocalTrajectoryParameters.h:89
DetId::subdetId
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
TrackerTopology::tibModule
unsigned int tibModule(const DetId &id) const
Definition: TrackerTopology.h:144
TrackerTopology::pxfModule
unsigned int pxfModule(const DetId &id) const
Definition: TrackerTopology.h:135
TrackerTopology::pxbLayer
unsigned int pxbLayer(const DetId &id) const
Definition: TrackerTopology.h:115
SiPixelErrorEstimation::tTrajectory
edm::EDGetTokenT< std::vector< Trajectory > > tTrajectory
Definition: SiPixelErrorEstimation.h:81
Local3DPoint
Point3DBase< float, LocalTag > Local3DPoint
Definition: LocalPoint.h:9
SiStripDetId
Detector identifier class for the strip tracker.
Definition: SiStripDetId.h:17
DetId
Definition: DetId.h:18
StripSubdetector::stereo
unsigned int stereo() const
stereo
Definition: StripSubdetector.h:43
TrackingRecHit::hit
virtual TrackingRecHit const * hit() const
Definition: TrackingRecHit.h:73
edm::Handle::product
T const * product() const
Definition: Handle.h:81
SiPixelErrorEstimation::computeAnglesFromDetPosition
void computeAnglesFromDetPosition(const SiPixelCluster &cl, const GeomDetUnit &det, float &alpha, float &beta)
Definition: SiPixelErrorEstimation.cc:1672
testEve_cfg.tracks
tuple tracks
Definition: testEve_cfg.py:39
SiPixelRecHit::cluster
ClusterRef cluster() const
Definition: SiPixelRecHit.h:49
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
TrackerTopology::tibIsStereo
bool tibIsStereo(const DetId &id) const
Definition: TrackerTopology.h:242
sysUtil.pid
tuple pid
Definition: sysUtil.py:22
PixelGeomDetUnit::specificTopology
virtual const PixelTopology & specificTopology() const
Returns a reference to the pixel proxy topology.
Definition: PixelGeomDetUnit.cc:18
edm::ESHandle::product
T const * product() const
Definition: ESHandle.h:86
PSimHit::processType
unsigned short processType() const
Definition: PSimHit.h:118
TrackerTopology::tobModule
unsigned int tobModule(const DetId &id) const
Definition: TrackerTopology.h:138
GetRecoTauVFromDQM_MC_cff.cl
tuple cl
Definition: GetRecoTauVFromDQM_MC_cff.py:37
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
SiPixelCluster
Pixel cluster – collection of neighboring pixels above threshold.
Definition: SiPixelCluster.h:33
PixelTopology::isItEdgePixelInY
virtual bool isItEdgePixelInY(int iybin) const =0
PSimHit::particleType
int particleType() const
Definition: PSimHit.h:85
edmNew::DetSet::end
iterator end()
Definition: DetSetNew.h:70
position
static int position[264][3]
Definition: ReadPGInfo.cc:509
SiPixelCluster::y
float y() const
Definition: SiPixelCluster.h:127
LocalPoint
Local3DPoint LocalPoint
Definition: LocalPoint.h:11
TrackerTopology::pxfSide
unsigned int pxfSide(const DetId &id) const
Definition: TrackerTopology.h:169
TrackerHitAssociator::associateHit
std::vector< PSimHit > associateHit(const TrackingRecHit &thit) const
Definition: TrackerHitAssociator.cc:121
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
SiPixelErrorEstimation::tPixRecHitCollection
edm::EDGetTokenT< SiPixelRecHitCollection > tPixRecHitCollection
Definition: SiPixelErrorEstimation.h:82
x
Definition: DDAxes.h:10
TrackerTopology::tobRod
unsigned int tobRod(const DetId &id) const
Definition: TrackerTopology.h:174
SiStripDetId::moduleGeometry
ModuleGeometry moduleGeometry() const
Definition: SiStripDetId.h:132
SiPixelErrorEstimation::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: SiPixelErrorEstimation.cc:371
TrackerTopology::tibIsInternalString
bool tibIsInternalString(const DetId &id) const
Definition: TrackerTopology.h:335
PV2DBase::x
T x() const
Definition: PV2DBase.h:45
PV3DBase::x
T x() const
Definition: PV3DBase.h:62
GloballyPositioned::position
const PositionType & position() const
Definition: GloballyPositioned.h:42
edm::SimTrackContainer
std::vector< SimTrack > SimTrackContainer
Definition: SimTrackContainer.h:12
findQualityFiles.size
tuple size
Write out results.
Definition: findQualityFiles.py:442
SiPixelCluster::x
float x() const
Definition: SiPixelCluster.h:119
TrackerTopology::pxfPanel
unsigned int pxfPanel(const DetId &id) const
Definition: TrackerTopology.h:365
TrackerTopology::pxfBlade
unsigned int pxfBlade(const DetId &id) const
Definition: TrackerTopology.h:362
TrackerTopology::tobLayer
unsigned int tobLayer(const DetId &id) const
Definition: TrackerTopology.h:118
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:23
TrackerTopology::tibOrder
unsigned int tibOrder(const DetId &id) const
Definition: TrackerTopology.h:190
edmNew::DetSet::begin
iterator begin()
Definition: DetSetNew.h:67
Generated for CMSSW Reference Manual by doxygen 1.8.5