File indexing completed on 2022-01-17 23:31:48
0001
0002
0003
0004
0005
0006
0007 #include <memory>
0008 #include <string>
0009 #include <iostream>
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "FWCore/Framework/interface/ESHandle.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014
0015 #include "FWCore/Framework/interface/Frameworkfwd.h"
0016 #include "FWCore/Framework/interface/EDAnalyzer.h"
0017
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020
0021
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025
0026
0027 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0028 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0029 #include "DataFormats/Common/interface/DetSetVector.h"
0030 #include "DataFormats/Common/interface/Ref.h"
0031 #include "DataFormats/DetId/interface/DetId.h"
0032
0033 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
0034 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
0035 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0036
0037 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0038 #include "Geometry/CommonDetUnit/interface/PixelGeomDetType.h"
0039 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0041 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0042 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0043
0044
0045
0046 #include "L1Trigger/GlobalTriggerAnalyzer/interface/L1GtUtils.h"
0047 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0048 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
0049 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0050 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0051
0052
0053 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
0054 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0055 #include "DataFormats/Common/interface/TriggerResults.h"
0056 #include "HLTrigger/HLTcore/interface/HLTConfigProvider.h"
0057 #include "FWCore/Common/interface/TriggerNames.h"
0058
0059
0060 #include "DataFormats/TrackReco/interface/Track.h"
0061 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0062
0063
0064
0065 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0066
0067 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0068
0069
0070 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0084 #include <DataFormats/VertexReco/interface/VertexFwd.h>
0085
0086
0087 #include "DataFormats/Luminosity/interface/LumiSummary.h"
0088 #include "DataFormats/Common/interface/ConditionsInEdm.h"
0089
0090
0091 #include "FWCore/ServiceRegistry/interface/Service.h"
0092 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0093
0094
0095 #include <TROOT.h>
0096 #include <TChain.h>
0097 #include <TFile.h>
0098 #include <TF1.h>
0099 #include <TH2F.h>
0100 #include <TH1F.h>
0101 #include <TProfile.h>
0102 #include <TVector3.h>
0103
0104 #define HISTOS
0105
0106
0107
0108 using namespace std;
0109
0110 class TestPixTracks : public edm::EDAnalyzer {
0111 public:
0112 explicit TestPixTracks(const edm::ParameterSet &conf);
0113 virtual ~TestPixTracks();
0114 virtual void analyze(const edm::Event &e, const edm::EventSetup &c) override;
0115 virtual void beginRun(const edm::EventSetup &iSetup) override;
0116 virtual void beginJob() override;
0117 virtual void endJob() override;
0118
0119 private:
0120 edm::EDGetTokenT<LumiSummary> lumiToken_;
0121 edm::EDGetTokenT<edm::ConditionsInLumiBlock> condToken_;
0122 edm::EDGetTokenT<L1GlobalTriggerReadoutRecord> l1gtrrToken_;
0123 edm::EDGetTokenT<edm::TriggerResults> hltToken_;
0124 edm::EDGetTokenT<reco::VertexCollection> vtxToken_;
0125 edm::EDGetTokenT<reco::TrackCollection> srcToken_;
0126 edm::EDGetTokenT<TrajTrackAssociationCollection> trackAssocToken_;
0127 edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerGeomToken_;
0128
0129 bool PRINT;
0130 float countTracks, countGoodTracks, countTracksInPix, countPVs, countEvents, countLumi;
0131
0132
0133
0134
0135
0136
0137
0138
0139 TH1D *hcharge1, *hcharge2, *hcharge3, *hcharge;
0140 TH1D *hpixcharge1, *hpixcharge2, *hpixcharge3, *hpixcharge;
0141 TH1D *hcols1, *hcols2, *hcols3, *hrows1, *hrows2, *hrows3;
0142 TH1D *hsize1, *hsize2, *hsize3, *hsizex1, *hsizex2, *hsizex3, *hsizey1, *hsizey2, *hsizey3;
0143
0144 TH1D *hclusPerTrk1, *hclusPerTrk2, *hclusPerTrk3;
0145 TH1D *hclusPerLay1, *hclusPerLay2, *hclusPerLay3;
0146 TH1D *hpixPerLay1, *hpixPerLay2, *hpixPerLay3;
0147
0148
0149
0150
0151 TH2F *hDetMap1, *hDetMap2, *hDetMap3;
0152
0153 TH2F *hcluDetMap1, *hcluDetMap2, *hcluDetMap3;
0154
0155 TH2F *hpvxy, *hclusMap1, *hclusMap2, *hclusMap3;
0156
0157 TH1D *hpvz, *hpvr, *hNumPv, *hNumPvClean;
0158 TH1D *hPt, *hEta, *hDz, *hD0, *hzdiff;
0159
0160
0161
0162
0163
0164 TH1D *hl1a, *hl1t, *hlt1;
0165
0166 TH1D *hclusBpix, *hpixBpix;
0167 TH1D *htracks, *htracksGood, *htracksGoodInPix;
0168
0169 TProfile *hclumult1, *hclumult2, *hclumult3;
0170 TProfile *hclumultx1, *hclumultx2, *hclumultx3;
0171 TProfile *hclumulty1, *hclumulty2, *hclumulty3;
0172 TProfile *hcluchar1, *hcluchar2, *hcluchar3;
0173 TProfile *hpixchar1, *hpixchar2, *hpixchar3;
0174
0175 TProfile *htracksls, *hpvsls, *htrackslsn, *hpvslsn, *hintgl, *hinstl, *hbeam1, *hbeam2;
0176
0177 TH1D *hlumi, *hlumi0, *hbx, *hbx0;
0178 };
0179
0180
0181 TestPixTracks::TestPixTracks(edm::ParameterSet const &conf) {
0182 PRINT = conf.getUntrackedParameter<bool>("Verbosity", false);
0183 lumiToken_ = consumes<LumiSummary>(edm::InputTag("lumiProducer"));
0184 condToken_ = consumes<edm::ConditionsInLumiBlock>(edm::InputTag("conditionsInEdm"));
0185 l1gtrrToken_ = consumes<L1GlobalTriggerReadoutRecord>(edm::InputTag("gtDigis"));
0186 hltToken_ = consumes<edm::TriggerResults>(edm::InputTag("TriggerResults", "", "HLT"));
0187 vtxToken_ = consumes<reco::VertexCollection>(edm::InputTag("offlinePrimaryVertices"));
0188 srcToken_ = consumes<reco::TrackCollection>(conf.getParameter<edm::InputTag>("src"));
0189 trackAssocToken_ =
0190 consumes<TrajTrackAssociationCollection>(edm::InputTag(conf.getParameter<std::string>("trajectoryInput")));
0191 trackerGeomToken_ = esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>();
0192
0193 }
0194
0195
0196 TestPixTracks::~TestPixTracks() {}
0197
0198
0199 void TestPixTracks::beginRun(const edm::EventSetup &iSetup) { cout << "BeginRun, Verbosity = " << PRINT << endl; }
0200
0201
0202 void TestPixTracks::beginJob() {
0203 cout << "BeginJob, Verbosity " << PRINT << endl;
0204
0205 countTracks = 0.;
0206 countGoodTracks = 0.;
0207 countTracksInPix = 0.;
0208 countPVs = 0.;
0209 countEvents = 0.;
0210 countLumi = 0.;
0211
0212 #ifdef HISTOS
0213
0214
0215 edm::Service<TFileService> fs;
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227 int sizeH = 20;
0228 float lowH = -0.5;
0229 float highH = 19.5;
0230
0231 hclusPerTrk1 = fs->make<TH1D>("hclusPerTrk1", "Clus per track l1", sizeH, lowH, highH);
0232 hclusPerTrk2 = fs->make<TH1D>("hclusPerTrk2", "Clus per track l2", sizeH, lowH, highH);
0233 hclusPerTrk3 = fs->make<TH1D>("hclusPerTrk3", "Clus per track l3", sizeH, lowH, highH);
0234
0235 sizeH = 2000;
0236 highH = 1999.5;
0237 hclusPerLay1 = fs->make<TH1D>("hclusPerLay1", "Clus per layer l1", sizeH, lowH, highH);
0238 hclusPerLay2 = fs->make<TH1D>("hclusPerLay2", "Clus per layer l2", sizeH, lowH, highH);
0239 hclusPerLay3 = fs->make<TH1D>("hclusPerLay3", "Clus per layer l3", sizeH, lowH, highH);
0240
0241 highH = 9999.5;
0242 hpixPerLay1 = fs->make<TH1D>("hpixPerLay1", "Pix per layer l1", sizeH, lowH, highH);
0243 hpixPerLay2 = fs->make<TH1D>("hpixPerLay2", "Pix per layer l2", sizeH, lowH, highH);
0244 hpixPerLay3 = fs->make<TH1D>("hpixPerLay3", "Pix per layer l3", sizeH, lowH, highH);
0245
0246
0247
0248
0249
0250
0251
0252
0253 hcharge1 = fs->make<TH1D>("hcharge1", "Clu charge l1", 400, 0., 200.);
0254 hcharge2 = fs->make<TH1D>("hcharge2", "Clu charge l2", 400, 0., 200.);
0255 hcharge3 = fs->make<TH1D>("hcharge3", "Clu charge l3", 400, 0., 200.);
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272
0273
0274
0275
0276 hcols1 = fs->make<TH1D>("hcols1", "Layer 1 cols", 500, -0.5, 499.5);
0277 hcols2 = fs->make<TH1D>("hcols2", "Layer 2 cols", 500, -0.5, 499.5);
0278 hcols3 = fs->make<TH1D>("hcols3", "Layer 3 cols", 500, -0.5, 499.5);
0279
0280 hrows1 = fs->make<TH1D>("hrows1", "Layer 1 rows", 200, -0.5, 199.5);
0281 hrows2 = fs->make<TH1D>("hrows2", "Layer 2 rows", 200, -0.5, 199.5);
0282 hrows3 = fs->make<TH1D>("hrows3", "layer 3 rows", 200, -0.5, 199.5);
0283
0284 hsize1 = fs->make<TH1D>("hsize1", "layer 1 clu size", 100, -0.5, 99.5);
0285 hsize2 = fs->make<TH1D>("hsize2", "layer 2 clu size", 100, -0.5, 99.5);
0286 hsize3 = fs->make<TH1D>("hsize3", "layer 3 clu size", 100, -0.5, 99.5);
0287 hsizex1 = fs->make<TH1D>("hsizex1", "lay1 clu size in x", 20, -0.5, 19.5);
0288 hsizex2 = fs->make<TH1D>("hsizex2", "lay2 clu size in x", 20, -0.5, 19.5);
0289 hsizex3 = fs->make<TH1D>("hsizex3", "lay3 clu size in x", 20, -0.5, 19.5);
0290 hsizey1 = fs->make<TH1D>("hsizey1", "lay1 clu size in y", 30, -0.5, 29.5);
0291 hsizey2 = fs->make<TH1D>("hsizey2", "lay2 clu size in y", 30, -0.5, 29.5);
0292 hsizey3 = fs->make<TH1D>("hsizey3", "lay3 clu size in y", 30, -0.5, 29.5);
0293
0294 hDetMap1 = fs->make<TH2F>("hDetMap1", "layer 1 clus map", 9, 0., 9., 23, 0., 23.);
0295 hDetMap2 = fs->make<TH2F>("hDetMap2", "layer 2 clus map", 9, 0., 9., 33, 0., 33.);
0296 hDetMap3 = fs->make<TH2F>("hDetMap3", "layer 3 clus map", 9, 0., 9., 45, 0., 45.);
0297
0298
0299
0300
0301
0302
0303 hcluDetMap1 = fs->make<TH2F>("hcluDetMap1", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0304 hcluDetMap2 = fs->make<TH2F>("hcluDetMap2", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0305 hcluDetMap3 = fs->make<TH2F>("hcluDetMap3", "clu det layer 1", 416, 0., 416., 160, 0., 160.);
0306
0307 htracksGoodInPix = fs->make<TH1D>("htracksGoodInPix", "count good tracks in pix", 2000, -0.5, 1999.5);
0308 htracksGood = fs->make<TH1D>("htracksGood", "count good tracks", 2000, -0.5, 1999.5);
0309 htracks = fs->make<TH1D>("htracks", "count tracks", 2000, -0.5, 1999.5);
0310 hclusBpix = fs->make<TH1D>("hclusBpix", "count clus in bpix", 200, -0.5, 1999.5);
0311 hpixBpix = fs->make<TH1D>("hpixBpix", "count pixels", 200, -0.5, 1999.5);
0312
0313 hpvxy = fs->make<TH2F>("hpvxy", "pv xy", 100, -1., 1., 100, -1., 1.);
0314 hpvz = fs->make<TH1D>("hpvz", "pv z", 1000, -50., 50.);
0315 hpvr = fs->make<TH1D>("hpvr", "pv r", 100, 0., 1.);
0316 hNumPv = fs->make<TH1D>("hNumPv", "num of pv", 100, 0., 100.);
0317 hNumPvClean = fs->make<TH1D>("hNumPvClean", "num of pv clean", 100, 0., 100.);
0318
0319 hPt = fs->make<TH1D>("hPt", "pt", 100, 0., 100.);
0320 hEta = fs->make<TH1D>("hEta", "eta", 50, -2.5, 2.5);
0321 hD0 = fs->make<TH1D>("hD0", "d0", 500, 0., 5.);
0322 hDz = fs->make<TH1D>("hDz", "pt", 250, -25., 25.);
0323 hzdiff = fs->make<TH1D>("hzdiff", "PVz-Trackz", 200, -10., 10.);
0324
0325 hl1a = fs->make<TH1D>("hl1a", "l1a", 128, -0.5, 127.5);
0326 hl1t = fs->make<TH1D>("hl1t", "l1t", 128, -0.5, 127.5);
0327 hlt1 = fs->make<TH1D>("hlt1", "hlt1", 256, -0.5, 255.5);
0328
0329 hclumult1 = fs->make<TProfile>("hclumult1", "cluster size layer 1", 60, -3., 3., 0.0, 100.);
0330 hclumult2 = fs->make<TProfile>("hclumult2", "cluster size layer 2", 60, -3., 3., 0.0, 100.);
0331 hclumult3 = fs->make<TProfile>("hclumult3", "cluster size layer 3", 60, -3., 3., 0.0, 100.);
0332
0333 hclumultx1 = fs->make<TProfile>("hclumultx1", "cluster x-size layer 1", 60, -3., 3., 0.0, 100.);
0334 hclumultx2 = fs->make<TProfile>("hclumultx2", "cluster x-size layer 2", 60, -3., 3., 0.0, 100.);
0335 hclumultx3 = fs->make<TProfile>("hclumultx3", "cluster x-size layer 3", 60, -3., 3., 0.0, 100.);
0336
0337 hclumulty1 = fs->make<TProfile>("hclumulty1", "cluster y-size layer 1", 60, -3., 3., 0.0, 100.);
0338 hclumulty2 = fs->make<TProfile>("hclumulty2", "cluster y-size layer 2", 60, -3., 3., 0.0, 100.);
0339 hclumulty3 = fs->make<TProfile>("hclumulty3", "cluster y-size layer 3", 60, -3., 3., 0.0, 100.);
0340
0341 hcluchar1 = fs->make<TProfile>("hcluchar1", "cluster char layer 1", 60, -3., 3., 0.0, 1000.);
0342 hcluchar2 = fs->make<TProfile>("hcluchar2", "cluster char layer 2", 60, -3., 3., 0.0, 1000.);
0343 hcluchar3 = fs->make<TProfile>("hcluchar3", "cluster char layer 3", 60, -3., 3., 0.0, 1000.);
0344
0345 hpixchar1 = fs->make<TProfile>("hpixchar1", "pix char layer 1", 60, -3., 3., 0.0, 1000.);
0346 hpixchar2 = fs->make<TProfile>("hpixchar2", "pix char layer 2", 60, -3., 3., 0.0, 1000.);
0347 hpixchar3 = fs->make<TProfile>("hpixchar3", "pix char layer 3", 60, -3., 3., 0.0, 1000.);
0348
0349 hintgl = fs->make<TProfile>("hintgl", "inst lumi vs ls ", 1000, 0., 3000., 0.0, 10000.);
0350 hinstl = fs->make<TProfile>("hinstl", "intg lumi vs ls ", 1000, 0., 3000., 0.0, 100.);
0351 hbeam1 = fs->make<TProfile>("hbeam1", "beam1 vs ls ", 1000, 0., 3000., 0.0, 1000.);
0352 hbeam2 = fs->make<TProfile>("hbeam2", "beam2 vs ls ", 1000, 0., 3000., 0.0, 1000.);
0353
0354 htracksls = fs->make<TProfile>("htracksls", "tracks with pix hits vs ls", 1000, 0., 3000., 0.0, 10000.);
0355 hpvsls = fs->make<TProfile>("hpvsls", "pvs vs ls", 1000, 0., 3000., 0.0, 1000.);
0356 htrackslsn = fs->make<TProfile>("htrackslsn", "tracks with pix hits/lumi vs ls", 1000, 0., 3000., 0.0, 10000.);
0357 hpvslsn = fs->make<TProfile>("hpvslsn", "pvs/lumi vs ls", 1000, 0., 3000., 0.0, 1000.);
0358
0359 hlumi0 = fs->make<TH1D>("hlumi0", "lumi", 2000, 0, 2000.);
0360 hlumi = fs->make<TH1D>("hlumi", "lumi", 2000, 0, 2000.);
0361 hbx0 = fs->make<TH1D>("hbx0", "bx", 4000, 0, 4000.);
0362 hbx = fs->make<TH1D>("hbx", "bx", 4000, 0, 4000.);
0363
0364 hclusMap1 = fs->make<TH2F>("hclusMap1", "clus - lay1", 260, -26., 26., 350, -3.5, 3.5);
0365 hclusMap2 = fs->make<TH2F>("hclusMap2", "clus - lay2", 260, -26., 26., 350, -3.5, 3.5);
0366 hclusMap3 = fs->make<TH2F>("hclusMap3", "clus - lay3", 260, -26., 26., 350, -3.5, 3.5);
0367
0368 #endif
0369 }
0370
0371 void TestPixTracks::endJob() {
0372 cout << " End PixelTracksTest " << endl;
0373
0374 if (countEvents > 0.) {
0375 countTracks /= countEvents;
0376 countGoodTracks /= countEvents;
0377 countTracksInPix /= countEvents;
0378 countPVs /= countEvents;
0379 countLumi /= 1000.;
0380 cout << " Average tracks/event " << countTracks << " good " << countGoodTracks << " in pix " << countTracksInPix
0381 << " PVs " << countPVs << " events " << countEvents << " lumi pb-1 " << countLumi << "/10, bug!" << endl;
0382 }
0383 }
0384
0385
0386 void TestPixTracks::analyze(const edm::Event &e, const edm::EventSetup &es) {
0387 using namespace edm;
0388 using namespace reco;
0389 static int lumiBlockOld = -9999;
0390
0391 const float CLU_SIZE_PT_CUT = 1.;
0392
0393 int trackNumber = 0;
0394 int countNiceTracks = 0;
0395 int countPixTracks = 0;
0396
0397 int numberOfDetUnits = 0;
0398 int numberOfClusters = 0;
0399 int numberOfPixels = 0;
0400
0401 int numberOfDetUnits1 = 0;
0402 int numOfClustersPerDet1 = 0;
0403 int numOfClustersPerLay1 = 0;
0404 int numOfPixelsPerLay1 = 0;
0405
0406 int numberOfDetUnits2 = 0;
0407 int numOfClustersPerDet2 = 0;
0408 int numOfClustersPerLay2 = 0;
0409 int numOfPixelsPerLay2 = 0;
0410
0411 int numberOfDetUnits3 = 0;
0412 int numOfClustersPerDet3 = 0;
0413 int numOfClustersPerLay3 = 0;
0414 int numOfPixelsPerLay3 = 0;
0415
0416 int run = e.id().run();
0417 int event = e.id().event();
0418 int lumiBlock = e.luminosityBlock();
0419 int bx = e.bunchCrossing();
0420 int orbit = e.orbitNumber();
0421
0422 if (PRINT)
0423 cout << "Run " << run << " Event " << event << " LS " << lumiBlock << endl;
0424
0425 hbx0->Fill(float(bx));
0426 hlumi0->Fill(float(lumiBlock));
0427
0428 edm::LuminosityBlock const &iLumi = e.getLuminosityBlock();
0429 edm::Handle<LumiSummary> lumi;
0430 iLumi.getByToken(lumiToken_, lumi);
0431 edm::Handle<edm::ConditionsInLumiBlock> cond;
0432 float intlumi = 0, instlumi = 0;
0433 int beamint1 = 0, beamint2 = 0;
0434 iLumi.getByToken(condToken_, cond);
0435
0436
0437
0438
0439 if (lumi.isValid()) {
0440 intlumi = (lumi->intgRecLumi()) / 1000.;
0441 instlumi = (lumi->avgInsDelLumi()) / 1000.;
0442 beamint1 = (cond->totalIntensityBeam1) / 1000;
0443 beamint2 = (cond->totalIntensityBeam2) / 1000;
0444 } else {
0445
0446 }
0447
0448 hinstl->Fill(float(lumiBlock), float(instlumi));
0449 hintgl->Fill(float(lumiBlock), float(intlumi));
0450 hbeam1->Fill(float(lumiBlock), float(beamint1));
0451 hbeam2->Fill(float(lumiBlock), float(beamint2));
0452
0453 #ifdef L1
0454
0455 Handle<L1GlobalTriggerReadoutRecord> L1GTRR;
0456 e.getByToken(l1gtrrToken_, L1GTRR);
0457
0458 if (L1GTRR.isValid()) {
0459
0460
0461 for (unsigned int i = 0; i < L1GTRR->decisionWord().size(); ++i) {
0462 int l1flag = L1GTRR->decisionWord()[i];
0463 int t1flag = L1GTRR->technicalTriggerWord()[i];
0464
0465 if (l1flag > 0)
0466 hl1a->Fill(float(i));
0467 if (t1flag > 0 && i < 64)
0468 hl1t->Fill(float(i));
0469 }
0470 }
0471 #endif
0472
0473 #ifdef HLT
0474
0475 bool hlt[256];
0476 for (int i = 0; i < 256; ++i)
0477 hlt[i] = false;
0478
0479 edm::TriggerNames TrigNames;
0480 edm::Handle<edm::TriggerResults> HLTResults;
0481
0482
0483 e.getByToken(hltToken_, HLTResults);
0484 if ((HLTResults.isValid() == true) && (HLTResults->size() > 0)) {
0485
0486 const edm::TriggerNames &TrigNames = e.triggerNames(*HLTResults);
0487
0488
0489
0490 for (unsigned int i = 0; i < TrigNames.triggerNames().size(); i++) {
0491
0492
0493 if ((HLTResults->wasrun(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0494 (HLTResults->accept(TrigNames.triggerIndex(TrigNames.triggerName(i))) == true) &&
0495 (HLTResults->error(TrigNames.triggerIndex(TrigNames.triggerName(i))) == false)) {
0496 hlt[i] = true;
0497 hlt1->Fill(float(i));
0498
0499 }
0500
0501 }
0502 }
0503 #endif
0504
0505
0506 edm::ESHandle<TrackerGeometry> geom = es.getHandle(trackerGeomToken_);
0507 const TrackerGeometry &theTracker(*geom);
0508
0509
0510
0511 edm::Handle<reco::VertexCollection> vertices;
0512 e.getByToken(vtxToken_, vertices);
0513
0514 if (PRINT)
0515 cout << " PV list " << vertices->size() << endl;
0516 int pvNotFake = 0, pvsTrue = 0;
0517 vector<float> pvzVector;
0518 for (reco::VertexCollection::const_iterator iv = vertices->begin(); iv != vertices->end(); ++iv) {
0519 if ((iv->isFake()) == 1)
0520 continue;
0521 pvNotFake++;
0522 float pvx = iv->x();
0523 float pvy = iv->y();
0524 float pvz = iv->z();
0525 int numTracksPerPV = iv->tracksSize();
0526
0527
0528
0529
0530
0531
0532
0533
0534 if (PRINT)
0535 cout << " PV " << pvNotFake << " pos = " << pvx << "/" << pvy << "/" << pvz << ", Num of tracks "
0536 << numTracksPerPV << endl;
0537
0538 hpvz->Fill(pvz);
0539 if (pvz > -22. && pvz < 22.) {
0540 float pvr = sqrt(pvx * pvx + pvy * pvy);
0541 hpvxy->Fill(pvx, pvy);
0542 hpvr->Fill(pvr);
0543 if (pvr < 0.3) {
0544 pvsTrue++;
0545 pvzVector.push_back(pvz);
0546
0547 }
0548 }
0549
0550
0551
0552 }
0553 hNumPv->Fill(float(pvNotFake));
0554 hNumPvClean->Fill(float(pvsTrue));
0555
0556 if (PRINT)
0557 cout << " Not fake PVs = " << pvNotFake << " good position " << pvsTrue << endl;
0558
0559 Handle<reco::TrackCollection> recTracks;
0560 e.getByToken(srcToken_, recTracks);
0561
0562 if (PRINT)
0563 cout << " Tracks " << recTracks->size() << endl;
0564 for (reco::TrackCollection::const_iterator t = recTracks->begin(); t != recTracks->end(); ++t) {
0565 trackNumber++;
0566 numOfClustersPerDet1 = 0;
0567 numOfClustersPerDet2 = 0;
0568 numOfClustersPerDet3 = 0;
0569 int pixelHits = 0;
0570
0571 int size = t->recHitsSize();
0572 float pt = t->pt();
0573 float eta = t->eta();
0574 float phi = t->phi();
0575 float trackCharge = t->charge();
0576 float d0 = t->d0();
0577 float dz = t->dz();
0578 float tkvx = t->vx();
0579 float tkvy = t->vy();
0580 float tkvz = t->vz();
0581
0582 if (PRINT)
0583 cout << "Track " << trackNumber << " Pt " << pt << " Eta " << eta << " d0/dz " << d0 << " " << dz << " Hits "
0584 << size << endl;
0585
0586 hEta->Fill(eta);
0587 hDz->Fill(dz);
0588 if (abs(eta) > 2.8 || abs(dz) > 25.)
0589 continue;
0590
0591 hD0->Fill(d0);
0592 if (d0 > 1.0)
0593 continue;
0594
0595 bool goodTrack = false;
0596 for (vector<float>::iterator m = pvzVector.begin(); m != pvzVector.end(); ++m) {
0597 float z = *m;
0598 float tmp = abs(z - dz);
0599 hzdiff->Fill(tmp);
0600 if (tmp < 1.)
0601 goodTrack = true;
0602 }
0603
0604 if (!goodTrack)
0605 continue;
0606 countNiceTracks++;
0607 hPt->Fill(pt);
0608
0609
0610 for (trackingRecHit_iterator recHit = t->recHitsBegin(); recHit != t->recHitsEnd(); ++recHit) {
0611 if (!((*recHit)->isValid()))
0612 continue;
0613
0614 if ((*recHit)->geographicalId().det() != DetId::Tracker)
0615 continue;
0616
0617 const DetId &hit_detId = (*recHit)->geographicalId();
0618 uint IntSubDetID = (hit_detId.subdetId());
0619
0620
0621 if (IntSubDetID == 0)
0622 continue;
0623 if (IntSubDetID != PixelSubdetector::PixelBarrel)
0624 continue;
0625
0626
0627 PXBDetId pdetId = PXBDetId(hit_detId);
0628
0629
0630
0631 int layer = pdetId.layer();
0632
0633 int ladder = pdetId.ladder();
0634
0635 int zindex = pdetId.module();
0636
0637 if (PRINT)
0638 cout << "barrel layer/ladder/module: " << layer << "/" << ladder << "/" << zindex << endl;
0639
0640
0641 const PixelGeomDetUnit *theGeomDet = dynamic_cast<const PixelGeomDetUnit *>(theTracker.idToDet(hit_detId));
0642 double detZ = theGeomDet->surface().position().z();
0643 double detR = theGeomDet->surface().position().perp();
0644 const PixelTopology *topol = &(theGeomDet->specificTopology());
0645
0646
0647
0648
0649 const SiPixelRecHit *hit = dynamic_cast<const SiPixelRecHit *>((*recHit).get());
0650
0651
0652
0653 if (hit) {
0654
0655
0656
0657
0658
0659
0660 edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = hit->cluster();
0661
0662 if (!clust.isNonnull())
0663 continue;
0664
0665 numberOfClusters++;
0666 pixelHits++;
0667 float charge = (clust->charge()) / 1000.0;
0668 int size = clust->size();
0669 int sizeX = clust->sizeX();
0670 int sizeY = clust->sizeY();
0671 float row = clust->x();
0672 float col = clust->y();
0673 numberOfPixels += size;
0674
0675
0676
0677 if (PRINT)
0678 cout << " cluster " << numberOfClusters << " charge = " << charge << " size = " << size << endl;
0679
0680 LocalPoint lp = topol->localPosition(MeasurementPoint(clust->x(), clust->y()));
0681
0682
0683
0684
0685 GlobalPoint clustgp = theGeomDet->surface().toGlobal(lp);
0686 double gX = clustgp.x();
0687 double gY = clustgp.y();
0688 double gZ = clustgp.z();
0689
0690
0691
0692 TVector3 v(gX, gY, gZ);
0693 float phi = v.Phi();
0694
0695
0696
0697
0698
0699
0700
0701
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715 if (layer == 1) {
0716 hDetMap1->Fill(float(zindex), float(ladder));
0717 hcluDetMap1->Fill(col, row);
0718 hcharge1->Fill(charge);
0719 hcols1->Fill(col);
0720 hrows1->Fill(row);
0721
0722 hclusMap1->Fill(gZ, phi);
0723
0724 if (pt > CLU_SIZE_PT_CUT) {
0725 hsize1->Fill(float(size));
0726 hsizex1->Fill(float(sizeX));
0727 hsizey1->Fill(float(sizeY));
0728
0729 hclumult1->Fill(eta, float(size));
0730 hclumultx1->Fill(eta, float(sizeX));
0731 hclumulty1->Fill(eta, float(sizeY));
0732 hcluchar1->Fill(eta, float(charge));
0733 }
0734
0735 numOfClustersPerDet1++;
0736 numOfClustersPerLay1++;
0737 numOfPixelsPerLay1 += size;
0738
0739 } else if (layer == 2) {
0740 hDetMap2->Fill(float(zindex), float(ladder));
0741 hcluDetMap2->Fill(col, row);
0742 hcharge2->Fill(charge);
0743 hcols2->Fill(col);
0744 hrows2->Fill(row);
0745
0746 hclusMap2->Fill(gZ, phi);
0747
0748 if (pt > CLU_SIZE_PT_CUT) {
0749 hsize2->Fill(float(size));
0750 hsizex2->Fill(float(sizeX));
0751 hsizey2->Fill(float(sizeY));
0752
0753 hclumult2->Fill(eta, float(size));
0754 hclumultx2->Fill(eta, float(sizeX));
0755 hclumulty2->Fill(eta, float(sizeY));
0756 hcluchar2->Fill(eta, float(charge));
0757 }
0758
0759 numOfClustersPerDet2++;
0760 numOfClustersPerLay2++;
0761 numOfPixelsPerLay2 += size;
0762
0763 } else if (layer == 3) {
0764 hDetMap3->Fill(float(zindex), float(ladder));
0765 hcluDetMap3->Fill(col, row);
0766 hcharge3->Fill(charge);
0767 hcols3->Fill(col);
0768 hrows3->Fill(row);
0769
0770 hclusMap3->Fill(gZ, phi);
0771
0772 if (pt > CLU_SIZE_PT_CUT) {
0773 hsize3->Fill(float(size));
0774 hsizex3->Fill(float(sizeX));
0775 hsizey3->Fill(float(sizeY));
0776 hclumult3->Fill(eta, float(size));
0777 hclumultx3->Fill(eta, float(sizeX));
0778 hclumulty3->Fill(eta, float(sizeY));
0779 hcluchar3->Fill(eta, float(charge));
0780 }
0781
0782 numOfClustersPerDet3++;
0783 numOfClustersPerLay3++;
0784 numOfPixelsPerLay3 += size;
0785 }
0786
0787 }
0788
0789 }
0790
0791 if (pixelHits > 0)
0792 countPixTracks++;
0793
0794 if (PRINT)
0795 cout << " Clusters for track " << trackNumber << " num of clusters " << numberOfClusters << " num of pixels "
0796 << pixelHits << endl;
0797
0798 #ifdef HISTOS
0799 if (numberOfClusters > 0) {
0800 hclusPerTrk1->Fill(float(numOfClustersPerDet1));
0801 if (PRINT)
0802 cout << "Lay1: number of clusters per track = " << numOfClustersPerDet1 << endl;
0803 hclusPerTrk2->Fill(float(numOfClustersPerDet2));
0804 if (PRINT)
0805 cout << "Lay2: number of clusters per track = " << numOfClustersPerDet1 << endl;
0806 hclusPerTrk3->Fill(float(numOfClustersPerDet3));
0807 if (PRINT)
0808 cout << "Lay3: number of clusters per track = " << numOfClustersPerDet1 << endl;
0809 }
0810 #endif
0811
0812 }
0813
0814 #ifdef HISTOS
0815 if (numberOfClusters > 0) {
0816 hclusPerLay1->Fill(float(numOfClustersPerLay1));
0817 hclusPerLay2->Fill(float(numOfClustersPerLay2));
0818 hclusPerLay3->Fill(float(numOfClustersPerLay3));
0819
0820
0821
0822
0823 hpixPerLay1->Fill(float(numOfPixelsPerLay1));
0824 hpixPerLay2->Fill(float(numOfPixelsPerLay2));
0825 hpixPerLay3->Fill(float(numOfPixelsPerLay3));
0826
0827 hclusBpix->Fill(float(numberOfClusters));
0828 hpixBpix->Fill(float(numberOfPixels));
0829 }
0830 htracksGood->Fill(float(countNiceTracks));
0831 htracksGoodInPix->Fill(float(countPixTracks));
0832 htracks->Fill(float(trackNumber));
0833
0834 hbx->Fill(float(bx));
0835 hlumi->Fill(float(lumiBlock));
0836
0837 htracksls->Fill(float(lumiBlock), float(countPixTracks));
0838 hpvsls->Fill(float(lumiBlock), float(pvsTrue));
0839 if (instlumi > 0.) {
0840 float tmp = float(countPixTracks) / instlumi;
0841 htrackslsn->Fill(float(lumiBlock), tmp);
0842 tmp = float(pvsTrue) / instlumi;
0843 hpvslsn->Fill(float(lumiBlock), tmp);
0844 }
0845
0846 #endif
0847
0848
0849 countTracks += float(trackNumber);
0850 countGoodTracks += float(countNiceTracks);
0851 countTracksInPix += float(countPixTracks);
0852 countPVs += float(pvsTrue);
0853 countEvents++;
0854 if (lumiBlock != lumiBlockOld) {
0855 countLumi += intlumi;
0856 lumiBlockOld = lumiBlock;
0857 }
0858
0859 if (PRINT)
0860 cout << " event with tracks = " << trackNumber << " " << countNiceTracks << endl;
0861
0862 return;
0863
0864 #ifdef USE_TRAJ
0865
0866
0867
0868
0869 edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
0870 e.getByToken(trackAssocToken_, trajTrackCollectionHandle);
0871
0872 TrajectoryStateCombiner tsoscomb;
0873
0874 int NbrTracks = trajTrackCollectionHandle->size();
0875 std::cout << " track measurements " << trajTrackCollectionHandle->size() << std::endl;
0876
0877 int trackNumber = 0;
0878 int numberOfClusters = 0;
0879
0880 for (TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(),
0881 itEnd = trajTrackCollectionHandle->end();
0882 it != itEnd;
0883 ++it) {
0884 int pixelHits = 0;
0885 int stripHits = 0;
0886 const Track &track = *it->val;
0887 const Trajectory &traj = *it->key;
0888
0889 std::vector<TrajectoryMeasurement> checkColl = traj.measurements();
0890 for (std::vector<TrajectoryMeasurement>::const_iterator checkTraj = checkColl.begin(); checkTraj != checkColl.end();
0891 ++checkTraj) {
0892 if (!checkTraj->updatedState().isValid())
0893 continue;
0894 TransientTrackingRecHit::ConstRecHitPointer testhit = checkTraj->recHit();
0895 if (!testhit->isValid() || testhit->geographicalId().det() != DetId::Tracker)
0896 continue;
0897 uint testSubDetID = (testhit->geographicalId().subdetId());
0898 if (testSubDetID == PixelSubdetector::PixelBarrel || testSubDetID == PixelSubdetector::PixelEndcap)
0899 pixelHits++;
0900 else if (testSubDetID == StripSubdetector::TIB || testSubDetID == StripSubdetector::TOB ||
0901 testSubDetID == StripSubdetector::TID || testSubDetID == StripSubdetector::TEC)
0902 stripHits++;
0903 }
0904
0905 if (pixelHits == 0)
0906 continue;
0907
0908 trackNumber++;
0909 std::cout << " track " << trackNumber << " has pixelhits " << pixelHits << std::endl;
0910 pixelHits = 0;
0911
0912
0913 for (std::vector<TrajectoryMeasurement>::const_iterator itTraj = checkColl.begin(); itTraj != checkColl.end();
0914 ++itTraj) {
0915 if (!itTraj->updatedState().isValid())
0916 continue;
0917
0918 TrajectoryStateOnSurface tsos = tsoscomb(itTraj->forwardPredictedState(), itTraj->backwardPredictedState());
0919 TransientTrackingRecHit::ConstRecHitPointer hit = itTraj->recHit();
0920 if (!hit->isValid() || hit->geographicalId().det() != DetId::Tracker)
0921 continue;
0922
0923 const DetId &hit_detId = hit->geographicalId();
0924 uint IntSubDetID = (hit_detId.subdetId());
0925
0926 if (IntSubDetID == 0)
0927 continue;
0928 if (IntSubDetID != PixelSubdetector::PixelBarrel)
0929 continue;
0930
0931
0932
0933
0934
0935
0936
0937
0938
0939
0940 const TrackingRecHit *persistentHit = hit->hit();
0941
0942 if ((persistentHit != 0) && (typeid(*persistentHit) == typeid(SiPixelRecHit))) {
0943
0944 const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(hit->hit());
0945
0946 edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const &clust = (*pixhit).cluster();
0947
0948 if (clust.isNonnull()) {
0949 numberOfClusters++;
0950 pixelHits++;
0951 float charge = (clust->charge()) / 1000.0;
0952 int size = clust->size();
0953 int size_x = clust->sizeX();
0954 int size_y = clust->sizeY();
0955 float row = clust->x();
0956 float col = clust->y();
0957
0958
0959
0960
0961
0962 int maxPixelCol = clust->maxPixelCol();
0963 int maxPixelRow = clust->maxPixelRow();
0964 int minPixelCol = clust->minPixelCol();
0965 int minPixelRow = clust->minPixelRow();
0966
0967
0968
0969
0970
0971
0972
0973
0974
0975
0976
0977
0978
0979
0980
0981
0982
0983
0984
0985
0986
0987 }
0988 }
0989
0990 }
0991
0992 if (PRINT)
0993 cout << " Cluster for track " << trackNumber << " cluaters " << numberOfClusters << " " << pixelHits << endl;
0994
0995 }
0996
0997 #endif
0998
0999 cout << " event with tracks = " << trackNumber << " " << countGoodTracks << endl;
1000
1001 }
1002
1003
1004 DEFINE_FWK_MODULE(TestPixTracks);