CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_1/src/CalibTracker/SiPixelErrorEstimation/src/SiPixelErrorEstimation.cc

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