Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:32 

0001 // Package:          SiPixelErrorEstimation
0002 // Class:            SiPixelErrorEstimation
0003 // Original Author:  Gavril Giurgiu (JHU)
0004 // Created:          Fri May  4 17:48:24 CDT 2007
0005 
0006 #include <iostream>
0007 #include "CalibTracker/SiPixelErrorEstimation/interface/SiPixelErrorEstimation.h"
0008 
0009 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0010 #include "DataFormats/GeometryVector/interface/LocalVector.h"
0011 
0012 #include "DataFormats/TrackReco/interface/Track.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0015 
0016 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0017 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0018 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0019 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0020 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0021 
0022 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0023 
0024 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
0025 
0026 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0027 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0028 
0029 #include "CalibTracker/SiStripCommon/interface/SiStripDetInfoFileReader.h"
0030 #include "FWCore/ParameterSet/interface/FileInPath.h"
0031 
0032 using namespace std;
0033 using namespace edm;
0034 
0035 SiPixelErrorEstimation::SiPixelErrorEstimation(const edm::ParameterSet& ps)
0036     : tfile_(nullptr),
0037       ttree_all_hits_(nullptr),
0038       ttree_track_hits_(nullptr),
0039       ttree_track_hits_strip_(nullptr),
0040       trackerHitAssociatorConfig_(consumesCollector()) {
0041   //Read config file
0042   outputFile_ = ps.getUntrackedParameter<string>("outputFile", "SiPixelErrorEstimation_Ntuple.root");
0043 
0044   // Replace  "ctfWithMaterialTracks" with "generalTracks"
0045   //src_ = ps.getUntrackedParameter<std::string>( "src", "ctfWithMaterialTracks" );
0046   src_ = ps.getUntrackedParameter<std::string>("src", "generalTracks");
0047 
0048   checkType_ = ps.getParameter<bool>("checkType");
0049   genType_ = ps.getParameter<int>("genType");
0050   include_trk_hits_ = ps.getParameter<bool>("include_trk_hits");
0051 
0052   tTrajectory = consumes<std::vector<Trajectory>>(src_);
0053   tPixRecHitCollection = consumes<SiPixelRecHitCollection>(edm::InputTag("siPixelRecHits"));
0054   tSimTrackContainer = consumes<edm::SimTrackContainer>(edm::InputTag("g4SimHits"));
0055   tTrackCollection = consumes<reco::TrackCollection>(src_);
0056   trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0057   trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0058   magneticFieldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0059 }
0060 
0061 SiPixelErrorEstimation::~SiPixelErrorEstimation() {}
0062 
0063 void SiPixelErrorEstimation::beginJob() {
0064   int bufsize = 64000;
0065 
0066   tfile_ = new TFile(outputFile_.c_str(), "RECREATE");
0067 
0068   ttree_track_hits_strip_ = new TTree("TrackHitNtupleStrip", "TrackHitNtupleStrip");
0069 
0070   ttree_track_hits_strip_->Branch("strip_rechitx", &strip_rechitx, "strip_rechitx/F", bufsize);
0071   ttree_track_hits_strip_->Branch("strip_rechity", &strip_rechity, "strip_rechity/F", bufsize);
0072   ttree_track_hits_strip_->Branch("strip_rechitz", &strip_rechitz, "strip_rechitz/F", bufsize);
0073 
0074   ttree_track_hits_strip_->Branch("strip_rechiterrx", &strip_rechiterrx, "strip_rechiterrx/F", bufsize);
0075   ttree_track_hits_strip_->Branch("strip_rechiterry", &strip_rechiterry, "strip_rechiterry/F", bufsize);
0076 
0077   ttree_track_hits_strip_->Branch("strip_rechitresx", &strip_rechitresx, "strip_rechitresx/F", bufsize);
0078 
0079   ttree_track_hits_strip_->Branch("strip_rechitresx2", &strip_rechitresx2, "strip_rechitresx2/F", bufsize);
0080 
0081   ttree_track_hits_strip_->Branch("strip_rechitresy", &strip_rechitresy, "strip_rechitresy/F", bufsize);
0082 
0083   ttree_track_hits_strip_->Branch("strip_rechitpullx", &strip_rechitpullx, "strip_rechitpullx/F", bufsize);
0084   ttree_track_hits_strip_->Branch("strip_rechitpully", &strip_rechitpully, "strip_rechitpully/F", bufsize);
0085 
0086   ttree_track_hits_strip_->Branch("strip_is_stereo", &strip_is_stereo, "strip_is_stereo/I", bufsize);
0087   ttree_track_hits_strip_->Branch("strip_hit_type", &strip_hit_type, "strip_hit_type/I", bufsize);
0088   ttree_track_hits_strip_->Branch("detector_type", &detector_type, "detector_type/I", bufsize);
0089 
0090   ttree_track_hits_strip_->Branch("strip_trk_pt", &strip_trk_pt, "strip_trk_pt/F", bufsize);
0091   ttree_track_hits_strip_->Branch("strip_cotalpha", &strip_cotalpha, "strip_cotalpha/F", bufsize);
0092   ttree_track_hits_strip_->Branch("strip_cotbeta", &strip_cotbeta, "strip_cotbeta/F", bufsize);
0093   ttree_track_hits_strip_->Branch("strip_locbx", &strip_locbx, "strip_locbx/F", bufsize);
0094   ttree_track_hits_strip_->Branch("strip_locby", &strip_locby, "strip_locby/F", bufsize);
0095   ttree_track_hits_strip_->Branch("strip_locbz", &strip_locbz, "strip_locbz/F", bufsize);
0096   ttree_track_hits_strip_->Branch("strip_charge", &strip_charge, "strip_charge/F", bufsize);
0097   ttree_track_hits_strip_->Branch("strip_size", &strip_size, "strip_size/I", bufsize);
0098 
0099   ttree_track_hits_strip_->Branch("strip_edge", &strip_edge, "strip_edge/I", bufsize);
0100   ttree_track_hits_strip_->Branch("strip_nsimhit", &strip_nsimhit, "strip_nsimhit/I", bufsize);
0101   ttree_track_hits_strip_->Branch("strip_pidhit", &strip_pidhit, "strip_pidhit/I", bufsize);
0102   ttree_track_hits_strip_->Branch("strip_simproc", &strip_simproc, "strip_simproc/I", bufsize);
0103 
0104   ttree_track_hits_strip_->Branch("strip_subdet_id", &strip_subdet_id, "strip_subdet_id/I", bufsize);
0105 
0106   ttree_track_hits_strip_->Branch("strip_tib_layer", &strip_tib_layer, "strip_tib_layer/I", bufsize);
0107   ttree_track_hits_strip_->Branch("strip_tib_module", &strip_tib_module, "strip_tib_module/I", bufsize);
0108   ttree_track_hits_strip_->Branch("strip_tib_order", &strip_tib_order, "strip_tib_order/I", bufsize);
0109   ttree_track_hits_strip_->Branch("strip_tib_side", &strip_tib_side, "strip_tib_side/I", bufsize);
0110   ttree_track_hits_strip_->Branch(
0111       "strip_tib_is_double_side", &strip_tib_is_double_side, "strip_tib_is_double_side/I", bufsize);
0112   ttree_track_hits_strip_->Branch(
0113       "strip_tib_is_z_plus_side", &strip_tib_is_z_plus_side, "strip_tib_is_z_plus_side/I", bufsize);
0114   ttree_track_hits_strip_->Branch(
0115       "strip_tib_is_z_minus_side", &strip_tib_is_z_minus_side, "strip_tib_is_z_minus_side/I", bufsize);
0116   ttree_track_hits_strip_->Branch(
0117       "strip_tib_layer_number", &strip_tib_layer_number, "strip_tib_layer_number/I", bufsize);
0118   ttree_track_hits_strip_->Branch(
0119       "strip_tib_string_number", &strip_tib_string_number, "strip_tib_string_number/I", bufsize);
0120   ttree_track_hits_strip_->Branch(
0121       "strip_tib_module_number", &strip_tib_module_number, "strip_tib_module_number/I", bufsize);
0122   ttree_track_hits_strip_->Branch(
0123       "strip_tib_is_internal_string", &strip_tib_is_internal_string, "strip_tib_is_internal_string/I", bufsize);
0124   ttree_track_hits_strip_->Branch(
0125       "strip_tib_is_external_string", &strip_tib_is_external_string, "strip_tib_is_external_string/I", bufsize);
0126   ttree_track_hits_strip_->Branch("strip_tib_is_rphi", &strip_tib_is_rphi, "strip_tib_is_rphi/I", bufsize);
0127   ttree_track_hits_strip_->Branch("strip_tib_is_stereo", &strip_tib_is_stereo, "strip_tib_is_stereo/I", bufsize);
0128   ttree_track_hits_strip_->Branch("strip_tob_layer", &strip_tob_layer, "strip_tob_layer/I", bufsize);
0129   ttree_track_hits_strip_->Branch("strip_tob_module", &strip_tob_module, "strip_tob_module/I", bufsize);
0130   ttree_track_hits_strip_->Branch("strip_tob_side", &strip_tob_side, "strip_tob_side/I", bufsize);
0131   ttree_track_hits_strip_->Branch(
0132       "strip_tob_is_double_side", &strip_tob_is_double_side, "strip_tob_is_double_side/I", bufsize);
0133   ttree_track_hits_strip_->Branch(
0134       "strip_tob_is_z_plus_side", &strip_tob_is_z_plus_side, "strip_tob_is_z_plus_side/I", bufsize);
0135   ttree_track_hits_strip_->Branch(
0136       "strip_tob_is_z_minus_side", &strip_tob_is_z_minus_side, "strip_tob_is_z_minus_side/I", bufsize);
0137   ttree_track_hits_strip_->Branch(
0138       "strip_tob_layer_number", &strip_tob_layer_number, "strip_tob_layer_number/I", bufsize);
0139   ttree_track_hits_strip_->Branch("strip_tob_rod_number", &strip_tob_rod_number, "strip_tob_rod_number/I", bufsize);
0140   ttree_track_hits_strip_->Branch(
0141       "strip_tob_module_number", &strip_tob_module_number, "strip_tob_module_number/I", bufsize);
0142 
0143   ttree_track_hits_strip_->Branch("strip_prob", &strip_prob, "strip_prob/F", bufsize);
0144   ttree_track_hits_strip_->Branch("strip_qbin", &strip_qbin, "strip_qbin/I", bufsize);
0145 
0146   ttree_track_hits_strip_->Branch("strip_nprm", &strip_nprm, "strip_nprm/I", bufsize);
0147 
0148   ttree_track_hits_strip_->Branch("strip_pidhit1", &strip_pidhit1, "strip_pidhit1/I", bufsize);
0149   ttree_track_hits_strip_->Branch("strip_simproc1", &strip_simproc1, "strip_simproc1/I", bufsize);
0150 
0151   ttree_track_hits_strip_->Branch("strip_pidhit2", &strip_pidhit2, "strip_pidhit2/I", bufsize);
0152   ttree_track_hits_strip_->Branch("strip_simproc2", &strip_simproc2, "strip_simproc2/I", bufsize);
0153 
0154   ttree_track_hits_strip_->Branch("strip_pidhit3", &strip_pidhit3, "strip_pidhit3/I", bufsize);
0155   ttree_track_hits_strip_->Branch("strip_simproc3", &strip_simproc3, "strip_simproc3/I", bufsize);
0156 
0157   ttree_track_hits_strip_->Branch("strip_pidhit4", &strip_pidhit4, "strip_pidhit4/I", bufsize);
0158   ttree_track_hits_strip_->Branch("strip_simproc4", &strip_simproc4, "strip_simproc4/I", bufsize);
0159 
0160   ttree_track_hits_strip_->Branch("strip_pidhit5", &strip_pidhit5, "strip_pidhit5/I", bufsize);
0161   ttree_track_hits_strip_->Branch("strip_simproc5", &strip_simproc5, "strip_simproc5/I", bufsize);
0162 
0163   ttree_track_hits_strip_->Branch("strip_split", &strip_split, "strip_split/I", bufsize);
0164 
0165   ttree_track_hits_strip_->Branch("strip_clst_err_x", &strip_clst_err_x, "strip_clst_err_x/F", bufsize);
0166   ttree_track_hits_strip_->Branch("strip_clst_err_y", &strip_clst_err_y, "strip_clst_err_y/F", bufsize);
0167 
0168   if (include_trk_hits_) {
0169     //tfile_ = new TFile ("SiPixelErrorEstimation_Ntuple.root" , "RECREATE");
0170     //const char* tmp_name = outputFile_.c_str();
0171 
0172     ttree_track_hits_ = new TTree("TrackHitNtuple", "TrackHitNtuple");
0173 
0174     ttree_track_hits_->Branch("evt", &evt, "evt/I", bufsize);
0175     ttree_track_hits_->Branch("run", &run, "run/I", bufsize);
0176 
0177     ttree_track_hits_->Branch("subdetId", &subdetId, "subdetId/I", bufsize);
0178 
0179     ttree_track_hits_->Branch("layer", &layer, "layer/I", bufsize);
0180     ttree_track_hits_->Branch("ladder", &ladder, "ladder/I", bufsize);
0181     ttree_track_hits_->Branch("mod", &mod, "mod/I", bufsize);
0182 
0183     ttree_track_hits_->Branch("side", &side, "side/I", bufsize);
0184     ttree_track_hits_->Branch("disk", &disk, "disk/I", bufsize);
0185     ttree_track_hits_->Branch("blade", &blade, "blade/I", bufsize);
0186     ttree_track_hits_->Branch("panel", &panel, "panel/I", bufsize);
0187     ttree_track_hits_->Branch("plaq", &plaq, "plaq/I", bufsize);
0188 
0189     ttree_track_hits_->Branch("half", &half, "half/I", bufsize);
0190     ttree_track_hits_->Branch("flipped", &flipped, "flipped/I", bufsize);
0191 
0192     ttree_track_hits_->Branch("rechitx", &rechitx, "rechitx/F", bufsize);
0193     ttree_track_hits_->Branch("rechity", &rechity, "rechity/F", bufsize);
0194     ttree_track_hits_->Branch("rechitz", &rechitz, "rechitz/F", bufsize);
0195 
0196     ttree_track_hits_->Branch("rechiterrx", &rechiterrx, "rechiterrx/F", bufsize);
0197     ttree_track_hits_->Branch("rechiterry", &rechiterry, "rechiterry/F", bufsize);
0198 
0199     ttree_track_hits_->Branch("rechitresx", &rechitresx, "rechitresx/F", bufsize);
0200     ttree_track_hits_->Branch("rechitresy", &rechitresy, "rechitresy/F", bufsize);
0201 
0202     ttree_track_hits_->Branch("rechitpullx", &rechitpullx, "rechitpullx/F", bufsize);
0203     ttree_track_hits_->Branch("rechitpully", &rechitpully, "rechitpully/F", bufsize);
0204 
0205     ttree_track_hits_->Branch("npix", &npix, "npix/I", bufsize);
0206     ttree_track_hits_->Branch("nxpix", &nxpix, "nxpix/I", bufsize);
0207     ttree_track_hits_->Branch("nypix", &nypix, "nypix/I", bufsize);
0208     ttree_track_hits_->Branch("charge", &charge, "charge/F", bufsize);
0209 
0210     ttree_track_hits_->Branch("edgex", &edgex, "edgex/I", bufsize);
0211     ttree_track_hits_->Branch("edgey", &edgey, "edgey/I", bufsize);
0212 
0213     ttree_track_hits_->Branch("bigx", &bigx, "bigx/I", bufsize);
0214     ttree_track_hits_->Branch("bigy", &bigy, "bigy/I", bufsize);
0215 
0216     ttree_track_hits_->Branch("alpha", &alpha, "alpha/F", bufsize);
0217     ttree_track_hits_->Branch("beta", &beta, "beta/F", bufsize);
0218 
0219     ttree_track_hits_->Branch("trk_alpha", &trk_alpha, "trk_alpha/F", bufsize);
0220     ttree_track_hits_->Branch("trk_beta", &trk_beta, "trk_beta/F", bufsize);
0221 
0222     ttree_track_hits_->Branch("phi", &phi, "phi/F", bufsize);
0223     ttree_track_hits_->Branch("eta", &eta, "eta/F", bufsize);
0224 
0225     ttree_track_hits_->Branch("simhitx", &simhitx, "simhitx/F", bufsize);
0226     ttree_track_hits_->Branch("simhity", &simhity, "simhity/F", bufsize);
0227 
0228     ttree_track_hits_->Branch("nsimhit", &nsimhit, "nsimhit/I", bufsize);
0229     ttree_track_hits_->Branch("pidhit", &pidhit, "pidhit/I", bufsize);
0230     ttree_track_hits_->Branch("simproc", &simproc, "simproc/I", bufsize);
0231 
0232     ttree_track_hits_->Branch("pixel_split", &pixel_split, "pixel_split/I", bufsize);
0233 
0234     ttree_track_hits_->Branch("pixel_clst_err_x", &pixel_clst_err_x, "pixel_clst_err_x/F", bufsize);
0235     ttree_track_hits_->Branch("pixel_clst_err_y", &pixel_clst_err_y, "pixel_clst_err_y/F", bufsize);
0236 
0237   }  // if ( include_trk_hits_ )
0238 
0239   // ----------------------------------------------------------------------
0240 
0241   ttree_all_hits_ = new TTree("AllHitNtuple", "AllHitNtuple");
0242 
0243   ttree_all_hits_->Branch("evt", &evt, "evt/I", bufsize);
0244   ttree_all_hits_->Branch("run", &run, "run/I", bufsize);
0245 
0246   ttree_all_hits_->Branch("subdetid", &all_subdetid, "subdetid/I", bufsize);
0247 
0248   ttree_all_hits_->Branch("layer", &all_layer, "layer/I", bufsize);
0249   ttree_all_hits_->Branch("ladder", &all_ladder, "ladder/I", bufsize);
0250   ttree_all_hits_->Branch("mod", &all_mod, "mod/I", bufsize);
0251 
0252   ttree_all_hits_->Branch("side", &all_side, "side/I", bufsize);
0253   ttree_all_hits_->Branch("disk", &all_disk, "disk/I", bufsize);
0254   ttree_all_hits_->Branch("blade", &all_blade, "blade/I", bufsize);
0255   ttree_all_hits_->Branch("panel", &all_panel, "panel/I", bufsize);
0256   ttree_all_hits_->Branch("plaq", &all_plaq, "plaq/I", bufsize);
0257 
0258   ttree_all_hits_->Branch("half", &all_half, "half/I", bufsize);
0259   ttree_all_hits_->Branch("flipped", &all_flipped, "flipped/I", bufsize);
0260 
0261   ttree_all_hits_->Branch("cols", &all_cols, "cols/I", bufsize);
0262   ttree_all_hits_->Branch("rows", &all_rows, "rows/I", bufsize);
0263 
0264   ttree_all_hits_->Branch("rechitx", &all_rechitx, "rechitx/F", bufsize);
0265   ttree_all_hits_->Branch("rechity", &all_rechity, "rechity/F", bufsize);
0266   ttree_all_hits_->Branch("rechitz", &all_rechitz, "rechitz/F", bufsize);
0267 
0268   ttree_all_hits_->Branch("rechiterrx", &all_rechiterrx, "rechiterrx/F", bufsize);
0269   ttree_all_hits_->Branch("rechiterry", &all_rechiterry, "rechiterry/F", bufsize);
0270 
0271   ttree_all_hits_->Branch("rechitresx", &all_rechitresx, "rechitresx/F", bufsize);
0272   ttree_all_hits_->Branch("rechitresy", &all_rechitresy, "rechitresy/F", bufsize);
0273 
0274   ttree_all_hits_->Branch("rechitpullx", &all_rechitpullx, "rechitpullx/F", bufsize);
0275   ttree_all_hits_->Branch("rechitpully", &all_rechitpully, "rechitpully/F", bufsize);
0276 
0277   ttree_all_hits_->Branch("npix", &all_npix, "npix/I", bufsize);
0278   ttree_all_hits_->Branch("nxpix", &all_nxpix, "nxpix/I", bufsize);
0279   ttree_all_hits_->Branch("nypix", &all_nypix, "nypix/I", bufsize);
0280 
0281   ttree_all_hits_->Branch("edgex", &all_edgex, "edgex/I", bufsize);
0282   ttree_all_hits_->Branch("edgey", &all_edgey, "edgey/I", bufsize);
0283 
0284   ttree_all_hits_->Branch("bigx", &all_bigx, "bigx/I", bufsize);
0285   ttree_all_hits_->Branch("bigy", &all_bigy, "bigy/I", bufsize);
0286 
0287   ttree_all_hits_->Branch("alpha", &all_alpha, "alpha/F", bufsize);
0288   ttree_all_hits_->Branch("beta", &all_beta, "beta/F", bufsize);
0289 
0290   ttree_all_hits_->Branch("simhitx", &all_simhitx, "simhitx/F", bufsize);
0291   ttree_all_hits_->Branch("simhity", &all_simhity, "simhity/F", bufsize);
0292 
0293   ttree_all_hits_->Branch("nsimhit", &all_nsimhit, "nsimhit/I", bufsize);
0294   ttree_all_hits_->Branch("pidhit", &all_pidhit, "pidhit/I", bufsize);
0295   ttree_all_hits_->Branch("simproc", &all_simproc, "simproc/I", bufsize);
0296 
0297   ttree_all_hits_->Branch("vtxr", &all_vtxr, "vtxr/F", bufsize);
0298   ttree_all_hits_->Branch("vtxz", &all_vtxz, "vtxz/F", bufsize);
0299 
0300   ttree_all_hits_->Branch("simpx", &all_simpx, "simpx/F", bufsize);
0301   ttree_all_hits_->Branch("simpy", &all_simpy, "simpy/F", bufsize);
0302   ttree_all_hits_->Branch("simpz", &all_simpz, "simpz/F", bufsize);
0303 
0304   ttree_all_hits_->Branch("eloss", &all_eloss, "eloss/F", bufsize);
0305 
0306   ttree_all_hits_->Branch("simphi", &all_simphi, "simphi/F", bufsize);
0307   ttree_all_hits_->Branch("simtheta", &all_simtheta, "simtheta/F", bufsize);
0308 
0309   ttree_all_hits_->Branch("trkid", &all_trkid, "trkid/I", bufsize);
0310 
0311   ttree_all_hits_->Branch("x1", &all_x1, "x1/F", bufsize);
0312   ttree_all_hits_->Branch("x2", &all_x2, "x2/F", bufsize);
0313   ttree_all_hits_->Branch("y1", &all_x1, "y1/F", bufsize);
0314   ttree_all_hits_->Branch("y2", &all_x2, "y2/F", bufsize);
0315   ttree_all_hits_->Branch("z1", &all_x1, "z1/F", bufsize);
0316   ttree_all_hits_->Branch("z2", &all_x2, "z2/F", bufsize);
0317 
0318   ttree_all_hits_->Branch("row1", &all_row1, "row1/F", bufsize);
0319   ttree_all_hits_->Branch("row2", &all_row2, "row2/F", bufsize);
0320   ttree_all_hits_->Branch("col1", &all_col1, "col1/F", bufsize);
0321   ttree_all_hits_->Branch("col2", &all_col2, "col2/F", bufsize);
0322 
0323   ttree_all_hits_->Branch("gx1", &all_gx1, "gx1/F", bufsize);
0324   ttree_all_hits_->Branch("gx2", &all_gx2, "gx2/F", bufsize);
0325   ttree_all_hits_->Branch("gy1", &all_gx1, "gy1/F", bufsize);
0326   ttree_all_hits_->Branch("gy2", &all_gx2, "gy2/F", bufsize);
0327   ttree_all_hits_->Branch("gz1", &all_gx1, "gz1/F", bufsize);
0328   ttree_all_hits_->Branch("gz2", &all_gx2, "gz2/F", bufsize);
0329 
0330   ttree_all_hits_->Branch("simtrketa", &all_simtrketa, "simtrketa/F", bufsize);
0331   ttree_all_hits_->Branch("simtrkphi", &all_simtrkphi, "simtrkphi/F", bufsize);
0332 
0333   ttree_all_hits_->Branch("clust_row", &all_clust_row, "clust_row/F", bufsize);
0334   ttree_all_hits_->Branch("clust_col", &all_clust_col, "clust_col/F", bufsize);
0335 
0336   ttree_all_hits_->Branch("clust_x", &all_clust_x, "clust_x/F", bufsize);
0337   ttree_all_hits_->Branch("clust_y", &all_clust_y, "clust_y/F", bufsize);
0338 
0339   ttree_all_hits_->Branch("clust_q", &all_clust_q, "clust_q/F", bufsize);
0340 
0341   ttree_all_hits_->Branch("clust_maxpixcol", &all_clust_maxpixcol, "clust_maxpixcol/I", bufsize);
0342   ttree_all_hits_->Branch("clust_maxpixrow", &all_clust_maxpixrow, "clust_maxpixrow/I", bufsize);
0343   ttree_all_hits_->Branch("clust_minpixcol", &all_clust_minpixcol, "clust_minpixcol/I", bufsize);
0344   ttree_all_hits_->Branch("clust_minpixrow", &all_clust_minpixrow, "clust_minpixrow/I", bufsize);
0345 
0346   ttree_all_hits_->Branch("clust_geoid", &all_clust_geoid, "clust_geoid/I", bufsize);
0347 
0348   ttree_all_hits_->Branch("clust_alpha", &all_clust_alpha, "clust_alpha/F", bufsize);
0349   ttree_all_hits_->Branch("clust_beta", &all_clust_beta, "clust_beta/F", bufsize);
0350 
0351   ttree_all_hits_->Branch("rowpix", all_pixrow, "row[npix]/F", bufsize);
0352   ttree_all_hits_->Branch("colpix", all_pixcol, "col[npix]/F", bufsize);
0353   ttree_all_hits_->Branch("adc", all_pixadc, "adc[npix]/F", bufsize);
0354   ttree_all_hits_->Branch("xpix", all_pixx, "x[npix]/F", bufsize);
0355   ttree_all_hits_->Branch("ypix", all_pixy, "y[npix]/F", bufsize);
0356   ttree_all_hits_->Branch("gxpix", all_pixgx, "gx[npix]/F", bufsize);
0357   ttree_all_hits_->Branch("gypix", all_pixgy, "gy[npix]/F", bufsize);
0358   ttree_all_hits_->Branch("gzpix", all_pixgz, "gz[npix]/F", bufsize);
0359 
0360   ttree_all_hits_->Branch("hit_probx", &all_hit_probx, "hit_probx/F", bufsize);
0361   ttree_all_hits_->Branch("hit_proby", &all_hit_proby, "hit_proby/F", bufsize);
0362 
0363   ttree_all_hits_->Branch("all_pixel_split", &all_pixel_split, "all_pixel_split/I", bufsize);
0364 
0365   ttree_all_hits_->Branch("all_pixel_clst_err_x", &all_pixel_clst_err_x, "all_pixel_clst_err_x/F", bufsize);
0366   ttree_all_hits_->Branch("all_pixel_clst_err_y", &all_pixel_clst_err_y, "all_pixel_clst_err_y/F", bufsize);
0367 }
0368 
0369 void SiPixelErrorEstimation::endJob() {
0370   tfile_->Write();
0371   tfile_->Close();
0372 }
0373 
0374 void SiPixelErrorEstimation::analyze(const edm::Event& e, const edm::EventSetup& es) {
0375   //Retrieve tracker topology from geometry
0376   edm::ESHandle<TrackerTopology> tTopoHandle = es.getHandle(trackerTopoToken_);
0377   const TrackerTopology* const tTopo = tTopoHandle.product();
0378 
0379   using namespace edm;
0380 
0381   run = e.id().run();
0382   evt = e.id().event();
0383 
0384   if (evt % 1000 == 0)
0385     cout << "evt = " << evt << endl;
0386 
0387   float math_pi = 3.14159265;
0388   float radtodeg = 180.0 / math_pi;
0389 
0390   LocalPoint position;
0391   LocalError error;
0392   float mindist = 999999.9;
0393 
0394   std::vector<PSimHit> matched;
0395   TrackerHitAssociator associate(e, trackerHitAssociatorConfig_);
0396 
0397   edm::ESHandle<TrackerGeometry> pDD = es.getHandle(trackerGeomToken_);
0398   const TrackerGeometry* tracker = &(*pDD);
0399 
0400   //cout << "...1..." << endl;
0401 
0402   edm::ESHandle<MagneticField> magneticField = es.getHandle(magneticFieldToken_);
0403   //const MagneticField* magField_ = magFieldHandle.product();
0404 
0405   edm::FileInPath FileInPath_;
0406 
0407   // Strip hits ==============================================================================================================
0408 
0409   edm::Handle<vector<Trajectory>> trajCollectionHandle;
0410 
0411   e.getByToken(tTrajectory, trajCollectionHandle);
0412   //e.getByLabel( "generalTracks", trajCollectionHandle);
0413 
0414   for (vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it != trajCollectionHandle->end(); ++it) {
0415     vector<TrajectoryMeasurement> tmColl = it->measurements();
0416     for (vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj != tmColl.end(); ++itTraj) {
0417       if (!itTraj->updatedState().isValid())
0418         continue;
0419 
0420       strip_rechitx = -9999999.9;
0421       strip_rechity = -9999999.9;
0422       strip_rechitz = -9999999.9;
0423       strip_rechiterrx = -9999999.9;
0424       strip_rechiterry = -9999999.9;
0425       strip_rechitresx = -9999999.9;
0426       strip_rechitresx2 = -9999999.9;
0427 
0428       strip_rechitresy = -9999999.9;
0429       strip_rechitpullx = -9999999.9;
0430       strip_rechitpully = -9999999.9;
0431       strip_is_stereo = -9999999;
0432       strip_hit_type = -9999999;
0433       detector_type = -9999999;
0434 
0435       strip_trk_pt = -9999999.9;
0436       strip_cotalpha = -9999999.9;
0437       strip_cotbeta = -9999999.9;
0438       strip_locbx = -9999999.9;
0439       strip_locby = -9999999.9;
0440       strip_locbz = -9999999.9;
0441       strip_charge = -9999999.9;
0442       strip_size = -9999999;
0443 
0444       strip_edge = -9999999;
0445       strip_nsimhit = -9999999;
0446       strip_pidhit = -9999999;
0447       strip_simproc = -9999999;
0448 
0449       strip_subdet_id = -9999999;
0450 
0451       strip_tib_layer = -9999999;
0452       strip_tib_module = -9999999;
0453       strip_tib_order = -9999999;
0454       strip_tib_side = -9999999;
0455       strip_tib_is_double_side = -9999999;
0456       strip_tib_is_z_plus_side = -9999999;
0457       strip_tib_is_z_minus_side = -9999999;
0458       strip_tib_layer_number = -9999999;
0459       strip_tib_string_number = -9999999;
0460       strip_tib_module_number = -9999999;
0461       strip_tib_is_internal_string = -9999999;
0462       strip_tib_is_external_string = -9999999;
0463       strip_tib_is_rphi = -9999999;
0464       strip_tib_is_stereo = -9999999;
0465 
0466       strip_tob_layer = -9999999;
0467       strip_tob_module = -9999999;
0468       strip_tob_side = -9999999;
0469       strip_tob_is_double_side = -9999999;
0470       strip_tob_is_z_plus_side = -9999999;
0471       strip_tob_is_z_minus_side = -9999999;
0472       strip_tob_layer_number = -9999999;
0473       strip_tob_rod_number = -9999999;
0474       strip_tob_module_number = -9999999;
0475 
0476       strip_tob_is_rphi = -9999999;
0477       strip_tob_is_stereo = -9999999;
0478 
0479       strip_prob = -9999999.9;
0480       strip_qbin = -9999999;
0481 
0482       strip_nprm = -9999999;
0483 
0484       strip_pidhit1 = -9999999;
0485       strip_simproc1 = -9999999;
0486 
0487       strip_pidhit2 = -9999999;
0488       strip_simproc2 = -9999999;
0489 
0490       strip_pidhit3 = -9999999;
0491       strip_simproc3 = -9999999;
0492 
0493       strip_pidhit4 = -9999999;
0494       strip_simproc4 = -9999999;
0495 
0496       strip_pidhit5 = -9999999;
0497       strip_simproc5 = -9999999;
0498 
0499       strip_split = -9999999;
0500       strip_clst_err_x = -9999999.9;
0501       strip_clst_err_y = -9999999.9;
0502 
0503       const TransientTrackingRecHit::ConstRecHitPointer trans_trk_rec_hit_point = itTraj->recHit();
0504 
0505       if (trans_trk_rec_hit_point == nullptr)
0506         continue;
0507 
0508       const TrackingRecHit* trk_rec_hit = (*trans_trk_rec_hit_point).hit();
0509 
0510       if (trk_rec_hit == nullptr)
0511         continue;
0512 
0513       DetId detid = (trk_rec_hit)->geographicalId();
0514 
0515       strip_subdet_id = (int)detid.subdetId();
0516 
0517       if ((int)detid.subdetId() != (int)(StripSubdetector::TIB) &&
0518           (int)detid.subdetId() != (int)(StripSubdetector::TOB))
0519         continue;
0520 
0521       const SiStripMatchedRecHit2D* matchedhit =
0522           dynamic_cast<const SiStripMatchedRecHit2D*>((*trans_trk_rec_hit_point).hit());
0523       const SiStripRecHit2D* hit2d = dynamic_cast<const SiStripRecHit2D*>((*trans_trk_rec_hit_point).hit());
0524       const SiStripRecHit1D* hit1d = dynamic_cast<const SiStripRecHit1D*>((*trans_trk_rec_hit_point).hit());
0525 
0526       if (!matchedhit && !hit2d && !hit1d)
0527         continue;
0528 
0529       position = (trk_rec_hit)->localPosition();
0530       error = (trk_rec_hit)->localPositionError();
0531 
0532       strip_rechitx = position.x();
0533       strip_rechity = position.y();
0534       strip_rechitz = position.z();
0535       strip_rechiterrx = sqrt(error.xx());
0536       strip_rechiterry = sqrt(error.yy());
0537 
0538       //cout << "strip_rechitx = " << strip_rechitx << endl;
0539       //cout << "strip_rechity = " << strip_rechity << endl;
0540       //cout << "strip_rechitz = " << strip_rechitz << endl;
0541 
0542       //cout << "strip_rechiterrx = " << strip_rechiterrx << endl;
0543       //cout << "strip_rechiterry = " << strip_rechiterry << endl;
0544 
0545       TrajectoryStateOnSurface tsos = itTraj->updatedState();
0546 
0547       strip_trk_pt = tsos.globalMomentum().perp();
0548 
0549       LocalTrajectoryParameters ltp = tsos.localParameters();
0550 
0551       LocalVector localDir = ltp.momentum() / ltp.momentum().mag();
0552 
0553       float locx = localDir.x();
0554       float locy = localDir.y();
0555       float locz = localDir.z();
0556 
0557       //alpha_ = atan2( locz, locx );
0558       //beta_  = atan2( locz, locy );
0559 
0560       strip_cotalpha = locx / locz;
0561       strip_cotbeta = locy / locz;
0562 
0563       StripSubdetector StripSubdet = (StripSubdetector)detid;
0564 
0565       if (StripSubdet.stereo() == 0)
0566         strip_is_stereo = 0;
0567       else
0568         strip_is_stereo = 1;
0569 
0570       SiStripDetId si_strip_det_id = SiStripDetId(detid);
0571 
0572       //const StripGeomDetUnit* strip_geom_det_unit = dynamic_cast<const StripGeomDetUnit*> ( tracker->idToDet( detid ) );
0573       const StripGeomDetUnit* strip_geom_det_unit = (const StripGeomDetUnit*)tracker->idToDetUnit(detid);
0574 
0575       if (strip_geom_det_unit != nullptr) {
0576         LocalVector lbfield = (strip_geom_det_unit->surface())
0577                                   .toLocal((*magneticField).inTesla(strip_geom_det_unit->surface().position()));
0578 
0579         strip_locbx = lbfield.x();
0580         strip_locby = lbfield.y();
0581         strip_locbz = lbfield.z();
0582       }
0583 
0584       //enum ModuleGeometry {UNKNOWNGEOMETRY, IB1, IB2, OB1, OB2, W1A, W2A, W3A, W1B, W2B, W3B, W4, W5, W6, W7};
0585 
0586       if (si_strip_det_id.moduleGeometry() == SiStripModuleGeometry::IB1) {
0587         detector_type = 1;
0588         //cout << "si_strip_det_id.moduleGeometry() = IB1" << endl;
0589         //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
0590       } else if (si_strip_det_id.moduleGeometry() == SiStripModuleGeometry::IB2) {
0591         detector_type = 2;
0592         //cout << "si_strip_det_id.moduleGeometry() = IB2" << endl;
0593         //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
0594       } else if (si_strip_det_id.moduleGeometry() == SiStripModuleGeometry::OB1) {
0595         detector_type = 3;
0596         //cout << "si_strip_det_id.moduleGeometry() = OB1" << endl;
0597         //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
0598       } else if (si_strip_det_id.moduleGeometry() == SiStripModuleGeometry::OB2) {
0599         detector_type = 4;
0600         //cout << "si_strip_det_id.moduleGeometry() = OB2" << endl;
0601         //cout << "si_strip_det_id.moduleGeometry() = " << si_strip_det_id.moduleGeometry() << endl;
0602       }
0603 
0604       // Store ntuple variables
0605 
0606       if ((int)detid.subdetId() == int(StripSubdetector::TIB)) {
0607         strip_tib_layer = (int)tTopo->tibLayer(detid);
0608         strip_tib_module = (int)tTopo->tibModule(detid);
0609         strip_tib_order = (int)tTopo->tibOrder(detid);
0610         strip_tib_side = (int)tTopo->tibSide(detid);
0611         strip_tib_is_double_side = (int)tTopo->tibIsDoubleSide(detid);
0612         strip_tib_is_z_plus_side = (int)tTopo->tibIsZPlusSide(detid);
0613         strip_tib_is_z_minus_side = (int)tTopo->tibIsZMinusSide(detid);
0614         strip_tib_layer_number = (int)tTopo->tibLayer(detid);
0615         strip_tib_string_number = (int)tTopo->tibString(detid);
0616         strip_tib_module_number = (int)tTopo->tibModule(detid);
0617         strip_tib_is_internal_string = (int)tTopo->tibIsInternalString(detid);
0618         strip_tib_is_external_string = (int)tTopo->tibIsExternalString(detid);
0619         strip_tib_is_rphi = (int)tTopo->tibIsRPhi(detid);
0620         strip_tib_is_stereo = (int)tTopo->tibIsStereo(detid);
0621       }
0622 
0623       if ((int)detid.subdetId() == int(StripSubdetector::TOB)) {
0624         strip_tob_layer = (int)tTopo->tobLayer(detid);
0625         strip_tob_module = (int)tTopo->tobModule(detid);
0626 
0627         strip_tob_side = (int)tTopo->tobSide(detid);
0628         strip_tob_is_double_side = (int)tTopo->tobIsDoubleSide(detid);
0629         strip_tob_is_z_plus_side = (int)tTopo->tobIsZPlusSide(detid);
0630         strip_tob_is_z_minus_side = (int)tTopo->tobIsZMinusSide(detid);
0631         strip_tob_layer_number = (int)tTopo->tobLayer(detid);
0632         strip_tob_rod_number = (int)tTopo->tobRod(detid);
0633         strip_tob_module_number = (int)tTopo->tobModule(detid);
0634 
0635         strip_tob_is_rphi = (int)tTopo->tobIsRPhi(detid);
0636         strip_tob_is_stereo = (int)tTopo->tobIsStereo(detid);
0637       }
0638 
0639       if (matchedhit) {
0640         //cout << endl << endl << endl;
0641         //cout << "Found a SiStripMatchedRecHit2D !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endl;
0642         //cout << endl << endl << endl;
0643         strip_hit_type = 0;
0644 
0645       }  //  if ( matchedhit )
0646 
0647       if (hit1d) {
0648         strip_hit_type = 1;
0649 
0650         const SiStripRecHit1D::ClusterRef& cluster = hit1d->cluster();
0651 
0652         if (cluster->getSplitClusterError() > 0.0)
0653           strip_split = 1;
0654         else
0655           strip_split = 0;
0656 
0657         strip_clst_err_x = cluster->getSplitClusterError();
0658         //strip_clst_err_y = ...
0659 
0660         // Get cluster total charge
0661         const auto& stripCharges = cluster->amplitudes();
0662         uint16_t charge = 0;
0663         for (unsigned int i = 0; i < stripCharges.size(); ++i) {
0664           charge += stripCharges[i];
0665         }
0666 
0667         strip_charge = (float)charge;
0668         strip_size = (int)((cluster->amplitudes()).size());
0669 
0670         // Association of the rechit to the simhit
0671         float mindist = 999999.9;
0672         matched.clear();
0673 
0674         matched = associate.associateHit(*hit1d);
0675 
0676         strip_nsimhit = (int)matched.size();
0677 
0678         if (!matched.empty()) {
0679           PSimHit closest;
0680 
0681           // Get the closest simhit
0682 
0683           int strip_nprimaries = 0;
0684           int current_index = 0;
0685 
0686           for (vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); ++m) {
0687             ++current_index;
0688 
0689             if ((*m).processType() == 2)
0690               ++strip_nprimaries;
0691 
0692             if (current_index == 1) {
0693               strip_pidhit1 = (*m).particleType();
0694               strip_simproc1 = (*m).processType();
0695             } else if (current_index == 2) {
0696               strip_pidhit2 = (*m).particleType();
0697               strip_simproc2 = (*m).processType();
0698             } else if (current_index == 3) {
0699               strip_pidhit3 = (*m).particleType();
0700               strip_simproc3 = (*m).processType();
0701             } else if (current_index == 4) {
0702               strip_pidhit4 = (*m).particleType();
0703               strip_simproc4 = (*m).processType();
0704             } else if (current_index == 5) {
0705               strip_pidhit5 = (*m).particleType();
0706               strip_simproc5 = (*m).processType();
0707             }
0708 
0709             float dist = abs((hit1d)->localPosition().x() - (*m).localPosition().x());
0710 
0711             if (dist < mindist) {
0712               mindist = dist;
0713               closest = (*m);
0714             }
0715 
0716           }  // for ( vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); ++m)
0717 
0718           strip_nprm = strip_nprimaries;
0719 
0720           strip_rechitresx = strip_rechitx - closest.localPosition().x();
0721           strip_rechitresy = strip_rechity - closest.localPosition().y();
0722 
0723           strip_rechitresx2 = strip_rechitx - matched[0].localPosition().x();
0724 
0725           strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
0726           strip_rechitpully = strip_rechitresy / strip_rechiterry;
0727 
0728           strip_pidhit = closest.particleType();
0729           strip_simproc = closest.processType();
0730 
0731         }  //   if( !matched.empty())
0732 
0733         //strip_prob = (hit1d)->getTemplProb();
0734         //strip_qbin = (hit1d)->getTemplQbin();
0735 
0736         //cout << endl;
0737         //cout << "SiPixelErrorEstimation 1d hit: " << endl;
0738         //cout << "prob 1d = " << strip_prob << endl;
0739         //cout << "qbin 1d = " << strip_qbin << endl;
0740         //cout << endl;
0741 
0742         // Is the cluster on edge ?
0743         /*
0744                 const auto stripDetInfo = SiStripDetInfoFileReader::read(edm::FileInPath{SiStripDetInfoFileReader::kDefaultFile}.fullPath());
0745         
0746         uint16_t firstStrip = cluster->firstStrip();
0747         uint16_t lastStrip = firstStrip + (cluster->amplitudes()).size() -1;
0748         unsigned short Nstrips;
0749         Nstrips = stripDetInfo.getNumberOfApvsAndStripLength(id1).first*128;
0750         
0751         if ( firstStrip == 0 || lastStrip == (Nstrips-1) ) 
0752         strip_edge = 1;
0753         else
0754         strip_edge = 0;
0755           */
0756 
0757       }  // if ( hit1d )
0758 
0759       if (hit2d) {
0760         strip_hit_type = 2;
0761 
0762         const SiStripRecHit1D::ClusterRef& cluster = hit2d->cluster();
0763 
0764         //if ( cluster->getSplitClusterError()  > 0.0 )
0765         //strip_split = 1;
0766         //else
0767         //strip_split = 0;
0768 
0769         //strip_clst_err_x = cluster->getSplitClusterError();
0770         //strip_clst_err_y = ...
0771 
0772         // Get cluster total charge
0773         auto& stripCharges = cluster->amplitudes();
0774         uint16_t charge = 0;
0775         for (unsigned int i = 0; i < stripCharges.size(); ++i) {
0776           charge += stripCharges[i];
0777         }
0778 
0779         strip_charge = (float)charge;
0780         strip_size = (int)((cluster->amplitudes()).size());
0781 
0782         // Association of the rechit to the simhit
0783         float mindist = 999999.9;
0784         matched.clear();
0785 
0786         matched = associate.associateHit(*hit2d);
0787 
0788         strip_nsimhit = (int)matched.size();
0789 
0790         if (!matched.empty()) {
0791           PSimHit closest;
0792 
0793           // Get the closest simhit
0794 
0795           for (vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); ++m) {
0796             float dist = abs((hit2d)->localPosition().x() - (*m).localPosition().x());
0797 
0798             if (dist < mindist) {
0799               mindist = dist;
0800               closest = (*m);
0801             }
0802           }
0803 
0804           strip_rechitresx = strip_rechitx - closest.localPosition().x();
0805           strip_rechitresy = strip_rechity - closest.localPosition().y();
0806 
0807           strip_rechitpullx = strip_rechitresx / strip_rechiterrx;
0808           strip_rechitpully = strip_rechitresy / strip_rechiterry;
0809 
0810           strip_pidhit = closest.particleType();
0811           strip_simproc = closest.processType();
0812 
0813         }  //   if( !matched.empty())
0814 
0815         //strip_prob = (hit2d)->getTemplProb();
0816         //strip_qbin = (hit2d)->getTemplQbin();
0817 
0818         //cout << endl;
0819         //cout << "SiPixelErrorEstimation 2d hit: " << endl;
0820         //cout << "prob 2d = " << strip_prob << endl;
0821         //cout << "qbin 2d = " << strip_qbin << endl;
0822         //cout << endl;
0823 
0824       }  // if ( hit2d )
0825 
0826       ttree_track_hits_strip_->Fill();
0827 
0828     }  // for ( vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj!=tmColl.end(); ++itTraj )
0829 
0830   }  // for( vector<Trajectory>::const_iterator it = trajCollectionHandle->begin(); it!=trajCollectionHandle->end(); ++it)
0831 
0832   //cout << "...2..." << endl;
0833 
0834   // --------------------------------------- all hits -----------------------------------------------------------
0835   edm::Handle<SiPixelRecHitCollection> recHitColl;
0836   e.getByToken(tPixRecHitCollection, recHitColl);
0837 
0838   Handle<edm::SimTrackContainer> simtracks;
0839   e.getByToken(tSimTrackContainer, simtracks);
0840 
0841   //-----Iterate over detunits
0842   for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin(); it != pDD->dets().end(); it++) {
0843     DetId detId = ((*it)->geographicalId());
0844 
0845     SiPixelRecHitCollection::const_iterator dsmatch = recHitColl->find(detId);
0846     if (dsmatch == recHitColl->end())
0847       continue;
0848 
0849     SiPixelRecHitCollection::DetSet pixelrechitRange = *dsmatch;
0850     SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorBegin = pixelrechitRange.begin();
0851     SiPixelRecHitCollection::DetSet::const_iterator pixelrechitRangeIteratorEnd = pixelrechitRange.end();
0852     SiPixelRecHitCollection::DetSet::const_iterator pixeliter = pixelrechitRangeIteratorBegin;
0853     std::vector<PSimHit> matched;
0854 
0855     //----Loop over rechits for this detId
0856     for (; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter) {
0857       matched.clear();
0858       matched = associate.associateHit(*pixeliter);
0859       // only consider rechits that have associated simhit
0860       // otherwise cannot determine residiual
0861       if (matched.empty()) {
0862         cout << "SiPixelErrorEstimation::analyze: rechits without associated simhit !!!!!!!" << endl;
0863         continue;
0864       }
0865 
0866       all_subdetid = -9999;
0867 
0868       all_layer = -9999;
0869       all_ladder = -9999;
0870       all_mod = -9999;
0871 
0872       all_side = -9999;
0873       all_disk = -9999;
0874       all_blade = -9999;
0875       all_panel = -9999;
0876       all_plaq = -9999;
0877 
0878       all_half = -9999;
0879       all_flipped = -9999;
0880 
0881       all_cols = -9999;
0882       all_rows = -9999;
0883 
0884       all_rechitx = -9999;
0885       all_rechity = -9999;
0886       all_rechitz = -9999;
0887 
0888       all_simhitx = -9999;
0889       all_simhity = -9999;
0890 
0891       all_rechiterrx = -9999;
0892       all_rechiterry = -9999;
0893 
0894       all_rechitresx = -9999;
0895       all_rechitresy = -9999;
0896 
0897       all_rechitpullx = -9999;
0898       all_rechitpully = -9999;
0899 
0900       all_npix = -9999;
0901       all_nxpix = -9999;
0902       all_nypix = -9999;
0903 
0904       all_edgex = -9999;
0905       all_edgey = -9999;
0906 
0907       all_bigx = -9999;
0908       all_bigy = -9999;
0909 
0910       all_alpha = -9999;
0911       all_beta = -9999;
0912 
0913       all_simphi = -9999;
0914       all_simtheta = -9999;
0915 
0916       all_simhitx = -9999;
0917       all_simhity = -9999;
0918 
0919       all_nsimhit = -9999;
0920       all_pidhit = -9999;
0921       all_simproc = -9999;
0922 
0923       all_vtxr = -9999;
0924       all_vtxz = -9999;
0925 
0926       all_simpx = -9999;
0927       all_simpy = -9999;
0928       all_simpz = -9999;
0929 
0930       all_eloss = -9999;
0931 
0932       all_trkid = -9999;
0933 
0934       all_x1 = -9999;
0935       all_x2 = -9999;
0936       all_y1 = -9999;
0937       all_y2 = -9999;
0938       all_z1 = -9999;
0939       all_z2 = -9999;
0940 
0941       all_row1 = -9999;
0942       all_row2 = -9999;
0943       all_col1 = -9999;
0944       all_col2 = -9999;
0945 
0946       all_gx1 = -9999;
0947       all_gx2 = -9999;
0948       all_gy1 = -9999;
0949       all_gy2 = -9999;
0950       all_gz1 = -9999;
0951       all_gz2 = -9999;
0952 
0953       all_simtrketa = -9999;
0954       all_simtrkphi = -9999;
0955 
0956       all_clust_row = -9999;
0957       all_clust_col = -9999;
0958 
0959       all_clust_x = -9999;
0960       all_clust_y = -9999;
0961 
0962       all_clust_q = -9999;
0963 
0964       all_clust_maxpixcol = -9999;
0965       all_clust_maxpixrow = -9999;
0966       all_clust_minpixcol = -9999;
0967       all_clust_minpixrow = -9999;
0968 
0969       all_clust_geoid = -9999;
0970 
0971       all_clust_alpha = -9999;
0972       all_clust_beta = -9999;
0973 
0974       /*
0975         for (int i=0; i<all_npix; ++i)
0976         {
0977         all_pixrow[i] = -9999;
0978         all_pixcol[i] = -9999;
0979         all_pixadc[i] = -9999;
0980         all_pixx[i] = -9999;
0981         all_pixy[i] = -9999;
0982         all_pixgx[i] = -9999;
0983         all_pixgy[i] = -9999;
0984         all_pixgz[i] = -9999;
0985         }
0986       */
0987 
0988       all_hit_probx = -9999;
0989       all_hit_proby = -9999;
0990       all_hit_cprob0 = -9999;
0991       all_hit_cprob1 = -9999;
0992       all_hit_cprob2 = -9999;
0993 
0994       all_pixel_split = -9999;
0995       all_pixel_clst_err_x = -9999.9;
0996       all_pixel_clst_err_y = -9999.9;
0997 
0998       all_nsimhit = (int)matched.size();
0999 
1000       all_subdetid = (int)detId.subdetId();
1001       // only consider rechits in pixel barrel and pixel forward
1002       if (!(all_subdetid == 1 || all_subdetid == 2)) {
1003         cout << "SiPixelErrorEstimation::analyze: Not in a pixel detector !!!!!" << endl;
1004         continue;
1005       }
1006 
1007       const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>(tracker->idToDet(detId));
1008 
1009       const PixelTopology* topol = &(theGeomDet->specificTopology());
1010 
1011       if (pixeliter->cluster()->getSplitClusterErrorX() > 0.0 && pixeliter->cluster()->getSplitClusterErrorY() > 0.0) {
1012         all_pixel_split = 1;
1013       } else {
1014         all_pixel_split = 0;
1015       }
1016 
1017       all_pixel_clst_err_x = pixeliter->cluster()->getSplitClusterErrorX();
1018       all_pixel_clst_err_y = pixeliter->cluster()->getSplitClusterErrorY();
1019 
1020       const int maxPixelCol = pixeliter->cluster()->maxPixelCol();
1021       const int maxPixelRow = pixeliter->cluster()->maxPixelRow();
1022       const int minPixelCol = pixeliter->cluster()->minPixelCol();
1023       const int minPixelRow = pixeliter->cluster()->minPixelRow();
1024 
1025       //all_hit_probx  = (float)pixeliter->probabilityX();
1026       //all_hit_proby  = (float)pixeliter->probabilityY();
1027       all_hit_cprob0 = (float)pixeliter->clusterProbability(0);
1028       all_hit_cprob1 = (float)pixeliter->clusterProbability(1);
1029       all_hit_cprob2 = (float)pixeliter->clusterProbability(2);
1030 
1031       // check whether the cluster is at the module edge
1032       if (topol->isItEdgePixelInX(minPixelRow) || topol->isItEdgePixelInX(maxPixelRow))
1033         all_edgex = 1;
1034       else
1035         all_edgex = 0;
1036 
1037       if (topol->isItEdgePixelInY(minPixelCol) || topol->isItEdgePixelInY(maxPixelCol))
1038         all_edgey = 1;
1039       else
1040         all_edgey = 0;
1041 
1042       // check whether this rechit contains big pixels
1043       if (topol->containsBigPixelInX(minPixelRow, maxPixelRow))
1044         all_bigx = 1;
1045       else
1046         all_bigx = 0;
1047 
1048       if (topol->containsBigPixelInY(minPixelCol, maxPixelCol))
1049         all_bigy = 1;
1050       else
1051         all_bigy = 0;
1052 
1053       if ((int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel) {
1054         all_layer = tTopo->pxbLayer(detId);
1055         all_ladder = tTopo->pxbLadder(detId);
1056         all_mod = tTopo->pxbModule(detId);
1057 
1058         int tmp_nrows = theGeomDet->specificTopology().nrows();
1059         if (tmp_nrows == 80)
1060           all_half = 1;
1061         else if (tmp_nrows == 160)
1062           all_half = 0;
1063         else
1064           cout << "-------------------------------------------------- Wrong module size !!!" << endl;
1065 
1066         float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 0.)).perp();
1067         float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 1.)).perp();
1068 
1069         if (tmp2 < tmp1)
1070           all_flipped = 1;
1071         else
1072           all_flipped = 0;
1073       } else if ((int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap) {
1074         all_side = tTopo->pxfSide(detId);
1075         all_disk = tTopo->pxfDisk(detId);
1076         all_blade = tTopo->pxfBlade(detId);
1077         all_panel = tTopo->pxfPanel(detId);
1078         all_plaq = tTopo->pxfModule(detId);  // also known as plaquette
1079 
1080       }  // else if ( detId.subdetId()==PixelSubdetector::PixelEndcap )
1081       else
1082         std::cout << "We are not in the pixel detector" << (int)detId.subdetId() << endl;
1083 
1084       all_cols = theGeomDet->specificTopology().ncolumns();
1085       all_rows = theGeomDet->specificTopology().nrows();
1086 
1087       LocalPoint lp = pixeliter->localPosition();
1088       // gavril: change this name
1089       all_rechitx = lp.x();
1090       all_rechity = lp.y();
1091       all_rechitz = lp.z();
1092 
1093       LocalError le = pixeliter->localPositionError();
1094       all_rechiterrx = sqrt(le.xx());
1095       all_rechiterry = sqrt(le.yy());
1096 
1097       bool found_hit_from_generated_particle = false;
1098 
1099       //---Loop over sim hits, fill closest
1100       float closest_dist = 99999.9;
1101       std::vector<PSimHit>::const_iterator closest_simhit = matched.begin();
1102 
1103       for (std::vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); m++) {
1104         if (checkType_) {
1105           int pid = (*m).particleType();
1106           if (abs(pid) != genType_)
1107             continue;
1108         }
1109 
1110         float simhitx = 0.5 * ((*m).entryPoint().x() + (*m).exitPoint().x());
1111         float simhity = 0.5 * ((*m).entryPoint().y() + (*m).exitPoint().y());
1112 
1113         float x_res = simhitx - rechitx;
1114         float y_res = simhity - rechity;
1115 
1116         float dist = sqrt(x_res * x_res + y_res * y_res);
1117 
1118         if (dist < closest_dist) {
1119           closest_dist = dist;
1120           closest_simhit = m;
1121           found_hit_from_generated_particle = true;
1122         }
1123       }  // end sim hit loop
1124 
1125       // If this recHit does not have any simHit with the same particleType as the particles generated
1126       // ignore it as most probably comes from delta rays.
1127       if (checkType_ && !found_hit_from_generated_particle)
1128         continue;
1129 
1130       all_x1 = (*closest_simhit).entryPoint().x();  // width (row index, in col direction)
1131       all_y1 = (*closest_simhit).entryPoint().y();  // length (col index, in row direction)
1132       all_z1 = (*closest_simhit).entryPoint().z();
1133       all_x2 = (*closest_simhit).exitPoint().x();
1134       all_y2 = (*closest_simhit).exitPoint().y();
1135       all_z2 = (*closest_simhit).exitPoint().z();
1136       GlobalPoint GP1 = theGeomDet->surface().toGlobal(Local3DPoint(
1137           (*closest_simhit).entryPoint().x(), (*closest_simhit).entryPoint().y(), (*closest_simhit).entryPoint().z()));
1138       GlobalPoint GP2 = theGeomDet->surface().toGlobal(Local3DPoint(
1139           (*closest_simhit).exitPoint().x(), (*closest_simhit).exitPoint().y(), (*closest_simhit).exitPoint().z()));
1140       all_gx1 = GP1.x();
1141       all_gx2 = GP2.x();
1142       all_gy1 = GP1.y();
1143       all_gy2 = GP2.y();
1144       all_gz1 = GP1.z();
1145       all_gz2 = GP2.z();
1146 
1147       MeasurementPoint mp1 = topol->measurementPosition(LocalPoint(
1148           (*closest_simhit).entryPoint().x(), (*closest_simhit).entryPoint().y(), (*closest_simhit).entryPoint().z()));
1149       MeasurementPoint mp2 = topol->measurementPosition(LocalPoint(
1150           (*closest_simhit).exitPoint().x(), (*closest_simhit).exitPoint().y(), (*closest_simhit).exitPoint().z()));
1151       all_row1 = mp1.x();
1152       all_col1 = mp1.y();
1153       all_row2 = mp2.x();
1154       all_col2 = mp2.y();
1155 
1156       all_simhitx = 0.5 * (all_x1 + all_x2);
1157       all_simhity = 0.5 * (all_y1 + all_y2);
1158 
1159       all_rechitresx = all_rechitx - all_simhitx;
1160       all_rechitresy = all_rechity - all_simhity;
1161 
1162       all_rechitpullx = all_rechitresx / all_rechiterrx;
1163       all_rechitpully = all_rechitresy / all_rechiterry;
1164 
1165       SiPixelRecHit::ClusterRef const& clust = pixeliter->cluster();
1166 
1167       all_npix = clust->size();
1168       all_nxpix = clust->sizeX();
1169       all_nypix = clust->sizeY();
1170 
1171       all_clust_row = clust->x();
1172       all_clust_col = clust->y();
1173 
1174       LocalPoint lp2 = topol->localPosition(MeasurementPoint(all_clust_row, all_clust_col));
1175       all_clust_x = lp2.x();
1176       all_clust_y = lp2.y();
1177 
1178       all_clust_q = clust->charge();
1179 
1180       all_clust_maxpixcol = clust->maxPixelCol();
1181       all_clust_maxpixrow = clust->maxPixelRow();
1182       all_clust_minpixcol = clust->minPixelCol();
1183       all_clust_minpixrow = clust->minPixelRow();
1184 
1185       all_clust_geoid = 0;  // never set!
1186 
1187       all_simpx = (*closest_simhit).momentumAtEntry().x();
1188       all_simpy = (*closest_simhit).momentumAtEntry().y();
1189       all_simpz = (*closest_simhit).momentumAtEntry().z();
1190       all_eloss = (*closest_simhit).energyLoss();
1191       all_simphi = (*closest_simhit).phiAtEntry();
1192       all_simtheta = (*closest_simhit).thetaAtEntry();
1193       all_pidhit = (*closest_simhit).particleType();
1194       all_trkid = (*closest_simhit).trackId();
1195 
1196       //--- Fill alpha and beta -- more useful for exploring the residuals...
1197       all_beta = atan2(all_simpz, all_simpy);
1198       all_alpha = atan2(all_simpz, all_simpx);
1199 
1200       all_simproc = (int)closest_simhit->processType();
1201 
1202       const edm::SimTrackContainer& trks = *(simtracks.product());
1203       SimTrackContainer::const_iterator trksiter;
1204       for (trksiter = trks.begin(); trksiter != trks.end(); trksiter++)
1205         if ((int)trksiter->trackId() == all_trkid) {
1206           all_simtrketa = trksiter->momentum().eta();
1207           all_simtrkphi = trksiter->momentum().phi();
1208         }
1209 
1210       all_vtxz = theGeomDet->surface().position().z();
1211       all_vtxr = theGeomDet->surface().position().perp();
1212 
1213       //computeAnglesFromDetPosition(clust,
1214       //               theGeomDet,
1215       //               all_clust_alpha, all_clust_beta )
1216 
1217       const std::vector<SiPixelCluster::Pixel>& pixvector = clust->pixels();
1218       for (int i = 0; i < (int)pixvector.size(); ++i) {
1219         SiPixelCluster::Pixel holdpix = pixvector[i];
1220         all_pixrow[i] = holdpix.x;
1221         all_pixcol[i] = holdpix.y;
1222         all_pixadc[i] = holdpix.adc;
1223         LocalPoint lp = topol->localPosition(MeasurementPoint(holdpix.x, holdpix.y));
1224         all_pixx[i] = lp.x();
1225         all_pixy[i] = lp.y();
1226         GlobalPoint GP = theGeomDet->surface().toGlobal(Local3DPoint(lp.x(), lp.y(), lp.z()));
1227         all_pixgx[i] = GP.x();
1228         all_pixgy[i] = GP.y();
1229         all_pixgz[i] = GP.z();
1230       }
1231 
1232       ttree_all_hits_->Fill();
1233 
1234     }  // for ( ; pixeliter != pixelrechitRangeIteratorEnd; ++pixeliter)
1235 
1236   }  // for (TrackerGeometry::DetContainer::const_iterator it = pDD->dets().begin(); it != pDD->dets().end(); it++)
1237 
1238   // ------------------------------------------------ all hits ---------------------------------------------------------------
1239 
1240   //cout << "...3..." << endl;
1241 
1242   // ------------------------------------------------ track hits only --------------------------------------------------------
1243 
1244   if (include_trk_hits_) {
1245     // Get tracks
1246     edm::Handle<reco::TrackCollection> trackCollection;
1247     e.getByToken(tTrackCollection, trackCollection);
1248     const reco::TrackCollection* tracks = trackCollection.product();
1249     reco::TrackCollection::const_iterator tciter;
1250 
1251     if (!tracks->empty()) {
1252       // Loop on tracks
1253       for (tciter = tracks->begin(); tciter != tracks->end(); ++tciter) {
1254         // First loop on hits: find matched hits
1255         for (auto const hit : tciter->recHits()) {
1256           // Is it a matched hit?
1257           const SiPixelRecHit* matchedhit = dynamic_cast<const SiPixelRecHit*>(hit);
1258 
1259           if (matchedhit) {
1260             rechitx = -9999.9;
1261             rechity = -9999.9;
1262             rechitz = -9999.9;
1263             rechiterrx = -9999.9;
1264             rechiterry = -9999.9;
1265             rechitresx = -9999.9;
1266             rechitresy = -9999.9;
1267             rechitpullx = -9999.9;
1268             rechitpully = -9999.9;
1269 
1270             npix = -9999;
1271             nxpix = -9999;
1272             nypix = -9999;
1273             charge = -9999.9;
1274 
1275             edgex = -9999;
1276             edgey = -9999;
1277 
1278             bigx = -9999;
1279             bigy = -9999;
1280 
1281             alpha = -9999.9;
1282             beta = -9999.9;
1283 
1284             phi = -9999.9;
1285             eta = -9999.9;
1286 
1287             subdetId = -9999;
1288 
1289             layer = -9999;
1290             ladder = -9999;
1291             mod = -9999;
1292             side = -9999;
1293             disk = -9999;
1294             blade = -9999;
1295             panel = -9999;
1296             plaq = -9999;
1297 
1298             half = -9999;
1299             flipped = -9999;
1300 
1301             nsimhit = -9999;
1302             pidhit = -9999;
1303             simproc = -9999;
1304 
1305             simhitx = -9999.9;
1306             simhity = -9999.9;
1307 
1308             hit_probx = -9999.9;
1309             hit_proby = -9999.9;
1310             hit_cprob0 = -9999.9;
1311             hit_cprob1 = -9999.9;
1312             hit_cprob2 = -9999.9;
1313 
1314             pixel_split = -9999;
1315 
1316             pixel_clst_err_x = -9999.9;
1317             pixel_clst_err_y = -9999.9;
1318 
1319             position = hit->localPosition();
1320             error = hit->localPositionError();
1321 
1322             rechitx = position.x();
1323             rechity = position.y();
1324             rechitz = position.z();
1325             rechiterrx = sqrt(error.xx());
1326             rechiterry = sqrt(error.yy());
1327 
1328             npix = matchedhit->cluster()->size();
1329             nxpix = matchedhit->cluster()->sizeX();
1330             nypix = matchedhit->cluster()->sizeY();
1331             charge = matchedhit->cluster()->charge();
1332 
1333             if (matchedhit->cluster()->getSplitClusterErrorX() > 0.0 &&
1334                 matchedhit->cluster()->getSplitClusterErrorY() > 0.0)
1335               pixel_split = 1;
1336             else
1337               pixel_split = 0;
1338 
1339             pixel_clst_err_x = matchedhit->cluster()->getSplitClusterErrorX();
1340             pixel_clst_err_y = matchedhit->cluster()->getSplitClusterErrorY();
1341 
1342             //hit_probx  = (float)matchedhit->probabilityX();
1343             //hit_proby  = (float)matchedhit->probabilityY();
1344             hit_cprob0 = (float)matchedhit->clusterProbability(0);
1345             hit_cprob1 = (float)matchedhit->clusterProbability(1);
1346             hit_cprob2 = (float)matchedhit->clusterProbability(2);
1347 
1348             //Association of the rechit to the simhit
1349             matched.clear();
1350             matched = associate.associateHit(*matchedhit);
1351 
1352             nsimhit = (int)matched.size();
1353 
1354             if (!matched.empty()) {
1355               mindist = 999999.9;
1356               float distx, disty, dist;
1357               bool found_hit_from_generated_particle = false;
1358 
1359               vector<PSimHit>::const_iterator closestit = matched.begin();
1360               for (vector<PSimHit>::const_iterator m = matched.begin(); m < matched.end(); ++m) {
1361                 if (checkType_) {  // only consider associated simhits with the generated pid (muons)
1362                   int pid = (*m).particleType();
1363                   if (abs(pid) != genType_)
1364                     continue;
1365                 }
1366 
1367                 float simhitx = 0.5 * ((*m).entryPoint().x() + (*m).exitPoint().x());
1368                 float simhity = 0.5 * ((*m).entryPoint().y() + (*m).exitPoint().y());
1369 
1370                 distx = fabs(rechitx - simhitx);
1371                 disty = fabs(rechity - simhity);
1372                 dist = sqrt(distx * distx + disty * disty);
1373 
1374                 if (dist < mindist) {
1375                   mindist = dist;
1376                   closestit = m;
1377                   found_hit_from_generated_particle = true;
1378                 }
1379               }  // for (vector<PSimHit>::const_iterator m=matched.begin(); m<matched.end(); m++)
1380 
1381               // This recHit does not have any simHit with the same particleType as the particles generated
1382               // Ignore it as most probably come from delta rays.
1383               if (checkType_ && !found_hit_from_generated_particle)
1384                 continue;
1385 
1386               DetId detId = hit->geographicalId();
1387 
1388               const PixelGeomDetUnit* theGeomDet = dynamic_cast<const PixelGeomDetUnit*>((*tracker).idToDet(detId));
1389 
1390               const PixelTopology* theTopol = &(theGeomDet->specificTopology());
1391 
1392               pidhit = (*closestit).particleType();
1393               simproc = (int)(*closestit).processType();
1394 
1395               simhitx = 0.5 * ((*closestit).entryPoint().x() + (*closestit).exitPoint().x());
1396               simhity = 0.5 * ((*closestit).entryPoint().y() + (*closestit).exitPoint().y());
1397 
1398               rechitresx = rechitx - simhitx;
1399               rechitresy = rechity - simhity;
1400               rechitpullx = (rechitx - simhitx) / sqrt(error.xx());
1401               rechitpully = (rechity - simhity) / sqrt(error.yy());
1402 
1403               float simhitpx = (*closestit).momentumAtEntry().x();
1404               float simhitpy = (*closestit).momentumAtEntry().y();
1405               float simhitpz = (*closestit).momentumAtEntry().z();
1406 
1407               beta = atan2(simhitpz, simhitpy) * radtodeg;
1408               alpha = atan2(simhitpz, simhitpx) * radtodeg;
1409 
1410               //beta  = fabs(atan2(simhitpz, simhitpy)) * radtodeg;
1411               //alpha = fabs(atan2(simhitpz, simhitpx)) * radtodeg;
1412 
1413               // calculate alpha and beta exactly as in PixelCPEBase.cc
1414               float locx = simhitpx;
1415               float locy = simhitpy;
1416               float locz = simhitpz;
1417 
1418               bool isFlipped = false;
1419               float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 0.)).perp();
1420               float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 1.)).perp();
1421               if (tmp2 < tmp1)
1422                 isFlipped = true;
1423               else
1424                 isFlipped = false;
1425 
1426               trk_alpha = acos(locx / sqrt(locx * locx + locz * locz)) * radtodeg;
1427               if (isFlipped)  // &&& check for FPIX !!!
1428                 trk_alpha = 180.0 - trk_alpha;
1429 
1430               trk_beta = acos(locy / sqrt(locy * locy + locz * locz)) * radtodeg;
1431 
1432               phi = tciter->momentum().phi() / math_pi * 180.0;
1433               eta = tciter->momentum().eta();
1434 
1435               const int maxPixelCol = (*matchedhit).cluster()->maxPixelCol();
1436               const int maxPixelRow = (*matchedhit).cluster()->maxPixelRow();
1437               const int minPixelCol = (*matchedhit).cluster()->minPixelCol();
1438               const int minPixelRow = (*matchedhit).cluster()->minPixelRow();
1439 
1440               // check whether the cluster is at the module edge
1441               if (theTopol->isItEdgePixelInX(minPixelRow) || theTopol->isItEdgePixelInX(maxPixelRow))
1442                 edgex = 1;
1443               else
1444                 edgex = 0;
1445 
1446               if (theTopol->isItEdgePixelInY(minPixelCol) || theTopol->isItEdgePixelInY(maxPixelCol))
1447                 edgey = 1;
1448               else
1449                 edgey = 0;
1450 
1451               // check whether this rechit contains big pixels
1452               if (theTopol->containsBigPixelInX(minPixelRow, maxPixelRow))
1453                 bigx = 1;
1454               else
1455                 bigx = 0;
1456 
1457               if (theTopol->containsBigPixelInY(minPixelCol, maxPixelCol))
1458                 bigy = 1;
1459               else
1460                 bigy = 0;
1461 
1462               subdetId = (int)detId.subdetId();
1463 
1464               if ((int)detId.subdetId() == (int)PixelSubdetector::PixelBarrel) {
1465                 int tmp_nrows = theGeomDet->specificTopology().nrows();
1466                 if (tmp_nrows == 80)
1467                   half = 1;
1468                 else if (tmp_nrows == 160)
1469                   half = 0;
1470                 else
1471                   cout << "-------------------------------------------------- Wrong module size !!!" << endl;
1472 
1473                 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 0.)).perp();
1474                 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 1.)).perp();
1475 
1476                 if (tmp2 < tmp1)
1477                   flipped = 1;
1478                 else
1479                   flipped = 0;
1480 
1481                 layer = tTopo->pxbLayer(detId);    // Layer: 1,2,3.
1482                 ladder = tTopo->pxbLadder(detId);  // Ladder: 1-20, 32, 44.
1483                 mod = tTopo->pxbModule(detId);     // Mod: 1-8.
1484               } else if ((int)detId.subdetId() == (int)PixelSubdetector::PixelEndcap) {
1485                 side = tTopo->pxfSide(detId);
1486                 disk = tTopo->pxfDisk(detId);
1487                 blade = tTopo->pxfBlade(detId);
1488                 panel = tTopo->pxfPanel(detId);
1489                 plaq = tTopo->pxfModule(detId);  // also known as plaquette
1490 
1491                 float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 0.)).perp();
1492                 float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0., 0., 1.)).perp();
1493 
1494                 if (tmp2 < tmp1)
1495                   flipped = 1;
1496                 else
1497                   flipped = 0;
1498 
1499               }  // else if ( detId.subdetId()==PixelSubdetector::PixelEndcap )
1500               //else std::// cout << "We are not in the pixel detector. detId.subdetId() = " << (int)detId.subdetId() << endl;
1501 
1502               ttree_track_hits_->Fill();
1503 
1504             }  // if ( !matched.empty() )
1505             else
1506               cout << "---------------- RecHit with no associated SimHit !!! -------------------------- " << endl;
1507 
1508           }  // if ( matchedhit )
1509 
1510         }  // end of loop on hits
1511 
1512       }  //end of loop on track
1513 
1514     }  // tracks > 0.
1515 
1516   }  // if ( include_trk_hits_ )
1517 
1518   // ----------------------------------------------- track hits only -----------------------------------------------------------
1519 }
1520 
1521 void SiPixelErrorEstimation::computeAnglesFromDetPosition(const SiPixelCluster& cl,
1522                                                           const GeomDetUnit& det,
1523                                                           float& alpha,
1524                                                           float& beta) {
1525   //--- This is a new det unit, so cache it
1526   const PixelGeomDetUnit* theDet = dynamic_cast<const PixelGeomDetUnit*>(&det);
1527   if (!theDet) {
1528     cout << "---------------------------------------------- Not a pixel detector !!!!!!!!!!!!!!" << endl;
1529     assert(0);
1530   }
1531 
1532   const PixelTopology* theTopol = &(theDet->specificTopology());
1533 
1534   // get cluster center of gravity (of charge)
1535   float xcenter = cl.x();
1536   float ycenter = cl.y();
1537 
1538   // get the cluster position in local coordinates (cm)
1539   LocalPoint lp = theTopol->localPosition(MeasurementPoint(xcenter, ycenter));
1540   //float lp_mod = sqrt( lp.x()*lp.x() + lp.y()*lp.y() + lp.z()*lp.z() );
1541 
1542   // get the cluster position in global coordinates (cm)
1543   GlobalPoint gp = theDet->surface().toGlobal(lp);
1544   float gp_mod = sqrt(gp.x() * gp.x() + gp.y() * gp.y() + gp.z() * gp.z());
1545 
1546   // normalize
1547   float gpx = gp.x() / gp_mod;
1548   float gpy = gp.y() / gp_mod;
1549   float gpz = gp.z() / gp_mod;
1550 
1551   // make a global vector out of the global point; this vector will point from the
1552   // origin of the detector to the cluster
1553   GlobalVector gv(gpx, gpy, gpz);
1554 
1555   // make local unit vector along local X axis
1556   const Local3DVector lvx(1.0, 0.0, 0.0);
1557 
1558   // get the unit X vector in global coordinates/
1559   GlobalVector gvx = theDet->surface().toGlobal(lvx);
1560 
1561   // make local unit vector along local Y axis
1562   const Local3DVector lvy(0.0, 1.0, 0.0);
1563 
1564   // get the unit Y vector in global coordinates
1565   GlobalVector gvy = theDet->surface().toGlobal(lvy);
1566 
1567   // make local unit vector along local Z axis
1568   const Local3DVector lvz(0.0, 0.0, 1.0);
1569 
1570   // get the unit Z vector in global coordinates
1571   GlobalVector gvz = theDet->surface().toGlobal(lvz);
1572 
1573   // calculate the components of gv (the unit vector pointing to the cluster)
1574   // in the local coordinate system given by the basis {gvx, gvy, gvz}
1575   // note that both gv and the basis {gvx, gvy, gvz} are given in global coordinates
1576   float gv_dot_gvx = gv.x() * gvx.x() + gv.y() * gvx.y() + gv.z() * gvx.z();
1577   float gv_dot_gvy = gv.x() * gvy.x() + gv.y() * gvy.y() + gv.z() * gvy.z();
1578   float gv_dot_gvz = gv.x() * gvz.x() + gv.y() * gvz.y() + gv.z() * gvz.z();
1579 
1580   // calculate angles
1581   alpha = atan2(gv_dot_gvz, gv_dot_gvx);
1582   beta = atan2(gv_dot_gvz, gv_dot_gvy);
1583 
1584   // calculate cotalpha and cotbeta
1585   //   cotalpha_ = 1.0/tan(alpha_);
1586   //   cotbeta_  = 1.0/tan(beta_ );
1587   // or like this
1588   //cotalpha_ = gv_dot_gvx / gv_dot_gvz;
1589   //cotbeta_  = gv_dot_gvy / gv_dot_gvz;
1590 }
1591 
1592 //define this as a plug-in
1593 DEFINE_FWK_MODULE(SiPixelErrorEstimation);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!