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