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