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/TrackerCommon/interface/TrackerTopology.h"
00019 #include "Geometry/Records/interface/IdealGeometryRecord.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
00367 edm::ESHandle<TrackerTopology> tTopoHandle;
00368 es.get<IdealGeometryRecord>().get(tTopoHandle);
00369 const TrackerTopology* const tTopo = tTopoHandle.product();
00370
00371 using namespace edm;
00372
00373 run = e.id().run();
00374 evt = e.id().event();
00375
00376 if ( evt%1000 == 0 )
00377 cout << "evt = " << evt << endl;
00378
00379 float math_pi = 3.14159265;
00380 float radtodeg = 180.0 / math_pi;
00381
00382 LocalPoint position;
00383 LocalError error;
00384 float mindist = 999999.9;
00385
00386 std::vector<PSimHit> matched;
00387 TrackerHitAssociator associate(e);
00388
00389 edm::ESHandle<TrackerGeometry> pDD;
00390 es.get<TrackerDigiGeometryRecord> ().get (pDD);
00391 const TrackerGeometry* tracker = &(* pDD);
00392
00393
00394
00395
00396
00397 edm::ESHandle<MagneticField> magneticField;
00398 es.get<IdealMagneticFieldRecord>().get( magneticField );
00399
00400
00401 edm::FileInPath FileInPath_;
00402
00403
00404
00405
00406
00407
00408 edm::Handle<vector<Trajectory> > trajCollectionHandle;
00409
00410 e.getByLabel( src_, trajCollectionHandle);
00411
00412
00413 for ( vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end(); ++it )
00414 {
00415
00416 vector<TrajectoryMeasurement> tmColl = it->measurements();
00417 for ( vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj!=tmColl.end(); ++itTraj )
00418 {
00419
00420 if ( !itTraj->updatedState().isValid() )
00421 continue;
00422
00423 strip_rechitx = -9999999.9;
00424 strip_rechity = -9999999.9;
00425 strip_rechitz = -9999999.9;
00426 strip_rechiterrx = -9999999.9;
00427 strip_rechiterry = -9999999.9;
00428 strip_rechitresx = -9999999.9;
00429 strip_rechitresx2 = -9999999.9;
00430
00431
00432 strip_rechitresy = -9999999.9;
00433 strip_rechitpullx = -9999999.9;
00434 strip_rechitpully = -9999999.9;
00435 strip_is_stereo = -9999999 ;
00436 strip_hit_type = -9999999 ;
00437 detector_type = -9999999 ;
00438
00439 strip_trk_pt = -9999999.9;
00440 strip_cotalpha = -9999999.9;
00441 strip_cotbeta = -9999999.9;
00442 strip_locbx = -9999999.9;
00443 strip_locby = -9999999.9;
00444 strip_locbz = -9999999.9;
00445 strip_charge = -9999999.9;
00446 strip_size = -9999999 ;
00447
00448 strip_edge = -9999999 ;
00449 strip_nsimhit = -9999999 ;
00450 strip_pidhit = -9999999 ;
00451 strip_simproc = -9999999 ;
00452
00453
00454 strip_subdet_id = -9999999;
00455
00456 strip_tib_layer = -9999999;
00457 strip_tib_module = -9999999;
00458 strip_tib_order = -9999999;
00459 strip_tib_side = -9999999;
00460 strip_tib_is_double_side = -9999999;
00461 strip_tib_is_z_plus_side = -9999999;
00462 strip_tib_is_z_minus_side = -9999999;
00463 strip_tib_layer_number = -9999999;
00464 strip_tib_string_number = -9999999;
00465 strip_tib_module_number = -9999999;
00466 strip_tib_is_internal_string= -9999999;
00467 strip_tib_is_external_string= -9999999;
00468 strip_tib_is_rphi = -9999999;
00469 strip_tib_is_stereo = -9999999;
00470
00471 strip_tob_layer = -9999999;
00472 strip_tob_module = -9999999;
00473 strip_tob_side = -9999999;
00474 strip_tob_is_double_side = -9999999;
00475 strip_tob_is_z_plus_side = -9999999;
00476 strip_tob_is_z_minus_side = -9999999;
00477 strip_tob_layer_number = -9999999;
00478 strip_tob_rod_number = -9999999;
00479 strip_tob_module_number = -9999999;
00480
00481 strip_tob_is_rphi = -9999999;
00482 strip_tob_is_stereo = -9999999;
00483
00484 strip_prob = -9999999.9;
00485 strip_qbin = -9999999;
00486
00487 strip_nprm = -9999999;
00488
00489 strip_pidhit1 = -9999999 ;
00490 strip_simproc1 = -9999999 ;
00491
00492 strip_pidhit2 = -9999999 ;
00493 strip_simproc2 = -9999999 ;
00494
00495 strip_pidhit3 = -9999999 ;
00496 strip_simproc3 = -9999999 ;
00497
00498 strip_pidhit4 = -9999999 ;
00499 strip_simproc4 = -9999999 ;
00500
00501 strip_pidhit5 = -9999999 ;
00502 strip_simproc5 = -9999999 ;
00503
00504 strip_split = -9999999;
00505 strip_clst_err_x = -9999999.9;
00506 strip_clst_err_y = -9999999.9;
00507
00508 const TransientTrackingRecHit::ConstRecHitPointer trans_trk_rec_hit_point = itTraj->recHit();
00509
00510 if ( trans_trk_rec_hit_point == NULL )
00511 continue;
00512
00513 const TrackingRecHit *trk_rec_hit = (*trans_trk_rec_hit_point).hit();
00514
00515 if ( trk_rec_hit == NULL )
00516 continue;
00517
00518 DetId detid = (trk_rec_hit)->geographicalId();
00519
00520 strip_subdet_id = (int)detid.subdetId();
00521
00522 if ( (int)detid.subdetId() != (int)(StripSubdetector::TIB) && (int)detid.subdetId() != (int)(StripSubdetector::TOB) )
00523 continue;
00524
00525 const SiStripMatchedRecHit2D* matchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>( (*trans_trk_rec_hit_point).hit() );
00526 const SiStripRecHit2D * hit2d = dynamic_cast<const SiStripRecHit2D *>( (*trans_trk_rec_hit_point).hit() );
00527 const SiStripRecHit1D * hit1d = dynamic_cast<const SiStripRecHit1D *>( (*trans_trk_rec_hit_point).hit() );
00528
00529 if ( !matchedhit && !hit2d && !hit1d )
00530 continue;
00531
00532 position = (trk_rec_hit)->localPosition();
00533 error = (trk_rec_hit)->localPositionError();
00534
00535 strip_rechitx = position.x();
00536 strip_rechity = position.y();
00537 strip_rechitz = position.z();
00538 strip_rechiterrx = sqrt( error.xx() );
00539 strip_rechiterry = sqrt( error.yy() );
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549 TrajectoryStateOnSurface tsos = itTraj->updatedState();
00550
00551 strip_trk_pt = tsos.globalMomentum().perp();
00552
00553 LocalTrajectoryParameters ltp = tsos.localParameters();
00554
00555 LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
00556
00557 float locx = localDir.x();
00558 float locy = localDir.y();
00559 float locz = localDir.z();
00560
00561
00562
00563
00564 strip_cotalpha = locx/locz;
00565 strip_cotbeta = locy/locz;
00566
00567
00568 StripSubdetector StripSubdet = (StripSubdetector)detid;
00569
00570 if ( StripSubdet.stereo() == 0 )
00571 strip_is_stereo = 0;
00572 else
00573 strip_is_stereo = 1;
00574
00575
00576 SiStripDetId si_strip_det_id = SiStripDetId( detid );
00577
00578
00579 const StripGeomDetUnit* strip_geom_det_unit = (const StripGeomDetUnit*)tracker->idToDetUnit( detid );
00580
00581 if ( strip_geom_det_unit != NULL )
00582 {
00583 LocalVector lbfield
00584 = (strip_geom_det_unit->surface()).toLocal( (*magneticField).inTesla( strip_geom_det_unit->surface().position() ) );
00585
00586 strip_locbx = lbfield.x();
00587 strip_locby = lbfield.y();
00588 strip_locbz = lbfield.z();
00589 }
00590
00591
00592
00593
00594 if ( si_strip_det_id.moduleGeometry() == 1 )
00595 {
00596 detector_type = 1;
00597
00598
00599 }
00600 else if ( si_strip_det_id.moduleGeometry() == 2 )
00601 {
00602 detector_type = 2;
00603
00604
00605 }
00606 else if ( si_strip_det_id.moduleGeometry() == 3 )
00607 {
00608 detector_type = 3;
00609
00610
00611 }
00612 else if ( si_strip_det_id.moduleGeometry() == 4 )
00613 {
00614 detector_type = 4;
00615
00616
00617 }
00618
00619
00620
00621
00622 if ( (int)detid.subdetId() == int(StripSubdetector::TIB) )
00623 {
00624
00625
00626
00627 strip_tib_layer = (int)tTopo->tibLayer( detid );
00628 strip_tib_module = (int)tTopo->tibModule( detid );
00629 strip_tib_order = (int)tTopo->tibOrder( detid );
00630 strip_tib_side = (int)tTopo->tibSide( detid );
00631 strip_tib_is_double_side = (int)tTopo->tibIsDoubleSide( detid );
00632 strip_tib_is_z_plus_side = (int)tTopo->tibIsZPlusSide( detid );
00633 strip_tib_is_z_minus_side = (int)tTopo->tibIsZMinusSide( detid );
00634 strip_tib_layer_number = (int)tTopo->tibLayer( detid );
00635 strip_tib_string_number = (int)tTopo->tibString( detid ) ;
00636 strip_tib_module_number = (int)tTopo->tibModule( detid );
00637 strip_tib_is_internal_string = (int)tTopo->tibIsInternalString( detid );
00638 strip_tib_is_external_string = (int)tTopo->tibIsExternalString( detid );
00639 strip_tib_is_rphi = (int)tTopo->tibIsRPhi( detid );
00640 strip_tib_is_stereo = (int)tTopo->tibIsStereo( detid );
00641 }
00642
00643
00644 if ( (int)detid.subdetId() == int(StripSubdetector::TOB) )
00645 {
00646
00647
00648
00649 strip_tob_layer = (int)tTopo->tobLayer( detid );
00650 strip_tob_module = (int)tTopo->tobModule( detid );
00651
00652 strip_tob_side = (int)tTopo->tobSide( detid );
00653 strip_tob_is_double_side = (int)tTopo->tobIsDoubleSide( detid );
00654 strip_tob_is_z_plus_side = (int)tTopo->tobIsZPlusSide( detid );
00655 strip_tob_is_z_minus_side = (int)tTopo->tobIsZMinusSide( detid );
00656 strip_tob_layer_number = (int)tTopo->tobLayer( detid );
00657 strip_tob_rod_number = (int)tTopo->tobRod( detid );
00658 strip_tob_module_number = (int)tTopo->tobModule( detid );
00659
00660
00661 strip_tob_is_rphi = (int)tTopo->tobIsRPhi( detid );
00662 strip_tob_is_stereo = (int)tTopo->tobIsStereo( detid );
00663
00664 }
00665
00666
00667 if ( matchedhit )
00668 {
00669
00670
00671
00672 strip_hit_type = 0;
00673
00674 }
00675
00676
00677 if ( hit1d )
00678 {
00679 strip_hit_type = 1;
00680
00681 const SiStripRecHit1D::ClusterRef & cluster = hit1d->cluster();
00682
00683 if ( cluster->getSplitClusterError() > 0.0 )
00684 strip_split = 1;
00685 else
00686 strip_split = 0;
00687
00688 strip_clst_err_x = cluster->getSplitClusterError();
00689
00690
00691
00692 const std::vector<uint8_t>& stripCharges = cluster->amplitudes();
00693 uint16_t charge = 0;
00694 for (unsigned int i = 0; i < stripCharges.size(); ++i)
00695 {
00696 charge += stripCharges.at(i);
00697 }
00698
00699 strip_charge = (float)charge;
00700 strip_size = (int)( (cluster->amplitudes()).size() );
00701
00702
00703
00704 float mindist = 999999.9;
00705 matched.clear();
00706
00707 matched = associate.associateHit(*hit1d);
00708
00709 strip_nsimhit = (int)matched.size();
00710
00711 if ( !matched.empty())
00712 {
00713 PSimHit closest;
00714
00715
00716
00717 int strip_nprimaries = 0;
00718 int current_index = 0;
00719
00720 for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
00721 {
00722 ++current_index;
00723
00724 if ( (*m).processType() == 2 )
00725 ++strip_nprimaries;
00726
00727 if ( current_index == 1 )
00728 {
00729 strip_pidhit1 = (*m).particleType();
00730 strip_simproc1 = (*m).processType();
00731 }
00732 else if ( current_index == 2 )
00733 {
00734 strip_pidhit2 = (*m).particleType();
00735 strip_simproc2 = (*m).processType();
00736 }
00737 else if ( current_index == 3 )
00738 {
00739 strip_pidhit3 = (*m).particleType();
00740 strip_simproc3 = (*m).processType();
00741 }
00742 else if ( current_index == 4 )
00743 {
00744 strip_pidhit4 = (*m).particleType();
00745 strip_simproc4 = (*m).processType();
00746 }
00747 else if ( current_index == 5 )
00748 {
00749 strip_pidhit5 = (*m).particleType();
00750 strip_simproc5 = (*m).processType();
00751 }
00752
00753
00754 float dist = abs( (hit1d)->localPosition().x() - (*m).localPosition().x() );
00755
00756 if ( dist<mindist )
00757 {
00758 mindist = dist;
00759 closest = (*m);
00760 }
00761
00762 }
00763
00764 strip_nprm = strip_nprimaries;
00765
00766 strip_rechitresx = strip_rechitx - closest.localPosition().x();
00767 strip_rechitresy = strip_rechity - closest.localPosition().y();
00768
00769 strip_rechitresx2 = strip_rechitx - matched[0].localPosition().x();
00770
00771 strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
00772 strip_rechitpully = strip_rechitresy / strip_rechiterry;
00773
00774 strip_pidhit = closest.particleType();
00775 strip_simproc = closest.processType();
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
00809
00810 }
00811
00812
00813 if ( hit2d )
00814 {
00815 strip_hit_type = 2;
00816
00817 const SiStripRecHit1D::ClusterRef & cluster = hit2d->cluster();
00818
00819
00820
00821
00822
00823
00824
00825
00826
00827
00828 const std::vector<uint8_t>& stripCharges = cluster->amplitudes();
00829 uint16_t charge = 0;
00830 for (unsigned int i = 0; i < stripCharges.size(); ++i)
00831 {
00832 charge += stripCharges.at(i);
00833 }
00834
00835 strip_charge = (float)charge;
00836 strip_size = (int)( (cluster->amplitudes()).size() );
00837
00838
00839 float mindist = 999999.9;
00840 matched.clear();
00841
00842 matched = associate.associateHit(*hit2d);
00843
00844 strip_nsimhit = (int)matched.size();
00845
00846 if ( !matched.empty())
00847 {
00848 PSimHit closest;
00849
00850
00851
00852 for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
00853 {
00854 float dist = abs( (hit2d)->localPosition().x() - (*m).localPosition().x() );
00855
00856 if ( dist<mindist )
00857 {
00858 mindist = dist;
00859 closest = (*m);
00860 }
00861 }
00862
00863 strip_rechitresx = strip_rechitx - closest.localPosition().x();
00864 strip_rechitresy = strip_rechity - closest.localPosition().y();
00865
00866 strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
00867 strip_rechitpully = strip_rechitresy / strip_rechiterry;
00868
00869 strip_pidhit = closest.particleType();
00870 strip_simproc = closest.processType();
00871
00872
00873 }
00874
00875
00876
00877
00878
00879
00880
00881
00882
00883
00884
00885
00886 }
00887
00888
00889 ttree_track_hits_strip_->Fill();
00890
00891
00892 }
00893
00894 }
00895
00896
00897
00898
00899
00900
00901
00902
00903
00904
00905
00906
00907
00908
00909
00910
00911 edm::Handle<SiPixelRecHitCollection> recHitColl;
00912 e.getByLabel( "siPixelRecHits", recHitColl);
00913
00914 Handle<edm::SimTrackContainer> simtracks;
00915 e.getByLabel("g4SimHits", simtracks);
00916
00917
00918 for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin(); it != pDD->dets().end(); it++)
00919 {
00920 DetId detId = ((*it)->geographicalId());
00921
00922 SiPixelRecHitCollection::const_iterator dsmatch = recHitColl->find(detId);
00923 if (dsmatch == recHitColl->end()) continue;
00924
00925 SiPixelRecHitCollection::DetSet pixelrechitRange = *dsmatch;
00926 SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
00927 SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
00928 SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
00929 std::vector<PSimHit> matched;
00930
00931
00932 for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter)
00933 {
00934 matched.clear();
00935 matched = associate.associateHit(*pixeliter);
00936
00937
00938 if ( matched.empty() )
00939 {
00940 cout << "SiPixelErrorEstimation::analyze: rechits without associated simhit !!!!!!!"
00941 << endl;
00942 continue;
00943 }
00944
00945 all_subdetid = -9999;
00946
00947 all_layer = -9999;
00948 all_ladder = -9999;
00949 all_mod = -9999;
00950
00951 all_side = -9999;
00952 all_disk = -9999;
00953 all_blade = -9999;
00954 all_panel = -9999;
00955 all_plaq = -9999;
00956
00957 all_half = -9999;
00958 all_flipped = -9999;
00959
00960 all_cols = -9999;
00961 all_rows = -9999;
00962
00963 all_rechitx = -9999;
00964 all_rechity = -9999;
00965 all_rechitz = -9999;
00966
00967 all_simhitx = -9999;
00968 all_simhity = -9999;
00969
00970 all_rechiterrx = -9999;
00971 all_rechiterry = -9999;
00972
00973 all_rechitresx = -9999;
00974 all_rechitresy = -9999;
00975
00976 all_rechitpullx = -9999;
00977 all_rechitpully = -9999;
00978
00979 all_npix = -9999;
00980 all_nxpix = -9999;
00981 all_nypix = -9999;
00982
00983 all_edgex = -9999;
00984 all_edgey = -9999;
00985
00986 all_bigx = -9999;
00987 all_bigy = -9999;
00988
00989 all_alpha = -9999;
00990 all_beta = -9999;
00991
00992 all_simphi = -9999;
00993 all_simtheta = -9999;
00994
00995 all_simhitx = -9999;
00996 all_simhity = -9999;
00997
00998 all_nsimhit = -9999;
00999 all_pidhit = -9999;
01000 all_simproc = -9999;
01001
01002 all_vtxr = -9999;
01003 all_vtxz = -9999;
01004
01005 all_simpx = -9999;
01006 all_simpy = -9999;
01007 all_simpz = -9999;
01008
01009 all_eloss = -9999;
01010
01011 all_trkid = -9999;
01012
01013 all_x1 = -9999;
01014 all_x2 = -9999;
01015 all_y1 = -9999;
01016 all_y2 = -9999;
01017 all_z1 = -9999;
01018 all_z2 = -9999;
01019
01020 all_row1 = -9999;
01021 all_row2 = -9999;
01022 all_col1 = -9999;
01023 all_col2 = -9999;
01024
01025 all_gx1 = -9999;
01026 all_gx2 = -9999;
01027 all_gy1 = -9999;
01028 all_gy2 = -9999;
01029 all_gz1 = -9999;
01030 all_gz2 = -9999;
01031
01032 all_simtrketa = -9999;
01033 all_simtrkphi = -9999;
01034
01035 all_clust_row = -9999;
01036 all_clust_col = -9999;
01037
01038 all_clust_x = -9999;
01039 all_clust_y = -9999;
01040
01041 all_clust_q = -9999;
01042
01043 all_clust_maxpixcol = -9999;
01044 all_clust_maxpixrow = -9999;
01045 all_clust_minpixcol = -9999;
01046 all_clust_minpixrow = -9999;
01047
01048 all_clust_geoid = -9999;
01049
01050 all_clust_alpha = -9999;
01051 all_clust_beta = -9999;
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067 all_hit_probx = -9999;
01068 all_hit_proby = -9999;
01069 all_hit_cprob0 = -9999;
01070 all_hit_cprob1 = -9999;
01071 all_hit_cprob2 = -9999;
01072
01073 all_pixel_split = -9999;
01074 all_pixel_clst_err_x = -9999.9;
01075 all_pixel_clst_err_y = -9999.9;
01076
01077 all_nsimhit = (int)matched.size();
01078
01079 all_subdetid = (int)detId.subdetId();
01080
01081 if ( !(all_subdetid==1 || all_subdetid==2) )
01082 {
01083 cout << "SiPixelErrorEstimation::analyze: Not in a pixel detector !!!!!" << endl;
01084 continue;
01085 }
01086
01087 const PixelGeomDetUnit* theGeomDet
01088 = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(detId) );
01089
01090 const PixelTopology* topol = &(theGeomDet->specificTopology());
01091
01092 if ( pixeliter->cluster()->getSplitClusterErrorX() > 0.0 &&
01093 pixeliter->cluster()->getSplitClusterErrorY() > 0.0 )
01094 {
01095 all_pixel_split = 1;
01096 }
01097 else
01098 {
01099 all_pixel_split = 0;
01100 }
01101
01102 all_pixel_clst_err_x = pixeliter->cluster()->getSplitClusterErrorX();
01103 all_pixel_clst_err_y = pixeliter->cluster()->getSplitClusterErrorY();
01104
01105
01106 const int maxPixelCol = pixeliter->cluster()->maxPixelCol();
01107 const int maxPixelRow = pixeliter->cluster()->maxPixelRow();
01108 const int minPixelCol = pixeliter->cluster()->minPixelCol();
01109 const int minPixelRow = pixeliter->cluster()->minPixelRow();
01110
01111
01112
01113 all_hit_cprob0 = (float)pixeliter->clusterProbability(0);
01114 all_hit_cprob1 = (float)pixeliter->clusterProbability(1);
01115 all_hit_cprob2 = (float)pixeliter->clusterProbability(2);
01116
01117
01118 if ( topol->isItEdgePixelInX( minPixelRow ) ||
01119 topol->isItEdgePixelInX( maxPixelRow ) )
01120 all_edgex = 1;
01121 else
01122 all_edgex = 0;
01123
01124 if ( topol->isItEdgePixelInY( minPixelCol ) ||
01125 topol->isItEdgePixelInY( maxPixelCol ) )
01126 all_edgey = 1;
01127 else
01128 all_edgey = 0;
01129
01130
01131 if ( topol->containsBigPixelInX(minPixelRow, maxPixelRow) )
01132 all_bigx = 1;
01133 else
01134 all_bigx = 0;
01135
01136 if ( topol->containsBigPixelInY(minPixelCol, maxPixelCol) )
01137 all_bigy = 1;
01138 else
01139 all_bigy = 0;
01140
01141 if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel )
01142 {
01143
01144 all_layer = tTopo->pxbLayer(detId);
01145 all_ladder = tTopo->pxbLadder(detId);
01146 all_mod = tTopo->pxbModule(detId);
01147
01148 int tmp_nrows = theGeomDet->specificTopology().nrows();
01149 if ( tmp_nrows == 80 )
01150 all_half = 1;
01151 else if ( tmp_nrows == 160 )
01152 all_half = 0;
01153 else
01154 cout << "-------------------------------------------------- Wrong module size !!!" << endl;
01155
01156 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01157 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01158
01159 if ( tmp2<tmp1 )
01160 all_flipped = 1;
01161 else
01162 all_flipped = 0;
01163 }
01164 else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
01165 {
01166
01167 all_side = tTopo->pxfSide(detId);
01168 all_disk = tTopo->pxfDisk(detId);
01169 all_blade = tTopo->pxfBlade(detId);
01170 all_panel = tTopo->pxfPanel(detId);
01171 all_plaq = tTopo->pxfModule(detId);
01172
01173 }
01174 else std::cout << "We are not in the pixel detector" << (int)detId.subdetId() << endl;
01175
01176 all_cols = theGeomDet->specificTopology().ncolumns();
01177 all_rows = theGeomDet->specificTopology().nrows();
01178
01179 LocalPoint lp = pixeliter->localPosition();
01180
01181 all_rechitx = lp.x();
01182 all_rechity = lp.y();
01183 all_rechitz = lp.z();
01184
01185 LocalError le = pixeliter->localPositionError();
01186 all_rechiterrx = sqrt( le.xx() );
01187 all_rechiterry = sqrt( le.yy() );
01188
01189 bool found_hit_from_generated_particle = false;
01190
01191
01192 float closest_dist = 99999.9;
01193 std::vector<PSimHit>::const_iterator closest_simhit = matched.begin();
01194
01195 for (std::vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); m++)
01196 {
01197 if ( checkType_ )
01198 {
01199 int pid = (*m).particleType();
01200 if ( abs(pid) != genType_ )
01201 continue;
01202 }
01203
01204 float simhitx = 0.5 * ( (*m).entryPoint().x() + (*m).exitPoint().x() );
01205 float simhity = 0.5 * ( (*m).entryPoint().y() + (*m).exitPoint().y() );
01206
01207 float x_res = simhitx - rechitx;
01208 float y_res = simhity - rechity;
01209
01210 float dist = sqrt( x_res*x_res + y_res*y_res );
01211
01212 if ( dist < closest_dist )
01213 {
01214 closest_dist = dist;
01215 closest_simhit = m;
01216 found_hit_from_generated_particle = true;
01217 }
01218 }
01219
01220
01221
01222 if ( checkType_ && !found_hit_from_generated_particle )
01223 continue;
01224
01225 all_x1 = (*closest_simhit).entryPoint().x();
01226 all_y1 = (*closest_simhit).entryPoint().y();
01227 all_z1 = (*closest_simhit).entryPoint().z();
01228 all_x2 = (*closest_simhit).exitPoint().x();
01229 all_y2 = (*closest_simhit).exitPoint().y();
01230 all_z2 = (*closest_simhit).exitPoint().z();
01231 GlobalPoint GP1 =
01232 theGeomDet->surface().toGlobal( Local3DPoint( (*closest_simhit).entryPoint().x(),
01233 (*closest_simhit).entryPoint().y(),
01234 (*closest_simhit).entryPoint().z() ) );
01235 GlobalPoint GP2 =
01236 theGeomDet->surface().toGlobal (Local3DPoint( (*closest_simhit).exitPoint().x(),
01237 (*closest_simhit).exitPoint().y(),
01238 (*closest_simhit).exitPoint().z() ) );
01239 all_gx1 = GP1.x();
01240 all_gx2 = GP2.x();
01241 all_gy1 = GP1.y();
01242 all_gy2 = GP2.y();
01243 all_gz1 = GP1.z();
01244 all_gz2 = GP2.z();
01245
01246 MeasurementPoint mp1 =
01247 topol->measurementPosition( LocalPoint( (*closest_simhit).entryPoint().x(),
01248 (*closest_simhit).entryPoint().y(),
01249 (*closest_simhit).entryPoint().z() ) );
01250 MeasurementPoint mp2 =
01251 topol->measurementPosition( LocalPoint( (*closest_simhit).exitPoint().x(),
01252 (*closest_simhit).exitPoint().y(),
01253 (*closest_simhit).exitPoint().z() ) );
01254 all_row1 = mp1.x();
01255 all_col1 = mp1.y();
01256 all_row2 = mp2.x();
01257 all_col2 = mp2.y();
01258
01259 all_simhitx = 0.5*(all_x1+all_x2);
01260 all_simhity = 0.5*(all_y1+all_y2);
01261
01262 all_rechitresx = all_rechitx - all_simhitx;
01263 all_rechitresy = all_rechity - all_simhity;
01264
01265 all_rechitpullx = all_rechitresx / all_rechiterrx;
01266 all_rechitpully = all_rechitresy / all_rechiterry;
01267
01268 SiPixelRecHit::ClusterRef const& clust = pixeliter->cluster();
01269
01270 all_npix = clust->size();
01271 all_nxpix = clust->sizeX();
01272 all_nypix = clust->sizeY();
01273
01274 all_clust_row = clust->x();
01275 all_clust_col = clust->y();
01276
01277 LocalPoint lp2 = topol->localPosition( MeasurementPoint( all_clust_row, all_clust_col ) );
01278 all_clust_x = lp2.x();
01279 all_clust_y = lp2.y();
01280
01281 all_clust_q = clust->charge();
01282
01283 all_clust_maxpixcol = clust->maxPixelCol();
01284 all_clust_maxpixrow = clust->maxPixelRow();
01285 all_clust_minpixcol = clust->minPixelCol();
01286 all_clust_minpixrow = clust->minPixelRow();
01287
01288 all_clust_geoid = 0;
01289
01290 all_simpx = (*closest_simhit).momentumAtEntry().x();
01291 all_simpy = (*closest_simhit).momentumAtEntry().y();
01292 all_simpz = (*closest_simhit).momentumAtEntry().z();
01293 all_eloss = (*closest_simhit).energyLoss();
01294 all_simphi = (*closest_simhit).phiAtEntry();
01295 all_simtheta = (*closest_simhit).thetaAtEntry();
01296 all_pidhit = (*closest_simhit).particleType();
01297 all_trkid = (*closest_simhit).trackId();
01298
01299
01300 all_beta = atan2(all_simpz, all_simpy);
01301 all_alpha = atan2(all_simpz, all_simpx);
01302
01303 all_simproc = (int)closest_simhit->processType();
01304
01305 const edm::SimTrackContainer& trks = *(simtracks.product());
01306 SimTrackContainer::const_iterator trksiter;
01307 for (trksiter = trks.begin(); trksiter != trks.end(); trksiter++)
01308 if ( (int)trksiter->trackId() == all_trkid )
01309 {
01310 all_simtrketa = trksiter->momentum().eta();
01311 all_simtrkphi = trksiter->momentum().phi();
01312 }
01313
01314 all_vtxz = theGeomDet->surface().position().z();
01315 all_vtxr = theGeomDet->surface().position().perp();
01316
01317
01318
01319
01320
01321 const std::vector<SiPixelCluster::Pixel>& pixvector = clust->pixels();
01322 for ( int i=0; i<(int)pixvector.size(); ++i)
01323 {
01324 SiPixelCluster::Pixel holdpix = pixvector[i];
01325 all_pixrow[i] = holdpix.x;
01326 all_pixcol[i] = holdpix.y;
01327 all_pixadc[i] = holdpix.adc;
01328 LocalPoint lp = topol->localPosition(MeasurementPoint(holdpix.x, holdpix.y));
01329 all_pixx[i] = lp.x();
01330 all_pixy[i]= lp.y();
01331 GlobalPoint GP = theGeomDet->surface().toGlobal(Local3DPoint(lp.x(),lp.y(),lp.z()));
01332 all_pixgx[i] = GP.x();
01333 all_pixgy[i]= GP.y();
01334 all_pixgz[i]= GP.z();
01335 }
01336
01337 ttree_all_hits_->Fill();
01338
01339 }
01340
01341 }
01342
01343
01344
01345
01346
01347
01348
01349
01350
01351
01352
01353 if ( include_trk_hits_ )
01354 {
01355
01356 edm::Handle<reco::TrackCollection> trackCollection;
01357 e.getByLabel(src_, trackCollection);
01358 const reco::TrackCollection *tracks = trackCollection.product();
01359 reco::TrackCollection::const_iterator tciter;
01360
01361 if ( tracks->size() > 0 )
01362 {
01363
01364 for ( tciter=tracks->begin(); tciter!=tracks->end(); ++tciter)
01365 {
01366
01367 for ( trackingRecHit_iterator it = tciter->recHitsBegin(); it != tciter->recHitsEnd(); ++it)
01368 {
01369 const TrackingRecHit &trk_rec_hit = **it;
01370
01371 const SiPixelRecHit* matchedhit = dynamic_cast<const SiPixelRecHit*>(&trk_rec_hit);
01372
01373 if ( matchedhit )
01374 {
01375 rechitx = -9999.9;
01376 rechity = -9999.9;
01377 rechitz = -9999.9;
01378 rechiterrx = -9999.9;
01379 rechiterry = -9999.9;
01380 rechitresx = -9999.9;
01381 rechitresy = -9999.9;
01382 rechitpullx = -9999.9;
01383 rechitpully = -9999.9;
01384
01385 npix = -9999;
01386 nxpix = -9999;
01387 nypix = -9999;
01388 charge = -9999.9;
01389
01390 edgex = -9999;
01391 edgey = -9999;
01392
01393 bigx = -9999;
01394 bigy = -9999;
01395
01396 alpha = -9999.9;
01397 beta = -9999.9;
01398
01399 phi = -9999.9;
01400 eta = -9999.9;
01401
01402 subdetId = -9999;
01403
01404 layer = -9999;
01405 ladder = -9999;
01406 mod = -9999;
01407 side = -9999;
01408 disk = -9999;
01409 blade = -9999;
01410 panel = -9999;
01411 plaq = -9999;
01412
01413 half = -9999;
01414 flipped = -9999;
01415
01416 nsimhit = -9999;
01417 pidhit = -9999;
01418 simproc = -9999;
01419
01420 simhitx = -9999.9;
01421 simhity = -9999.9;
01422
01423 hit_probx = -9999.9;
01424 hit_proby = -9999.9;
01425 hit_cprob0 = -9999.9;
01426 hit_cprob1 = -9999.9;
01427 hit_cprob2 = -9999.9;
01428
01429 pixel_split = -9999;
01430
01431 pixel_clst_err_x = -9999.9;
01432 pixel_clst_err_y = -9999.9;
01433
01434 position = (*it)->localPosition();
01435 error = (*it)->localPositionError();
01436
01437 rechitx = position.x();
01438 rechity = position.y();
01439 rechitz = position.z();
01440 rechiterrx = sqrt(error.xx());
01441 rechiterry = sqrt(error.yy());
01442
01443 npix = matchedhit->cluster()->size();
01444 nxpix = matchedhit->cluster()->sizeX();
01445 nypix = matchedhit->cluster()->sizeY();
01446 charge = matchedhit->cluster()->charge();
01447
01448 if ( matchedhit->cluster()->getSplitClusterErrorX() > 0.0 &&
01449 matchedhit->cluster()->getSplitClusterErrorY() > 0.0 )
01450 pixel_split = 1;
01451 else
01452 pixel_split = 0;
01453
01454 pixel_clst_err_x = matchedhit->cluster()->getSplitClusterErrorX();
01455 pixel_clst_err_y = matchedhit->cluster()->getSplitClusterErrorY();
01456
01457
01458
01459 hit_cprob0 = (float)matchedhit->clusterProbability(0);
01460 hit_cprob1 = (float)matchedhit->clusterProbability(1);
01461 hit_cprob2 = (float)matchedhit->clusterProbability(2);
01462
01463
01464
01465 matched.clear();
01466 matched = associate.associateHit(*matchedhit);
01467
01468 nsimhit = (int)matched.size();
01469
01470 if ( !matched.empty() )
01471 {
01472 mindist = 999999.9;
01473 float distx, disty, dist;
01474 bool found_hit_from_generated_particle = false;
01475
01476 int n_assoc_muon = 0;
01477
01478 vector<PSimHit>::const_iterator closestit = matched.begin();
01479 for (vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
01480 {
01481 if ( checkType_ )
01482 {
01483 int pid = (*m).particleType();
01484 if ( abs(pid) != genType_ )
01485 continue;
01486 }
01487
01488 float simhitx = 0.5 * ( (*m).entryPoint().x() + (*m).exitPoint().x() );
01489 float simhity = 0.5 * ( (*m).entryPoint().y() + (*m).exitPoint().y() );
01490
01491 distx = fabs(rechitx - simhitx);
01492 disty = fabs(rechity - simhity);
01493 dist = sqrt( distx*distx + disty*disty );
01494
01495 if ( dist < mindist )
01496 {
01497 n_assoc_muon++;
01498
01499 mindist = dist;
01500 closestit = m;
01501 found_hit_from_generated_particle = true;
01502 }
01503 }
01504
01505
01506
01507 if ( checkType_ && !found_hit_from_generated_particle )
01508 continue;
01509
01510
01511
01512
01513
01514
01515
01516 DetId detId = (*it)->geographicalId();
01517
01518 const PixelGeomDetUnit* theGeomDet =
01519 dynamic_cast<const PixelGeomDetUnit*> ((*tracker).idToDet(detId) );
01520
01521 const PixelTopology* theTopol = &(theGeomDet->specificTopology());
01522
01523 pidhit = (*closestit).particleType();
01524 simproc = (int)(*closestit).processType();
01525
01526 simhitx = 0.5*( (*closestit).entryPoint().x() + (*closestit).exitPoint().x() );
01527 simhity = 0.5*( (*closestit).entryPoint().y() + (*closestit).exitPoint().y() );
01528
01529 rechitresx = rechitx - simhitx;
01530 rechitresy = rechity - simhity;
01531 rechitpullx = ( rechitx - simhitx ) / sqrt(error.xx());
01532 rechitpully = ( rechity - simhity ) / sqrt(error.yy());
01533
01534 float simhitpx = (*closestit).momentumAtEntry().x();
01535 float simhitpy = (*closestit).momentumAtEntry().y();
01536 float simhitpz = (*closestit).momentumAtEntry().z();
01537
01538 beta = atan2(simhitpz, simhitpy) * radtodeg;
01539 alpha = atan2(simhitpz, simhitpx) * radtodeg;
01540
01541
01542
01543
01544
01545 float locx = simhitpx;
01546 float locy = simhitpy;
01547 float locz = simhitpz;
01548
01549 bool isFlipped = false;
01550 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01551 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01552 if ( tmp2<tmp1 )
01553 isFlipped = true;
01554 else
01555 isFlipped = false;
01556
01557 trk_alpha = acos(locx/sqrt(locx*locx+locz*locz)) * radtodeg;
01558 if ( isFlipped )
01559 trk_alpha = 180.0 - trk_alpha ;
01560
01561 trk_beta = acos(locy/sqrt(locy*locy+locz*locz)) * radtodeg;
01562
01563
01564 phi = tciter->momentum().phi() / math_pi*180.0;
01565 eta = tciter->momentum().eta();
01566
01567 const int maxPixelCol = (*matchedhit).cluster()->maxPixelCol();
01568 const int maxPixelRow = (*matchedhit).cluster()->maxPixelRow();
01569 const int minPixelCol = (*matchedhit).cluster()->minPixelCol();
01570 const int minPixelRow = (*matchedhit).cluster()->minPixelRow();
01571
01572
01573 if ( theTopol->isItEdgePixelInX( minPixelRow ) ||
01574 theTopol->isItEdgePixelInX( maxPixelRow ) )
01575 edgex = 1;
01576 else
01577 edgex = 0;
01578
01579 if ( theTopol->isItEdgePixelInY( minPixelCol ) ||
01580 theTopol->isItEdgePixelInY( maxPixelCol ) )
01581 edgey = 1;
01582 else
01583 edgey = 0;
01584
01585
01586 if ( theTopol->containsBigPixelInX(minPixelRow, maxPixelRow) )
01587 bigx = 1;
01588 else
01589 bigx = 0;
01590
01591 if ( theTopol->containsBigPixelInY(minPixelCol, maxPixelCol) )
01592 bigy = 1;
01593 else
01594 bigy = 0;
01595
01596 subdetId = (int)detId.subdetId();
01597
01598 if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel )
01599 {
01600
01601 int tmp_nrows = theGeomDet->specificTopology().nrows();
01602 if ( tmp_nrows == 80 )
01603 half = 1;
01604 else if ( tmp_nrows == 160 )
01605 half = 0;
01606 else
01607 cout << "-------------------------------------------------- Wrong module size !!!" << endl;
01608
01609 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01610 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01611
01612 if ( tmp2<tmp1 )
01613 flipped = 1;
01614 else
01615 flipped = 0;
01616
01617
01618 layer = tTopo->pxbLayer(detId);
01619 ladder = tTopo->pxbLadder(detId);
01620 mod = tTopo->pxbModule(detId);
01621 }
01622 else if ( (int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap )
01623 {
01624
01625 side = tTopo->pxfSide(detId);
01626 disk = tTopo->pxfDisk(detId);
01627 blade = tTopo->pxfBlade(detId);
01628 panel = tTopo->pxfPanel(detId);
01629 plaq = tTopo->pxfModule(detId);
01630
01631 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
01632 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
01633
01634 if ( tmp2<tmp1 )
01635 flipped = 1;
01636 else
01637 flipped = 0;
01638
01639 }
01640
01641
01642 ttree_track_hits_->Fill();
01643
01644 }
01645 else
01646 cout << "---------------- RecHit with no associated SimHit !!! -------------------------- " << endl;
01647
01648 }
01649
01650 }
01651
01652 }
01653
01654 }
01655
01656 }
01657
01658
01659
01660
01661
01662 }
01663
01664 void SiPixelErrorEstimation::
01665 computeAnglesFromDetPosition(const SiPixelCluster & cl,
01666 const GeomDetUnit & det,
01667 float& alpha, float& beta )
01668 {
01669
01670 const PixelGeomDetUnit* theDet = dynamic_cast<const PixelGeomDetUnit*>( &det );
01671 if (! theDet)
01672 {
01673 cout << "---------------------------------------------- Not a pixel detector !!!!!!!!!!!!!!" << endl;
01674 assert(0);
01675 }
01676
01677 const PixelTopology* theTopol = &(theDet->specificTopology());
01678
01679
01680 float xcenter = cl.x();
01681 float ycenter = cl.y();
01682
01683
01684 LocalPoint lp = theTopol->localPosition( MeasurementPoint(xcenter, ycenter) );
01685
01686
01687
01688 GlobalPoint gp = theDet->surface().toGlobal( lp );
01689 float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() );
01690
01691
01692 float gpx = gp.x()/gp_mod;
01693 float gpy = gp.y()/gp_mod;
01694 float gpz = gp.z()/gp_mod;
01695
01696
01697
01698 GlobalVector gv(gpx, gpy, gpz);
01699
01700
01701 const Local3DVector lvx(1.0, 0.0, 0.0);
01702
01703
01704 GlobalVector gvx = theDet->surface().toGlobal( lvx );
01705
01706
01707 const Local3DVector lvy(0.0, 1.0, 0.0);
01708
01709
01710 GlobalVector gvy = theDet->surface().toGlobal( lvy );
01711
01712
01713 const Local3DVector lvz(0.0, 0.0, 1.0);
01714
01715
01716 GlobalVector gvz = theDet->surface().toGlobal( lvz );
01717
01718
01719
01720
01721 float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z();
01722 float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z();
01723 float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z();
01724
01725
01726 alpha = atan2( gv_dot_gvz, gv_dot_gvx );
01727 beta = atan2( gv_dot_gvz, gv_dot_gvy );
01728
01729
01730
01731
01732
01733
01734
01735 }
01736
01737
01738 DEFINE_FWK_MODULE(SiPixelErrorEstimation);
1.7.3