00001
00002
00003
00004
00005
00006
00007 #include <iostream>
00008 #include "CalibTracker/SiPixelErrorEstimation/interface/SiPixelErrorEstimation.h"
00009
00010 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00011 #include "DataFormats/GeometryVector/interface/LocalVector.h"
00012
00013 #include "DataFormats/TrackReco/interface/Track.h"
00014 #include "DataFormats/DetId/interface/DetId.h"
00015 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
00016
00017 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00018 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00019 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00020 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00021 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00022
00023 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
00024
00025 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00026
00027 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
00028
00029 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
00030 #include "FWCore/ParameterSet/interface/FileInPath.h"
00031
00032 using namespace std;
00033 using namespace edm;
00034
00035 SiPixelErrorEstimation::SiPixelErrorEstimation(const edm::ParameterSet& ps):tfile_(0), ttree_all_hits_(0),
00036 ttree_track_hits_(0), ttree_track_hits_strip_(0)
00037 {
00038
00039 outputFile_ = ps.getUntrackedParameter<string>( "outputFile", "SiPixelErrorEstimation_Ntuple.root" );
00040
00041
00042
00043 src_ = ps.getUntrackedParameter<std::string>( "src", "generalTracks" );
00044
00045 checkType_ = ps.getParameter<bool>( "checkType" );
00046 genType_ = ps.getParameter<int>( "genType" );
00047 include_trk_hits_ = ps.getParameter<bool>( "include_trk_hits" );
00048 }
00049
00050 SiPixelErrorEstimation::~SiPixelErrorEstimation()
00051 {}
00052
00053 void SiPixelErrorEstimation::beginJob()
00054 {
00055 int bufsize = 64000;
00056
00057
00058 tfile_ = new TFile ( outputFile_.c_str() , "RECREATE");
00059
00060
00061 ttree_track_hits_strip_ = new TTree("TrackHitNtupleStrip", "TrackHitNtupleStrip");
00062
00063 ttree_track_hits_strip_->Branch("strip_rechitx", &strip_rechitx, "strip_rechitx/F" , bufsize);
00064 ttree_track_hits_strip_->Branch("strip_rechity", &strip_rechity, "strip_rechity/F" , bufsize);
00065 ttree_track_hits_strip_->Branch("strip_rechitz", &strip_rechitz, "strip_rechitz/F" , bufsize);
00066
00067 ttree_track_hits_strip_->Branch("strip_rechiterrx", &strip_rechiterrx, "strip_rechiterrx/F" , bufsize);
00068 ttree_track_hits_strip_->Branch("strip_rechiterry", &strip_rechiterry, "strip_rechiterry/F" , bufsize);
00069
00070 ttree_track_hits_strip_->Branch("strip_rechitresx", &strip_rechitresx, "strip_rechitresx/F" , bufsize);
00071
00072 ttree_track_hits_strip_->Branch("strip_rechitresx2", &strip_rechitresx2, "strip_rechitresx2/F" , bufsize);
00073
00074
00075 ttree_track_hits_strip_->Branch("strip_rechitresy", &strip_rechitresy, "strip_rechitresy/F" , bufsize);
00076
00077 ttree_track_hits_strip_->Branch("strip_rechitpullx", &strip_rechitpullx, "strip_rechitpullx/F", bufsize);
00078 ttree_track_hits_strip_->Branch("strip_rechitpully", &strip_rechitpully, "strip_rechitpully/F", bufsize);
00079
00080 ttree_track_hits_strip_->Branch("strip_is_stereo", &strip_is_stereo, "strip_is_stereo/I", bufsize);
00081 ttree_track_hits_strip_->Branch("strip_hit_type" , &strip_hit_type , "strip_hit_type/I" , bufsize);
00082 ttree_track_hits_strip_->Branch("detector_type" , &detector_type , "detector_type/I" , bufsize);
00083
00084 ttree_track_hits_strip_->Branch("strip_trk_pt" , &strip_trk_pt , "strip_trk_pt/F" , bufsize);
00085 ttree_track_hits_strip_->Branch("strip_cotalpha" , &strip_cotalpha , "strip_cotalpha/F" , bufsize);
00086 ttree_track_hits_strip_->Branch("strip_cotbeta" , &strip_cotbeta , "strip_cotbeta/F" , bufsize);
00087 ttree_track_hits_strip_->Branch("strip_locbx" , &strip_locbx , "strip_locbx/F" , bufsize);
00088 ttree_track_hits_strip_->Branch("strip_locby" , &strip_locby , "strip_locby/F" , bufsize);
00089 ttree_track_hits_strip_->Branch("strip_locbz" , &strip_locbz , "strip_locbz/F" , bufsize);
00090 ttree_track_hits_strip_->Branch("strip_charge" , &strip_charge , "strip_charge/F" , bufsize);
00091 ttree_track_hits_strip_->Branch("strip_size" , &strip_size , "strip_size/I" , bufsize);
00092
00093
00094 ttree_track_hits_strip_->Branch("strip_edge" , &strip_edge , "strip_edge/I" , bufsize);
00095 ttree_track_hits_strip_->Branch("strip_nsimhit", &strip_nsimhit, "strip_nsimhit/I" , bufsize);
00096 ttree_track_hits_strip_->Branch("strip_pidhit" , &strip_pidhit , "strip_pidhit/I" , bufsize);
00097 ttree_track_hits_strip_->Branch("strip_simproc", &strip_simproc, "strip_simproc/I" , bufsize);
00098
00099
00100 ttree_track_hits_strip_->Branch("strip_subdet_id" , &strip_subdet_id , "strip_subdet_id/I" , bufsize);
00101
00102 ttree_track_hits_strip_->Branch("strip_tib_layer" , &strip_tib_layer , "strip_tib_layer/I" , bufsize);
00103 ttree_track_hits_strip_->Branch("strip_tib_module" , &strip_tib_module , "strip_tib_module/I" , bufsize);
00104 ttree_track_hits_strip_->Branch("strip_tib_order" , &strip_tib_order , "strip_tib_order/I" , bufsize);
00105 ttree_track_hits_strip_->Branch("strip_tib_side" , &strip_tib_side , "strip_tib_side/I" , bufsize);
00106 ttree_track_hits_strip_->Branch("strip_tib_is_double_side" , &strip_tib_is_double_side , "strip_tib_is_double_side/I" , bufsize);
00107 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);
00108 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);
00109 ttree_track_hits_strip_->Branch("strip_tib_layer_number" , &strip_tib_layer_number , "strip_tib_layer_number/I" , bufsize);
00110 ttree_track_hits_strip_->Branch("strip_tib_string_number" , &strip_tib_string_number , "strip_tib_string_number/I" , bufsize);
00111 ttree_track_hits_strip_->Branch("strip_tib_module_number" , &strip_tib_module_number ,"strip_tib_module_number/I" , bufsize);
00112 ttree_track_hits_strip_->Branch("strip_tib_is_internal_string", &strip_tib_is_internal_string,"strip_tib_is_internal_string/I", bufsize);
00113 ttree_track_hits_strip_->Branch("strip_tib_is_external_string", &strip_tib_is_external_string,"strip_tib_is_external_string/I", bufsize);
00114 ttree_track_hits_strip_->Branch("strip_tib_is_rphi" , &strip_tib_is_rphi , "strip_tib_is_rphi/I" , bufsize);
00115 ttree_track_hits_strip_->Branch("strip_tib_is_stereo" , &strip_tib_is_stereo , "strip_tib_is_stereo/I" , bufsize);
00116 ttree_track_hits_strip_->Branch("strip_tob_layer" , &strip_tob_layer , "strip_tob_layer/I" , bufsize);
00117 ttree_track_hits_strip_->Branch("strip_tob_module" , &strip_tob_module , "strip_tob_module/I" , bufsize);
00118 ttree_track_hits_strip_->Branch("strip_tob_side" , &strip_tob_side , "strip_tob_side/I" , bufsize);
00119 ttree_track_hits_strip_->Branch("strip_tob_is_double_side" , &strip_tob_is_double_side , "strip_tob_is_double_side/I" , bufsize);
00120 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);
00121 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);
00122 ttree_track_hits_strip_->Branch("strip_tob_layer_number" , &strip_tob_layer_number , "strip_tob_layer_number/I" , bufsize);
00123 ttree_track_hits_strip_->Branch("strip_tob_rod_number" , &strip_tob_rod_number , "strip_tob_rod_number/I" , bufsize);
00124 ttree_track_hits_strip_->Branch("strip_tob_module_number" , &strip_tob_module_number , "strip_tob_module_number/I" , bufsize);
00125
00126
00127 ttree_track_hits_strip_->Branch("strip_prob", &strip_prob, "strip_prob/F" , bufsize);
00128 ttree_track_hits_strip_->Branch("strip_qbin", &strip_qbin, "strip_qbin/I", bufsize);
00129
00130 ttree_track_hits_strip_->Branch("strip_nprm", &strip_nprm, "strip_nprm/I", bufsize);
00131
00132 ttree_track_hits_strip_->Branch("strip_pidhit1" , &strip_pidhit1 , "strip_pidhit1/I" , bufsize);
00133 ttree_track_hits_strip_->Branch("strip_simproc1", &strip_simproc1, "strip_simproc1/I" , bufsize);
00134
00135 ttree_track_hits_strip_->Branch("strip_pidhit2" , &strip_pidhit2 , "strip_pidhit2/I" , bufsize);
00136 ttree_track_hits_strip_->Branch("strip_simproc2", &strip_simproc2, "strip_simproc2/I" , bufsize);
00137
00138 ttree_track_hits_strip_->Branch("strip_pidhit3" , &strip_pidhit3 , "strip_pidhit3/I" , bufsize);
00139 ttree_track_hits_strip_->Branch("strip_simproc3", &strip_simproc3, "strip_simproc3/I" , bufsize);
00140
00141 ttree_track_hits_strip_->Branch("strip_pidhit4" , &strip_pidhit4 , "strip_pidhit4/I" , bufsize);
00142 ttree_track_hits_strip_->Branch("strip_simproc4", &strip_simproc4, "strip_simproc4/I" , bufsize);
00143
00144 ttree_track_hits_strip_->Branch("strip_pidhit5" , &strip_pidhit5 , "strip_pidhit5/I" , bufsize);
00145 ttree_track_hits_strip_->Branch("strip_simproc5", &strip_simproc5, "strip_simproc5/I" , bufsize);
00146
00147 ttree_track_hits_strip_->Branch("strip_split", &strip_split, "strip_split/I" , bufsize);
00148
00149 ttree_track_hits_strip_->Branch("strip_clst_err_x", &strip_clst_err_x, "strip_clst_err_x/F" , bufsize);
00150 ttree_track_hits_strip_->Branch("strip_clst_err_y", &strip_clst_err_y, "strip_clst_err_y/F" , bufsize);
00151
00152 if ( include_trk_hits_ )
00153 {
00154
00155
00156
00157
00158 ttree_track_hits_ = new TTree("TrackHitNtuple", "TrackHitNtuple");
00159
00160 ttree_track_hits_->Branch("evt", &evt, "evt/I", bufsize);
00161 ttree_track_hits_->Branch("run", &run, "run/I", bufsize);
00162
00163 ttree_track_hits_->Branch("subdetId", &subdetId, "subdetId/I", bufsize);
00164
00165 ttree_track_hits_->Branch("layer" , &layer , "layer/I" , bufsize);
00166 ttree_track_hits_->Branch("ladder", &ladder, "ladder/I", bufsize);
00167 ttree_track_hits_->Branch("mod" , &mod , "mod/I" , bufsize);
00168
00169 ttree_track_hits_->Branch("side" , &side , "side/I" , bufsize);
00170 ttree_track_hits_->Branch("disk" , &disk , "disk/I" , bufsize);
00171 ttree_track_hits_->Branch("blade" , &blade , "blade/I" , bufsize);
00172 ttree_track_hits_->Branch("panel" , &panel , "panel/I" , bufsize);
00173 ttree_track_hits_->Branch("plaq" , &plaq , "plaq/I" , bufsize);
00174
00175 ttree_track_hits_->Branch("half" , &half , "half/I" , bufsize);
00176 ttree_track_hits_->Branch("flipped", &flipped, "flipped/I", bufsize);
00177
00178 ttree_track_hits_->Branch("rechitx", &rechitx, "rechitx/F" , bufsize);
00179 ttree_track_hits_->Branch("rechity", &rechity, "rechity/F" , bufsize);
00180 ttree_track_hits_->Branch("rechitz", &rechitz, "rechitz/F" , bufsize);
00181
00182 ttree_track_hits_->Branch("rechiterrx", &rechiterrx, "rechiterrx/F" , bufsize);
00183 ttree_track_hits_->Branch("rechiterry", &rechiterry, "rechiterry/F" , bufsize);
00184
00185 ttree_track_hits_->Branch("rechitresx", &rechitresx, "rechitresx/F" , bufsize);
00186 ttree_track_hits_->Branch("rechitresy", &rechitresy, "rechitresy/F" , bufsize);
00187
00188 ttree_track_hits_->Branch("rechitpullx", &rechitpullx, "rechitpullx/F", bufsize);
00189 ttree_track_hits_->Branch("rechitpully", &rechitpully, "rechitpully/F", bufsize);
00190
00191
00192 ttree_track_hits_->Branch("npix" , &npix , "npix/I" , bufsize);
00193 ttree_track_hits_->Branch("nxpix" , &nxpix , "nxpix/I" , bufsize);
00194 ttree_track_hits_->Branch("nypix" , &nypix , "nypix/I" , bufsize);
00195 ttree_track_hits_->Branch("charge", &charge, "charge/F", bufsize);
00196
00197 ttree_track_hits_->Branch("edgex", &edgex, "edgex/I", bufsize);
00198 ttree_track_hits_->Branch("edgey", &edgey, "edgey/I", bufsize);
00199
00200 ttree_track_hits_->Branch("bigx", &bigx, "bigx/I", bufsize);
00201 ttree_track_hits_->Branch("bigy", &bigy, "bigy/I", bufsize);
00202
00203 ttree_track_hits_->Branch("alpha", &alpha, "alpha/F", bufsize);
00204 ttree_track_hits_->Branch("beta" , &beta , "beta/F" , bufsize);
00205
00206 ttree_track_hits_->Branch("trk_alpha", &trk_alpha, "trk_alpha/F", bufsize);
00207 ttree_track_hits_->Branch("trk_beta" , &trk_beta , "trk_beta/F" , bufsize);
00208
00209 ttree_track_hits_->Branch("phi", &phi, "phi/F", bufsize);
00210 ttree_track_hits_->Branch("eta", &eta, "eta/F", bufsize);
00211
00212 ttree_track_hits_->Branch("simhitx", &simhitx, "simhitx/F", bufsize);
00213 ttree_track_hits_->Branch("simhity", &simhity, "simhity/F", bufsize);
00214
00215 ttree_track_hits_->Branch("nsimhit", &nsimhit, "nsimhit/I", bufsize);
00216 ttree_track_hits_->Branch("pidhit" , &pidhit , "pidhit/I" , bufsize);
00217 ttree_track_hits_->Branch("simproc", &simproc, "simproc/I", bufsize);
00218
00219 ttree_track_hits_->Branch("pixel_split", &pixel_split, "pixel_split/I", bufsize);
00220
00221 ttree_track_hits_->Branch("pixel_clst_err_x", &pixel_clst_err_x, "pixel_clst_err_x/F" , bufsize);
00222 ttree_track_hits_->Branch("pixel_clst_err_y", &pixel_clst_err_y, "pixel_clst_err_y/F" , bufsize);
00223
00224 }
00225
00226
00227
00228 ttree_all_hits_ = new TTree("AllHitNtuple", "AllHitNtuple");
00229
00230 ttree_all_hits_->Branch("evt", &evt, "evt/I", bufsize);
00231 ttree_all_hits_->Branch("run", &run, "run/I", bufsize);
00232
00233 ttree_all_hits_->Branch("subdetid", &all_subdetid, "subdetid/I", bufsize);
00234
00235 ttree_all_hits_->Branch("layer" , &all_layer , "layer/I" , bufsize);
00236 ttree_all_hits_->Branch("ladder", &all_ladder, "ladder/I", bufsize);
00237 ttree_all_hits_->Branch("mod" , &all_mod , "mod/I" , bufsize);
00238
00239 ttree_all_hits_->Branch("side" , &all_side , "side/I" , bufsize);
00240 ttree_all_hits_->Branch("disk" , &all_disk , "disk/I" , bufsize);
00241 ttree_all_hits_->Branch("blade" , &all_blade , "blade/I" , bufsize);
00242 ttree_all_hits_->Branch("panel" , &all_panel , "panel/I" , bufsize);
00243 ttree_all_hits_->Branch("plaq" , &all_plaq , "plaq/I" , bufsize);
00244
00245 ttree_all_hits_->Branch("half" , &all_half , "half/I" , bufsize);
00246 ttree_all_hits_->Branch("flipped", &all_flipped, "flipped/I", bufsize);
00247
00248 ttree_all_hits_->Branch("cols", &all_cols, "cols/I", bufsize);
00249 ttree_all_hits_->Branch("rows", &all_rows, "rows/I", bufsize);
00250
00251 ttree_all_hits_->Branch("rechitx" , &all_rechitx , "rechitx/F" , bufsize);
00252 ttree_all_hits_->Branch("rechity" , &all_rechity , "rechity/F" , bufsize);
00253 ttree_all_hits_->Branch("rechitz" , &all_rechitz , "rechitz/F" , bufsize);
00254
00255 ttree_all_hits_->Branch("rechiterrx" , &all_rechiterrx , "rechiterrx/F" , bufsize);
00256 ttree_all_hits_->Branch("rechiterry" , &all_rechiterry , "rechiterry/F" , bufsize);
00257
00258 ttree_all_hits_->Branch("rechitresx" , &all_rechitresx , "rechitresx/F" , bufsize);
00259 ttree_all_hits_->Branch("rechitresy" , &all_rechitresy , "rechitresy/F" , bufsize);
00260
00261 ttree_all_hits_->Branch("rechitpullx", &all_rechitpullx, "rechitpullx/F", bufsize);
00262 ttree_all_hits_->Branch("rechitpully", &all_rechitpully, "rechitpully/F", bufsize);
00263
00264 ttree_all_hits_->Branch("npix" , &all_npix , "npix/I" , bufsize);
00265 ttree_all_hits_->Branch("nxpix" , &all_nxpix , "nxpix/I" , bufsize);
00266 ttree_all_hits_->Branch("nypix" , &all_nypix , "nypix/I" , bufsize);
00267
00268 ttree_all_hits_->Branch("edgex", &all_edgex, "edgex/I", bufsize);
00269 ttree_all_hits_->Branch("edgey", &all_edgey, "edgey/I", bufsize);
00270
00271 ttree_all_hits_->Branch("bigx", &all_bigx, "bigx/I", bufsize);
00272 ttree_all_hits_->Branch("bigy", &all_bigy, "bigy/I", bufsize);
00273
00274 ttree_all_hits_->Branch("alpha", &all_alpha, "alpha/F", bufsize);
00275 ttree_all_hits_->Branch("beta" , &all_beta , "beta/F" , bufsize);
00276
00277 ttree_all_hits_->Branch("simhitx", &all_simhitx, "simhitx/F", bufsize);
00278 ttree_all_hits_->Branch("simhity", &all_simhity, "simhity/F", bufsize);
00279
00280 ttree_all_hits_->Branch("nsimhit", &all_nsimhit, "nsimhit/I", bufsize);
00281 ttree_all_hits_->Branch("pidhit" , &all_pidhit , "pidhit/I" , bufsize);
00282 ttree_all_hits_->Branch("simproc", &all_simproc, "simproc/I", bufsize);
00283
00284 ttree_all_hits_->Branch("vtxr", &all_vtxr, "vtxr/F", bufsize);
00285 ttree_all_hits_->Branch("vtxz", &all_vtxz, "vtxz/F", bufsize);
00286
00287 ttree_all_hits_->Branch("simpx", &all_simpx, "simpx/F", bufsize);
00288 ttree_all_hits_->Branch("simpy", &all_simpy, "simpy/F", bufsize);
00289 ttree_all_hits_->Branch("simpz", &all_simpz, "simpz/F", bufsize);
00290
00291 ttree_all_hits_->Branch("eloss", &all_eloss, "eloss/F", bufsize);
00292
00293 ttree_all_hits_->Branch("simphi", &all_simphi, "simphi/F", bufsize);
00294 ttree_all_hits_->Branch("simtheta", &all_simtheta, "simtheta/F", bufsize);
00295
00296 ttree_all_hits_->Branch("trkid", &all_trkid, "trkid/I", bufsize);
00297
00298 ttree_all_hits_->Branch("x1", &all_x1, "x1/F", bufsize);
00299 ttree_all_hits_->Branch("x2", &all_x2, "x2/F", bufsize);
00300 ttree_all_hits_->Branch("y1", &all_x1, "y1/F", bufsize);
00301 ttree_all_hits_->Branch("y2", &all_x2, "y2/F", bufsize);
00302 ttree_all_hits_->Branch("z1", &all_x1, "z1/F", bufsize);
00303 ttree_all_hits_->Branch("z2", &all_x2, "z2/F", bufsize);
00304
00305 ttree_all_hits_->Branch("row1", &all_row1, "row1/F", bufsize);
00306 ttree_all_hits_->Branch("row2", &all_row2, "row2/F", bufsize);
00307 ttree_all_hits_->Branch("col1", &all_col1, "col1/F", bufsize);
00308 ttree_all_hits_->Branch("col2", &all_col2, "col2/F", bufsize);
00309
00310 ttree_all_hits_->Branch("gx1", &all_gx1, "gx1/F", bufsize);
00311 ttree_all_hits_->Branch("gx2", &all_gx2, "gx2/F", bufsize);
00312 ttree_all_hits_->Branch("gy1", &all_gx1, "gy1/F", bufsize);
00313 ttree_all_hits_->Branch("gy2", &all_gx2, "gy2/F", bufsize);
00314 ttree_all_hits_->Branch("gz1", &all_gx1, "gz1/F", bufsize);
00315 ttree_all_hits_->Branch("gz2", &all_gx2, "gz2/F", bufsize);
00316
00317 ttree_all_hits_->Branch("simtrketa", &all_simtrketa, "simtrketa/F", bufsize);
00318 ttree_all_hits_->Branch("simtrkphi", &all_simtrkphi, "simtrkphi/F", bufsize);
00319
00320 ttree_all_hits_->Branch("clust_row", &all_clust_row, "clust_row/F", bufsize);
00321 ttree_all_hits_->Branch("clust_col", &all_clust_col, "clust_col/F", bufsize);
00322
00323 ttree_all_hits_->Branch("clust_x", &all_clust_x, "clust_x/F", bufsize);
00324 ttree_all_hits_->Branch("clust_y", &all_clust_y, "clust_y/F", bufsize);
00325
00326 ttree_all_hits_->Branch("clust_q", &all_clust_q, "clust_q/F", bufsize);
00327
00328 ttree_all_hits_->Branch("clust_maxpixcol", &all_clust_maxpixcol, "clust_maxpixcol/I", bufsize);
00329 ttree_all_hits_->Branch("clust_maxpixrow", &all_clust_maxpixrow, "clust_maxpixrow/I", bufsize);
00330 ttree_all_hits_->Branch("clust_minpixcol", &all_clust_minpixcol, "clust_minpixcol/I", bufsize);
00331 ttree_all_hits_->Branch("clust_minpixrow", &all_clust_minpixrow, "clust_minpixrow/I", bufsize);
00332
00333 ttree_all_hits_->Branch("clust_geoid", &all_clust_geoid, "clust_geoid/I", bufsize);
00334
00335 ttree_all_hits_->Branch("clust_alpha", &all_clust_alpha, "clust_alpha/F", bufsize);
00336 ttree_all_hits_->Branch("clust_beta" , &all_clust_beta , "clust_beta/F" , bufsize);
00337
00338 ttree_all_hits_->Branch("rowpix", all_pixrow, "row[npix]/F", bufsize);
00339 ttree_all_hits_->Branch("colpix", all_pixcol, "col[npix]/F", bufsize);
00340 ttree_all_hits_->Branch("adc", all_pixadc, "adc[npix]/F", bufsize);
00341 ttree_all_hits_->Branch("xpix", all_pixx, "x[npix]/F", bufsize);
00342 ttree_all_hits_->Branch("ypix", all_pixy, "y[npix]/F", bufsize);
00343 ttree_all_hits_->Branch("gxpix", all_pixgx, "gx[npix]/F", bufsize);
00344 ttree_all_hits_->Branch("gypix", all_pixgy, "gy[npix]/F", bufsize);
00345 ttree_all_hits_->Branch("gzpix", all_pixgz, "gz[npix]/F", bufsize);
00346
00347 ttree_all_hits_->Branch("hit_probx", &all_hit_probx, "hit_probx/F" , bufsize);
00348 ttree_all_hits_->Branch("hit_proby", &all_hit_proby, "hit_proby/F" , bufsize);
00349
00350 ttree_all_hits_->Branch("all_pixel_split", &all_pixel_split, "all_pixel_split/I" , bufsize);
00351
00352 ttree_all_hits_->Branch("all_pixel_clst_err_x", &all_pixel_clst_err_x, "all_pixel_clst_err_x/F" , bufsize);
00353 ttree_all_hits_->Branch("all_pixel_clst_err_y", &all_pixel_clst_err_y, "all_pixel_clst_err_y/F" , bufsize);
00354
00355 }
00356
00357 void SiPixelErrorEstimation::endJob()
00358 {
00359 tfile_->Write();
00360 tfile_->Close();
00361 }
00362
00363 void
00364 SiPixelErrorEstimation::analyze(const edm::Event& e, const edm::EventSetup& es)
00365 {
00366 using namespace edm;
00367
00368 run = e.id().run();
00369 evt = e.id().event();
00370
00371 if ( evt%1000 == 0 )
00372 cout << "evt = " << evt << endl;
00373
00374 float math_pi = 3.14159265;
00375 float radtodeg = 180.0 / math_pi;
00376
00377 LocalPoint position;
00378 LocalError error;
00379 float mindist = 999999.9;
00380
00381 std::vector<PSimHit> matched;
00382 TrackerHitAssociator associate(e);
00383
00384 edm::ESHandle<TrackerGeometry> pDD;
00385 es.get<TrackerDigiGeometryRecord> ().get (pDD);
00386 const TrackerGeometry* tracker = &(* pDD);
00387
00388
00389
00390
00391
00392 edm::ESHandle<MagneticField> magneticField;
00393 es.get<IdealMagneticFieldRecord>().get( magneticField );
00394
00395
00396 edm::FileInPath FileInPath_;
00397
00398
00399
00400
00401
00402
00403 edm::Handle<vector<Trajectory> > trajCollectionHandle;
00404
00405 e.getByLabel( src_, trajCollectionHandle);
00406
00407
00408 for ( vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end(); ++it )
00409 {
00410
00411 vector<TrajectoryMeasurement> tmColl = it->measurements();
00412 for ( vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj!=tmColl.end(); ++itTraj )
00413 {
00414
00415 if ( !itTraj->updatedState().isValid() )
00416 continue;
00417
00418 strip_rechitx = -9999999.9;
00419 strip_rechity = -9999999.9;
00420 strip_rechitz = -9999999.9;
00421 strip_rechiterrx = -9999999.9;
00422 strip_rechiterry = -9999999.9;
00423 strip_rechitresx = -9999999.9;
00424 strip_rechitresx2 = -9999999.9;
00425
00426
00427 strip_rechitresy = -9999999.9;
00428 strip_rechitpullx = -9999999.9;
00429 strip_rechitpully = -9999999.9;
00430 strip_is_stereo = -9999999 ;
00431 strip_hit_type = -9999999 ;
00432 detector_type = -9999999 ;
00433
00434 strip_trk_pt = -9999999.9;
00435 strip_cotalpha = -9999999.9;
00436 strip_cotbeta = -9999999.9;
00437 strip_locbx = -9999999.9;
00438 strip_locby = -9999999.9;
00439 strip_locbz = -9999999.9;
00440 strip_charge = -9999999.9;
00441 strip_size = -9999999 ;
00442
00443 strip_edge = -9999999 ;
00444 strip_nsimhit = -9999999 ;
00445 strip_pidhit = -9999999 ;
00446 strip_simproc = -9999999 ;
00447
00448
00449 strip_subdet_id = -9999999;
00450
00451 strip_tib_layer = -9999999;
00452 strip_tib_module = -9999999;
00453 strip_tib_order = -9999999;
00454 strip_tib_side = -9999999;
00455 strip_tib_is_double_side = -9999999;
00456 strip_tib_is_z_plus_side = -9999999;
00457 strip_tib_is_z_minus_side = -9999999;
00458 strip_tib_layer_number = -9999999;
00459 strip_tib_string_number = -9999999;
00460 strip_tib_module_number = -9999999;
00461 strip_tib_is_internal_string= -9999999;
00462 strip_tib_is_external_string= -9999999;
00463 strip_tib_is_rphi = -9999999;
00464 strip_tib_is_stereo = -9999999;
00465
00466 strip_tob_layer = -9999999;
00467 strip_tob_module = -9999999;
00468 strip_tob_side = -9999999;
00469 strip_tob_is_double_side = -9999999;
00470 strip_tob_is_z_plus_side = -9999999;
00471 strip_tob_is_z_minus_side = -9999999;
00472 strip_tob_layer_number = -9999999;
00473 strip_tob_rod_number = -9999999;
00474 strip_tob_module_number = -9999999;
00475
00476 strip_tob_is_rphi = -9999999;
00477 strip_tob_is_stereo = -9999999;
00478
00479 strip_prob = -9999999.9;
00480 strip_qbin = -9999999;
00481
00482 strip_nprm = -9999999;
00483
00484 strip_pidhit1 = -9999999 ;
00485 strip_simproc1 = -9999999 ;
00486
00487 strip_pidhit2 = -9999999 ;
00488 strip_simproc2 = -9999999 ;
00489
00490 strip_pidhit3 = -9999999 ;
00491 strip_simproc3 = -9999999 ;
00492
00493 strip_pidhit4 = -9999999 ;
00494 strip_simproc4 = -9999999 ;
00495
00496 strip_pidhit5 = -9999999 ;
00497 strip_simproc5 = -9999999 ;
00498
00499 strip_split = -9999999;
00500 strip_clst_err_x = -9999999.9;
00501 strip_clst_err_y = -9999999.9;
00502
00503 const TransientTrackingRecHit::ConstRecHitPointer trans_trk_rec_hit_point = itTraj->recHit();
00504
00505 if ( trans_trk_rec_hit_point == NULL )
00506 continue;
00507
00508 const TrackingRecHit *trk_rec_hit = (*trans_trk_rec_hit_point).hit();
00509
00510 if ( trk_rec_hit == NULL )
00511 continue;
00512
00513 DetId detid = (trk_rec_hit)->geographicalId();
00514
00515 strip_subdet_id = (int)detid.subdetId();
00516
00517 if ( (int)detid.subdetId() != (int)(StripSubdetector::TIB) && (int)detid.subdetId() != (int)(StripSubdetector::TOB) )
00518 continue;
00519
00520 const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>( (*trans_trk_rec_hit_point).hit() );
00521 const SiStripRecHit2D * hit2d = dynamic_cast<const SiStripRecHit2D *>( (*trans_trk_rec_hit_point).hit() );
00522 const SiStripRecHit1D * hit1d = dynamic_cast<const SiStripRecHit1D *>( (*trans_trk_rec_hit_point).hit() );
00523
00524 if ( !matchedhit && !hit2d && !hit1d )
00525 continue;
00526
00527 position = (trk_rec_hit)->localPosition();
00528 error = (trk_rec_hit)->localPositionError();
00529
00530 strip_rechitx = position.x();
00531 strip_rechity = position.y();
00532 strip_rechitz = position.z();
00533 strip_rechiterrx = sqrt( error.xx() );
00534 strip_rechiterry = sqrt( error.yy() );
00535
00536
00537
00538
00539
00540
00541
00542
00543
00544 TrajectoryStateOnSurface tsos = itTraj->updatedState();
00545
00546 strip_trk_pt = tsos.globalMomentum().perp();
00547
00548 LocalTrajectoryParameters ltp = tsos.localParameters();
00549
00550 LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
00551
00552 float locx = localDir.x();
00553 float locy = localDir.y();
00554 float locz = localDir.z();
00555
00556
00557
00558
00559 strip_cotalpha = locx/locz;
00560 strip_cotbeta = locy/locz;
00561
00562
00563 StripSubdetector StripSubdet = (StripSubdetector)detid;
00564
00565 if ( StripSubdet.stereo() == 0 )
00566 strip_is_stereo = 0;
00567 else
00568 strip_is_stereo = 1;
00569
00570
00571 SiStripDetId si_strip_det_id = SiStripDetId( detid );
00572
00573
00574 const StripGeomDetUnit* strip_geom_det_unit = (const StripGeomDetUnit*)tracker->idToDetUnit( detid );
00575
00576 if ( strip_geom_det_unit != NULL )
00577 {
00578 LocalVector lbfield
00579 = (strip_geom_det_unit->surface()).toLocal( (*magneticField).inTesla( strip_geom_det_unit->surface().position() ) );
00580
00581 strip_locbx = lbfield.x();
00582 strip_locby = lbfield.y();
00583 strip_locbz = lbfield.z();
00584 }
00585
00586
00587
00588
00589 if ( si_strip_det_id.moduleGeometry() == 1 )
00590 {
00591 detector_type = 1;
00592
00593
00594 }
00595 else if ( si_strip_det_id.moduleGeometry() == 2 )
00596 {
00597 detector_type = 2;
00598
00599
00600 }
00601 else if ( si_strip_det_id.moduleGeometry() == 3 )
00602 {
00603 detector_type = 3;
00604
00605
00606 }
00607 else if ( si_strip_det_id.moduleGeometry() == 4 )
00608 {
00609 detector_type = 4;
00610
00611
00612 }
00613
00614
00615
00616
00617 if ( (int)detid.subdetId() == int(StripSubdetector::TIB) )
00618 {
00619
00620 TIBDetId tib_detid( detid );
00621
00622 strip_tib_layer = (int)tib_detid.layer();
00623 strip_tib_module = (int)tib_detid.module();
00624 strip_tib_order = (int)tib_detid.order();
00625 strip_tib_side = (int)tib_detid.side();
00626 strip_tib_is_double_side = (int)tib_detid.isDoubleSide();
00627 strip_tib_is_z_plus_side = (int)tib_detid.isZPlusSide();
00628 strip_tib_is_z_minus_side = (int)tib_detid.isZMinusSide();
00629 strip_tib_layer_number = (int)tib_detid.layerNumber();
00630 strip_tib_string_number = (int)tib_detid.stringNumber() ;
00631 strip_tib_module_number = (int)tib_detid.moduleNumber();
00632 strip_tib_is_internal_string = (int)tib_detid.isInternalString();
00633 strip_tib_is_external_string = (int)tib_detid.isExternalString();
00634 strip_tib_is_rphi = (int)tib_detid.isRPhi();
00635 strip_tib_is_stereo = (int)tib_detid.isStereo();
00636 }
00637
00638
00639 if ( (int)detid.subdetId() == int(StripSubdetector::TOB) )
00640 {
00641
00642 TOBDetId tob_detid( detid );
00643
00644 strip_tob_layer = (int)tob_detid.layer();
00645 strip_tob_module = (int)tob_detid.module();
00646
00647 strip_tob_side = (int)tob_detid.side();
00648 strip_tob_is_double_side = (int)tob_detid.isDoubleSide();
00649 strip_tob_is_z_plus_side = (int)tob_detid.isZPlusSide();
00650 strip_tob_is_z_minus_side = (int)tob_detid.isZMinusSide();
00651 strip_tob_layer_number = (int)tob_detid.layerNumber();
00652 strip_tob_rod_number = (int)tob_detid.rodNumber();
00653 strip_tob_module_number = (int)tob_detid.moduleNumber();
00654
00655
00656 strip_tob_is_rphi = (int)tob_detid.isRPhi();
00657 strip_tob_is_stereo = (int)tob_detid.isStereo();
00658
00659 }
00660
00661
00662 if ( matchedhit )
00663 {
00664
00665
00666
00667 strip_hit_type = 0;
00668
00669 }
00670
00671
00672 if ( hit1d )
00673 {
00674 strip_hit_type = 1;
00675
00676 const SiStripRecHit1D::ClusterRef & cluster = hit1d->cluster();
00677
00678 if ( cluster->getSplitClusterError() > 0.0 )
00679 strip_split = 1;
00680 else
00681 strip_split = 0;
00682
00683 strip_clst_err_x = cluster->getSplitClusterError();
00684
00685
00686
00687 const std::vector<uint8_t>& stripCharges = cluster->amplitudes();
00688 uint16_t charge = 0;
00689 for (unsigned int i = 0; i < stripCharges.size(); ++i)
00690 {
00691 charge += stripCharges.at(i);
00692 }
00693
00694 strip_charge = (float)charge;
00695 strip_size = (int)( (cluster->amplitudes()).size() );
00696
00697
00698
00699 float mindist = 999999.9;
00700 matched.clear();
00701
00702 matched = associate.associateHit(*hit1d);
00703
00704 strip_nsimhit = (int)matched.size();
00705
00706 if ( !matched.empty())
00707 {
00708 PSimHit closest;
00709
00710
00711
00712 int strip_nprimaries = 0;
00713 int current_index = 0;
00714
00715 for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
00716 {
00717 ++current_index;
00718
00719 if ( (*m).processType() == 2 )
00720 ++strip_nprimaries;
00721
00722 if ( current_index == 1 )
00723 {
00724 strip_pidhit1 = (*m).particleType();
00725 strip_simproc1 = (*m).processType();
00726 }
00727 else if ( current_index == 2 )
00728 {
00729 strip_pidhit2 = (*m).particleType();
00730 strip_simproc2 = (*m).processType();
00731 }
00732 else if ( current_index == 3 )
00733 {
00734 strip_pidhit3 = (*m).particleType();
00735 strip_simproc3 = (*m).processType();
00736 }
00737 else if ( current_index == 4 )
00738 {
00739 strip_pidhit4 = (*m).particleType();
00740 strip_simproc4 = (*m).processType();
00741 }
00742 else if ( current_index == 5 )
00743 {
00744 strip_pidhit5 = (*m).particleType();
00745 strip_simproc5 = (*m).processType();
00746 }
00747
00748
00749 float dist = abs( (hit1d)->localPosition().x() - (*m).localPosition().x() );
00750
00751 if ( dist<mindist )
00752 {
00753 mindist = dist;
00754 closest = (*m);
00755 }
00756
00757 }
00758
00759 strip_nprm = strip_nprimaries;
00760
00761 strip_rechitresx = strip_rechitx - closest.localPosition().x();
00762 strip_rechitresy = strip_rechity - closest.localPosition().y();
00763
00764 strip_rechitresx2 = strip_rechitx - matched[0].localPosition().x();
00765
00766 strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
00767 strip_rechitpully = strip_rechitresy / strip_rechiterry;
00768
00769 strip_pidhit = closest.particleType();
00770 strip_simproc = closest.processType();
00771
00772
00773 }
00774
00775
00776
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790
00791
00792
00793
00794
00795
00796
00797
00798
00799
00800
00801
00802
00803
00804
00805 }
00806
00807
00808 if ( hit2d )
00809 {
00810 strip_hit_type = 2;
00811
00812 const SiStripRecHit1D::ClusterRef & cluster = hit2d->cluster();
00813
00814
00815
00816
00817
00818
00819
00820
00821
00822
00823 const std::vector<uint8_t>& stripCharges = cluster->amplitudes();
00824 uint16_t charge = 0;
00825 for (unsigned int i = 0; i < stripCharges.size(); ++i)
00826 {
00827 charge += stripCharges.at(i);
00828 }
00829
00830 strip_charge = (float)charge;
00831 strip_size = (int)( (cluster->amplitudes()).size() );
00832
00833
00834 float mindist = 999999.9;
00835 matched.clear();
00836
00837 matched = associate.associateHit(*hit2d);
00838
00839 strip_nsimhit = (int)matched.size();
00840
00841 if ( !matched.empty())
00842 {
00843 PSimHit closest;
00844
00845
00846
00847 for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
00848 {
00849 float dist = abs( (hit2d)->localPosition().x() - (*m).localPosition().x() );
00850
00851 if ( dist<mindist )
00852 {
00853 mindist = dist;
00854 closest = (*m);
00855 }
00856 }
00857
00858 strip_rechitresx = strip_rechitx - closest.localPosition().x();
00859 strip_rechitresy = strip_rechity - closest.localPosition().y();
00860
00861 strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
00862 strip_rechitpully = strip_rechitresy / strip_rechiterry;
00863
00864 strip_pidhit = closest.particleType();
00865 strip_simproc = closest.processType();
00866
00867
00868 }
00869
00870
00871
00872
00873
00874
00875
00876
00877
00878
00879
00880
00881 }
00882
00883
00884 ttree_track_hits_strip_->Fill();
00885
00886
00887 }
00888
00889 }
00890
00891
00892
00893
00894
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906 edm::Handle<SiPixelRecHitCollection> recHitColl;
00907 e.getByLabel( "siPixelRecHits", recHitColl);
00908
00909 Handle<edm::SimTrackContainer> simtracks;
00910 e.getByLabel("g4SimHits", simtracks);
00911
00912
00913 for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin(); it != pDD->dets().end(); it++)
00914 {
00915 DetId detId = ((*it)->geographicalId());
00916
00917 SiPixelRecHitCollection::const_iterator dsmatch = recHitColl->find(detId);
00918 if (dsmatch == recHitColl->end()) continue;
00919
00920 SiPixelRecHitCollection::DetSet pixelrechitRange = *dsmatch;
00921 SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
00922 SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
00923 SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
00924 std::vector<PSimHit> matched;
00925
00926
00927 for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter)
00928 {
00929 matched.clear();
00930 matched = associate.associateHit(*pixeliter);
00931
00932
00933 if ( matched.empty() )
00934 {
00935 cout << "SiPixelErrorEstimation::analyze: rechits without associated simhit !!!!!!!"
00936 << endl;
00937 continue;
00938 }
00939
00940 all_subdetid = -9999;
00941
00942 all_layer = -9999;
00943 all_ladder = -9999;
00944 all_mod = -9999;
00945
00946 all_side = -9999;
00947 all_disk = -9999;
00948 all_blade = -9999;
00949 all_panel = -9999;
00950 all_plaq = -9999;
00951
00952 all_half = -9999;
00953 all_flipped = -9999;
00954
00955 all_cols = -9999;
00956 all_rows = -9999;
00957
00958 all_rechitx = -9999;
00959 all_rechity = -9999;
00960 all_rechitz = -9999;
00961
00962 all_simhitx = -9999;
00963 all_simhity = -9999;
00964
00965 all_rechiterrx = -9999;
00966 all_rechiterry = -9999;
00967
00968 all_rechitresx = -9999;
00969 all_rechitresy = -9999;
00970
00971 all_rechitpullx = -9999;
00972 all_rechitpully = -9999;
00973
00974 all_npix = -9999;
00975 all_nxpix = -9999;
00976 all_nypix = -9999;
00977
00978 all_edgex = -9999;
00979 all_edgey = -9999;
00980
00981 all_bigx = -9999;
00982 all_bigy = -9999;
00983
00984 all_alpha = -9999;
00985 all_beta = -9999;
00986
00987 all_simphi = -9999;
00988 all_simtheta = -9999;
00989
00990 all_simhitx = -9999;
00991 all_simhity = -9999;
00992
00993 all_nsimhit = -9999;
00994 all_pidhit = -9999;
00995 all_simproc = -9999;
00996
00997 all_vtxr = -9999;
00998 all_vtxz = -9999;
00999
01000 all_simpx = -9999;
01001 all_simpy = -9999;
01002 all_simpz = -9999;
01003
01004 all_eloss = -9999;
01005
01006 all_trkid = -9999;
01007
01008 all_x1 = -9999;
01009 all_x2 = -9999;
01010 all_y1 = -9999;
01011 all_y2 = -9999;
01012 all_z1 = -9999;
01013 all_z2 = -9999;
01014
01015 all_row1 = -9999;
01016 all_row2 = -9999;
01017 all_col1 = -9999;
01018 all_col2 = -9999;
01019
01020 all_gx1 = -9999;
01021 all_gx2 = -9999;
01022 all_gy1 = -9999;
01023 all_gy2 = -9999;
01024 all_gz1 = -9999;
01025 all_gz2 = -9999;
01026
01027 all_simtrketa = -9999;
01028 all_simtrkphi = -9999;
01029
01030 all_clust_row = -9999;
01031 all_clust_col = -9999;
01032
01033 all_clust_x = -9999;
01034 all_clust_y = -9999;
01035
01036 all_clust_q = -9999;
01037
01038 all_clust_maxpixcol = -9999;
01039 all_clust_maxpixrow = -9999;
01040 all_clust_minpixcol = -9999;
01041 all_clust_minpixrow = -9999;
01042
01043 all_clust_geoid = -9999;
01044
01045 all_clust_alpha = -9999;
01046 all_clust_beta = -9999;
01047
01048
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062 all_hit_probx = -9999;
01063 all_hit_proby = -9999;
01064 all_hit_cprob0 = -9999;
01065 all_hit_cprob1 = -9999;
01066 all_hit_cprob2 = -9999;
01067
01068 all_pixel_split = -9999;
01069 all_pixel_clst_err_x = -9999.9;
01070 all_pixel_clst_err_y = -9999.9;
01071
01072 all_nsimhit = (int)matched.size();
01073
01074 all_subdetid = (int)detId.subdetId();
01075
01076 if ( !(all_subdetid==1 || all_subdetid==2) )
01077 {
01078 cout << "SiPixelErrorEstimation::analyze: Not in a pixel detector !!!!!" << endl;
01079 continue;
01080 }
01081
01082 const PixelGeomDetUnit* theGeomDet
01083 = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(detId) );
01084
01085 const PixelTopology* topol = &(theGeomDet->specificTopology());
01086
01087 if ( pixeliter->cluster()->getSplitClusterErrorX() > 0.0 &&
01088 pixeliter->cluster()->getSplitClusterErrorY() > 0.0 )
01089 {
01090 all_pixel_split = 1;
01091 }
01092 else
01093 {
01094 all_pixel_split = 0;
01095 }
01096
01097 all_pixel_clst_err_x = pixeliter->cluster()->getSplitClusterErrorX();
01098 all_pixel_clst_err_y = pixeliter->cluster()->getSplitClusterErrorY();
01099
01100
01101 const int maxPixelCol = pixeliter->cluster()->maxPixelCol();
01102 const int maxPixelRow = pixeliter->cluster()->maxPixelRow();
01103 const int minPixelCol = pixeliter->cluster()->minPixelCol();
01104 const int minPixelRow = pixeliter->cluster()->minPixelRow();
01105
01106
01107
01108 all_hit_cprob0 = (float)pixeliter->clusterProbability(0);
01109 all_hit_cprob1 = (float)pixeliter->clusterProbability(1);
01110 all_hit_cprob2 = (float)pixeliter->clusterProbability(2);
01111
01112
01113 if ( topol->isItEdgePixelInX( minPixelRow ) ||
01114 topol->isItEdgePixelInX( maxPixelRow ) )
01115 all_edgex = 1;
01116 else
01117 all_edgex = 0;
01118
01119 if ( topol->isItEdgePixelInY( minPixelCol ) ||
01120 topol->isItEdgePixelInY( maxPixelCol ) )
01121 all_edgey = 1;
01122 else
01123 all_edgey = 0;
01124
01125
01126 if ( topol->containsBigPixelInX(minPixelRow, maxPixelRow) )
01127 all_bigx = 1;
01128 else
01129 all_bigx = 0;
01130
01131 if ( topol->containsBigPixelInY(minPixelCol, maxPixelCol) )
01132 all_bigy = 1;
01133 else
01134 all_bigy = 0;
01135
01136 if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel )
01137 {
01138 PXBDetId bdetid(detId);
01139 all_layer = bdetid.layer();
01140 all_ladder = bdetid.ladder();
01141 all_mod = bdetid.module();
01142
01143 int tmp_nrows = theGeomDet->specificTopology().nrows();
01144 if ( tmp_nrows == 80 )
01145 all_half = 1;
01146 else if ( tmp_nrows == 160 )
01147 all_half = 0;
01148 else
01149 cout << "-------------------------------------------------- Wrong module size !!!" << endl;
01150
01151 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01152 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01153
01154 if ( tmp2<tmp1 )
01155 all_flipped = 1;
01156 else
01157 all_flipped = 0;
01158 }
01159 else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
01160 {
01161 PXFDetId fdetid(detId);
01162 all_side = fdetid.side();
01163 all_disk = fdetid.disk();
01164 all_blade = fdetid.blade();
01165 all_panel = fdetid.panel();
01166 all_plaq = fdetid.module();
01167
01168 }
01169 else std::cout << "We are not in the pixel detector" << (int)detId.subdetId() << endl;
01170
01171 all_cols = theGeomDet->specificTopology().ncolumns();
01172 all_rows = theGeomDet->specificTopology().nrows();
01173
01174 LocalPoint lp = pixeliter->localPosition();
01175
01176 all_rechitx = lp.x();
01177 all_rechity = lp.y();
01178 all_rechitz = lp.z();
01179
01180 LocalError le = pixeliter->localPositionError();
01181 all_rechiterrx = sqrt( le.xx() );
01182 all_rechiterry = sqrt( le.yy() );
01183
01184 bool found_hit_from_generated_particle = false;
01185
01186
01187 float closest_dist = 99999.9;
01188 std::vector<PSimHit>::const_iterator closest_simhit = matched.begin();
01189
01190 for (std::vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); m++)
01191 {
01192 if ( checkType_ )
01193 {
01194 int pid = (*m).particleType();
01195 if ( abs(pid) != genType_ )
01196 continue;
01197 }
01198
01199 float simhitx = 0.5 * ( (*m).entryPoint().x() + (*m).exitPoint().x() );
01200 float simhity = 0.5 * ( (*m).entryPoint().y() + (*m).exitPoint().y() );
01201
01202 float x_res = simhitx - rechitx;
01203 float y_res = simhity - rechity;
01204
01205 float dist = sqrt( x_res*x_res + y_res*y_res );
01206
01207 if ( dist < closest_dist )
01208 {
01209 closest_dist = dist;
01210 closest_simhit = m;
01211 found_hit_from_generated_particle = true;
01212 }
01213 }
01214
01215
01216
01217 if ( checkType_ && !found_hit_from_generated_particle )
01218 continue;
01219
01220 all_x1 = (*closest_simhit).entryPoint().x();
01221 all_y1 = (*closest_simhit).entryPoint().y();
01222 all_z1 = (*closest_simhit).entryPoint().z();
01223 all_x2 = (*closest_simhit).exitPoint().x();
01224 all_y2 = (*closest_simhit).exitPoint().y();
01225 all_z2 = (*closest_simhit).exitPoint().z();
01226 GlobalPoint GP1 =
01227 theGeomDet->surface().toGlobal( Local3DPoint( (*closest_simhit).entryPoint().x(),
01228 (*closest_simhit).entryPoint().y(),
01229 (*closest_simhit).entryPoint().z() ) );
01230 GlobalPoint GP2 =
01231 theGeomDet->surface().toGlobal (Local3DPoint( (*closest_simhit).exitPoint().x(),
01232 (*closest_simhit).exitPoint().y(),
01233 (*closest_simhit).exitPoint().z() ) );
01234 all_gx1 = GP1.x();
01235 all_gx2 = GP2.x();
01236 all_gy1 = GP1.y();
01237 all_gy2 = GP2.y();
01238 all_gz1 = GP1.z();
01239 all_gz2 = GP2.z();
01240
01241 MeasurementPoint mp1 =
01242 topol->measurementPosition( LocalPoint( (*closest_simhit).entryPoint().x(),
01243 (*closest_simhit).entryPoint().y(),
01244 (*closest_simhit).entryPoint().z() ) );
01245 MeasurementPoint mp2 =
01246 topol->measurementPosition( LocalPoint( (*closest_simhit).exitPoint().x(),
01247 (*closest_simhit).exitPoint().y(),
01248 (*closest_simhit).exitPoint().z() ) );
01249 all_row1 = mp1.x();
01250 all_col1 = mp1.y();
01251 all_row2 = mp2.x();
01252 all_col2 = mp2.y();
01253
01254 all_simhitx = 0.5*(all_x1+all_x2);
01255 all_simhity = 0.5*(all_y1+all_y2);
01256
01257 all_rechitresx = all_rechitx - all_simhitx;
01258 all_rechitresy = all_rechity - all_simhity;
01259
01260 all_rechitpullx = all_rechitresx / all_rechiterrx;
01261 all_rechitpully = all_rechitresy / all_rechiterry;
01262
01263 SiPixelRecHit::ClusterRef const& clust = pixeliter->cluster();
01264
01265 all_npix = clust->size();
01266 all_nxpix = clust->sizeX();
01267 all_nypix = clust->sizeY();
01268
01269 all_clust_row = clust->x();
01270 all_clust_col = clust->y();
01271
01272 LocalPoint lp2 = topol->localPosition( MeasurementPoint( all_clust_row, all_clust_col ) );
01273 all_clust_x = lp2.x();
01274 all_clust_y = lp2.y();
01275
01276 all_clust_q = clust->charge();
01277
01278 all_clust_maxpixcol = clust->maxPixelCol();
01279 all_clust_maxpixrow = clust->maxPixelRow();
01280 all_clust_minpixcol = clust->minPixelCol();
01281 all_clust_minpixrow = clust->minPixelRow();
01282
01283 all_clust_geoid = 0;
01284
01285 all_simpx = (*closest_simhit).momentumAtEntry().x();
01286 all_simpy = (*closest_simhit).momentumAtEntry().y();
01287 all_simpz = (*closest_simhit).momentumAtEntry().z();
01288 all_eloss = (*closest_simhit).energyLoss();
01289 all_simphi = (*closest_simhit).phiAtEntry();
01290 all_simtheta = (*closest_simhit).thetaAtEntry();
01291 all_pidhit = (*closest_simhit).particleType();
01292 all_trkid = (*closest_simhit).trackId();
01293
01294
01295 all_beta = atan2(all_simpz, all_simpy);
01296 all_alpha = atan2(all_simpz, all_simpx);
01297
01298 all_simproc = (int)closest_simhit->processType();
01299
01300 const edm::SimTrackContainer& trks = *(simtracks.product());
01301 SimTrackContainer::const_iterator trksiter;
01302 for (trksiter = trks.begin(); trksiter != trks.end(); trksiter++)
01303 if ( (int)trksiter->trackId() == all_trkid )
01304 {
01305 all_simtrketa = trksiter->momentum().eta();
01306 all_simtrkphi = trksiter->momentum().phi();
01307 }
01308
01309 all_vtxz = theGeomDet->surface().position().z();
01310 all_vtxr = theGeomDet->surface().position().perp();
01311
01312
01313
01314
01315
01316 const std::vector<SiPixelCluster::Pixel>& pixvector = clust->pixels();
01317 for ( int i=0; i<(int)pixvector.size(); ++i)
01318 {
01319 SiPixelCluster::Pixel holdpix = pixvector[i];
01320 all_pixrow[i] = holdpix.x;
01321 all_pixcol[i] = holdpix.y;
01322 all_pixadc[i] = holdpix.adc;
01323 LocalPoint lp = topol->localPosition(MeasurementPoint(holdpix.x, holdpix.y));
01324 all_pixx[i] = lp.x();
01325 all_pixy[i]= lp.y();
01326 GlobalPoint GP = theGeomDet->surface().toGlobal(Local3DPoint(lp.x(),lp.y(),lp.z()));
01327 all_pixgx[i] = GP.x();
01328 all_pixgy[i]= GP.y();
01329 all_pixgz[i]= GP.z();
01330 }
01331
01332 ttree_all_hits_->Fill();
01333
01334 }
01335
01336 }
01337
01338
01339
01340
01341
01342
01343
01344
01345
01346
01347
01348 if ( include_trk_hits_ )
01349 {
01350
01351 edm::Handle<reco::TrackCollection> trackCollection;
01352 e.getByLabel(src_, trackCollection);
01353 const reco::TrackCollection *tracks = trackCollection.product();
01354 reco::TrackCollection::const_iterator tciter;
01355
01356 if ( tracks->size() > 0 )
01357 {
01358
01359 for ( tciter=tracks->begin(); tciter!=tracks->end(); ++tciter)
01360 {
01361
01362 for ( trackingRecHit_iterator it = tciter->recHitsBegin(); it != tciter->recHitsEnd(); ++it)
01363 {
01364 const TrackingRecHit &trk_rec_hit = **it;
01365
01366 const SiPixelRecHit* matchedhit = dynamic_cast<const SiPixelRecHit*>(&trk_rec_hit);
01367
01368 if ( matchedhit )
01369 {
01370 rechitx = -9999.9;
01371 rechity = -9999.9;
01372 rechitz = -9999.9;
01373 rechiterrx = -9999.9;
01374 rechiterry = -9999.9;
01375 rechitresx = -9999.9;
01376 rechitresy = -9999.9;
01377 rechitpullx = -9999.9;
01378 rechitpully = -9999.9;
01379
01380 npix = -9999;
01381 nxpix = -9999;
01382 nypix = -9999;
01383 charge = -9999.9;
01384
01385 edgex = -9999;
01386 edgey = -9999;
01387
01388 bigx = -9999;
01389 bigy = -9999;
01390
01391 alpha = -9999.9;
01392 beta = -9999.9;
01393
01394 phi = -9999.9;
01395 eta = -9999.9;
01396
01397 subdetId = -9999;
01398
01399 layer = -9999;
01400 ladder = -9999;
01401 mod = -9999;
01402 side = -9999;
01403 disk = -9999;
01404 blade = -9999;
01405 panel = -9999;
01406 plaq = -9999;
01407
01408 half = -9999;
01409 flipped = -9999;
01410
01411 nsimhit = -9999;
01412 pidhit = -9999;
01413 simproc = -9999;
01414
01415 simhitx = -9999.9;
01416 simhity = -9999.9;
01417
01418 hit_probx = -9999.9;
01419 hit_proby = -9999.9;
01420 hit_cprob0 = -9999.9;
01421 hit_cprob1 = -9999.9;
01422 hit_cprob2 = -9999.9;
01423
01424 pixel_split = -9999;
01425
01426 pixel_clst_err_x = -9999.9;
01427 pixel_clst_err_y = -9999.9;
01428
01429 position = (*it)->localPosition();
01430 error = (*it)->localPositionError();
01431
01432 rechitx = position.x();
01433 rechity = position.y();
01434 rechitz = position.z();
01435 rechiterrx = sqrt(error.xx());
01436 rechiterry = sqrt(error.yy());
01437
01438 npix = matchedhit->cluster()->size();
01439 nxpix = matchedhit->cluster()->sizeX();
01440 nypix = matchedhit->cluster()->sizeY();
01441 charge = matchedhit->cluster()->charge();
01442
01443 if ( matchedhit->cluster()->getSplitClusterErrorX() > 0.0 &&
01444 matchedhit->cluster()->getSplitClusterErrorY() > 0.0 )
01445 pixel_split = 1;
01446 else
01447 pixel_split = 0;
01448
01449 pixel_clst_err_x = matchedhit->cluster()->getSplitClusterErrorX();
01450 pixel_clst_err_y = matchedhit->cluster()->getSplitClusterErrorY();
01451
01452
01453
01454 hit_cprob0 = (float)matchedhit->clusterProbability(0);
01455 hit_cprob1 = (float)matchedhit->clusterProbability(1);
01456 hit_cprob2 = (float)matchedhit->clusterProbability(2);
01457
01458
01459
01460 matched.clear();
01461 matched = associate.associateHit(*matchedhit);
01462
01463 nsimhit = (int)matched.size();
01464
01465 if ( !matched.empty() )
01466 {
01467 mindist = 999999.9;
01468 float distx, disty, dist;
01469 bool found_hit_from_generated_particle = false;
01470
01471 int n_assoc_muon = 0;
01472
01473 vector<PSimHit>::const_iterator closestit = matched.begin();
01474 for (vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
01475 {
01476 if ( checkType_ )
01477 {
01478 int pid = (*m).particleType();
01479 if ( abs(pid) != genType_ )
01480 continue;
01481 }
01482
01483 float simhitx = 0.5 * ( (*m).entryPoint().x() + (*m).exitPoint().x() );
01484 float simhity = 0.5 * ( (*m).entryPoint().y() + (*m).exitPoint().y() );
01485
01486 distx = fabs(rechitx - simhitx);
01487 disty = fabs(rechity - simhity);
01488 dist = sqrt( distx*distx + disty*disty );
01489
01490 if ( dist < mindist )
01491 {
01492 n_assoc_muon++;
01493
01494 mindist = dist;
01495 closestit = m;
01496 found_hit_from_generated_particle = true;
01497 }
01498 }
01499
01500
01501
01502 if ( checkType_ && !found_hit_from_generated_particle )
01503 continue;
01504
01505
01506
01507
01508
01509
01510
01511 DetId detId = (*it)->geographicalId();
01512
01513 const PixelGeomDetUnit* theGeomDet =
01514 dynamic_cast<const PixelGeomDetUnit*> ((*tracker).idToDet(detId) );
01515
01516 const PixelTopology* theTopol = &(theGeomDet->specificTopology());
01517
01518 pidhit = (*closestit).particleType();
01519 simproc = (int)(*closestit).processType();
01520
01521 simhitx = 0.5*( (*closestit).entryPoint().x() + (*closestit).exitPoint().x() );
01522 simhity = 0.5*( (*closestit).entryPoint().y() + (*closestit).exitPoint().y() );
01523
01524 rechitresx = rechitx - simhitx;
01525 rechitresy = rechity - simhity;
01526 rechitpullx = ( rechitx - simhitx ) / sqrt(error.xx());
01527 rechitpully = ( rechity - simhity ) / sqrt(error.yy());
01528
01529 float simhitpx = (*closestit).momentumAtEntry().x();
01530 float simhitpy = (*closestit).momentumAtEntry().y();
01531 float simhitpz = (*closestit).momentumAtEntry().z();
01532
01533 beta = atan2(simhitpz, simhitpy) * radtodeg;
01534 alpha = atan2(simhitpz, simhitpx) * radtodeg;
01535
01536
01537
01538
01539
01540 float locx = simhitpx;
01541 float locy = simhitpy;
01542 float locz = simhitpz;
01543
01544 bool isFlipped = false;
01545 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01546 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01547 if ( tmp2<tmp1 )
01548 isFlipped = true;
01549 else
01550 isFlipped = false;
01551
01552 trk_alpha = acos(locx/sqrt(locx*locx+locz*locz)) * radtodeg;
01553 if ( isFlipped )
01554 trk_alpha = 180.0 - trk_alpha ;
01555
01556 trk_beta = acos(locy/sqrt(locy*locy+locz*locz)) * radtodeg;
01557
01558
01559 phi = tciter->momentum().phi() / math_pi*180.0;
01560 eta = tciter->momentum().eta();
01561
01562 const int maxPixelCol = (*matchedhit).cluster()->maxPixelCol();
01563 const int maxPixelRow = (*matchedhit).cluster()->maxPixelRow();
01564 const int minPixelCol = (*matchedhit).cluster()->minPixelCol();
01565 const int minPixelRow = (*matchedhit).cluster()->minPixelRow();
01566
01567
01568 if ( theTopol->isItEdgePixelInX( minPixelRow ) ||
01569 theTopol->isItEdgePixelInX( maxPixelRow ) )
01570 edgex = 1;
01571 else
01572 edgex = 0;
01573
01574 if ( theTopol->isItEdgePixelInY( minPixelCol ) ||
01575 theTopol->isItEdgePixelInY( maxPixelCol ) )
01576 edgey = 1;
01577 else
01578 edgey = 0;
01579
01580
01581 if ( theTopol->containsBigPixelInX(minPixelRow, maxPixelRow) )
01582 bigx = 1;
01583 else
01584 bigx = 0;
01585
01586 if ( theTopol->containsBigPixelInY(minPixelCol, maxPixelCol) )
01587 bigy = 1;
01588 else
01589 bigy = 0;
01590
01591 subdetId = (int)detId.subdetId();
01592
01593 if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel )
01594 {
01595
01596 int tmp_nrows = theGeomDet->specificTopology().nrows();
01597 if ( tmp_nrows == 80 )
01598 half = 1;
01599 else if ( tmp_nrows == 160 )
01600 half = 0;
01601 else
01602 cout << "-------------------------------------------------- Wrong module size !!!" << endl;
01603
01604 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01605 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01606
01607 if ( tmp2<tmp1 )
01608 flipped = 1;
01609 else
01610 flipped = 0;
01611
01612 PXBDetId bdetid(detId);
01613 layer = bdetid.layer();
01614 ladder = bdetid.ladder();
01615 mod = bdetid.module();
01616 }
01617 else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
01618 {
01619 PXFDetId fdetid(detId);
01620 side = fdetid.side();
01621 disk = fdetid.disk();
01622 blade = fdetid.blade();
01623 panel = fdetid.panel();
01624 plaq = fdetid.module();
01625
01626 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01627 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01628
01629 if ( tmp2<tmp1 )
01630 flipped = 1;
01631 else
01632 flipped = 0;
01633
01634 }
01635
01636
01637 ttree_track_hits_->Fill();
01638
01639 }
01640 else
01641 cout << "---------------- RecHit with no associated SimHit !!! -------------------------- " << endl;
01642
01643 }
01644
01645 }
01646
01647 }
01648
01649 }
01650
01651 }
01652
01653
01654
01655
01656
01657 }
01658
01659 void SiPixelErrorEstimation::
01660 computeAnglesFromDetPosition(const SiPixelCluster & cl,
01661 const GeomDetUnit & det,
01662 float& alpha, float& beta )
01663 {
01664
01665 const PixelGeomDetUnit* theDet = dynamic_cast<const PixelGeomDetUnit*>( &det );
01666 if (! theDet)
01667 {
01668 cout << "---------------------------------------------- Not a pixel detector !!!!!!!!!!!!!!" << endl;
01669 assert(0);
01670 }
01671
01672 const PixelTopology* theTopol = &(theDet->specificTopology());
01673
01674
01675 float xcenter = cl.x();
01676 float ycenter = cl.y();
01677
01678
01679 LocalPoint lp = theTopol->localPosition( MeasurementPoint(xcenter, ycenter) );
01680
01681
01682
01683 GlobalPoint gp = theDet->surface().toGlobal( lp );
01684 float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() );
01685
01686
01687 float gpx = gp.x()/gp_mod;
01688 float gpy = gp.y()/gp_mod;
01689 float gpz = gp.z()/gp_mod;
01690
01691
01692
01693 GlobalVector gv(gpx, gpy, gpz);
01694
01695
01696 const Local3DVector lvx(1.0, 0.0, 0.0);
01697
01698
01699 GlobalVector gvx = theDet->surface().toGlobal( lvx );
01700
01701
01702 const Local3DVector lvy(0.0, 1.0, 0.0);
01703
01704
01705 GlobalVector gvy = theDet->surface().toGlobal( lvy );
01706
01707
01708 const Local3DVector lvz(0.0, 0.0, 1.0);
01709
01710
01711 GlobalVector gvz = theDet->surface().toGlobal( lvz );
01712
01713
01714
01715
01716 float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z();
01717 float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z();
01718 float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z();
01719
01720
01721 alpha = atan2( gv_dot_gvz, gv_dot_gvx );
01722 beta = atan2( gv_dot_gvz, gv_dot_gvy );
01723
01724
01725
01726
01727
01728
01729
01730 }
01731
01732
01733 DEFINE_FWK_MODULE(SiPixelErrorEstimation);