File indexing completed on 2023-03-17 11:19:38
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #include <memory>
0019 #include <iostream>
0020 #include <cmath>
0021
0022
0023 #include "TH1.h"
0024 #include "TH2.h"
0025 #include "TProfile.h"
0026
0027
0028 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0031
0032 #include "FWCore/Framework/interface/Event.h"
0033
0034 #include "FWCore/Framework/interface/MakerMacros.h"
0035
0036 #include <DataFormats/BeamSpot/interface/BeamSpot.h>
0037
0038 #include "DataFormats/VertexReco/interface/Vertex.h"
0039 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0040
0041
0042
0043
0044 #include "DataFormats/TrackReco/interface/Track.h"
0045 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0046 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0047
0048 #include <DataFormats/TrackReco/interface/HitPattern.h>
0049
0050
0051 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0052 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0053 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0054 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0055
0056 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0057
0058
0059 #include "FWCore/Framework/interface/ESHandle.h"
0060 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0061 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0062 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0063
0064 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0065 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0066
0067 #include <TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h>
0068 #include <MagneticField/Engine/interface/MagneticField.h>
0069
0070 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0071 #include <TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h>
0072 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0073
0074
0075
0076
0077
0078
0079 class myCounters {
0080 public:
0081 static int neve;
0082 static unsigned int prevrun;
0083 };
0084
0085 int myCounters::neve = 0;
0086 unsigned int myCounters::prevrun = 0;
0087
0088
0089
0090 class Triplet : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0091 public:
0092 explicit Triplet(const edm::ParameterSet &);
0093 ~Triplet();
0094
0095 private:
0096 virtual void beginJob();
0097 virtual void analyze(const edm::Event &, const edm::EventSetup &);
0098 virtual void endJob();
0099 void triplets(double x1,
0100 double y1,
0101 double z1,
0102 double x2,
0103 double y2,
0104 double z2,
0105 double x3,
0106 double y3,
0107 double z3,
0108 double ptsig,
0109 double &dc,
0110 double &dz);
0111
0112
0113 edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0114 edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0115 edm::EDGetTokenT<reco::TrackCollection> trackToken_;
0116
0117 edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> trackerTopoToken_;
0118 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0119 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackBuilderToken_;
0120 edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> transientTrackingRecHitBuilderToken_;
0121
0122 TH1D *h000, *h001, *h002, *h003, *h004, *h005, *h006, *h007, *h008, *h009;
0123 TH1D *h010, *h011, *h012, *h013, *h014, *h015, *h016, *h017, *h018, *h019;
0124 TH1D *h020, *h021, *h022, *h023, *h024, *h025, *h028, *h029;
0125 TH1D *h030, *h031, *h032, *h033, *h034, *h035, *h036, *h037, *h038, *h039;
0126 TH1D *h040, *h041, *h042, *h043, *h044, *h045, *h046, *h047, *h048, *h049;
0127 TH1D *h050, *h051, *h052, *h053, *h054, *h055, *h056, *h057, *h058, *h059;
0128 TH1D *h060, *h061, *h062, *h063, *h064, *h065, *h067, *h068, *h069;
0129 TH2D *h066;
0130 TH1D *h070, *h071, *h072, *h073, *h074, *h075, *h076, *h077, *h078, *h079;
0131 TH1D *h080, *h081, *h082, *h083, *h084, *h085, *h086, *h087, *h088, *h089;
0132 TH1D *h090, *h091, *h092, *h093, *h094;
0133 TH2D *h095, *h096, *h097, *h098, *h099;
0134
0135 TH1D *h100, *h101, *h102, *h103, *h104, *h105, *h108;
0136 TH2D *h106, *h107, *h109;
0137 TH1D *h110, *h112, *h113, *h114, *h115, *h116, *h118, *h119;
0138 TH2D *h111, *h117;
0139 TH1D *h120, *h121, *h122, *h123, *h124, *h125, *h126, *h127, *h128, *h129;
0140 TH1D *h130, *h131, *h132, *h133, *h134, *h135, *h136, *h137, *h138, *h139;
0141 TH1D *h140, *h141, *h142, *h143, *h144, *h145, *h146, *h147, *h148, *h149;
0142 TH1D *h150, *h151, *h152, *h153, *h154, *h155, *h156, *h157, *h158, *h159;
0143 TH1D *h160, *h161, *h162, *h163, *h164, *h165, *h166, *h167, *h168, *h169;
0144 TH1D *h170, *h171, *h172, *h173, *h174, *h175, *h176, *h177, *h178, *h179;
0145 TH1D *h180, *h181, *h182, *h183, *h184, *h185, *h186, *h187, *h188, *h189;
0146 TH1D *h190, *h191, *h192, *h193, *h194, *h195, *h196, *h197, *h198, *h199;
0147
0148 TH1D *h200, *h201, *h202, *h203, *h204, *h205, *h208;
0149 TH2D *h206, *h207, *h209;
0150
0151 TH1D *h300, *h301, *h302, *h303, *h304, *h305, *h308;
0152 TH2D *h306, *h307, *h309;
0153
0154 TH1D *h400, *h401, *h402, *h403, *h404, *h405, *h406, *h407, *h408;
0155 TProfile *h409;
0156 TH1D *h410, *h411;
0157 TProfile *h412, *h413, *h414, *h415, *h416, *h417, *h418, *h419;
0158 TH1D *h420, *h421;
0159 TProfile *h422, *h423, *h424, *h425, *h426, *h427, *h428, *h429;
0160 TH1D *h430, *h431, *h432, *h435, *h436, *h437, *h438, *h439;
0161 TProfile *h433, *h434;
0162 TH1D *h440, *h441;
0163 TProfile *h442, *h443, *h444, *h445, *h446, *h447, *h448, *h449;
0164 TH1D *h450, *h451;
0165 };
0166
0167
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178 Triplet::Triplet(const edm::ParameterSet &iConfig) {
0179 usesResource(TFileService::kSharedResource);
0180 bsToken_ = consumes<reco::BeamSpot>(edm::InputTag("offlineBeamSpot"));
0181 vtxToken_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0182 trackToken_ = consumes<reco::TrackCollection>(edm::InputTag("generalTracks"));
0183
0184 trackerTopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0185 trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0186 transientTrackBuilderToken_ =
0187 esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"));
0188 transientTrackingRecHitBuilderToken_ =
0189 esConsumes<TransientTrackingRecHitBuilder, TransientRecHitRecord>(edm::ESInputTag("", "WithTrackAngle"));
0190
0191 edm::LogPrint("Triplet") << "Triplet constructed\n";
0192 }
0193
0194
0195
0196 Triplet::~Triplet() = default;
0197
0198
0199
0200
0201
0202 void Triplet::beginJob() {
0203 edm::Service<TFileService> fs;
0204
0205 h000 = fs->make<TH1D>("h000", "beta star; beta star [cm]", 100, 0, 500);
0206 h001 = fs->make<TH1D>("h001", "emittance; emittance x [cm]", 100, 0, 1e-5);
0207 h002 = fs->make<TH1D>("h002", "beam width x; beam width x [um]", 100, 0, 200);
0208 h003 = fs->make<TH1D>("h003", "beam width y; beam width y [um]", 100, 0, 200);
0209 h004 = fs->make<TH1D>("h004", "beam spot sigma z; beam spot sigma z [cm]", 100, 0, 20);
0210 h005 = fs->make<TH1D>("h005", "beam spot x; beam spot x [cm]", 100, -1, 1);
0211 h006 = fs->make<TH1D>("h006", "beam spot y; beam spot y [cm]", 100, -1, 1);
0212 h007 = fs->make<TH1D>("h007", "beam spot z; beam spot z [cm]", 100, -5, 5);
0213 h008 = fs->make<TH1D>("h008", "beam slope dxdz; beam slope dxdz [mrad]", 100, -5, 5);
0214 h009 = fs->make<TH1D>("h009", "beam slope dydz; beam slope dydz [mrad]", 100, -5, 5);
0215
0216 h010 = fs->make<TH1D>("h010", "number of primary vertices; vertices; events", 31, -0.5, 30.5);
0217 h011 = fs->make<TH1D>("h011", "invalid z-vertex;z [cm]", 100, -50, 50);
0218 h012 = fs->make<TH1D>("h012", "fake z-vertex;z [cm]", 100, -50, 50);
0219 h013 = fs->make<TH1D>("h013", "non-fake z-vertex;z [cm]", 100, -50, 50);
0220 h014 = fs->make<TH1D>("h014", "vertex x;x [cm]", 100, -0.5, 0.5);
0221 h015 = fs->make<TH1D>("h015", "vertex y;y [cm]", 100, -0.5, 0.5);
0222 h016 = fs->make<TH1D>("h016", "tracks per vertex; tracks; vertices", 101, -0.5, 100.5);
0223 h017 = fs->make<TH1D>("h017", "tracks per vertex; tracks; vertices", 100, 5, 505);
0224 h018 = fs->make<TH1D>("h018", "z-vertex with refitted tracks;z [cm]", 100, -50, 50);
0225 h019 = fs->make<TH1D>("h019", "z-vertex without refitted tracks;z [cm]", 100, -50, 50);
0226
0227 h021 = fs->make<TH1D>("h021", "vertex sum pt; sum pt [GeV]", 100, 0, 100);
0228 h022 = fs->make<TH1D>("h022", "vertex max sum pt; max sum pt [GeV]", 100, 0, 100);
0229
0230 h023 = fs->make<TH1D>("h023", "best vertex x;x [cm]", 100, -0.25, 0.25);
0231 h024 = fs->make<TH1D>("h024", "best vertex y;y [cm]", 100, -0.25, 0.25);
0232 h025 = fs->make<TH1D>("h025", "best vertex z;z [cm]", 100, -25, 25);
0233
0234 h028 = fs->make<TH1D>("h028", "sum track pt; sum track Pt [GeV]", 100, 0, 500);
0235 h029 = fs->make<TH1D>("h029", "sum primary track charge; sum track charge", 41, -20.5, 20.5);
0236
0237 h040 = fs->make<TH1D>("h040", "number of tracks; tracks", 101, -5, 1005);
0238 h041 = fs->make<TH1D>("h041", "track charge; charge", 11, -5.5, 5.5);
0239 h042 = fs->make<TH1D>("h042", "pt; pt [GeV]", 100, 0, 5);
0240 h043 = fs->make<TH1D>("h043", "pt use logy, pt [GeV]", 100, 0, 100);
0241
0242 h044 = fs->make<TH1D>("h044", "number of rec hits; hits; tracks", 51, -0.5, 50.5);
0243 h045 = fs->make<TH1D>("h045", "valid tracker hits; tracker hits; tracks", 51, -0.5, 50.5);
0244 h046 = fs->make<TH1D>("h046", "valid pixel barrel hits; valid pixel barrel hits; tracks", 11, -0.5, 10.5);
0245 h047 = fs->make<TH1D>("h047", "tracker layers; tracker layers; tracks", 31, -0.5, 30.5);
0246 h048 = fs->make<TH1D>("h048", "pixel barrel layers; pixel barrel layers; tracks", 11, -0.5, 10.5);
0247
0248 h051 = fs->make<TH1D>("h051", "kap-kap; dkap; tracks", 100, -0.01, 0.01);
0249 h052 = fs->make<TH1D>("h052", "phi-phi; dphi; tracks", 100, -0.1, 0.1);
0250 h053 = fs->make<TH1D>("h053", "dca-dca; ddca; tracks", 100, -0.1, 0.1);
0251 h054 = fs->make<TH1D>("h054", "dip-dip; ddip; tracks", 100, -0.1, 0.1);
0252 h055 = fs->make<TH1D>("h055", "z0-z0; dz0; tracks", 100, -0.1, 0.1);
0253
0254 h056 = fs->make<TH1D>("h056", "tracks", 100, 0.01142, 0.01143);
0255 h049 = fs->make<TH1D>("h049", "tracks", 1000, -0.000001, 0.000001);
0256
0257 h057 = fs->make<TH1D>("h057", "tscp ref x; x [cm]; hits", 100, -1, 1);
0258 h058 = fs->make<TH1D>("h058", "tscp ref y; y [cm]; hits", 100, -1, 1);
0259 h059 = fs->make<TH1D>("h059", "tscp ref z; z [cm]; hits", 100, -10, 10);
0260
0261 h060 = fs->make<TH1D>("h060", "rec hit tracker subdet; subdet ID; tracks", 11, -0.5, 10.5);
0262 h061 = fs->make<TH1D>("h061", "rec hits local x; x [cm]; hits", 120, -6, 6);
0263 h062 = fs->make<TH1D>("h062", "rec hits local y; y [cm]; hits", 80, -4, 4);
0264
0265 h063 = fs->make<TH1D>("h063", "rec hits global R; R [cm]; hits", 120, 0, 120);
0266 h064 = fs->make<TH1D>("h064", "rec hits global phi; phi [deg]; hits", 180, -180, 180);
0267 h065 = fs->make<TH1D>("h065", "rec hits global z; z [cm]; hits", 300, -300, 300);
0268
0269 h066 = fs->make<TH2D>("h066", "rec hits barrel x-y; x [cm]; y [cm]", 240, -120, 120, 240, -120, 120);
0270
0271 h100 = fs->make<TH1D>("h100", "hits on tracks PXB layer; PXB layer; hits", 6, -0.5, 5.5);
0272
0273 h101 = fs->make<TH1D>("h101", "hits on tracks PXB1 ladder; ladder; hits", 22, -0.5, 21.5);
0274 h102 = fs->make<TH1D>("h102", "hits on tracks PXB1 module; module; hits", 10, -0.5, 9.5);
0275 h103 = fs->make<TH1D>("h103", "hits on tracks PXB1 R; R [cm]; hits", 150, 0, 15);
0276 h104 = fs->make<TH1D>("h104", "hits on tracks PXB1 phi; phi [deg]; hits", 180, -180, 180);
0277 h105 = fs->make<TH1D>("h105", "hits on tracks PXB1 z; z [cm]; hits", 600, -30, 30);
0278 h106 = fs->make<TH2D>("h106", "hits on tracks PXB1 phi-z; phi [deg]; z [cm]", 180, -180, 180, 600, -30, 30);
0279 h107 = fs->make<TH2D>("h107", "hits local x-y PXB1; x [cm]; y [cm]", 180, -0.9, 0.9, 440, -3.3, 3.3);
0280
0281 h111 = fs->make<TH2D>("h111", "hits on tracks PXB1 x-y; x [cm]; y [cm]", 240, -6, 6, 240, -6, 6);
0282 h112 = fs->make<TH1D>("h112", "residuals PXB1 dU; dU [um]; hits", 100, -250, 250);
0283 h113 = fs->make<TH1D>("h113", "residuals PXB1 dZ; dZ [um]; hits", 100, -250, 250);
0284 h114 = fs->make<TH1D>("h114", "residuals PXB1 dU; dU [um]; hits", 100, -250, 250);
0285 h115 = fs->make<TH1D>("h115", "residuals PXB1 dZ; dZ [um]; hits", 100, -250, 250);
0286
0287 h201 = fs->make<TH1D>("h201", "hits on tracks PXB2 ladder; ladder; hits", 34, -0.5, 33.5);
0288 h202 = fs->make<TH1D>("h202", "hits on tracks PXB2 module; module; hits", 10, -0.5, 9.5);
0289 h203 = fs->make<TH1D>("h203", "hits on tracks PXB2 R; R [cm]; hits", 150, 0, 15);
0290 h204 = fs->make<TH1D>("h204", "hits on tracks PXB2 phi; phi [deg]; hits", 180, -180, 180);
0291 h205 = fs->make<TH1D>("h205", "hits on tracks PXB2 z; z [cm]; hits", 600, -30, 30);
0292 h206 = fs->make<TH2D>("h206", "hits on tracks PXB2 phi-z; phi [deg]; z [cm]", 180, -180, 180, 600, -30, 30);
0293 h207 = fs->make<TH2D>("h207", "hits local x-y PXB2; x [cm]; y [cm]", 180, -0.9, 0.9, 440, -3.3, 3.3);
0294
0295 h301 = fs->make<TH1D>("h301", "hits on tracks PXB3 ladder; ladder; hits", 46, -0.5, 45.5);
0296 h302 = fs->make<TH1D>("h302", "hits on tracks PXB3 module; module; hits", 10, -0.5, 9.5);
0297 h303 = fs->make<TH1D>("h303", "hits on tracks PXB3 R; R [cm]; hits", 150, 0, 15);
0298 h304 = fs->make<TH1D>("h304", "hits on tracks PXB3 phi; phi [deg]; hits", 180, -180, 180);
0299 h305 = fs->make<TH1D>("h305", "hits on tracks PXB3 z; z [cm]; hits", 600, -30, 30);
0300 h306 = fs->make<TH2D>("h306", "hits on tracks PXB3 phi-z; phi [deg]; z [cm]", 180, -180, 180, 600, -30, 30);
0301 h307 = fs->make<TH2D>("h307", "hits local x-y PXB3; x [cm]; y [cm]", 180, -0.9, 0.9, 440, -3.3, 3.3);
0302
0303
0304
0305 h401 = fs->make<TH1D>("h401", "triplets z2; z [cm]; hits", 600, -30, 30);
0306 h402 = fs->make<TH1D>("h402", "uphi-phi; dphi; tracks", 100, -0.1, 0.1);
0307 h403 = fs->make<TH1D>("h403", "udca-dca; ddca; tracks", 100, -0.1, 0.1);
0308 h404 = fs->make<TH1D>("h404", "udip-dip; ddip; tracks", 100, -0.1, 0.1);
0309 h405 = fs->make<TH1D>("h405", "uz0-z0; dz0; tracks", 100, -0.1, 0.1);
0310
0311 h406 = fs->make<TH1D>("h406", "valid tracker hits; tracker hits; tracks", 51, -0.5, 50.5);
0312 h407 = fs->make<TH1D>("h407", "valid pixel barrel hits; valid pixel barrel hits; tracks", 11, -0.5, 10.5);
0313 h408 = fs->make<TH1D>("h408", "tracker layers; tracker layers; tracks", 31, -0.5, 30.5);
0314 h409 = fs->make<TProfile>("h409", "angle of incidence; phi2 [deg]; phi inc 2 [deg]", 180, -180, 180, -90, 90);
0315
0316 h410 = fs->make<TH1D>("h410", "residuals PXB2 dca2; dca2 [um]; hits", 100, -150, 150);
0317 h411 = fs->make<TH1D>("h411", "residuals PXB2 dz2 ; dz2 [um]; hits", 100, -300, 300);
0318
0319
0320
0321 h412 = fs->make<TProfile>("h412", "PXB2 dxy vs phi; phi2 [deg]; <dxy2> [um]", 180, -180, 180, -99, 99);
0322 h413 = fs->make<TProfile>("h413", "PXB2 dz vs phi; phi2 [deg]; <dz2> [um]", 180, -180, 180, -99, 199);
0323
0324 h414 = fs->make<TProfile>("h414", "PXB2 dxy vs z; z2 [cm]; <dxy2> [um]", 80, -20, 20, -99, 99);
0325 h415 = fs->make<TProfile>("h415", "PXB2 dz vs z; z2 [cm]; <dz2> [um]", 80, -20, 20, -199, 199);
0326
0327 h416 = fs->make<TProfile>("h416", "PXB2 dxy vs pt; log(pt [GeV]); <dxy2> [um]", 20, 0, 2, -99, 99);
0328 h417 = fs->make<TProfile>("h417", "PXB2 dz vs pt; log(pt [GeV]); <dz2> [um]", 20, 0, 2, -199, 199);
0329
0330 h418 = fs->make<TProfile>("h418", "PXB2 dxy vs pt +; log(pt [GeV]); <dxy2> [um]", 20, 0, 2, -99, 99);
0331 h419 = fs->make<TProfile>("h419", "PXB2 dxy vs pt -; log(pt [GeV]); <dxy2> [um]", 20, 0, 2, -99, 99);
0332
0333 h420 = fs->make<TH1D>("h420", "residuals PXB2 dca2, pt > 12; dca2 [um]; hits", 100, -150, 150);
0334 h421 = fs->make<TH1D>("h421", "residuals PXB2 dz2, pt > 12; dz2 [um]; hits", 100, -300, 300);
0335
0336
0337
0338 h422 = fs->make<TProfile>("h422", "PXB2 sxy vs phi; phi2 [deg]; sxy [um]", 360, -180, 180, 0, 99);
0339 h423 = fs->make<TProfile>("h423", "PXB2 sz vs phi; phi2 [deg]; sz [um]", 360, -180, 180, 0, 199);
0340
0341 h424 = fs->make<TProfile>("h424", "PXB2 sxy vs z; z2 [cm]; sxy [um]", 80, -20, 20, 0, 99);
0342 h425 = fs->make<TProfile>("h425", "PXB2 sz vs z; z2 [cm]; sz [um]", 80, -20, 20, 0, 199);
0343
0344 h426 = fs->make<TProfile>("h426", "PXB2 sxy vs pt; log(pt [GeV]); sxy [um]", 20, 0, 2, 0, 99);
0345 h427 = fs->make<TProfile>("h427", "PXB2 sz vs pt; log(pt [GeV]); sz [um]", 20, 0, 2, 0, 199);
0346
0347 h428 = fs->make<TProfile>("h428", "PXB2 sxy vs dip; dip [deg]; sxy [um]", 70, -70, 70, 0, 99);
0348 h429 = fs->make<TProfile>("h429", "PXB2 sz vs dip; dip [deg]; sz [um]", 70, -70, 70, 0, 199);
0349
0350 h430 = fs->make<TH1D>("h430", "residuals PXB2 dca2; dca2 [um]; hits", 100, -150, 150);
0351 h431 = fs->make<TH1D>("h431", "residuals PXB2 dz2; dz2 [um]; hits", 100, -300, 300);
0352
0353 h432 = fs->make<TH1D>("h432", "residuals PXB2 dz2, 18 < |dip| < 50; dz2 [um]; hits", 100, -300, 300);
0354 h433 = fs->make<TProfile>("h433", "PXB2 sz vs pt, best dip; log(pt [GeV]); sz [um]", 20, 0, 2, 0, 199);
0355
0356 h434 = fs->make<TProfile>("h434", "PXB2 sxy vs inc; phi inc 2 [deg]; sxy [um]", 40, -10, 10, 0, 99);
0357
0358 h435 = fs->make<TH1D>("h435", "ukap-kap; dkap; tracks", 100, -0.01, 0.01);
0359 h436 = fs->make<TH1D>("h436", "uphi-phi; dphi; tracks", 100, -0.1, 0.1);
0360 h437 = fs->make<TH1D>("h437", "udca-dca; ddca; tracks", 100, -0.1, 0.1);
0361
0362 h440 = fs->make<TH1D>("h440", "pixel track dcap, pt > 2; dcap [um]; hits", 100, -1000, 1000);
0363 h441 = fs->make<TH1D>("h441", "pixel track dcap, pt > 4; dcap [um]; hits", 100, -1000, 1000);
0364
0365 h442 = fs->make<TProfile>("h442", "pixel track dcap vs phi; phi2 [deg]; <dcap> [um]", 180, -180, 180, -500, 499);
0366 h443 = fs->make<TProfile>("h443", "pixel tracks dcap vs pt; log(pt [GeV]); <dcap> [um]", 20, 0, 2, -500, 499);
0367
0368 h444 = fs->make<TProfile>("h444", "pixel track sdcap vs phi; phi2 [deg]; sdcap [um]", 180, -180, 180, 0, 499);
0369 h445 = fs->make<TProfile>("h445", "pixel tracks sdcap vs pt; log(pt [GeV]); sdcap [um]", 20, 0, 2, 0, 499);
0370
0371 h450 = fs->make<TH1D>("h450", "residuals PXB2 dca2; dca2 [um]; hits", 100, -150, 150);
0372 h451 = fs->make<TH1D>("h451", "residuals PXB2 dz2; dz2 [um]; hits", 100, -300, 300);
0373 }
0374
0375
0376
0377
0378 void Triplet::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0379
0380 edm::ESHandle<TrackerTopology> tTopo = iSetup.getHandle(trackerTopoToken_);
0381
0382 using namespace std;
0383 using namespace edm;
0384 using namespace reco;
0385 using namespace math;
0386
0387 const double pi = 4 * atan(1);
0388 const double wt = 180 / pi;
0389 const double twopi = 2 * pi;
0390 const double pihalf = 2 * atan(1);
0391
0392
0393 myCounters::neve++;
0394
0395 if (myCounters::prevrun != iEvent.run()) {
0396 time_t unixZeit = iEvent.time().unixTime();
0397
0398 edm::LogPrint("Triplet") << "new run " << iEvent.run();
0399 edm::LogPrint("Triplet") << ", LumiBlock " << iEvent.luminosityBlock();
0400 edm::LogPrint("Triplet") << " taken " << ctime(&unixZeit);
0401
0402 myCounters::prevrun = iEvent.run();
0403
0404 }
0405
0406 int idbg = 0;
0407 if (myCounters::neve < 1)
0408 idbg = 1;
0409
0410
0411
0412
0413 edm::Handle<reco::BeamSpot> rbs;
0414 iEvent.getByToken(bsToken_, rbs);
0415
0416 XYZPoint bsP = XYZPoint(0, 0, 0);
0417
0418
0419 if (!rbs.failedToGet() && rbs.isValid()) {
0420
0421 h000->Fill(rbs->betaStar());
0422 h001->Fill(rbs->emittanceX());
0423 h002->Fill(rbs->BeamWidthX() * 1e4);
0424 h003->Fill(rbs->BeamWidthY() * 1e4);
0425 h004->Fill(rbs->sigmaZ());
0426 h005->Fill(rbs->x0());
0427 h006->Fill(rbs->y0());
0428 h007->Fill(rbs->z0());
0429 h008->Fill(rbs->dxdz() * 1e3);
0430 h009->Fill(rbs->dydz() * 1e3);
0431 bsP = XYZPoint(rbs->x0(), rbs->y0(), rbs->z0());
0432
0433 if (idbg) {
0434 edm::LogPrint("Triplet") << "beam spot x " << rbs->x0();
0435 edm::LogPrint("Triplet") << ", y " << rbs->y0();
0436 edm::LogPrint("Triplet") << ", z " << rbs->z0();
0437 edm::LogPrint("Triplet");
0438 }
0439
0440 }
0441
0442
0443
0444
0445 Handle<VertexCollection> vertices;
0446 iEvent.getByToken(vtxToken_, vertices);
0447
0448 if (vertices.failedToGet())
0449 return;
0450 if (!vertices.isValid())
0451 return;
0452
0453 h010->Fill(vertices->size());
0454
0455
0456
0457
0458
0459
0460
0461 XYZPoint vtxN = XYZPoint(0, 0, 0);
0462 XYZPoint vtxP = XYZPoint(0, 0, 0);
0463
0464 double bestNdof = 0;
0465 double maxSumPt = 0;
0466 Vertex bestPvx;
0467
0468 for (VertexCollection::const_iterator iVertex = vertices->begin(); iVertex != vertices->end(); ++iVertex) {
0469 if (!iVertex->isValid())
0470 h011->Fill(iVertex->z());
0471 else {
0472 if (iVertex->isFake())
0473 h012->Fill(iVertex->z());
0474
0475 else {
0476 h013->Fill(iVertex->z());
0477 h014->Fill(iVertex->x());
0478 h015->Fill(iVertex->y());
0479 h016->Fill(iVertex->ndof());
0480 h017->Fill(iVertex->ndof());
0481
0482 if (idbg) {
0483 edm::LogPrint("Triplet") << "vertex";
0484 edm::LogPrint("Triplet") << ": x " << iVertex->x();
0485 edm::LogPrint("Triplet") << ", y " << iVertex->y();
0486 edm::LogPrint("Triplet") << ", z " << iVertex->z();
0487 edm::LogPrint("Triplet") << ", ndof " << iVertex->ndof();
0488 edm::LogPrint("Triplet") << ", sumpt " << iVertex->p4().pt();
0489 edm::LogPrint("Triplet");
0490 }
0491
0492 if (iVertex->hasRefittedTracks())
0493 h018->Fill(iVertex->z());
0494 else
0495 h019->Fill(iVertex->z());
0496
0497 if (iVertex->ndof() > bestNdof) {
0498 bestNdof = iVertex->ndof();
0499 vtxN = XYZPoint(iVertex->x(), iVertex->y(), iVertex->z());
0500 }
0501
0502 h021->Fill(iVertex->p4().pt());
0503
0504 if (iVertex->p4().pt() > maxSumPt) {
0505 maxSumPt = iVertex->p4().pt();
0506 vtxP = XYZPoint(iVertex->x(), iVertex->y(), iVertex->z());
0507 bestPvx = *iVertex;
0508 }
0509 }
0510 }
0511 }
0512
0513 h022->Fill(maxSumPt);
0514
0515 #ifndef SIM
0516 if (maxSumPt < 1)
0517 return;
0518 #endif
0519 if (maxSumPt < 1)
0520 vtxP = vtxN;
0521
0522
0523
0524
0525
0526
0527
0528
0529 h023->Fill(vtxP.x());
0530 h024->Fill(vtxP.y());
0531 h025->Fill(vtxP.z());
0532
0533
0534
0535
0536
0537
0538
0539
0540
0541
0542
0543
0544
0545
0546
0547 Handle<TrackCollection> tracks;
0548
0549 iEvent.getByToken(trackToken_, tracks);
0550
0551 if (tracks.failedToGet())
0552 return;
0553 if (!tracks.isValid())
0554 return;
0555
0556 h040->Fill(tracks->size());
0557
0558
0559
0560 edm::ESHandle<TrackerGeometry> pDD = iSetup.getHandle(trackerGeomToken_);
0561
0562 if (!pDD.isValid()) {
0563 edm::LogPrint("Triplet") << "Unable to find TrackerDigiGeometry. Return\n";
0564 return;
0565 }
0566
0567
0568
0569 for (TrackerGeometry::DetContainer::const_iterator idet = pDD->dets().begin(); idet != pDD->dets().end(); ++idet) {
0570 DetId mydetId = (*idet)->geographicalId();
0571 uint32_t mysubDet = mydetId.subdetId();
0572
0573 if (idbg) {
0574 edm::LogPrint("Triplet") << "Det " << mydetId.det();
0575 edm::LogPrint("Triplet") << ", subDet " << mydetId.subdetId();
0576
0577 if (mysubDet == PixelSubdetector::PixelBarrel) {
0578 edm::LogPrint("Triplet") << ": PXB layer " << tTopo->pxbLayer(mydetId);
0579 edm::LogPrint("Triplet") << ", ladder " << tTopo->pxbLadder(mydetId);
0580 edm::LogPrint("Triplet") << ", module " << tTopo->pxbModule(mydetId);
0581 edm::LogPrint("Triplet") << ", at R " << (*idet)->position().perp();
0582 edm::LogPrint("Triplet") << ", F " << (*idet)->position().barePhi() * wt;
0583 edm::LogPrint("Triplet") << ", z " << (*idet)->position().z();
0584 edm::LogPrint("Triplet");
0585 edm::LogPrint("Triplet") << "rot x";
0586 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().xx();
0587 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().xy();
0588 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().xz();
0589 edm::LogPrint("Triplet");
0590 edm::LogPrint("Triplet") << "rot y";
0591 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().yx();
0592 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().yy();
0593 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().yz();
0594 edm::LogPrint("Triplet");
0595 edm::LogPrint("Triplet") << "rot z";
0596 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().zx();
0597 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().zy();
0598 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().zz();
0599 edm::LogPrint("Triplet");
0600
0601
0602
0603
0604 edm::LogPrint("Triplet") << "normal";
0605 edm::LogPrint("Triplet") << ": x " << (*idet)->surface().normalVector().x();
0606 edm::LogPrint("Triplet") << ", y " << (*idet)->surface().normalVector().y();
0607 edm::LogPrint("Triplet") << ", z " << (*idet)->surface().normalVector().z();
0608 edm::LogPrint("Triplet") << ", f " << (*idet)->surface().normalVector().barePhi() * wt;
0609
0610 }
0611
0612 if (mysubDet == PixelSubdetector::PixelEndcap) {
0613 edm::LogPrint("Triplet") << ": PXD side " << tTopo->pxfSide(mydetId);
0614 edm::LogPrint("Triplet") << ", disk " << tTopo->pxfDisk(mydetId);
0615 edm::LogPrint("Triplet") << ", blade " << tTopo->pxfBlade(mydetId);
0616 edm::LogPrint("Triplet") << ", panel " << tTopo->pxfPanel(mydetId);
0617 edm::LogPrint("Triplet") << ", module " << tTopo->pxfModule(mydetId);
0618 edm::LogPrint("Triplet") << ", at R " << (*idet)->position().perp();
0619 edm::LogPrint("Triplet") << ", F " << (*idet)->position().barePhi() * wt;
0620 edm::LogPrint("Triplet") << ", z " << (*idet)->position().z();
0621 edm::LogPrint("Triplet");
0622 edm::LogPrint("Triplet") << "rot x";
0623 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().xx();
0624 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().xy();
0625 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().xz();
0626 edm::LogPrint("Triplet");
0627 edm::LogPrint("Triplet") << "rot y";
0628 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().yx();
0629 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().yy();
0630 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().yz();
0631 edm::LogPrint("Triplet");
0632 edm::LogPrint("Triplet") << "rot z";
0633 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().zx();
0634 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().zy();
0635 edm::LogPrint("Triplet") << "\t" << (*idet)->rotation().zz();
0636 edm::LogPrint("Triplet");
0637 edm::LogPrint("Triplet") << "normal";
0638 edm::LogPrint("Triplet") << ": x " << (*idet)->surface().normalVector().x();
0639 edm::LogPrint("Triplet") << ", y " << (*idet)->surface().normalVector().y();
0640 edm::LogPrint("Triplet") << ", z " << (*idet)->surface().normalVector().z();
0641 edm::LogPrint("Triplet") << ", f " << (*idet)->surface().normalVector().barePhi() * wt;
0642
0643 }
0644
0645 edm::LogPrint("Triplet");
0646
0647 }
0648
0649 }
0650
0651
0652
0653
0654 edm::ESHandle<TransientTrackBuilder> theB = iSetup.getHandle(transientTrackBuilderToken_);
0655
0656
0657
0658 edm::ESHandle<TransientTrackingRecHitBuilder> hitBuilder = iSetup.getHandle(transientTrackingRecHitBuilderToken_);
0659
0660
0661
0662 double sumpt = 0;
0663 double sumq = 0;
0664 int kk = -1;
0665 Surface::GlobalPoint origin = Surface::GlobalPoint(0, 0, 0);
0666
0667
0668
0669
0670 for (TrackCollection::const_iterator iTrack = tracks->begin(); iTrack != tracks->end(); ++iTrack) {
0671 kk++;
0672
0673
0674
0675
0676
0677 double pt = iTrack->pt();
0678
0679 if (pt < 0.75)
0680 continue;
0681
0682
0683 if (abs(iTrack->dxy(vtxP)) > 5 * iTrack->dxyError())
0684 continue;
0685
0686 double logpt = log(pt) / log(10);
0687 double charge = iTrack->charge();
0688 h041->Fill(charge);
0689 h042->Fill(pt);
0690 h043->Fill(pt);
0691
0692 if (idbg) {
0693 edm::LogPrint("Triplet") << "Track " << kk;
0694 edm::LogPrint("Triplet") << ": pt " << iTrack->pt();
0695 edm::LogPrint("Triplet") << ", eta " << iTrack->eta();
0696 edm::LogPrint("Triplet") << ", phi " << iTrack->phi() * wt;
0697 edm::LogPrint("Triplet") << ", dxyv " << iTrack->dxy(vtxP);
0698 edm::LogPrint("Triplet") << ", dzv " << iTrack->dz(vtxP);
0699 edm::LogPrint("Triplet");
0700 }
0701
0702 const reco::HitPattern &hp = iTrack->hitPattern();
0703
0704 h045->Fill(hp.numberOfValidTrackerHits());
0705 h046->Fill(hp.numberOfValidPixelBarrelHits());
0706 h047->Fill(hp.trackerLayersWithMeasurement());
0707 h048->Fill(hp.pixelBarrelLayersWithMeasurement());
0708
0709 double phi = iTrack->phi();
0710 double dca = iTrack->d0();
0711
0712 double dip = iTrack->lambda();
0713 double z0 = iTrack->dz();
0714 double tet = pihalf - dip;
0715
0716
0717
0718
0719 TransientTrack tTrack = theB->build(*iTrack);
0720
0721 double kap = tTrack.initialFreeState().transverseCurvature();
0722
0723 TrajectoryStateClosestToPoint tscp = tTrack.trajectoryStateClosestToPoint(origin);
0724
0725
0726
0727 if (tscp.isValid()) {
0728 h057->Fill(tscp.referencePoint().x());
0729 h058->Fill(tscp.referencePoint().y());
0730 h059->Fill(tscp.referencePoint().z());
0731 kap = tscp.perigeeParameters().transverseCurvature();
0732 phi = tscp.perigeeParameters().phi();
0733 dca = tscp.perigeeParameters().transverseImpactParameter();
0734 tet = tscp.perigeeParameters().theta();
0735 z0 = tscp.perigeeParameters().longitudinalImpactParameter();
0736 dip = pihalf - tet;
0737
0738 h051->Fill(kap - tTrack.initialFreeState().transverseCurvature());
0739 h052->Fill(phi - iTrack->phi());
0740 h053->Fill(dca - iTrack->d0());
0741 h054->Fill(dip - iTrack->lambda());
0742 h055->Fill(z0 - iTrack->dz());
0743 }
0744
0745 double tmp = abs(kap * pt);
0746 h056->Fill(tmp);
0747 double rho1 = pt / 0.0114257;
0748 if (charge > 0)
0749 rho1 = -rho1;
0750 double kap1 = 1. / rho1;
0751 double tmp1 = (kap1 - kap);
0752 h049->Fill(tmp1);
0753
0754
0755
0756 double cf = cos(phi);
0757 double sf = sin(phi);
0758
0759
0760
0761
0762
0763
0764
0765 double erd = 1.0 - kap * dca;
0766 double drd = dca * (0.5 * kap * dca - 1.0);
0767 double hkk = 0.5 * kap * kap;
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780 double R1 = 4.4;
0781
0782 double s1 = 0;
0783 double fpos1 = phi - pihalf;
0784
0785 if (R1 >= abs(dca)) {
0786
0787
0788
0789 double sindp = (0.5 * kap * (R1 * R1 + dca * dca) - dca) / (R1 * erd);
0790 fpos1 = phi + asin(sindp);
0791
0792
0793
0794 double sina = R1 * kap * sqrt(1.0 - sindp * sindp);
0795
0796
0797
0798 if (sina >= 1.0)
0799 s1 = pi / kap;
0800 else {
0801 if (sina <= -1.0)
0802 s1 = -pi / kap;
0803 else
0804 s1 = asin(sina) / kap;
0805 }
0806
0807
0808
0809 if (hkk * (R1 * R1 - dca * dca) > erd)
0810 s1 = pi / abs(kap) - s1;
0811
0812 }
0813
0814 if (fpos1 > pi)
0815 fpos1 -= twopi;
0816 else if (fpos1 < -pi)
0817 fpos1 += twopi;
0818
0819 double zR1 = z0 + s1 * tan(dip);
0820
0821
0822
0823
0824 double xcrss[99];
0825 double ycrss[99];
0826 double zcrss[99];
0827 int ncrss = 0;
0828
0829 for (TrackerGeometry::DetContainer::const_iterator idet = pDD->dets().begin(); idet != pDD->dets().end(); ++idet) {
0830 DetId mydetId = (*idet)->geographicalId();
0831 uint32_t mysubDet = mydetId.subdetId();
0832
0833 if (mysubDet != PixelSubdetector::PixelBarrel)
0834 continue;
0835
0836
0837
0838
0839
0840
0841
0842
0843
0844
0845 if (tTopo->pxbLayer(mydetId) == 1) {
0846 double dz = zR1 - (*idet)->position().z();
0847
0848 if (abs(dz) > 4.0)
0849 continue;
0850
0851 double df = fpos1 - (*idet)->position().barePhi();
0852
0853 if (df > pi)
0854 df -= twopi;
0855 else if (df < -pi)
0856 df += twopi;
0857
0858 if (abs(df) * wt > 36.0)
0859 continue;
0860
0861
0862
0863
0864
0865
0866
0867
0868
0869
0870
0871
0872
0873 double phiN = (*idet)->surface().normalVector().barePhi();
0874
0875 double phidet = phiN - pihalf;
0876 double ux = cos(phidet);
0877 double uy = sin(phidet);
0878 double x = (*idet)->position().x();
0879 double y = (*idet)->position().y();
0880
0881
0882
0883
0884 double a = kap * (ux * ux + uy * uy) * 0.5;
0885 double b = erd * (ux * sf - uy * cf) + kap * (ux * x + uy * y);
0886 double c = drd + erd * (x * sf - y * cf) + kap * (x * x + y * y) * 0.5;
0887 double dis = b * b - 4.0 * a * c;
0888 double f = 0;
0889
0890 if (dis > 0) {
0891 dis = sqrt(dis);
0892 double f1 = 0;
0893 double f2 = 0;
0894
0895 if (b < 0) {
0896 f1 = 0.5 * (dis - b) / a;
0897 f2 = 2.0 * c / (dis - b);
0898 } else {
0899 f1 = -0.5 * (dis + b) / a;
0900 f2 = -2.0 * c / (dis + b);
0901 }
0902
0903 f = f1;
0904 if (abs(f2) < abs(f1))
0905 f = f2;
0906
0907 }
0908
0909 xcrss[ncrss] = x + f * ux;
0910 ycrss[ncrss] = y + f * uy;
0911 double r = sqrt(xcrss[ncrss] * xcrss[ncrss] + ycrss[ncrss] * ycrss[ncrss]);
0912
0913 double s = 0;
0914
0915 if (r >= abs(dca)) {
0916 double sindp = (0.5 * kap * (r * r + dca * dca) - dca) / (r * erd);
0917 double sina = r * kap * sqrt(1.0 - sindp * sindp);
0918 if (sina >= 1.0)
0919 s = pi / kap;
0920 else {
0921 if (sina <= -1.0)
0922 s = -pi / kap;
0923 else
0924 s = asin(sina) / kap;
0925 }
0926 if (hkk * (r * r - dca * dca) > erd)
0927 s = pi / abs(kap) - s;
0928 }
0929
0930 zcrss[ncrss] = z0 + s * tan(dip);
0931
0932 ncrss++;
0933
0934 }
0935
0936 }
0937
0938
0939
0940
0941 if (iTrack->extra().isNonnull() && iTrack->extra().isAvailable()) {
0942 h044->Fill(tTrack.recHitsSize());
0943
0944 double x1 = 0;
0945 double y1 = 0;
0946 double z1 = 0;
0947 double x2 = 0;
0948 double y2 = 0;
0949 double z2 = 0;
0950 double x3 = 0;
0951 double y3 = 0;
0952 double z3 = 0;
0953 int n1 = 0;
0954 int n2 = 0;
0955 int n3 = 0;
0956
0957
0958 for (trackingRecHit_iterator irecHit = iTrack->recHitsBegin(); irecHit != iTrack->recHitsEnd(); ++irecHit) {
0959 if ((*irecHit)->isValid()) {
0960 DetId detId = (*irecHit)->geographicalId();
0961
0962
0963
0964 if (detId.det() != 1) {
0965 edm::LogPrint("Triplet") << "rec hit ID = " << detId.det() << " not in tracker!?!?\n";
0966 continue;
0967 }
0968
0969 uint32_t subDet = detId.subdetId();
0970
0971
0972
0973
0974 h060->Fill(subDet);
0975
0976
0977
0978 TransientTrackingRecHit::RecHitPointer trecHit = hitBuilder->build(&*(*irecHit));
0979
0980 double xloc = trecHit->localPosition().x();
0981 double yloc = trecHit->localPosition().y();
0982
0983 h061->Fill(xloc);
0984 h062->Fill(yloc);
0985
0986 double gX = trecHit->globalPosition().x();
0987 double gY = trecHit->globalPosition().y();
0988 double gZ = trecHit->globalPosition().z();
0989 double gF = atan2(gY, gX);
0990 double gR = sqrt(gX * gX + gY * gY);
0991
0992 h063->Fill(gR);
0993 h064->Fill(gF * wt);
0994 h065->Fill(gZ);
0995
0996
0997
0998
0999 if (subDet == PixelSubdetector::PixelBarrel || subDet == StripSubdetector::TIB ||
1000 subDet == StripSubdetector::TOB) {
1001
1002 h066->Fill(gX, gY);
1003
1004 }
1005
1006 if (subDet == PixelSubdetector::PixelBarrel) {
1007 int ilay = tTopo->pxbLayer(detId);
1008 int ilad = tTopo->pxbLadder(detId);
1009 int imod = tTopo->pxbModule(detId);
1010 bool halfmod = 0;
1011
1012 h100->Fill(ilay);
1013
1014 if (ilay == 1) {
1015 n1++;
1016 x1 = gX;
1017 y1 = gY;
1018 z1 = gZ;
1019
1020 h101->Fill(ilad);
1021 h102->Fill(imod);
1022
1023 h103->Fill(gR);
1024 h104->Fill(gF * wt);
1025 h105->Fill(gZ);
1026
1027 h106->Fill(gF * wt, gZ);
1028
1029 if (ilad == 5)
1030 halfmod = 1;
1031 else if (ilad == 6)
1032 halfmod = 1;
1033 else if (ilad == 15)
1034 halfmod = 1;
1035 else if (ilad == 16)
1036 halfmod = 1;
1037
1038 if (!halfmod) {
1039 h107->Fill(xloc, yloc);
1040 }
1041
1042
1043
1044 for (int icrss = 0; icrss < ncrss; ++icrss) {
1045 double fcrss = atan2(ycrss[icrss], xcrss[icrss]);
1046 double df = gF - fcrss;
1047 if (df > pi)
1048 df -= twopi;
1049 else if (df < -pi)
1050 df += twopi;
1051 double du = gR * df;
1052 double dz = gZ - zcrss[icrss];
1053
1054 if (abs(du) < 0.01 && abs(dz) < 0.05)
1055 h111->Fill(gX, gY);
1056 h112->Fill(du * 1E4);
1057 h113->Fill(dz * 1E4);
1058
1059 if (abs(dz) < 0.02)
1060 h114->Fill(du * 1E4);
1061 if (abs(du) < 0.01)
1062 h115->Fill(dz * 1E4);
1063
1064 }
1065
1066 }
1067
1068 if (ilay == 2) {
1069 n2++;
1070 x2 = gX;
1071 y2 = gY;
1072 z2 = gZ;
1073
1074
1075 h201->Fill(ilad);
1076 h202->Fill(imod);
1077
1078 h203->Fill(gR);
1079 h204->Fill(gF * wt);
1080 h205->Fill(gZ);
1081
1082 h206->Fill(gF * wt, gZ);
1083
1084 if (ilad == 8)
1085 halfmod = 1;
1086 else if (ilad == 9)
1087 halfmod = 1;
1088 else if (ilad == 24)
1089 halfmod = 1;
1090 else if (ilad == 25)
1091 halfmod = 1;
1092
1093 if (!halfmod) {
1094 h207->Fill(xloc, yloc);
1095 }
1096
1097 }
1098
1099 if (ilay == 3) {
1100 n3++;
1101 x3 = gX;
1102 y3 = gY;
1103 z3 = gZ;
1104
1105 h301->Fill(ilad);
1106 h302->Fill(imod);
1107
1108 h303->Fill(gR);
1109 h304->Fill(gF * wt);
1110 h305->Fill(gZ);
1111
1112 h306->Fill(gF * wt, gZ);
1113
1114 if (ilad == 11)
1115 halfmod = 1;
1116 if (ilad == 12)
1117 halfmod = 1;
1118 if (ilad == 33)
1119 halfmod = 1;
1120 if (ilad == 34)
1121 halfmod = 1;
1122
1123 if (!halfmod) {
1124 h307->Fill(xloc, yloc);
1125 }
1126
1127 }
1128
1129 }
1130
1131 }
1132
1133 }
1134
1135
1136
1137
1138 if (n1 * n2 * n3 > 0) {
1139 double dca2 = 0., dz2 = 0.;
1140
1141 double ptsig = pt;
1142 if (charge < 0.)
1143 ptsig = -pt;
1144 triplets(x1, y1, z1, x2, y2, z2, x3, y3, z3, ptsig, dca2, dz2);
1145
1146 if (pt > 4) {
1147 h410->Fill(dca2 * 1E4);
1148 h411->Fill(dz2 * 1E4);
1149 }
1150 if (pt > 12) {
1151 h420->Fill(dca2 * 1E4);
1152 h421->Fill(dz2 * 1E4);
1153 if (hp.trackerLayersWithMeasurement() > 8) {
1154 h430->Fill(dca2 * 1E4);
1155 h431->Fill(dz2 * 1E4);
1156 }
1157
1158
1159
1160
1161 }
1162
1163
1164
1165
1166 if (pt > 4) {
1167
1168
1169
1170 h414->Fill(z2, dca2 * 1E4);
1171 h415->Fill(z2, dz2 * 1E4);
1172 }
1173
1174 h416->Fill(logpt, dca2 * 1E4);
1175 h417->Fill(logpt, dz2 * 1E4);
1176 if (iTrack->charge() > 0)
1177 h418->Fill(logpt, dca2 * 1E4);
1178 else
1179 h419->Fill(logpt, dca2 * 1E4);
1180
1181
1182
1183
1184
1185
1186
1187 if (pt > 4) {
1188
1189
1190
1191 h424->Fill(z2, abs(dca2) * 1E4);
1192 h425->Fill(z2, abs(dz2) * 1E4);
1193
1194
1195
1196
1197
1198
1199
1200 }
1201
1202 h426->Fill(logpt, abs(dca2) * 1E4);
1203 h427->Fill(logpt, abs(dz2) * 1E4);
1204
1205
1206 }
1207
1208 }
1209
1210 sumpt += iTrack->pt();
1211 sumq += iTrack->charge();
1212
1213 }
1214
1215 h028->Fill(sumpt);
1216 h029->Fill(sumq);
1217 }
1218
1219 void Triplet::triplets(double x1,
1220 double y1,
1221 double z1,
1222 double x2,
1223 double y2,
1224 double z2,
1225 double x3,
1226 double y3,
1227 double z3,
1228 double ptsig,
1229 double &dca2,
1230 double &dz2) {
1231 using namespace std;
1232 const double pi = 4 * atan(1);
1233
1234 const double twopi = 2 * pi;
1235
1236
1237
1238 double pt = abs(ptsig);
1239
1240 double kap = 0.0114257 / pt;
1241 if (ptsig > 0)
1242 kap = -kap;
1243
1244 double rho = 1 / kap;
1245 double rinv = -kap;
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268 double r1 = sqrt(x1 * x1 + y1 * y1);
1269 double r3 = sqrt(x3 * x3 + y3 * y3);
1270
1271 if (r3 - r1 < 2.0)
1272 edm::LogPrint("Triplet") << "warn r1 = " << r1 << ", r3 = " << r3;
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283 double lam = sqrt(-0.25 + rho * rho / ((x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3)));
1284
1285
1286
1287
1288 if (kap > 0)
1289 lam = -lam;
1290
1291
1292
1293 double x0 = 0.5 * (x1 + x3) + lam * (-y1 + y3);
1294 double y0 = 0.5 * (y1 + y3) + lam * (x1 - x3);
1295
1296
1297
1298 double num = (y3 - y0) * (x1 - x0) - (x3 - x0) * (y1 - y0);
1299 double den = (x1 - x0) * (x3 - x0) + (y1 - y0) * (y3 - y0);
1300 double tandip = kap * (z3 - z1) / atan(num / den);
1301
1302
1303
1304
1305
1306
1307 double uphi;
1308 if (kap > 0)
1309 uphi = atan2(-x0, y0);
1310 else
1311 uphi = atan2(x0, -y0);
1312
1313
1314
1315
1316 double udca;
1317 if (kap > 0)
1318 udca = rho - sqrt(x0 * x0 + y0 * y0);
1319 else
1320 udca = rho + sqrt(x0 * x0 + y0 * y0);
1321
1322
1323
1324 double dphi = atan(((x1 - x0) * y0 - (y1 - y0) * x0) / ((x1 - x0) * x0 + (y1 - y0) * y0));
1325
1326 double uz0 = z1 + tandip * dphi * rho;
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338 double cosphi = cos(uphi);
1339 double sinphi = sin(uphi);
1340 double dp = -x2 * sinphi + y2 * cosphi + udca;
1341 double dl = -x2 * cosphi - y2 * sinphi;
1342 double sa = 2 * dp + rinv * (dp * dp + dl * dl);
1343 dca2 = sa / (1 + sqrt(1 + rinv * sa));
1344
1345 double ud = 1 + rinv * udca;
1346 double phi2 = atan2(-rinv * x2 + ud * sinphi, rinv * y2 + ud * cosphi);
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358 double xx = x2 + dca2 * sin(phi2);
1359 double yy = y2 - dca2 * cos(phi2);
1360
1361 double f0 = uphi;
1362 double kx = kap * xx;
1363 double ky = kap * yy;
1364 double kd = kap * udca;
1365
1366
1367
1368 double dx = kx - (kd - 1) * sin(f0);
1369 double dy = ky + (kd - 1) * cos(f0);
1370 double ks = atan2(dx, -dy) - f0;
1371
1372
1373
1374 if (ks > pi)
1375 ks = ks - twopi;
1376 else if (ks < -pi)
1377 ks = ks + twopi;
1378
1379 double s = ks * rho;
1380 double uz2 = uz0 + s * tandip;
1381 dz2 = z2 - uz2;
1382 }
1383
1384
1385
1386 void Triplet::endJob() { edm::LogPrint("Triplet") << "end of job after " << myCounters::neve << " events.\n"; }
1387
1388
1389 DEFINE_FWK_MODULE(Triplet);