Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:19:38 

0001 // -*- C++ -*-

0002 //

0003 // Package:  Triplet

0004 // Class:    Triplet

0005 //

0006 // my/Triplet/src/Triplet.cc

0007 //

0008 // plot hits on CMS tracks on RECO

0009 //

0010 // Original Author:  Daniel Pitzl, DESY,,

0011 //         Created:  Sat Feb 12 12:12:42 CET 2011

0012 // $Id$

0013 // d.k.

0014 // Split into a sperate call.

0015 //

0016 
0017 // system include files:

0018 #include <memory>
0019 #include <iostream>
0020 #include <cmath>
0021 
0022 // ROOT:

0023 #include "TH1.h"
0024 #include "TH2.h"
0025 #include "TProfile.h"
0026 
0027 // CMS and user include files:

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 //#include <FWCore/Framework/interface/EventSetup.h>

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 //#include "DataFormats/METReco/interface/PFMET.h"

0042 //#include "DataFormats/METReco/interface/PFMETFwd.h"

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 // To convert detId to subdet/layer number:

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 //#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

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 //#include "TrackingTools/TransientTrack/interface/TransientTrack.h"

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 //#define SIM // use for single muon simulations

0075 
0076 //

0077 // class declaration:

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   // ----------member data:

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 // constants, enums and typedefs:

0169 //

0170 
0171 //

0172 // static data member definitions:

0173 //

0174 
0175 //

0176 // constructor:

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 // destructor:

0195 //

0196 Triplet::~Triplet() = default;
0197 
0198 //

0199 // member functions:

0200 // method called once each job just before starting event loop

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   // triplets:

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   // mean resid profiles:

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   // width profiles:

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 // method called for each event:

0377 //

0378 void Triplet::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0379   //Retrieve tracker topology from geometry

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   //const double sqrtpihalf = sqrt(pihalf);

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);  // ctime has endline

0401 
0402     myCounters::prevrun = iEvent.run();
0403 
0404   }  // new run

0405 
0406   int idbg = 0;
0407   if (myCounters::neve < 1)
0408     idbg = 1;
0409   //

0410   //--------------------------------------------------------------------

0411   // beam spot:

0412   //

0413   edm::Handle<reco::BeamSpot> rbs;
0414   iEvent.getByToken(bsToken_, rbs);
0415 
0416   XYZPoint bsP = XYZPoint(0, 0, 0);
0417   //int ibs = 0;

0418 
0419   if (!rbs.failedToGet() && rbs.isValid()) {
0420     //ibs = 1;

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   }  // bs valid

0441   //

0442   //--------------------------------------------------------------------

0443   // primary vertices:

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   // need vertex global point for tracks

0456   // from #include "DataFormats/GeometryVector/interface/GlobalPoint.h"

0457   // Global points are three-dimensional by default

0458   // typedef Global3DPoint  GlobalPoint;

0459   // typedef Point3DBase< float, GlobalTag> Global3DPoint;

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());  // empty in Zmumu sample

0494         else
0495           h019->Fill(iVertex->z());  // all here: why?

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       }  // non-fake

0510     }    //valid

0511   }      // loop over vertices

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   if( ibs ) {

0524     vtxP.SetX( bsP.x() ); // beam spot. should take tilt into account!

0525     vtxP.SetY( bsP.y() ); // beam spot. should take tilt into account!

0526   }

0527   */
0528 
0529   h023->Fill(vtxP.x());
0530   h024->Fill(vtxP.y());
0531   h025->Fill(vtxP.z());
0532 
0533   //double xbs = 0;

0534   //double ybs = 0;

0535   //if( ibs ) {

0536   //xbs = bsP.x();

0537   //ybs = bsP.y();

0538   //}

0539   //else {

0540   //xbs = vtxP.x();

0541   //ybs = vtxP.y();

0542   //}

0543   //

0544   //--------------------------------------------------------------------

0545   // tracks:

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   // get tracker geometry:

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   // loop over tracker detectors:

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         // normal vector: includes alignment (varies from module to module along z on one ladder)

0602         // neighbouring ladders alternate with inward/outward orientation

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       }  //PXB

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       }  //PXD

0644 
0645       edm::LogPrint("Triplet");
0646 
0647     }  //idbg

0648 
0649   }  //idet

0650 
0651   //

0652   // transient track builder, needs B-field from data base (global tag in .py)

0653   //

0654   edm::ESHandle<TransientTrackBuilder> theB = iSetup.getHandle(transientTrackBuilderToken_);
0655   //

0656   // transient rec hits:

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   // Tracks:

0669   //

0670   for (TrackCollection::const_iterator iTrack = tracks->begin(); iTrack != tracks->end(); ++iTrack) {
0671     kk++;
0672 
0673     // cpt = cqRB = 0.3*R[m]*B[T] = 1.14*R[m] for B=3.8T

0674     // D = 2R = 2*pt/1.14

0675     // calo: D = 1.3 m => pt = 0.74 GeV/c

0676 
0677     double pt = iTrack->pt();
0678 
0679     if (pt < 0.75)
0680       continue;  // curls up

0681     //if( pt < 1.75 ) continue;// want sharper image

0682 
0683     if (abs(iTrack->dxy(vtxP)) > 5 * iTrack->dxyError())
0684       continue;  // not prompt

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();  // w.r.t. origin

0711     //double dca = -iTrack->dxy(); // dxy = -d0

0712     double dip = iTrack->lambda();
0713     double z0 = iTrack->dz();
0714     double tet = pihalf - dip;
0715     //double eta = iTrack->eta();

0716     //

0717     // transient track:

0718     //

0719     TransientTrack tTrack = theB->build(*iTrack);
0720 
0721     double kap = tTrack.initialFreeState().transverseCurvature();
0722 
0723     TrajectoryStateClosestToPoint tscp = tTrack.trajectoryStateClosestToPoint(origin);
0724 
0725     //edm::LogPrint("Triplet")<<pt<<" "<<kap<<" ";

0726 
0727     if (tscp.isValid()) {
0728       h057->Fill(tscp.referencePoint().x());  // 0.0

0729       h058->Fill(tscp.referencePoint().y());  // 0.0

0730       h059->Fill(tscp.referencePoint().z());  // 0.0

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     //edm::LogPrint("Triplet")<<pt<<" "<<kap<<" "<<tmp<<" "<<charge<<" "<<kap1<<endl;

0755 
0756     double cf = cos(phi);
0757     double sf = sin(phi);
0758     //double xdca =  dca * sf;

0759     //double ydca = -dca * cf;

0760 
0761     //double tt = tan(tet);

0762 
0763     //double rinv = -kap; // Karimaki

0764     //double rho = 1/kap;

0765     double erd = 1.0 - kap * dca;
0766     double drd = dca * (0.5 * kap * dca - 1.0);  // 0.5 * kap * dca**2 - dca;

0767     double hkk = 0.5 * kap * kap;
0768     //

0769     // track w.r.t. beam (cirmov):

0770     //

0771     //double dp = -xbs*sf + ybs*cf + dca;

0772     //double dl = -xbs*cf - ybs*sf;

0773     //double sa = 2*dp + rinv*(dp*dp+dl*dl);

0774     //double dcap = sa / ( 1 + sqrt(1 + rinv*sa) );// distance to beam

0775     //double ud = 1 + rinv*dca;

0776     //double phip = atan2( -rinv*xbs + ud*sf, rinv*ybs + ud*cf );//direction

0777     //

0778     // track at R(PXB1), from FUNPHI, FUNLEN:

0779     //

0780     double R1 = 4.4;  // PXB1

0781 
0782     double s1 = 0;
0783     double fpos1 = phi - pihalf;
0784 
0785     if (R1 >= abs(dca)) {
0786       //

0787       // sin(delta phi):

0788       //

0789       double sindp = (0.5 * kap * (R1 * R1 + dca * dca) - dca) / (R1 * erd);
0790       fpos1 = phi + asin(sindp);  // phi position at R1

0791       //

0792       // sin(alpha):

0793       //

0794       double sina = R1 * kap * sqrt(1.0 - sindp * sindp);
0795       //

0796       // s = alpha / kappa:

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;  //always positive

0805       }
0806       //

0807       // Check direction: limit to half-turn

0808       //

0809       if (hkk * (R1 * R1 - dca * dca) > erd)
0810         s1 = pi / abs(kap) - s1;  // always positive

0811 
0812     }  // R1 > dca

0813 
0814     if (fpos1 > pi)
0815       fpos1 -= twopi;
0816     else if (fpos1 < -pi)
0817       fpos1 += twopi;
0818 
0819     double zR1 = z0 + s1 * tan(dip);  // z at R1

0820     //

0821     //--------------------------------------------------------------------------

0822     // loop over tracker detectors:

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     edm::LogPrint("Triplet") << ": PXB layer " << tTopo->pxbLayer(mydetId);

0837     edm::LogPrint("Triplet") << ", ladder " << tTopo->pxbLadder(mydetId);

0838     edm::LogPrint("Triplet") << ", module " << tTopo->pxbModule(mydetId);

0839     edm::LogPrint("Triplet") << ", at R1 " << (*idet)->position().perp();

0840     edm::LogPrint("Triplet") << ", F " << (*idet)->position().barePhi()*wt;

0841     edm::LogPrint("Triplet") << ", z " << (*idet)->position().z();

0842     edm::LogPrint("Triplet") ;

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         // normal vector: includes alignment (varies from module to module along z on one ladder)

0862         // neighbouring ladders alternate with inward/outward orientation

0863         //

0864         /*

0865       edm::LogPrint("Triplet") << "normal";

0866       edm::LogPrint("Triplet") << ": x " << (*idet)->surface().normalVector().x();

0867       edm::LogPrint("Triplet") << ", y " << (*idet)->surface().normalVector().y();

0868       edm::LogPrint("Triplet") << ", z " << (*idet)->surface().normalVector().z();

0869       edm::LogPrint("Triplet") << ", f " << (*idet)->surface().normalVector().barePhi()*wt;

0870       edm::LogPrint("Triplet") ;

0871     */
0872 
0873         double phiN = (*idet)->surface().normalVector().barePhi();  //normal vector

0874 
0875         double phidet = phiN - pihalf;  // orientation of sensor plane in x-y

0876         double ux = cos(phidet);        // vector in sensor plane

0877         double uy = sin(phidet);
0878         double x = (*idet)->position().x();
0879         double y = (*idet)->position().y();
0880         //

0881         // intersect helix with line: FUNRXY (in FUNPHI) from V. Blobel

0882         // factor f for intersect point (x + f*ux, y + f*uy)

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         }  //dis

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);  // z at r

0931 
0932         ncrss++;
0933 
0934       }  //PXB1

0935 
0936     }  //idet

0937     //

0938     //--------------------------------------------------------------------------

0939     // rec hits from track extra:

0940     //

0941     if (iTrack->extra().isNonnull() && iTrack->extra().isAvailable()) {
0942       h044->Fill(tTrack.recHitsSize());  // tTrack

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       //double phiN2 = 0;

0957 
0958       for (trackingRecHit_iterator irecHit = iTrack->recHitsBegin(); irecHit != iTrack->recHitsEnd(); ++irecHit) {
0959         if ((*irecHit)->isValid()) {
0960           DetId detId = (*irecHit)->geographicalId();
0961 
0962           // enum Detector { Tracker=1, Muon=2, Ecal=3, Hcal=4, Calo=5 };

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           // enum SubDetector{ PixelBarrel=1, PixelEndcap=2 };

0972           // enum SubDetector{ TIB=3, TID=4, TOB=5, TEC=6 };

0973 
0974           h060->Fill(subDet);
0975           //

0976           // build hit: (from what?)

0977           //

0978           TransientTrackingRecHit::RecHitPointer trecHit = hitBuilder->build(&*(*irecHit));
0979 
0980           double xloc = trecHit->localPosition().x();  // 1st meas coord

0981           double yloc = trecHit->localPosition().y();  // 2nd meas coord or zero

0982           //double zloc = trecHit->localPosition().z();// up, always zero

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           //      GeomDet* igeomdet = trecHit->det();

0997           //double phiN = trecHit->det()->surface().normalVector().barePhi();//normal vector

0998 
0999           if (subDet == PixelSubdetector::PixelBarrel || subDet == StripSubdetector::TIB ||
1000               subDet == StripSubdetector::TOB) {  // barrel

1001 
1002             h066->Fill(gX, gY);
1003 
1004           }  //barrel

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);  // 1,2,3

1013 
1014             if (ilay == 1) {
1015               n1++;
1016               x1 = gX;
1017               y1 = gY;
1018               z1 = gZ;
1019 
1020               h101->Fill(ilad);  // 1..20

1021               h102->Fill(imod);  // 1..8

1022 
1023               h103->Fill(gR);
1024               h104->Fill(gF * wt);
1025               h105->Fill(gZ);
1026 
1027               h106->Fill(gF * wt, gZ);  // phi-z of hit

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);  // hit within one module

1040               }
1041               //

1042               // my crossings:

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               }  //crss

1065 
1066             }  //PXB1

1067 
1068             if (ilay == 2) {
1069               n2++;
1070               x2 = gX;
1071               y2 = gY;
1072               z2 = gZ;
1073               //phiN2 = phiN;

1074 
1075               h201->Fill(ilad);  // 1..32

1076               h202->Fill(imod);  //1..8

1077 
1078               h203->Fill(gR);
1079               h204->Fill(gF * wt);
1080               h205->Fill(gZ);
1081 
1082               h206->Fill(gF * wt, gZ);  // phi-z of hit

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);  // hit within one module

1095               }
1096 
1097             }  //PXB2

1098 
1099             if (ilay == 3) {
1100               n3++;
1101               x3 = gX;
1102               y3 = gY;
1103               z3 = gZ;
1104 
1105               h301->Fill(ilad);  //1..44

1106               h302->Fill(imod);  //1..8

1107 
1108               h303->Fill(gR);
1109               h304->Fill(gF * wt);
1110               h305->Fill(gZ);
1111 
1112               h306->Fill(gF * wt, gZ);  // phi-z of hit

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);  // hit within one module

1125               }
1126 
1127             }  //PXB3

1128 
1129           }  //PXB

1130 
1131         }  //valid

1132 
1133       }  //loop rechits

1134       //

1135       // 1-2-3 triplet:

1136       //

1137       //if( hp.pixelBarrelLayersWithMeasurement() == 3 ){

1138       if (n1 * n2 * n3 > 0) {
1139         double dca2 = 0., dz2 = 0.;
1140         //triplets(x1,y1,z1,x2,y2,z2,x3,y3,z3,kap,dca2,dz2);

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           //if( phiinc*wt > -1 && phiinc*wt < 7 ){

1158           //h450->Fill( dca2*1E4 );

1159           //h451->Fill( dz2*1E4 );

1160           //}

1161         }
1162 
1163         //

1164         // residual profiles: alignment check

1165         //

1166         if (pt > 4) {
1167           //h412->Fill( f2*wt, dca2*1E4 );

1168           //h413->Fill( f2*wt, dz2*1E4 );

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         // profile of abs(dca) gives mean abs(dca):

1182         // mean of abs(Gauss) = 0.7979 * RMS = 1/sqrt(pi/2)

1183         // => rms = sqrt(pi/2) * mean of abs (sqrt(pi/2) = 1.2533)

1184         // point resolution = 1/sqrt(3/2) * triplet middle residual width

1185         // => sqrt(pi/2)*sqrt(2/3) = sqrt(pi/3) = 1.0233, almost one

1186 
1187         if (pt > 4) {
1188           //h422->Fill( f2*wt, abs(dca2)*1E4 );

1189           //h423->Fill( f2*wt, abs(dz2)*1E4 );

1190 
1191           h424->Fill(z2, abs(dca2) * 1E4);
1192           h425->Fill(z2, abs(dz2) * 1E4);
1193 
1194           //h428->Fill( udip*wt, abs(dca2)*1E4 );

1195           //h429->Fill( udip*wt, abs(dz2)*1E4 );

1196           //if( abs(dip)*wt > 18 && abs(dip)*wt < 50 ) h432->Fill( dz2*1E4 );

1197 
1198           //h434->Fill( phiinc*wt, abs(dca2)*1E4 );

1199 
1200         }  //pt

1201 
1202         h426->Fill(logpt, abs(dca2) * 1E4);
1203         h427->Fill(logpt, abs(dz2) * 1E4);
1204         //if( abs(dip)*wt > 18 && abs(dip)*wt < 50 ) h433->Fill( logpt, abs(dz2)*1E4 );

1205 
1206       }  //3 PXB layers

1207 
1208     }  //extra

1209 
1210     sumpt += iTrack->pt();
1211     sumq += iTrack->charge();
1212 
1213   }  // loop over tracks

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   //const double wt = 180/pi;

1234   const double twopi = 2 * pi;
1235   //const double pihalf = 2*atan(1);

1236   //const double sqrtpihalf = sqrt(pihalf);

1237 
1238   double pt = abs(ptsig);
1239   //double rho = pt/0.0114257;

1240   double kap = 0.0114257 / pt;
1241   if (ptsig > 0)
1242     kap = -kap;  // kap i snegative for positive charge

1243 
1244   double rho = 1 / kap;
1245   double rinv = -kap;  // Karimaki

1246 
1247   //double f2 = atan2( y2, x2 );//position angle

1248 
1249   //h406->Fill( hp.numberOfValidTrackerHits() );

1250   //h407->Fill( hp.numberOfValidPixelBarrelHits() );

1251   //h408->Fill( hp.trackerLayersWithMeasurement() );

1252 
1253   // Author: Johannes Gassner (15.11.1996)

1254   // Make track from 2 space points and kappa (cmz98/ftn/csmktr.f)

1255   // Definition of the Helix :

1256   //

1257   // x( t ) = X0 + KAPPA^-1 * SIN( PHI0 + t )

1258   // y( t ) = Y0 - KAPPA^-1 * COS( PHI0 + t )          t > 0

1259   // z( t ) = Z0 + KAPPA^-1 * TAN( DIP ) * t

1260 
1261   // Center of the helix in the xy-projection:

1262 
1263   // X0 = + ( DCA - KAPPA^-1 ) * SIN( PHI0 )

1264   // Y0 = - ( DCA - KAPPA^-1 ) * COS( PHI0 )

1265   //

1266   // Point 1 must be in the inner layer, 3 in the outer:

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   // Calculate the centre of the helix in xy-projection that

1275   // transverses the two spacepoints. The points with the same

1276   // distance from the two points are lying on a line.

1277   // LAMBDA is the distance between the point in the middle of

1278   // the two spacepoints and the centre of the helix.

1279   //

1280   // we already have kap and rho = 1/kap

1281   //

1282 
1283   double lam = sqrt(-0.25 + rho * rho / ((x1 - x3) * (x1 - x3) + (y1 - y3) * (y1 - y3)));
1284   //

1285   // There are two solutions, the sign of kap gives the information

1286   // which of them is correct.

1287   //

1288   if (kap > 0)
1289     lam = -lam;
1290   //

1291   // ( X0, Y0 ) is the centre of the circle that describes the helix in xy-projection.

1292   //

1293   double x0 = 0.5 * (x1 + x3) + lam * (-y1 + y3);
1294   double y0 = 0.5 * (y1 + y3) + lam * (x1 - x3);
1295   //

1296   // Calculate theta :

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   //double udip = atan(tandip);

1302   //double utet = pihalf - udip;

1303   //

1304   // To get phi0 in the right intervall one must differ two cases

1305   // with positve and negative kap:

1306   //

1307   double uphi;
1308   if (kap > 0)
1309     uphi = atan2(-x0, y0);
1310   else
1311     uphi = atan2(x0, -y0);
1312   //

1313   // The distance of the closest approach DCA depends on the sign

1314   // of kappa :

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   // angle from first hit to dca point:

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   //h401->Fill( z2 );

1329   //h402->Fill( uphi - iTrack->phi() );

1330   //h403->Fill( udca - iTrack->d0() );

1331   //h404->Fill( udip - iTrack->lambda() );

1332   //h405->Fill( uz0  - iTrack->dz() );

1333   //

1334   // interpolate to middle hit:

1335   // cirmov

1336   // we already have rinv = -kap

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));  // distance to hit 2

1344 
1345   double ud = 1 + rinv * udca;
1346   double phi2 = atan2(-rinv * x2 + ud * sinphi, rinv * y2 + ud * cosphi);  //direction

1347 
1348   //double phiinc = phi2 - phiN2;//angle of incidence in rphi w.r.t. normal vector

1349   //

1350   // phiN alternates inward/outward

1351   // reduce phiinc

1352   //if( phiinc > pihalf ) phiinc -= pi;

1353   //else if( phiinc < -pihalf ) phiinc += pi;

1354   //h409->Fill( f2*wt, phiinc*wt );

1355   //

1356   // arc length:

1357   //

1358   double xx = x2 + dca2 * sin(phi2);  // point on track

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   // Solve track equation for s:

1367   //

1368   double dx = kx - (kd - 1) * sin(f0);
1369   double dy = ky + (kd - 1) * cos(f0);
1370   double ks = atan2(dx, -dy) - f0;  // turn angle

1371   //

1372   //---  Limit to half-turn:

1373   //

1374   if (ks > pi)
1375     ks = ks - twopi;
1376   else if (ks < -pi)
1377     ks = ks + twopi;
1378 
1379   double s = ks * rho;            // signed

1380   double uz2 = uz0 + s * tandip;  //track z at R2

1381   dz2 = z2 - uz2;
1382 }
1383 //----------------------------------------------------------------------

1384 // method called just after ending the event loop:

1385 //

1386 void Triplet::endJob() { edm::LogPrint("Triplet") << "end of job after " << myCounters::neve << " events.\n"; }
1387 
1388 //define this as a plug-in

1389 DEFINE_FWK_MODULE(Triplet);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!