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