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