File indexing completed on 2021-11-07 23:57:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <memory>
0020 #include <vector>
0021 #include <regex>
0022 #include <cassert>
0023 #include <chrono>
0024 #include <iomanip>
0025 #include <boost/range/adaptor/indexed.hpp>
0026
0027
0028 #include "Alignment/OfflineValidation/plugins/PrimaryVertexValidation.h"
0029
0030
0031 #include "TH1F.h"
0032 #include "TH2F.h"
0033 #include "TF1.h"
0034 #include "TVector3.h"
0035 #include "TFile.h"
0036 #include "TMath.h"
0037 #include "TROOT.h"
0038 #include "TChain.h"
0039 #include "TNtuple.h"
0040 #include "TMatrixD.h"
0041 #include "TVectorD.h"
0042
0043
0044 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0045 #include "DataFormats/GeometryCommonDetAlgo/interface/Measurement1D.h"
0046 #include "DataFormats/TrackReco/interface/Track.h"
0047 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0048 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0049 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0050 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0051 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0052 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0053 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0054 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0055 #include "Geometry/TrackerGeometryBuilder/interface/PixelTopologyMap.h"
0056 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0057 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
0058 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ_vect.h"
0059 #include "RecoVertex/PrimaryVertexProducer/interface/DAClusterizerInZ.h"
0060 #include "RecoVertex/PrimaryVertexProducer/interface/GapClusterizerInZ.h"
0061 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
0062 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0063 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateClosestToPoint.h"
0064
0065 const int PrimaryVertexValidation::nMaxtracks_;
0066
0067
0068 PrimaryVertexValidation::PrimaryVertexValidation(const edm::ParameterSet& iConfig)
0069 : magFieldToken_(esConsumes<MagneticField, IdealMagneticFieldRecord>()),
0070 trackingGeomToken_(esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord>()),
0071 ttkToken_(esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", "TransientTrackBuilder"))),
0072 topoToken_(esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0073 topoTokenBR_(esConsumes<TrackerTopology, TrackerTopologyRcd, edm::Transition::BeginRun>()),
0074 geomTokenBR_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord, edm::Transition::BeginRun>()),
0075 runInfoTokenBR_(esConsumes<RunInfo, RunInfoRcd, edm::Transition::BeginRun>()),
0076 compressionSettings_(iConfig.getUntrackedParameter<int>("compressionSettings", -1)),
0077 storeNtuple_(iConfig.getParameter<bool>("storeNtuple")),
0078 lightNtupleSwitch_(iConfig.getParameter<bool>("isLightNtuple")),
0079 useTracksFromRecoVtx_(iConfig.getParameter<bool>("useTracksFromRecoVtx")),
0080 vertexZMax_(iConfig.getUntrackedParameter<double>("vertexZMax", 99.)),
0081 intLumi_(iConfig.getUntrackedParameter<double>("intLumi", 0.)),
0082 askFirstLayerHit_(iConfig.getParameter<bool>("askFirstLayerHit")),
0083 doBPix_(iConfig.getUntrackedParameter<bool>("doBPix", true)),
0084 doFPix_(iConfig.getUntrackedParameter<bool>("doFPix", true)),
0085 ptOfProbe_(iConfig.getUntrackedParameter<double>("probePt", 0.)),
0086 pOfProbe_(iConfig.getUntrackedParameter<double>("probeP", 0.)),
0087 etaOfProbe_(iConfig.getUntrackedParameter<double>("probeEta", 2.4)),
0088 nHitsOfProbe_(iConfig.getUntrackedParameter<double>("probeNHits", 0.)),
0089 nBins_(iConfig.getUntrackedParameter<int>("numberOfBins", 24)),
0090 minPt_(iConfig.getUntrackedParameter<double>("minPt", 1.)),
0091 maxPt_(iConfig.getUntrackedParameter<double>("maxPt", 20.)),
0092 debug_(iConfig.getParameter<bool>("Debug")),
0093 runControl_(iConfig.getUntrackedParameter<bool>("runControl", false)),
0094 forceBeamSpotContraint_(iConfig.getUntrackedParameter<bool>("forceBeamSpot", false)) {
0095
0096
0097
0098 usesResource(TFileService::kSharedResource);
0099
0100 std::vector<unsigned int> defaultRuns;
0101 defaultRuns.push_back(0);
0102 runControlNumbers_ = iConfig.getUntrackedParameter<std::vector<unsigned int>>("runControlNumber", defaultRuns);
0103
0104 edm::InputTag TrackCollectionTag_ = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");
0105 theTrackCollectionToken_ = consumes<reco::TrackCollection>(TrackCollectionTag_);
0106
0107 edm::InputTag VertexCollectionTag_ = iConfig.getParameter<edm::InputTag>("VertexCollectionTag");
0108 theVertexCollectionToken_ = consumes<reco::VertexCollection>(VertexCollectionTag_);
0109
0110 edm::InputTag BeamspotTag_ = iConfig.getParameter<edm::InputTag>("BeamSpotTag");
0111 theBeamspotToken_ = consumes<reco::BeamSpot>(BeamspotTag_);
0112
0113
0114 theTrackFilter_ =
0115 std::make_unique<TrackFilterForPVFinding>(iConfig.getParameter<edm::ParameterSet>("TkFilterParameters"));
0116
0117 std::string clusteringAlgorithm =
0118 iConfig.getParameter<edm::ParameterSet>("TkClusParameters").getParameter<std::string>("algorithm");
0119 if (clusteringAlgorithm == "gap") {
0120 theTrackClusterizer_ =
0121 std::make_unique<GapClusterizerInZ>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
0122 .getParameter<edm::ParameterSet>("TkGapClusParameters"));
0123 } else if (clusteringAlgorithm == "DA") {
0124 theTrackClusterizer_ =
0125 std::make_unique<DAClusterizerInZ>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
0126 .getParameter<edm::ParameterSet>("TkDAClusParameters"));
0127
0128 } else if (clusteringAlgorithm == "DA_vect") {
0129 theTrackClusterizer_ =
0130 std::make_unique<DAClusterizerInZ_vect>(iConfig.getParameter<edm::ParameterSet>("TkClusParameters")
0131 .getParameter<edm::ParameterSet>("TkDAClusParameters"));
0132 } else {
0133 throw VertexException("PrimaryVertexProducerAlgorithm: unknown clustering algorithm: " + clusteringAlgorithm);
0134 }
0135
0136 theDetails_.histobins = 500;
0137 theDetails_.setMap(PVValHelper::dxy, PVValHelper::phi, -2000., 2000.);
0138 theDetails_.setMap(PVValHelper::dxy, PVValHelper::eta, -3000., 3000.);
0139 theDetails_.setMap(PVValHelper::dxy, PVValHelper::pT, -1000., 1000.);
0140 theDetails_.setMap(PVValHelper::dxy, PVValHelper::pTCentral, -1000., 1000.);
0141 theDetails_.setMap(PVValHelper::dxy, PVValHelper::ladder, -1000., 1000.);
0142 theDetails_.setMap(PVValHelper::dxy, PVValHelper::modZ, -1000., 1000.);
0143
0144 for (int i = PVValHelper::phi; i < PVValHelper::END_OF_PLOTS; i++) {
0145 for (int j = PVValHelper::dx; j < PVValHelper::END_OF_TYPES; j++) {
0146 auto plot_index = static_cast<PVValHelper::plotVariable>(i);
0147 auto res_index = static_cast<PVValHelper::residualType>(j);
0148
0149 if (debug_) {
0150 edm::LogInfo("PrimaryVertexValidation")
0151 << "==> " << std::get<0>(PVValHelper::getTypeString(res_index)) << " " << std::setw(10)
0152 << std::get<0>(PVValHelper::getVarString(plot_index)) << std::endl;
0153 }
0154 if (res_index != PVValHelper::d3D && res_index != PVValHelper::norm_d3D)
0155 theDetails_.setMap(res_index,
0156 plot_index,
0157 theDetails_.getLow(PVValHelper::dxy, plot_index),
0158 theDetails_.getHigh(PVValHelper::dxy, plot_index));
0159 else
0160 theDetails_.setMap(res_index, plot_index, 0., theDetails_.getHigh(PVValHelper::dxy, plot_index));
0161 }
0162 }
0163
0164 edm::LogVerbatim("PrimaryVertexValidation") << "######################################";
0165 for (const auto& it : theDetails_.range) {
0166 edm::LogVerbatim("PrimaryVertexValidation")
0167 << "|" << std::setw(10) << std::get<0>(PVValHelper::getTypeString(it.first.first)) << "|" << std::setw(10)
0168 << std::get<0>(PVValHelper::getVarString(it.first.second)) << "| (" << std::setw(5) << it.second.first << ";"
0169 << std::setw(5) << it.second.second << ") |" << std::endl;
0170 }
0171
0172 theDetails_.trendbins[PVValHelper::phi] = PVValHelper::generateBins(nBins_ + 1, -180., 360.);
0173 theDetails_.trendbins[PVValHelper::eta] = PVValHelper::generateBins(nBins_ + 1, -etaOfProbe_, 2 * etaOfProbe_);
0174
0175 if (debug_) {
0176 edm::LogVerbatim("PrimaryVertexValidation") << "etaBins: ";
0177 for (auto ieta : theDetails_.trendbins[PVValHelper::eta]) {
0178 edm::LogVerbatim("PrimaryVertexValidation") << ieta << " ";
0179 }
0180 edm::LogVerbatim("PrimaryVertexValidation") << "\n";
0181
0182 edm::LogVerbatim("PrimaryVertexValidation") << "phiBins: ";
0183 for (auto iphi : theDetails_.trendbins[PVValHelper::phi]) {
0184 edm::LogVerbatim("PrimaryVertexValidation") << iphi << " ";
0185 }
0186 edm::LogVerbatim("PrimaryVertexValidation") << "\n";
0187 }
0188
0189
0190
0191 mypT_bins_ = PVValHelper::makeLogBins<float, nPtBins_>(minPt_, maxPt_);
0192
0193 std::string toOutput = "";
0194 for (auto ptbin : mypT_bins_) {
0195 toOutput += " ";
0196 toOutput += std::to_string(ptbin);
0197 toOutput += ",";
0198 }
0199
0200 edm::LogVerbatim("PrimaryVertexValidation") << "######################################\n";
0201 edm::LogVerbatim("PrimaryVertexValidation") << "The pT binning is: [" << toOutput << "] \n";
0202 }
0203
0204
0205 PrimaryVertexValidation::~PrimaryVertexValidation() = default;
0206
0207
0208
0209
0210
0211 void PrimaryVertexValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0212 using namespace std;
0213 using namespace reco;
0214 using namespace IPTools;
0215
0216 if (nBins_ != 24 && debug_) {
0217 edm::LogInfo("PrimaryVertexValidation") << "Using: " << nBins_ << " bins plots";
0218 }
0219
0220 bool passesRunControl = false;
0221
0222 if (runControl_) {
0223 for (const auto& runControlNumber : runControlNumbers_) {
0224 if (iEvent.eventAuxiliary().run() == runControlNumber) {
0225 if (debug_) {
0226 edm::LogInfo("PrimaryVertexValidation")
0227 << " run number: " << iEvent.eventAuxiliary().run() << " keeping run:" << runControlNumber;
0228 }
0229 passesRunControl = true;
0230 break;
0231 }
0232 }
0233 if (!passesRunControl)
0234 return;
0235 }
0236
0237 Nevt_++;
0238
0239
0240
0241
0242
0243 SetVarToZero();
0244
0245
0246
0247
0248
0249 const TrackerTopology* const tTopo = &iSetup.getData(topoToken_);
0250
0251
0252
0253
0254
0255 edm::ESHandle<MagneticField> theMGField = iSetup.getHandle(magFieldToken_);
0256
0257
0258
0259
0260
0261 edm::ESHandle<GlobalTrackingGeometry> theTrackingGeometry = iSetup.getHandle(trackingGeomToken_);
0262
0263
0264
0265
0266
0267 edm::ESHandle<TransientTrackBuilder> theB_ = iSetup.getHandle(ttkToken_);
0268 double fBfield_ = ((*theB_).field()->inTesla(GlobalPoint(0., 0., 0.))).z();
0269
0270
0271
0272
0273
0274 edm::Handle<TrackCollection> trackCollectionHandle = iEvent.getHandle(theTrackCollectionToken_);
0275 if (!trackCollectionHandle.isValid())
0276 return;
0277 auto const& tracks = *trackCollectionHandle;
0278
0279
0280
0281
0282
0283
0284 edm::Handle<std::vector<Vertex>> vertices;
0285
0286 try {
0287 vertices = iEvent.getHandle(theVertexCollectionToken_);
0288 } catch (cms::Exception& er) {
0289 LogTrace("PrimaryVertexValidation") << "caught std::exception " << er.what() << std::endl;
0290 }
0291
0292 std::vector<Vertex> vsorted = *(vertices);
0293
0294
0295 std::sort(vsorted.begin(), vsorted.end(), PrimaryVertexValidation::vtxSort);
0296
0297
0298
0299 if (vsorted.empty())
0300 return;
0301
0302
0303 if (std::abs(vsorted[0].z()) > vertexZMax_)
0304 return;
0305
0306 if (vsorted[0].isValid()) {
0307 xOfflineVertex_ = (vsorted)[0].x();
0308 yOfflineVertex_ = (vsorted)[0].y();
0309 zOfflineVertex_ = (vsorted)[0].z();
0310
0311 xErrOfflineVertex_ = (vsorted)[0].xError();
0312 yErrOfflineVertex_ = (vsorted)[0].yError();
0313 zErrOfflineVertex_ = (vsorted)[0].zError();
0314 }
0315
0316 h_xOfflineVertex->Fill(xOfflineVertex_);
0317 h_yOfflineVertex->Fill(yOfflineVertex_);
0318 h_zOfflineVertex->Fill(zOfflineVertex_);
0319 h_xErrOfflineVertex->Fill(xErrOfflineVertex_);
0320 h_yErrOfflineVertex->Fill(yErrOfflineVertex_);
0321 h_zErrOfflineVertex->Fill(zErrOfflineVertex_);
0322
0323 unsigned int vertexCollectionSize = vsorted.size();
0324 int nvvertex = 0;
0325
0326 for (unsigned int i = 0; i < vertexCollectionSize; i++) {
0327 const Vertex& vertex = vsorted.at(i);
0328 if (vertex.isValid())
0329 nvvertex++;
0330 }
0331
0332 nOfflineVertices_ = nvvertex;
0333 h_nOfflineVertices->Fill(nvvertex);
0334
0335 if (!vsorted.empty() && useTracksFromRecoVtx_) {
0336 double sumpt = 0;
0337 size_t ntracks = 0;
0338 double chi2ndf = 0.;
0339 double chi2prob = 0.;
0340
0341 if (!vsorted.at(0).isFake()) {
0342 Vertex pv = vsorted.at(0);
0343
0344 ntracks = pv.tracksSize();
0345 chi2ndf = pv.normalizedChi2();
0346 chi2prob = TMath::Prob(pv.chi2(), (int)pv.ndof());
0347
0348 h_recoVtxNtracks_->Fill(ntracks);
0349 h_recoVtxChi2ndf_->Fill(chi2ndf);
0350 h_recoVtxChi2Prob_->Fill(chi2prob);
0351
0352 for (Vertex::trackRef_iterator itrk = pv.tracks_begin(); itrk != pv.tracks_end(); ++itrk) {
0353 double pt = (**itrk).pt();
0354 sumpt += pt * pt;
0355
0356 const math::XYZPoint myVertex(pv.position().x(), pv.position().y(), pv.position().z());
0357
0358 double dxyRes = (**itrk).dxy(myVertex);
0359 double dzRes = (**itrk).dz(myVertex);
0360
0361 double dxy_err = (**itrk).dxyError();
0362 double dz_err = (**itrk).dzError();
0363
0364 float trackphi = ((**itrk).phi()) * (180 / M_PI);
0365 float tracketa = (**itrk).eta();
0366
0367 for (int i = 0; i < nBins_; i++) {
0368 float phiF = theDetails_.trendbins[PVValHelper::phi][i];
0369 float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
0370
0371 float etaF = theDetails_.trendbins[PVValHelper::eta][i];
0372 float etaL = theDetails_.trendbins[PVValHelper::eta][i + 1];
0373
0374 if (tracketa >= etaF && tracketa < etaL) {
0375 PVValHelper::fillByIndex(a_dxyEtaBiasResiduals, i, dxyRes * cmToum, "1");
0376 PVValHelper::fillByIndex(a_dzEtaBiasResiduals, i, dzRes * cmToum, "2");
0377 PVValHelper::fillByIndex(n_dxyEtaBiasResiduals, i, (dxyRes) / dxy_err, "3");
0378 PVValHelper::fillByIndex(n_dzEtaBiasResiduals, i, (dzRes) / dz_err, "4");
0379 }
0380
0381 if (trackphi >= phiF && trackphi < phiL) {
0382 PVValHelper::fillByIndex(a_dxyPhiBiasResiduals, i, dxyRes * cmToum, "5");
0383 PVValHelper::fillByIndex(a_dzPhiBiasResiduals, i, dzRes * cmToum, "6");
0384 PVValHelper::fillByIndex(n_dxyPhiBiasResiduals, i, (dxyRes) / dxy_err, "7");
0385 PVValHelper::fillByIndex(n_dzPhiBiasResiduals, i, (dzRes) / dz_err, "8");
0386
0387 for (int j = 0; j < nBins_; j++) {
0388 float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
0389 float etaK = theDetails_.trendbins[PVValHelper::eta][j + 1];
0390
0391 if (tracketa >= etaJ && tracketa < etaK) {
0392 a_dxyBiasResidualsMap[i][j]->Fill(dxyRes * cmToum);
0393 a_dzBiasResidualsMap[i][j]->Fill(dzRes * cmToum);
0394
0395 n_dxyBiasResidualsMap[i][j]->Fill((dxyRes) / dxy_err);
0396 n_dzBiasResidualsMap[i][j]->Fill((dzRes) / dz_err);
0397 }
0398 }
0399 }
0400 }
0401 }
0402
0403 h_recoVtxSumPt_->Fill(sumpt);
0404 }
0405 }
0406
0407
0408
0409
0410
0411 BeamSpot beamSpot;
0412 edm::Handle<BeamSpot> beamSpotHandle = iEvent.getHandle(theBeamspotToken_);
0413
0414 if (beamSpotHandle.isValid()) {
0415 beamSpot = *beamSpotHandle;
0416 BSx0_ = beamSpot.x0();
0417 BSy0_ = beamSpot.y0();
0418 BSz0_ = beamSpot.z0();
0419 Beamsigmaz_ = beamSpot.sigmaZ();
0420 Beamdxdz_ = beamSpot.dxdz();
0421 BeamWidthX_ = beamSpot.BeamWidthX();
0422 BeamWidthY_ = beamSpot.BeamWidthY();
0423
0424 wxy2_ = TMath::Power(BeamWidthX_, 2) + TMath::Power(BeamWidthY_, 2);
0425
0426 } else {
0427 edm::LogWarning("PrimaryVertexValidation") << "No BeamSpot found!";
0428 }
0429
0430 h_BSx0->Fill(BSx0_);
0431 h_BSy0->Fill(BSy0_);
0432 h_BSz0->Fill(BSz0_);
0433 h_Beamsigmaz->Fill(Beamsigmaz_);
0434 h_BeamWidthX->Fill(BeamWidthX_);
0435 h_BeamWidthY->Fill(BeamWidthY_);
0436
0437 if (debug_)
0438 edm::LogInfo("PrimaryVertexValidation") << "Beamspot x:" << BSx0_ << " y:" << BSy0_ << " z:" << BSz0_;
0439
0440
0441
0442
0443
0444 RunNumber_ = iEvent.eventAuxiliary().run();
0445 LuminosityBlockNumber_ = iEvent.eventAuxiliary().luminosityBlock();
0446 EventNumber_ = iEvent.eventAuxiliary().id().event();
0447
0448 if (debug_)
0449 edm::LogInfo("PrimaryVertexValidation") << " looping over " << trackCollectionHandle->size() << "tracks";
0450
0451 h_nTracks->Fill(trackCollectionHandle->size());
0452
0453
0454
0455
0456
0457 std::vector<TransientTrack> t_tks;
0458 for (const auto& track : tracks) {
0459 TransientTrack tt = theB_->build(&(track));
0460 tt.setBeamSpot(beamSpot);
0461 t_tks.push_back(tt);
0462 }
0463
0464 if (debug_) {
0465 edm::LogInfo("PrimaryVertexValidation") << "Found: " << t_tks.size() << " reconstructed tracks";
0466 }
0467
0468
0469
0470
0471
0472 std::vector<TransientTrack> seltks = theTrackFilter_->select(t_tks);
0473
0474
0475
0476
0477
0478 vector<vector<TransientTrack>> clusters = theTrackClusterizer_->clusterize(seltks);
0479
0480 if (debug_) {
0481 edm::LogInfo("PrimaryVertexValidation")
0482 << " looping over: " << clusters.size() << " clusters from " << t_tks.size() << " selected tracks";
0483 }
0484
0485 nClus_ = clusters.size();
0486 h_nClus->Fill(nClus_);
0487
0488
0489
0490
0491 for (const auto& iclus : clusters) {
0492 nTracksPerClus_ = 0;
0493
0494 unsigned int i = 0;
0495 for (const auto& theTTrack : iclus) {
0496 i++;
0497
0498 if (nTracks_ >= nMaxtracks_) {
0499 edm::LogError("PrimaryVertexValidation")
0500 << " Warning - Number of tracks: " << nTracks_ << " , greater than " << nMaxtracks_;
0501 continue;
0502 }
0503
0504 const Track& theTrack = theTTrack.track();
0505
0506 pt_[nTracks_] = theTrack.pt();
0507 p_[nTracks_] = theTrack.p();
0508 nhits_[nTracks_] = theTrack.numberOfValidHits();
0509 eta_[nTracks_] = theTrack.eta();
0510 theta_[nTracks_] = theTrack.theta();
0511 phi_[nTracks_] = theTrack.phi();
0512 chi2_[nTracks_] = theTrack.chi2();
0513 chi2ndof_[nTracks_] = theTrack.normalizedChi2();
0514 charge_[nTracks_] = theTrack.charge();
0515 qoverp_[nTracks_] = theTrack.qoverp();
0516 dz_[nTracks_] = theTrack.dz();
0517 dxy_[nTracks_] = theTrack.dxy();
0518
0519 TrackBase::TrackQuality _trackQuality = TrackBase::qualityByName("highPurity");
0520 isHighPurity_[nTracks_] = theTrack.quality(_trackQuality);
0521
0522 math::XYZPoint point(BSx0_, BSy0_, BSz0_);
0523 dxyBs_[nTracks_] = theTrack.dxy(point);
0524 dzBs_[nTracks_] = theTrack.dz(point);
0525
0526 xPCA_[nTracks_] = theTrack.vertex().x();
0527 yPCA_[nTracks_] = theTrack.vertex().y();
0528 zPCA_[nTracks_] = theTrack.vertex().z();
0529
0530
0531
0532
0533
0534 const reco::HitPattern& hits = theTrack.hitPattern();
0535
0536 int nRecHit1D = 0;
0537 int nRecHit2D = 0;
0538 int nhitinTIB = hits.numberOfValidStripTIBHits();
0539 int nhitinTOB = hits.numberOfValidStripTOBHits();
0540 int nhitinTID = hits.numberOfValidStripTIDHits();
0541 int nhitinTEC = hits.numberOfValidStripTECHits();
0542 int nhitinBPIX = hits.numberOfValidPixelBarrelHits();
0543 int nhitinFPIX = hits.numberOfValidPixelEndcapHits();
0544 for (trackingRecHit_iterator iHit = theTTrack.recHitsBegin(); iHit != theTTrack.recHitsEnd(); ++iHit) {
0545 if ((*iHit)->isValid()) {
0546 if (this->isHit2D(**iHit, phase_)) {
0547 ++nRecHit2D;
0548 } else {
0549 ++nRecHit1D;
0550 }
0551 }
0552 }
0553
0554 nhits1D_[nTracks_] = nRecHit1D;
0555 nhits2D_[nTracks_] = nRecHit2D;
0556 nhitsBPIX_[nTracks_] = nhitinBPIX;
0557 nhitsFPIX_[nTracks_] = nhitinFPIX;
0558 nhitsTIB_[nTracks_] = nhitinTIB;
0559 nhitsTID_[nTracks_] = nhitinTID;
0560 nhitsTOB_[nTracks_] = nhitinTOB;
0561 nhitsTEC_[nTracks_] = nhitinTEC;
0562
0563
0564
0565
0566
0567 bool pass = true;
0568 if (askFirstLayerHit_)
0569 pass = this->hasFirstLayerPixelHits(theTTrack);
0570 if (pass && (theTrack.pt() >= ptOfProbe_) && std::abs(theTrack.eta()) <= etaOfProbe_ &&
0571 (theTrack.numberOfValidHits()) >= nHitsOfProbe_ && (theTrack.p()) >= pOfProbe_) {
0572 isGoodTrack_[nTracks_] = 1;
0573 }
0574
0575
0576
0577
0578
0579 vector<TransientTrack> theFinalTracks;
0580 theFinalTracks.clear();
0581
0582 for (const auto& tk : iclus) {
0583 pass = this->hasFirstLayerPixelHits(tk);
0584 if (pass) {
0585 if (tk == theTTrack)
0586 continue;
0587 else {
0588 theFinalTracks.push_back(tk);
0589 }
0590 }
0591 }
0592
0593 if (theFinalTracks.size() > 1) {
0594 if (debug_)
0595 edm::LogInfo("PrimaryVertexValidation") << "Transient Track Collection size: " << theFinalTracks.size();
0596 try {
0597
0598 auto theFitter = std::unique_ptr<VertexFitter<5>>(new AdaptiveVertexFitter());
0599 TransientVertex theFittedVertex;
0600
0601 if (forceBeamSpotContraint_) {
0602 theFittedVertex = theFitter->vertex(theFinalTracks, beamSpot);
0603 } else {
0604 theFittedVertex = theFitter->vertex(theFinalTracks);
0605 }
0606
0607 double totalTrackWeights = 0;
0608 if (theFittedVertex.isValid()) {
0609 if (theFittedVertex.hasTrackWeight()) {
0610 for (const auto& theFinalTrack : theFinalTracks) {
0611 sumOfWeightsUnbiasedVertex_[nTracks_] += theFittedVertex.trackWeight(theFinalTrack);
0612 totalTrackWeights += theFittedVertex.trackWeight(theFinalTrack);
0613 h_fitVtxTrackWeights_->Fill(theFittedVertex.trackWeight(theFinalTrack));
0614 }
0615 }
0616
0617 h_fitVtxTrackAverageWeight_->Fill(totalTrackWeights / theFinalTracks.size());
0618
0619 const math::XYZPoint theRecoVertex(xOfflineVertex_, yOfflineVertex_, zOfflineVertex_);
0620 const math::XYZPoint myVertex(
0621 theFittedVertex.position().x(), theFittedVertex.position().y(), theFittedVertex.position().z());
0622
0623 const Vertex vertex = theFittedVertex;
0624 fillTrackHistos(hDA, "all", &theTTrack, vertex, beamSpot, fBfield_);
0625
0626 hasRecVertex_[nTracks_] = 1;
0627 xUnbiasedVertex_[nTracks_] = theFittedVertex.position().x();
0628 yUnbiasedVertex_[nTracks_] = theFittedVertex.position().y();
0629 zUnbiasedVertex_[nTracks_] = theFittedVertex.position().z();
0630
0631 chi2normUnbiasedVertex_[nTracks_] = theFittedVertex.normalisedChiSquared();
0632 chi2UnbiasedVertex_[nTracks_] = theFittedVertex.totalChiSquared();
0633 DOFUnbiasedVertex_[nTracks_] = theFittedVertex.degreesOfFreedom();
0634 chi2ProbUnbiasedVertex_[nTracks_] =
0635 TMath::Prob(theFittedVertex.totalChiSquared(), (int)theFittedVertex.degreesOfFreedom());
0636 tracksUsedForVertexing_[nTracks_] = theFinalTracks.size();
0637
0638 h_fitVtxNtracks_->Fill(theFinalTracks.size());
0639 h_fitVtxChi2_->Fill(theFittedVertex.totalChiSquared());
0640 h_fitVtxNdof_->Fill(theFittedVertex.degreesOfFreedom());
0641 h_fitVtxChi2ndf_->Fill(theFittedVertex.normalisedChiSquared());
0642 h_fitVtxChi2Prob_->Fill(
0643 TMath::Prob(theFittedVertex.totalChiSquared(), (int)theFittedVertex.degreesOfFreedom()));
0644
0645
0646 double dxyFromMyVertex = theTrack.dxy(myVertex);
0647 double dzFromMyVertex = theTrack.dz(myVertex);
0648
0649 GlobalPoint vert(
0650 theFittedVertex.position().x(), theFittedVertex.position().y(), theFittedVertex.position().z());
0651
0652
0653
0654
0655
0656 double dz_err = sqrt(std::pow(theTrack.dzError(), 2) + theFittedVertex.positionError().czz());
0657
0658
0659 std::pair<bool, Measurement1D> s_ip2dpv = signedTransverseImpactParameter(
0660 theTTrack, GlobalVector(theTrack.px(), theTrack.py(), theTrack.pz()), theFittedVertex);
0661
0662 double s_ip2dpv_corr = s_ip2dpv.second.value();
0663 double s_ip2dpv_err = s_ip2dpv.second.error();
0664
0665
0666 std::pair<bool, Measurement1D> s_ip3dpv = signedImpactParameter3D(
0667 theTTrack, GlobalVector(theTrack.px(), theTrack.py(), theTrack.pz()), theFittedVertex);
0668
0669 double s_ip3dpv_corr = s_ip3dpv.second.value();
0670 double s_ip3dpv_err = s_ip3dpv.second.error();
0671
0672
0673 std::pair<bool, Measurement1D> ip3dpv = absoluteImpactParameter3D(theTTrack, theFittedVertex);
0674 double ip3d_corr = ip3dpv.second.value();
0675 double ip3d_err = ip3dpv.second.error();
0676
0677
0678 TrajectoryStateClosestToPoint traj = (theTTrack).trajectoryStateClosestToPoint(vert);
0679
0680 GlobalPoint refPoint = traj.position();
0681 GlobalPoint cPToVtx = traj.theState().position();
0682
0683 float my_dx = refPoint.x() - myVertex.x();
0684 float my_dy = refPoint.y() - myVertex.y();
0685
0686 float my_dx2 = cPToVtx.x() - myVertex.x();
0687 float my_dy2 = cPToVtx.y() - myVertex.y();
0688
0689 float my_dxy = std::sqrt(my_dx * my_dx + my_dy * my_dy);
0690
0691 double d0 = traj.perigeeParameters().transverseImpactParameter();
0692
0693 double z0 = traj.perigeeParameters().longitudinalImpactParameter();
0694 double z0_error = traj.perigeeError().longitudinalImpactParameterError();
0695
0696 if (debug_) {
0697 edm::LogInfo("PrimaryVertexValidation")
0698 << "my_dx:" << my_dx << " my_dy:" << my_dy << " my_dxy:" << my_dxy << " my_dx2:" << my_dx2
0699 << " my_dy2:" << my_dy2 << " d0: " << d0 << " dxyFromVtx:" << dxyFromMyVertex << "\n"
0700 << " ============================== "
0701 << "\n"
0702 << "diff1:" << std::abs(d0) - std::abs(my_dxy) << "\n"
0703 << "diff2:" << std::abs(d0) - std::abs(dxyFromMyVertex) << "\n"
0704 << "diff3:" << (my_dx - my_dx2) << " " << (my_dy - my_dy2) << "\n"
0705 << std::endl;
0706 }
0707
0708
0709
0710 dxyFromMyVertex_[nTracks_] = dxyFromMyVertex;
0711 dxyErrorFromMyVertex_[nTracks_] = s_ip2dpv_err;
0712 IPTsigFromMyVertex_[nTracks_] = dxyFromMyVertex / s_ip2dpv_err;
0713
0714 dzFromMyVertex_[nTracks_] = dzFromMyVertex;
0715 dzErrorFromMyVertex_[nTracks_] = dz_err;
0716 IPLsigFromMyVertex_[nTracks_] = dzFromMyVertex / dz_err;
0717
0718 d3DFromMyVertex_[nTracks_] = ip3d_corr;
0719 d3DErrorFromMyVertex_[nTracks_] = ip3d_err;
0720 IP3DsigFromMyVertex_[nTracks_] = (ip3d_corr / ip3d_err);
0721
0722
0723
0724 float trackphi = (theTrack.phi()) * (180. / M_PI);
0725 float tracketa = theTrack.eta();
0726 float trackpt = theTrack.pt();
0727 float trackp = theTrack.p();
0728 float tracknhits = theTrack.numberOfValidHits();
0729
0730
0731
0732 int ladder_num = -1.;
0733 int module_num = -1.;
0734 int L1BPixHitCount = 0;
0735
0736 for (auto const& hit : theTrack.recHits()) {
0737 const DetId& detId = hit->geographicalId();
0738 unsigned int subid = detId.subdetId();
0739
0740 if (hit->isValid() && (subid == PixelSubdetector::PixelBarrel)) {
0741 int layer = tTopo->pxbLayer(detId);
0742 if (layer == 1) {
0743 const SiPixelRecHit* prechit = dynamic_cast<const SiPixelRecHit*>(
0744 hit);
0745 double clusterProbability = prechit->clusterProbability(0);
0746 if (clusterProbability > 0) {
0747 h_probeL1ClusterProb_->Fill(log10(clusterProbability));
0748 }
0749
0750 L1BPixHitCount += 1;
0751 ladder_num = tTopo->pxbLadder(detId);
0752 module_num = tTopo->pxbModule(detId);
0753 }
0754 }
0755 }
0756
0757 h_probeL1Ladder_->Fill(ladder_num);
0758 h_probeL1Module_->Fill(module_num);
0759 h2_probeLayer1Map_->Fill(module_num, ladder_num);
0760 h_probeHasBPixL1Overlap_->Fill(L1BPixHitCount);
0761
0762
0763 if (module_num > 0 && ladder_num > 0) {
0764 a_dxyL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dxyFromMyVertex * cmToum);
0765 a_dzL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dzFromMyVertex * cmToum);
0766 n_dxyL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dxyFromMyVertex / s_ip2dpv_err);
0767 n_dzL1ResidualsMap[ladder_num - 1][module_num - 1]->Fill(dzFromMyVertex / dz_err);
0768 }
0769
0770
0771
0772 for (int ipTBin = 0; ipTBin < nPtBins_; ipTBin++) {
0773 float pTF = mypT_bins_[ipTBin];
0774 float pTL = mypT_bins_[ipTBin + 1];
0775
0776 if (debug_)
0777 edm::LogInfo("PrimaryVertexValidation") << "ipTBin:" << ipTBin << " " << mypT_bins_[ipTBin]
0778 << " < pT < " << mypT_bins_[ipTBin + 1] << std::endl;
0779
0780 if (std::abs(tracketa) < 1.5 && (trackpt >= pTF && trackpt < pTL)) {
0781 if (debug_)
0782 edm::LogInfo("PrimaryVertexValidation") << "passes this cut: " << mypT_bins_[ipTBin] << std::endl;
0783 PVValHelper::fillByIndex(h_dxy_pT_, ipTBin, dxyFromMyVertex * cmToum, "9");
0784 PVValHelper::fillByIndex(h_dz_pT_, ipTBin, dzFromMyVertex * cmToum, "10");
0785 PVValHelper::fillByIndex(h_norm_dxy_pT_, ipTBin, dxyFromMyVertex / s_ip2dpv_err, "11");
0786 PVValHelper::fillByIndex(h_norm_dz_pT_, ipTBin, dzFromMyVertex / dz_err, "12");
0787
0788 if (std::abs(tracketa) < 1.) {
0789 if (debug_)
0790 edm::LogInfo("PrimaryVertexValidation")
0791 << "passes tight eta cut: " << mypT_bins_[ipTBin] << std::endl;
0792 PVValHelper::fillByIndex(h_dxy_Central_pT_, ipTBin, dxyFromMyVertex * cmToum, "13");
0793 PVValHelper::fillByIndex(h_dz_Central_pT_, ipTBin, dzFromMyVertex * cmToum, "14");
0794 PVValHelper::fillByIndex(h_norm_dxy_Central_pT_, ipTBin, dxyFromMyVertex / s_ip2dpv_err, "15");
0795 PVValHelper::fillByIndex(h_norm_dz_Central_pT_, ipTBin, dzFromMyVertex / dz_err, "16");
0796 }
0797 }
0798 }
0799
0800
0801 if (trackpt >= ptOfProbe_ && std::abs(tracketa) <= etaOfProbe_ && tracknhits >= nHitsOfProbe_ &&
0802 trackp >= pOfProbe_) {
0803 std::pair<bool, bool> pixelOcc = pixelHitsCheck((theTTrack));
0804
0805 if (debug_) {
0806 if (pixelOcc.first == true)
0807 edm::LogInfo("PrimaryVertexValidation") << "has BPIx hits" << std::endl;
0808 if (pixelOcc.second == true)
0809 edm::LogInfo("PrimaryVertexValidation") << "has FPix hits" << std::endl;
0810 }
0811
0812 if (!doBPix_ && (pixelOcc.first == true))
0813 continue;
0814 if (!doFPix_ && (pixelOcc.second == true))
0815 continue;
0816
0817 fillTrackHistos(hDA, "sel", &(theTTrack), vertex, beamSpot, fBfield_);
0818
0819
0820 h_probePt_->Fill(theTrack.pt());
0821 h_probePtRebin_->Fill(theTrack.pt());
0822 h_probeP_->Fill(theTrack.p());
0823 h_probeEta_->Fill(theTrack.eta());
0824 h_probePhi_->Fill(theTrack.phi());
0825 h2_probeEtaPhi_->Fill(theTrack.eta(), theTrack.phi());
0826 h2_probeEtaPt_->Fill(theTrack.eta(), theTrack.pt());
0827
0828 h_probeChi2_->Fill(theTrack.chi2());
0829 h_probeNormChi2_->Fill(theTrack.normalizedChi2());
0830 h_probeCharge_->Fill(theTrack.charge());
0831 h_probeQoverP_->Fill(theTrack.qoverp());
0832 h_probeHits_->Fill(theTrack.numberOfValidHits());
0833 h_probeHits1D_->Fill(nRecHit1D);
0834 h_probeHits2D_->Fill(nRecHit2D);
0835 h_probeHitsInTIB_->Fill(nhitinTIB);
0836 h_probeHitsInTOB_->Fill(nhitinTOB);
0837 h_probeHitsInTID_->Fill(nhitinTID);
0838 h_probeHitsInTEC_->Fill(nhitinTEC);
0839 h_probeHitsInBPIX_->Fill(nhitinBPIX);
0840 h_probeHitsInFPIX_->Fill(nhitinFPIX);
0841
0842 float dxyRecoV = theTrack.dz(theRecoVertex);
0843 float dzRecoV = theTrack.dxy(theRecoVertex);
0844 float dxysigmaRecoV =
0845 TMath::Sqrt(theTrack.d0Error() * theTrack.d0Error() + xErrOfflineVertex_ * yErrOfflineVertex_);
0846 float dzsigmaRecoV =
0847 TMath::Sqrt(theTrack.dzError() * theTrack.dzError() + zErrOfflineVertex_ * zErrOfflineVertex_);
0848
0849 double zTrack = (theTTrack.stateAtBeamLine().trackStateAtPCA()).position().z();
0850 double zVertex = theFittedVertex.position().z();
0851 double tantheta = tan((theTTrack.stateAtBeamLine().trackStateAtPCA()).momentum().theta());
0852
0853 double dz2 = pow(theTrack.dzError(), 2) + wxy2_ / pow(tantheta, 2);
0854 double restrkz = zTrack - zVertex;
0855 double pulltrkz = (zTrack - zVertex) / TMath::Sqrt(dz2);
0856
0857 h_probedxyRecoV_->Fill(dxyRecoV);
0858 h_probedzRecoV_->Fill(dzRecoV);
0859
0860 h_probedzRefitV_->Fill(dxyFromMyVertex);
0861 h_probedxyRefitV_->Fill(dzFromMyVertex);
0862
0863 h_probed0RefitV_->Fill(d0);
0864 h_probez0RefitV_->Fill(z0);
0865
0866 h_probesignIP2DRefitV_->Fill(s_ip2dpv_corr);
0867 h_probed3DRefitV_->Fill(ip3d_corr);
0868 h_probereszRefitV_->Fill(restrkz);
0869
0870 h_probeRecoVSigZ_->Fill(dzRecoV / dzsigmaRecoV);
0871 h_probeRecoVSigXY_->Fill(dxyRecoV / dxysigmaRecoV);
0872 h_probeRefitVSigZ_->Fill(dzFromMyVertex / dz_err);
0873 h_probeRefitVSigXY_->Fill(dxyFromMyVertex / s_ip2dpv_err);
0874 h_probeRefitVSig3D_->Fill(ip3d_corr / ip3d_err);
0875 h_probeRefitVLogSig3D_->Fill(log10(ip3d_corr / ip3d_err));
0876 h_probeRefitVSigResZ_->Fill(pulltrkz);
0877
0878 a_dxyVsPhi->Fill(trackphi, dxyFromMyVertex * cmToum);
0879 a_dzVsPhi->Fill(trackphi, z0 * cmToum);
0880 n_dxyVsPhi->Fill(trackphi, dxyFromMyVertex / s_ip2dpv_err);
0881 n_dzVsPhi->Fill(trackphi, z0 / z0_error);
0882
0883 a_dxyVsEta->Fill(tracketa, dxyFromMyVertex * cmToum);
0884 a_dzVsEta->Fill(tracketa, z0 * cmToum);
0885 n_dxyVsEta->Fill(tracketa, dxyFromMyVertex / s_ip2dpv_err);
0886 n_dzVsEta->Fill(tracketa, z0 / z0_error);
0887
0888 if (ladder_num > 0 && module_num > 0) {
0889 LogDebug("PrimaryVertexValidation")
0890 << " ladder_num: " << ladder_num << " module_num: " << module_num << std::endl;
0891
0892 PVValHelper::fillByIndex(h_dxy_modZ_, module_num - 1, dxyFromMyVertex * cmToum, "17");
0893 PVValHelper::fillByIndex(h_dz_modZ_, module_num - 1, dzFromMyVertex * cmToum, "18");
0894 PVValHelper::fillByIndex(h_norm_dxy_modZ_, module_num - 1, dxyFromMyVertex / s_ip2dpv_err, "19");
0895 PVValHelper::fillByIndex(h_norm_dz_modZ_, module_num - 1, dzFromMyVertex / dz_err, "20");
0896
0897 PVValHelper::fillByIndex(h_dxy_ladder_, ladder_num - 1, dxyFromMyVertex * cmToum, "21");
0898
0899 LogDebug("PrimaryVertexValidation") << "h_dxy_ladder size:" << h_dxy_ladder_.size() << std::endl;
0900
0901 if (L1BPixHitCount == 1) {
0902 PVValHelper::fillByIndex(h_dxy_ladderNoOverlap_, ladder_num - 1, dxyFromMyVertex * cmToum);
0903 } else {
0904 PVValHelper::fillByIndex(h_dxy_ladderOverlap_, ladder_num - 1, dxyFromMyVertex * cmToum);
0905 }
0906
0907 h2_probePassingLayer1Map_->Fill(module_num, ladder_num);
0908
0909 PVValHelper::fillByIndex(h_dz_ladder_, ladder_num - 1, dzFromMyVertex * cmToum, "22");
0910 PVValHelper::fillByIndex(h_norm_dxy_ladder_, ladder_num - 1, dxyFromMyVertex / s_ip2dpv_err, "23");
0911 PVValHelper::fillByIndex(h_norm_dz_ladder_, ladder_num - 1, dzFromMyVertex / dz_err, "24");
0912 }
0913
0914
0915 for (int i = 0; i < nBins_; i++) {
0916 float phiF = theDetails_.trendbins[PVValHelper::phi][i];
0917 float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
0918
0919 float etaF = theDetails_.trendbins[PVValHelper::eta][i];
0920 float etaL = theDetails_.trendbins[PVValHelper::eta][i + 1];
0921
0922 if (tracketa >= etaF && tracketa < etaL) {
0923 PVValHelper::fillByIndex(a_dxyEtaResiduals, i, dxyFromMyVertex * cmToum, "25");
0924 PVValHelper::fillByIndex(a_dxEtaResiduals, i, my_dx * cmToum, "26");
0925 PVValHelper::fillByIndex(a_dyEtaResiduals, i, my_dy * cmToum, "27");
0926 PVValHelper::fillByIndex(a_dzEtaResiduals, i, dzFromMyVertex * cmToum, "28");
0927 PVValHelper::fillByIndex(n_dxyEtaResiduals, i, dxyFromMyVertex / s_ip2dpv_err, "29");
0928 PVValHelper::fillByIndex(n_dzEtaResiduals, i, dzFromMyVertex / dz_err, "30");
0929 PVValHelper::fillByIndex(a_IP2DEtaResiduals, i, s_ip2dpv_corr * cmToum, "31");
0930 PVValHelper::fillByIndex(n_IP2DEtaResiduals, i, s_ip2dpv_corr / s_ip2dpv_err, "32");
0931 PVValHelper::fillByIndex(a_reszEtaResiduals, i, restrkz * cmToum, "33");
0932 PVValHelper::fillByIndex(n_reszEtaResiduals, i, pulltrkz, "34");
0933 PVValHelper::fillByIndex(a_d3DEtaResiduals, i, ip3d_corr * cmToum, "35");
0934 PVValHelper::fillByIndex(n_d3DEtaResiduals, i, ip3d_corr / ip3d_err, "36");
0935 PVValHelper::fillByIndex(a_IP3DEtaResiduals, i, s_ip3dpv_corr * cmToum, "37");
0936 PVValHelper::fillByIndex(n_IP3DEtaResiduals, i, s_ip3dpv_corr / s_ip3dpv_err, "38");
0937 }
0938
0939 if (trackphi >= phiF && trackphi < phiL) {
0940 PVValHelper::fillByIndex(a_dxyPhiResiduals, i, dxyFromMyVertex * cmToum, "39");
0941 PVValHelper::fillByIndex(a_dxPhiResiduals, i, my_dx * cmToum, "40");
0942 PVValHelper::fillByIndex(a_dyPhiResiduals, i, my_dy * cmToum, "41");
0943 PVValHelper::fillByIndex(a_dzPhiResiduals, i, dzFromMyVertex * cmToum, "42");
0944 PVValHelper::fillByIndex(n_dxyPhiResiduals, i, dxyFromMyVertex / s_ip2dpv_err, "43");
0945 PVValHelper::fillByIndex(n_dzPhiResiduals, i, dzFromMyVertex / dz_err, "44");
0946 PVValHelper::fillByIndex(a_IP2DPhiResiduals, i, s_ip2dpv_corr * cmToum, "45");
0947 PVValHelper::fillByIndex(n_IP2DPhiResiduals, i, s_ip2dpv_corr / s_ip2dpv_err, "46");
0948 PVValHelper::fillByIndex(a_reszPhiResiduals, i, restrkz * cmToum, "47");
0949 PVValHelper::fillByIndex(n_reszPhiResiduals, i, pulltrkz, "48");
0950 PVValHelper::fillByIndex(a_d3DPhiResiduals, i, ip3d_corr * cmToum, "49");
0951 PVValHelper::fillByIndex(n_d3DPhiResiduals, i, ip3d_corr / ip3d_err, "50");
0952 PVValHelper::fillByIndex(a_IP3DPhiResiduals, i, s_ip3dpv_corr * cmToum, "51");
0953 PVValHelper::fillByIndex(n_IP3DPhiResiduals, i, s_ip3dpv_corr / s_ip3dpv_err, "52");
0954
0955 for (int j = 0; j < nBins_; j++) {
0956 float etaJ = theDetails_.trendbins[PVValHelper::eta][j];
0957 float etaK = theDetails_.trendbins[PVValHelper::eta][j + 1];
0958
0959 if (tracketa >= etaJ && tracketa < etaK) {
0960 a_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex * cmToum);
0961 a_dzResidualsMap[i][j]->Fill(dzFromMyVertex * cmToum);
0962 n_dxyResidualsMap[i][j]->Fill(dxyFromMyVertex / s_ip2dpv_err);
0963 n_dzResidualsMap[i][j]->Fill(dzFromMyVertex / dz_err);
0964 a_d3DResidualsMap[i][j]->Fill(ip3d_corr * cmToum);
0965 n_d3DResidualsMap[i][j]->Fill(ip3d_corr / ip3d_err);
0966 }
0967 }
0968 }
0969 }
0970 }
0971
0972 if (debug_) {
0973 edm::LogInfo("PrimaryVertexValidation")
0974 << " myVertex.x()= " << myVertex.x() << "\n"
0975 << " myVertex.y()= " << myVertex.y() << " \n"
0976 << " myVertex.z()= " << myVertex.z() << " \n"
0977 << " theTrack.dz(myVertex)= " << theTrack.dz(myVertex) << " \n"
0978 << " zPCA -myVertex.z() = " << (theTrack.vertex().z() - myVertex.z());
0979
0980 }
0981 }
0982
0983
0984
0985 } catch (cms::Exception& er) {
0986 LogTrace("PrimaryVertexValidation") << "caught std::exception " << er.what() << std::endl;
0987 }
0988
0989 }
0990
0991 else {
0992 if (debug_)
0993 edm::LogInfo("PrimaryVertexValidation") << "Not enough tracks to make a vertex. Returns no vertex info";
0994 }
0995
0996 ++nTracks_;
0997 ++nTracksPerClus_;
0998
0999 if (debug_)
1000 edm::LogInfo("PrimaryVertexValidation") << "Track " << i << " : pT = " << theTrack.pt();
1001
1002 }
1003 }
1004
1005
1006
1007 if (storeNtuple_) {
1008 rootTree_->Fill();
1009 }
1010 }
1011
1012
1013 bool PrimaryVertexValidation::isHit2D(const TrackingRecHit& hit, const PVValHelper::detectorPhase& thePhase) const {
1014 if (hit.dimension() < 2) {
1015 return false;
1016 } else {
1017 const DetId detId(hit.geographicalId());
1018 if (detId.det() == DetId::Tracker) {
1019 if (detId.subdetId() == PixelSubdetector::PixelBarrel || detId.subdetId() == PixelSubdetector::PixelEndcap) {
1020 return true;
1021 } else if (thePhase != PVValHelper::phase2) {
1022 if (dynamic_cast<const SiStripRecHit2D*>(&hit))
1023 return false;
1024 else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
1025 return true;
1026 else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit))
1027 return false;
1028 else {
1029 edm::LogError("UnknownType") << "@SUB=PrimaryVertexValidation::isHit2D"
1030 << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
1031 << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
1032 return false;
1033 }
1034 } else {
1035 return false;
1036 }
1037 } else {
1038 edm::LogWarning("DetectorMismatch") << "@SUB=PrimaryVertexValidation::isHit2D"
1039 << "Hit not in tracker with 'official' dimension >=2.";
1040 return true;
1041 }
1042 }
1043
1044 }
1045
1046
1047 std::pair<bool, bool> PrimaryVertexValidation::pixelHitsCheck(const reco::TransientTrack& track) {
1048 bool hasBPixHits = false;
1049 bool hasFPixHits = false;
1050
1051 const reco::HitPattern& p = track.hitPattern();
1052 if (p.numberOfValidPixelEndcapHits() != 0) {
1053 hasFPixHits = true;
1054 }
1055 if (p.numberOfValidPixelBarrelHits() != 0) {
1056 hasBPixHits = true;
1057 }
1058
1059 return std::make_pair(hasBPixHits, hasFPixHits);
1060 }
1061
1062
1063 bool PrimaryVertexValidation::hasFirstLayerPixelHits(const reco::TransientTrack& track) {
1064 using namespace reco;
1065 const HitPattern& p = track.hitPattern();
1066 for (int i = 0; i < p.numberOfAllHits(HitPattern::TRACK_HITS); i++) {
1067 uint32_t pattern = p.getHitPattern(HitPattern::TRACK_HITS, i);
1068 if (p.pixelBarrelHitFilter(pattern) || p.pixelEndcapHitFilter(pattern)) {
1069 if (p.getLayer(pattern) == 1) {
1070 if (p.validHitFilter(pattern)) {
1071 return true;
1072 }
1073 }
1074 }
1075 }
1076 return false;
1077 }
1078
1079
1080 void PrimaryVertexValidation::beginJob() {
1081 edm::LogInfo("PrimaryVertexValidation") << "######################################\n"
1082 << "Begin Job \n"
1083 << "######################################";
1084
1085
1086 Nevt_ = 0;
1087 if (compressionSettings_ > 0) {
1088 fs->file().SetCompressionSettings(compressionSettings_);
1089 }
1090
1091
1092 rootTree_ = fs->make<TTree>("tree", "PV Validation tree");
1093
1094
1095
1096 if (lightNtupleSwitch_) {
1097 rootTree_->Branch("EventNumber", &EventNumber_, "EventNumber/i");
1098 rootTree_->Branch("RunNumber", &RunNumber_, "RunNumber/i");
1099 rootTree_->Branch("LuminosityBlockNumber", &LuminosityBlockNumber_, "LuminosityBlockNumber/i");
1100 rootTree_->Branch("nOfflineVertices", &nOfflineVertices_, "nOfflineVertices/I");
1101 rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
1102 rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
1103 rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
1104 rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
1105 rootTree_->Branch("dxyFromMyVertex", &dxyFromMyVertex_, "dxyFromMyVertex[nTracks]/D");
1106 rootTree_->Branch("dzFromMyVertex", &dzFromMyVertex_, "dzFromMyVertex[nTracks]/D");
1107 rootTree_->Branch("d3DFromMyVertex", &d3DFromMyVertex_, "d3DFromMyVertex[nTracks]/D");
1108 rootTree_->Branch("IPTsigFromMyVertex", &IPTsigFromMyVertex_, "IPTsigFromMyVertex_[nTracks]/D");
1109 rootTree_->Branch("IPLsigFromMyVertex", &IPLsigFromMyVertex_, "IPLsigFromMyVertex_[nTracks]/D");
1110 rootTree_->Branch("IP3DsigFromMyVertex", &IP3DsigFromMyVertex_, "IP3DsigFromMyVertex_[nTracks]/D");
1111 rootTree_->Branch("hasRecVertex", &hasRecVertex_, "hasRecVertex[nTracks]/I");
1112 rootTree_->Branch("isGoodTrack", &isGoodTrack_, "isGoodTrack[nTracks]/I");
1113 rootTree_->Branch("isHighPurity", &isHighPurity_, "isHighPurity_[nTracks]/I");
1114
1115 } else {
1116 rootTree_->Branch("nTracks", &nTracks_, "nTracks/I");
1117 rootTree_->Branch("nTracksPerClus", &nTracksPerClus_, "nTracksPerClus/I");
1118 rootTree_->Branch("nClus", &nClus_, "nClus/I");
1119 rootTree_->Branch("xOfflineVertex", &xOfflineVertex_, "xOfflineVertex/D");
1120 rootTree_->Branch("yOfflineVertex", &yOfflineVertex_, "yOfflineVertex/D");
1121 rootTree_->Branch("zOfflineVertex", &zOfflineVertex_, "zOfflineVertex/D");
1122 rootTree_->Branch("BSx0", &BSx0_, "BSx0/D");
1123 rootTree_->Branch("BSy0", &BSy0_, "BSy0/D");
1124 rootTree_->Branch("BSz0", &BSz0_, "BSz0/D");
1125 rootTree_->Branch("Beamsigmaz", &Beamsigmaz_, "Beamsigmaz/D");
1126 rootTree_->Branch("Beamdxdz", &Beamdxdz_, "Beamdxdz/D");
1127 rootTree_->Branch("BeamWidthX", &BeamWidthX_, "BeamWidthX/D");
1128 rootTree_->Branch("BeamWidthY", &BeamWidthY_, "BeamWidthY/D");
1129 rootTree_->Branch("pt", &pt_, "pt[nTracks]/D");
1130 rootTree_->Branch("p", &p_, "p[nTracks]/D");
1131 rootTree_->Branch("nhits", &nhits_, "nhits[nTracks]/I");
1132 rootTree_->Branch("nhits1D", &nhits1D_, "nhits1D[nTracks]/I");
1133 rootTree_->Branch("nhits2D", &nhits2D_, "nhits2D[nTracks]/I");
1134 rootTree_->Branch("nhitsBPIX", &nhitsBPIX_, "nhitsBPIX[nTracks]/I");
1135 rootTree_->Branch("nhitsFPIX", &nhitsFPIX_, "nhitsFPIX[nTracks]/I");
1136 rootTree_->Branch("nhitsTIB", &nhitsTIB_, "nhitsTIB[nTracks]/I");
1137 rootTree_->Branch("nhitsTID", &nhitsTID_, "nhitsTID[nTracks]/I");
1138 rootTree_->Branch("nhitsTOB", &nhitsTOB_, "nhitsTOB[nTracks]/I");
1139 rootTree_->Branch("nhitsTEC", &nhitsTEC_, "nhitsTEC[nTracks]/I");
1140 rootTree_->Branch("eta", &eta_, "eta[nTracks]/D");
1141 rootTree_->Branch("theta", &theta_, "theta[nTracks]/D");
1142 rootTree_->Branch("phi", &phi_, "phi[nTracks]/D");
1143 rootTree_->Branch("chi2", &chi2_, "chi2[nTracks]/D");
1144 rootTree_->Branch("chi2ndof", &chi2ndof_, "chi2ndof[nTracks]/D");
1145 rootTree_->Branch("charge", &charge_, "charge[nTracks]/I");
1146 rootTree_->Branch("qoverp", &qoverp_, "qoverp[nTracks]/D");
1147 rootTree_->Branch("dz", &dz_, "dz[nTracks]/D");
1148 rootTree_->Branch("dxy", &dxy_, "dxy[nTracks]/D");
1149 rootTree_->Branch("dzBs", &dzBs_, "dzBs[nTracks]/D");
1150 rootTree_->Branch("dxyBs", &dxyBs_, "dxyBs[nTracks]/D");
1151 rootTree_->Branch("xPCA", &xPCA_, "xPCA[nTracks]/D");
1152 rootTree_->Branch("yPCA", &yPCA_, "yPCA[nTracks]/D");
1153 rootTree_->Branch("zPCA", &zPCA_, "zPCA[nTracks]/D");
1154 rootTree_->Branch("xUnbiasedVertex", &xUnbiasedVertex_, "xUnbiasedVertex[nTracks]/D");
1155 rootTree_->Branch("yUnbiasedVertex", &yUnbiasedVertex_, "yUnbiasedVertex[nTracks]/D");
1156 rootTree_->Branch("zUnbiasedVertex", &zUnbiasedVertex_, "zUnbiasedVertex[nTracks]/D");
1157 rootTree_->Branch("chi2normUnbiasedVertex", &chi2normUnbiasedVertex_, "chi2normUnbiasedVertex[nTracks]/F");
1158 rootTree_->Branch("chi2UnbiasedVertex", &chi2UnbiasedVertex_, "chi2UnbiasedVertex[nTracks]/F");
1159 rootTree_->Branch("DOFUnbiasedVertex", &DOFUnbiasedVertex_, " DOFUnbiasedVertex[nTracks]/F");
1160 rootTree_->Branch("chi2ProbUnbiasedVertex", &chi2ProbUnbiasedVertex_, "chi2ProbUnbiasedVertex[nTracks]/F");
1161 rootTree_->Branch(
1162 "sumOfWeightsUnbiasedVertex", &sumOfWeightsUnbiasedVertex_, "sumOfWeightsUnbiasedVertex[nTracks]/F");
1163 rootTree_->Branch("tracksUsedForVertexing", &tracksUsedForVertexing_, "tracksUsedForVertexing[nTracks]/I");
1164 rootTree_->Branch("dxyFromMyVertex", &dxyFromMyVertex_, "dxyFromMyVertex[nTracks]/D");
1165 rootTree_->Branch("dzFromMyVertex", &dzFromMyVertex_, "dzFromMyVertex[nTracks]/D");
1166 rootTree_->Branch("dxyErrorFromMyVertex", &dxyErrorFromMyVertex_, "dxyErrorFromMyVertex_[nTracks]/D");
1167 rootTree_->Branch("dzErrorFromMyVertex", &dzErrorFromMyVertex_, "dzErrorFromMyVertex_[nTracks]/D");
1168 rootTree_->Branch("IPTsigFromMyVertex", &IPTsigFromMyVertex_, "IPTsigFromMyVertex_[nTracks]/D");
1169 rootTree_->Branch("IPLsigFromMyVertex", &IPLsigFromMyVertex_, "IPLsigFromMyVertex_[nTracks]/D");
1170 rootTree_->Branch("hasRecVertex", &hasRecVertex_, "hasRecVertex[nTracks]/I");
1171 rootTree_->Branch("isGoodTrack", &isGoodTrack_, "isGoodTrack[nTracks]/I");
1172 }
1173
1174
1175 TFileDirectory EventFeatures = fs->mkdir("EventFeatures");
1176
1177 TH1F::SetDefaultSumw2(kTRUE);
1178
1179 h_lumiFromConfig =
1180 EventFeatures.make<TH1F>("h_lumiFromConfig", "luminosity from config;;luminosity of present run", 1, -0.5, 0.5);
1181 h_lumiFromConfig->SetBinContent(1, intLumi_);
1182
1183 h_runFromConfig = EventFeatures.make<TH1I>("h_runFromConfig",
1184 "run number from config;;run number (from configuration)",
1185 runControlNumbers_.size(),
1186 0.,
1187 runControlNumbers_.size());
1188
1189 for (const auto& run : runControlNumbers_ | boost::adaptors::indexed(1)) {
1190 h_runFromConfig->SetBinContent(run.index(), run.value());
1191 }
1192
1193 h_runFromEvent =
1194 EventFeatures.make<TH1I>("h_runFromEvent", "run number from event;;run number (from event)", 1, -0.5, 0.5);
1195 h_nTracks =
1196 EventFeatures.make<TH1F>("h_nTracks", "number of tracks per event;n_{tracks}/event;n_{events}", 300, -0.5, 299.5);
1197 h_nClus =
1198 EventFeatures.make<TH1F>("h_nClus", "number of track clusters;n_{clusters}/event;n_{events}", 50, -0.5, 49.5);
1199 h_nOfflineVertices = EventFeatures.make<TH1F>(
1200 "h_nOfflineVertices", "number of offline reconstructed vertices;n_{vertices}/event;n_{events}", 50, -0.5, 49.5);
1201 h_runNumber = EventFeatures.make<TH1F>("h_runNumber", "run number;run number;n_{events}", 100000, 250000., 350000.);
1202 h_xOfflineVertex = EventFeatures.make<TH1F>(
1203 "h_xOfflineVertex", "x-coordinate of offline vertex;x_{vertex};n_{events}", 100, -0.1, 0.1);
1204 h_yOfflineVertex = EventFeatures.make<TH1F>(
1205 "h_yOfflineVertex", "y-coordinate of offline vertex;y_{vertex};n_{events}", 100, -0.1, 0.1);
1206 h_zOfflineVertex = EventFeatures.make<TH1F>(
1207 "h_zOfflineVertex", "z-coordinate of offline vertex;z_{vertex};n_{events}", 100, -30., 30.);
1208 h_xErrOfflineVertex = EventFeatures.make<TH1F>(
1209 "h_xErrOfflineVertex", "x-coordinate error of offline vertex;err_{x}^{vtx};n_{events}", 100, 0., 0.01);
1210 h_yErrOfflineVertex = EventFeatures.make<TH1F>(
1211 "h_yErrOfflineVertex", "y-coordinate error of offline vertex;err_{y}^{vtx};n_{events}", 100, 0., 0.01);
1212 h_zErrOfflineVertex = EventFeatures.make<TH1F>(
1213 "h_zErrOfflineVertex", "z-coordinate error of offline vertex;err_{z}^{vtx};n_{events}", 100, 0., 10.);
1214 h_BSx0 = EventFeatures.make<TH1F>("h_BSx0", "x-coordinate of reco beamspot;x^{BS}_{0};n_{events}", 100, -0.1, 0.1);
1215 h_BSy0 = EventFeatures.make<TH1F>("h_BSy0", "y-coordinate of reco beamspot;y^{BS}_{0};n_{events}", 100, -0.1, 0.1);
1216 h_BSz0 = EventFeatures.make<TH1F>("h_BSz0", "z-coordinate of reco beamspot;z^{BS}_{0};n_{events}", 100, -1., 1.);
1217 h_Beamsigmaz =
1218 EventFeatures.make<TH1F>("h_Beamsigmaz", "z-coordinate beam width;#sigma_{Z}^{beam};n_{events}", 100, 0., 1.);
1219 h_BeamWidthX =
1220 EventFeatures.make<TH1F>("h_BeamWidthX", "x-coordinate beam width;#sigma_{X}^{beam};n_{events}", 100, 0., 0.01);
1221 h_BeamWidthY =
1222 EventFeatures.make<TH1F>("h_BeamWidthY", "y-coordinate beam width;#sigma_{Y}^{beam};n_{events}", 100, 0., 0.01);
1223
1224 h_etaMax = EventFeatures.make<TH1F>("etaMax", "etaMax", 1, -0.5, 0.5);
1225 h_pTinfo = EventFeatures.make<TH1F>("pTinfo", "pTinfo", 3, -1.5, 1.5);
1226 h_pTinfo->GetXaxis()->SetBinLabel(1, "n. bins");
1227 h_pTinfo->GetXaxis()->SetBinLabel(2, "pT min");
1228 h_pTinfo->GetXaxis()->SetBinLabel(3, "pT max");
1229
1230 h_nbins = EventFeatures.make<TH1F>("nbins", "nbins", 1, -0.5, 0.5);
1231 h_nLadders = EventFeatures.make<TH1F>("nladders", "n. ladders", 1, -0.5, 0.5);
1232 h_nModZ = EventFeatures.make<TH1F>("nModZ", "n. modules along z", 1, -0.5, 0.5);
1233
1234
1235 TFileDirectory ProbeFeatures = fs->mkdir("ProbeTrackFeatures");
1236
1237 h_probePt_ = ProbeFeatures.make<TH1F>("h_probePt", "p_{T} of probe track;track p_{T} (GeV); tracks", 100, 0., 50.);
1238 h_probePtRebin_ = ProbeFeatures.make<TH1F>(
1239 "h_probePtRebin", "p_{T} of probe track;track p_{T} (GeV); tracks", mypT_bins_.size() - 1, mypT_bins_.data());
1240 h_probeP_ = ProbeFeatures.make<TH1F>("h_probeP", "momentum of probe track;track p (GeV); tracks", 100, 0., 100.);
1241 h_probeEta_ = ProbeFeatures.make<TH1F>(
1242 "h_probeEta", "#eta of the probe track;track #eta;tracks", 54, -etaOfProbe_, etaOfProbe_);
1243 h_probePhi_ = ProbeFeatures.make<TH1F>("h_probePhi", "#phi of probe track;track #phi (rad);tracks", 100, -3.15, 3.15);
1244
1245 h2_probeEtaPhi_ =
1246 ProbeFeatures.make<TH2F>("h2_probeEtaPhi",
1247 "probe track #phi vs #eta;#eta of probe track;track #phi of probe track (rad); tracks",
1248 54,
1249 -etaOfProbe_,
1250 etaOfProbe_,
1251 100,
1252 -M_PI,
1253 M_PI);
1254 h2_probeEtaPt_ = ProbeFeatures.make<TH2F>("h2_probeEtaPt",
1255 "probe track p_{T} vs #eta;#eta of probe track;track p_{T} (GeV); tracks",
1256 54,
1257 -etaOfProbe_,
1258 etaOfProbe_,
1259 100,
1260 0.,
1261 50.);
1262
1263 h_probeChi2_ =
1264 ProbeFeatures.make<TH1F>("h_probeChi2", "#chi^{2} of probe track;track #chi^{2}; tracks", 100, 0., 100.);
1265 h_probeNormChi2_ = ProbeFeatures.make<TH1F>(
1266 "h_probeNormChi2", " normalized #chi^{2} of probe track;track #chi^{2}/ndof; tracks", 100, 0., 10.);
1267 h_probeCharge_ =
1268 ProbeFeatures.make<TH1F>("h_probeCharge", "charge of probe track;track charge Q;tracks", 3, -1.5, 1.5);
1269 h_probeQoverP_ =
1270 ProbeFeatures.make<TH1F>("h_probeQoverP", "q/p of probe track; track Q/p (GeV^{-1});tracks", 200, -1., 1.);
1271 h_probedzRecoV_ = ProbeFeatures.make<TH1F>(
1272 "h_probedzRecoV", "d_{z}(V_{offline}) of probe track;track d_{z}(V_{off}) (cm);tracks", 200, -1., 1.);
1273 h_probedxyRecoV_ = ProbeFeatures.make<TH1F>(
1274 "h_probedxyRecoV", "d_{xy}(V_{offline}) of probe track;track d_{xy}(V_{off}) (cm);tracks", 200, -1., 1.);
1275 h_probedzRefitV_ = ProbeFeatures.make<TH1F>(
1276 "h_probedzRefitV", "d_{z}(V_{refit}) of probe track;track d_{z}(V_{fit}) (cm);tracks", 200, -0.5, 0.5);
1277 h_probesignIP2DRefitV_ = ProbeFeatures.make<TH1F>(
1278 "h_probesignIPRefitV", "ip_{2D}(V_{refit}) of probe track;track ip_{2D}(V_{fit}) (cm);tracks", 200, -1., 1.);
1279 h_probedxyRefitV_ = ProbeFeatures.make<TH1F>(
1280 "h_probedxyRefitV", "d_{xy}(V_{refit}) of probe track;track d_{xy}(V_{fit}) (cm);tracks", 200, -0.5, 0.5);
1281
1282 h_probez0RefitV_ = ProbeFeatures.make<TH1F>(
1283 "h_probez0RefitV", "z_{0}(V_{refit}) of probe track;track z_{0}(V_{fit}) (cm);tracks", 200, -1., 1.);
1284 h_probed0RefitV_ = ProbeFeatures.make<TH1F>(
1285 "h_probed0RefitV", "d_{0}(V_{refit}) of probe track;track d_{0}(V_{fit}) (cm);tracks", 200, -1., 1.);
1286
1287 h_probed3DRefitV_ = ProbeFeatures.make<TH1F>(
1288 "h_probed3DRefitV", "d_{3D}(V_{refit}) of probe track;track d_{3D}(V_{fit}) (cm);tracks", 200, 0., 1.);
1289 h_probereszRefitV_ = ProbeFeatures.make<TH1F>(
1290 "h_probeReszRefitV", "z_{track} -z_{V_{refit}};track res_{z}(V_{refit}) (cm);tracks", 200, -1., 1.);
1291
1292 h_probeRecoVSigZ_ = ProbeFeatures.make<TH1F>(
1293 "h_probeRecoVSigZ", "Longitudinal DCA Significance (reco);d_{z}(V_{off})/#sigma_{dz};tracks", 100, -8, 8);
1294 h_probeRecoVSigXY_ = ProbeFeatures.make<TH1F>(
1295 "h_probeRecoVSigXY", "Transverse DCA Significance (reco);d_{xy}(V_{off})/#sigma_{dxy};tracks", 100, -8, 8);
1296 h_probeRefitVSigZ_ = ProbeFeatures.make<TH1F>(
1297 "h_probeRefitVSigZ", "Longitudinal DCA Significance (refit);d_{z}(V_{fit})/#sigma_{dz};tracks", 100, -8, 8);
1298 h_probeRefitVSigXY_ = ProbeFeatures.make<TH1F>(
1299 "h_probeRefitVSigXY", "Transverse DCA Significance (refit);d_{xy}(V_{fit})/#sigma_{dxy};tracks", 100, -8, 8);
1300 h_probeRefitVSig3D_ = ProbeFeatures.make<TH1F>(
1301 "h_probeRefitVSig3D", "3D DCA Significance (refit);d_{3D}/#sigma_{3D};tracks", 100, 0., 20.);
1302 h_probeRefitVLogSig3D_ =
1303 ProbeFeatures.make<TH1F>("h_probeRefitVLogSig3D",
1304 "log_{10}(3D DCA-Significance) (refit);log_{10}(d_{3D}/#sigma_{3D});tracks",
1305 100,
1306 -5.,
1307 4.);
1308 h_probeRefitVSigResZ_ = ProbeFeatures.make<TH1F>(
1309 "h_probeRefitVSigResZ",
1310 "Longitudinal residual significance (refit);(z_{track} -z_{V_{fit}})/#sigma_{res_{z}};tracks",
1311 100,
1312 -8,
1313 8);
1314
1315 h_probeHits_ = ProbeFeatures.make<TH1F>("h_probeNRechits", "N_{hits} ;N_{hits} ;tracks", 40, -0.5, 39.5);
1316 h_probeHits1D_ = ProbeFeatures.make<TH1F>("h_probeNRechits1D", "N_{hits} 1D ;N_{hits} 1D ;tracks", 40, -0.5, 39.5);
1317 h_probeHits2D_ = ProbeFeatures.make<TH1F>("h_probeNRechits2D", "N_{hits} 2D ;N_{hits} 2D ;tracks", 40, -0.5, 39.5);
1318 h_probeHitsInTIB_ =
1319 ProbeFeatures.make<TH1F>("h_probeNRechitsTIB", "N_{hits} TIB ;N_{hits} TIB;tracks", 40, -0.5, 39.5);
1320 h_probeHitsInTOB_ =
1321 ProbeFeatures.make<TH1F>("h_probeNRechitsTOB", "N_{hits} TOB ;N_{hits} TOB;tracks", 40, -0.5, 39.5);
1322 h_probeHitsInTID_ =
1323 ProbeFeatures.make<TH1F>("h_probeNRechitsTID", "N_{hits} TID ;N_{hits} TID;tracks", 40, -0.5, 39.5);
1324 h_probeHitsInTEC_ =
1325 ProbeFeatures.make<TH1F>("h_probeNRechitsTEC", "N_{hits} TEC ;N_{hits} TEC;tracks", 40, -0.5, 39.5);
1326 h_probeHitsInBPIX_ =
1327 ProbeFeatures.make<TH1F>("h_probeNRechitsBPIX", "N_{hits} BPIX;N_{hits} BPIX;tracks", 40, -0.5, 39.5);
1328 h_probeHitsInFPIX_ =
1329 ProbeFeatures.make<TH1F>("h_probeNRechitsFPIX", "N_{hits} FPIX;N_{hits} FPIX;tracks", 40, -0.5, 39.5);
1330
1331 h_probeL1Ladder_ = ProbeFeatures.make<TH1F>(
1332 "h_probeL1Ladder", "Ladder number (L1 hit); ladder number", nLadders_ + 2, -1.5, nLadders_ + 0.5);
1333 h_probeL1Module_ = ProbeFeatures.make<TH1F>(
1334 "h_probeL1Module", "Module number (L1 hit); module number", nModZ_ + 2, -1.5, nModZ_ + 0.5);
1335
1336 h2_probeLayer1Map_ = ProbeFeatures.make<TH2F>("h2_probeLayer1Map",
1337 "Position in Layer 1 of first hit;module number;ladder number",
1338 nModZ_,
1339 0.5,
1340 nModZ_ + 0.5,
1341 nLadders_,
1342 0.5,
1343 nLadders_ + 0.5);
1344
1345 h2_probePassingLayer1Map_ = ProbeFeatures.make<TH2F>("h2_probePassingLayer1Map",
1346 "Position in Layer 1 of first hit;module number;ladder number",
1347 nModZ_,
1348 0.5,
1349 nModZ_ + 0.5,
1350 nLadders_,
1351 0.5,
1352 nLadders_ + 0.5);
1353 h_probeHasBPixL1Overlap_ =
1354 ProbeFeatures.make<TH1I>("h_probeHasBPixL1Overlap", "n. hits in L1;n. L1-BPix hits;tracks", 5, -0.5, 4.5);
1355 h_probeL1ClusterProb_ = ProbeFeatures.make<TH1F>(
1356 "h_probeL1ClusterProb",
1357 "log_{10}(Cluster Probability) for Layer1 hits;log_{10}(cluster probability); n. Layer1 hits",
1358 100,
1359 -10.,
1360 0.);
1361
1362
1363 TFileDirectory RefitVertexFeatures = fs->mkdir("RefitVertexFeatures");
1364 h_fitVtxNtracks_ = RefitVertexFeatures.make<TH1F>(
1365 "h_fitVtxNtracks", "N_{trks} used in vertex fit;N^{fit}_{tracks};vertices", 100, -0.5, 99.5);
1366 h_fitVtxNdof_ = RefitVertexFeatures.make<TH1F>(
1367 "h_fitVtxNdof", "N_{DOF} of vertex fit;N_{DOF} of refit vertex;vertices", 100, -0.5, 99.5);
1368 h_fitVtxChi2_ = RefitVertexFeatures.make<TH1F>(
1369 "h_fitVtxChi2", "#chi^{2} of vertex fit;vertex #chi^{2};vertices", 100, -0.5, 99.5);
1370 h_fitVtxChi2ndf_ = RefitVertexFeatures.make<TH1F>(
1371 "h_fitVtxChi2ndf", "#chi^{2}/ndf of vertex fit;vertex #chi^{2}/ndf;vertices", 100, -0.5, 9.5);
1372 h_fitVtxChi2Prob_ = RefitVertexFeatures.make<TH1F>(
1373 "h_fitVtxChi2Prob", "Prob(#chi^{2},ndf) of vertex fit;Prob(#chi^{2},ndf);vertices", 40, 0., 1.);
1374 h_fitVtxTrackWeights_ = RefitVertexFeatures.make<TH1F>(
1375 "h_fitVtxTrackWeights", "track weights associated to track;track weights;tracks", 40, 0., 1.);
1376 h_fitVtxTrackAverageWeight_ = RefitVertexFeatures.make<TH1F>(
1377 "h_fitVtxTrackAverageWeight_", "average track weight per vertex;#LT track weight #GT;vertices", 40, 0., 1.);
1378
1379 if (useTracksFromRecoVtx_) {
1380 TFileDirectory RecoVertexFeatures = fs->mkdir("RecoVertexFeatures");
1381 h_recoVtxNtracks_ =
1382 RecoVertexFeatures.make<TH1F>("h_recoVtxNtracks", "N^{vtx}_{trks};N^{vtx}_{trks};vertices", 100, -0.5, 99.5);
1383 h_recoVtxChi2ndf_ =
1384 RecoVertexFeatures.make<TH1F>("h_recoVtxChi2ndf", "#chi^{2}/ndf vtx;#chi^{2}/ndf vtx;vertices", 10, -0.5, 9.5);
1385 h_recoVtxChi2Prob_ = RecoVertexFeatures.make<TH1F>(
1386 "h_recoVtxChi2Prob", "Prob(#chi^{2},ndf);Prob(#chi^{2},ndf);vertices", 40, 0., 1.);
1387 h_recoVtxSumPt_ =
1388 RecoVertexFeatures.make<TH1F>("h_recoVtxSumPt", "Sum(p^{trks}_{T});Sum(p^{trks}_{T});vertices", 100, 0., 200.);
1389 }
1390
1391 TFileDirectory DA = fs->mkdir("DA");
1392
1393 hDA = bookVertexHistograms(DA);
1394
1395
1396
1397
1398
1399
1400
1401 const float dxymax_phi = theDetails_.getHigh(PVValHelper::dxy, PVValHelper::phi);
1402 const float dzmax_phi = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1403 const float dxymax_eta = theDetails_.getHigh(PVValHelper::dxy, PVValHelper::phi);
1404 const float dzmax_eta = theDetails_.getHigh(PVValHelper::dz, PVValHelper::eta);
1405
1406 const float d3Dmax_eta = theDetails_.getHigh(PVValHelper::d3D, PVValHelper::eta);
1407
1408
1409
1410
1411
1412
1413
1414
1415
1416
1417
1418
1419
1420
1421 TFileDirectory AbsTransPhiRes = fs->mkdir("Abs_Transv_Phi_Residuals");
1422 a_dxyPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::dxy, PVValHelper::phi);
1423 a_dxPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::dx, PVValHelper::phi);
1424 a_dyPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::dy, PVValHelper::phi);
1425 a_IP2DPhiResiduals = bookResidualsHistogram(AbsTransPhiRes, nBins_, PVValHelper::IP2D, PVValHelper::phi);
1426
1427 TFileDirectory AbsTransEtaRes = fs->mkdir("Abs_Transv_Eta_Residuals");
1428 a_dxyEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::dxy, PVValHelper::eta);
1429 a_dxEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::dx, PVValHelper::eta);
1430 a_dyEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::dy, PVValHelper::eta);
1431 a_IP2DEtaResiduals = bookResidualsHistogram(AbsTransEtaRes, nBins_, PVValHelper::IP2D, PVValHelper::eta);
1432
1433 TFileDirectory AbsLongPhiRes = fs->mkdir("Abs_Long_Phi_Residuals");
1434 a_dzPhiResiduals = bookResidualsHistogram(AbsLongPhiRes, nBins_, PVValHelper::dz, PVValHelper::phi);
1435 a_reszPhiResiduals = bookResidualsHistogram(AbsLongPhiRes, nBins_, PVValHelper::resz, PVValHelper::phi);
1436
1437 TFileDirectory AbsLongEtaRes = fs->mkdir("Abs_Long_Eta_Residuals");
1438 a_dzEtaResiduals = bookResidualsHistogram(AbsLongEtaRes, nBins_, PVValHelper::dz, PVValHelper::eta);
1439 a_reszEtaResiduals = bookResidualsHistogram(AbsLongEtaRes, nBins_, PVValHelper::resz, PVValHelper::eta);
1440
1441 TFileDirectory Abs3DPhiRes = fs->mkdir("Abs_3D_Phi_Residuals");
1442 a_d3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes, nBins_, PVValHelper::d3D, PVValHelper::phi);
1443 a_IP3DPhiResiduals = bookResidualsHistogram(Abs3DPhiRes, nBins_, PVValHelper::IP3D, PVValHelper::phi);
1444
1445 TFileDirectory Abs3DEtaRes = fs->mkdir("Abs_3D_Eta_Residuals");
1446 a_d3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes, nBins_, PVValHelper::d3D, PVValHelper::eta);
1447 a_IP3DEtaResiduals = bookResidualsHistogram(Abs3DEtaRes, nBins_, PVValHelper::IP3D, PVValHelper::eta);
1448
1449 TFileDirectory NormTransPhiRes = fs->mkdir("Norm_Transv_Phi_Residuals");
1450 n_dxyPhiResiduals = bookResidualsHistogram(NormTransPhiRes, nBins_, PVValHelper::norm_dxy, PVValHelper::phi, true);
1451 n_IP2DPhiResiduals = bookResidualsHistogram(NormTransPhiRes, nBins_, PVValHelper::norm_IP2D, PVValHelper::phi, true);
1452
1453 TFileDirectory NormTransEtaRes = fs->mkdir("Norm_Transv_Eta_Residuals");
1454 n_dxyEtaResiduals = bookResidualsHistogram(NormTransEtaRes, nBins_, PVValHelper::norm_dxy, PVValHelper::eta, true);
1455 n_IP2DEtaResiduals = bookResidualsHistogram(NormTransEtaRes, nBins_, PVValHelper::norm_IP2D, PVValHelper::eta, true);
1456
1457 TFileDirectory NormLongPhiRes = fs->mkdir("Norm_Long_Phi_Residuals");
1458 n_dzPhiResiduals = bookResidualsHistogram(NormLongPhiRes, nBins_, PVValHelper::norm_dz, PVValHelper::phi, true);
1459 n_reszPhiResiduals = bookResidualsHistogram(NormLongPhiRes, nBins_, PVValHelper::norm_resz, PVValHelper::phi, true);
1460
1461 TFileDirectory NormLongEtaRes = fs->mkdir("Norm_Long_Eta_Residuals");
1462 n_dzEtaResiduals = bookResidualsHistogram(NormLongEtaRes, nBins_, PVValHelper::norm_dz, PVValHelper::eta, true);
1463 n_reszEtaResiduals = bookResidualsHistogram(NormLongEtaRes, nBins_, PVValHelper::norm_resz, PVValHelper::eta, true);
1464
1465 TFileDirectory Norm3DPhiRes = fs->mkdir("Norm_3D_Phi_Residuals");
1466 n_d3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes, nBins_, PVValHelper::norm_d3D, PVValHelper::phi, true);
1467 n_IP3DPhiResiduals = bookResidualsHistogram(Norm3DPhiRes, nBins_, PVValHelper::norm_IP3D, PVValHelper::phi, true);
1468
1469 TFileDirectory Norm3DEtaRes = fs->mkdir("Norm_3D_Eta_Residuals");
1470 n_d3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes, nBins_, PVValHelper::norm_d3D, PVValHelper::eta, true);
1471 n_IP3DEtaResiduals = bookResidualsHistogram(Norm3DEtaRes, nBins_, PVValHelper::norm_IP3D, PVValHelper::eta, true);
1472
1473 TFileDirectory AbsDoubleDiffRes = fs->mkdir("Abs_DoubleDiffResiduals");
1474 TFileDirectory NormDoubleDiffRes = fs->mkdir("Norm_DoubleDiffResiduals");
1475
1476 TFileDirectory AbsL1Map = fs->mkdir("Abs_L1Residuals");
1477 TFileDirectory NormL1Map = fs->mkdir("Norm_L1Residuals");
1478
1479
1480
1481 TFileDirectory AbsTranspTRes = fs->mkdir("Abs_Transv_pT_Residuals");
1482 h_dxy_pT_ = bookResidualsHistogram(AbsTranspTRes, nPtBins_, PVValHelper::dxy, PVValHelper::pT);
1483
1484 TFileDirectory AbsLongpTRes = fs->mkdir("Abs_Long_pT_Residuals");
1485 h_dz_pT_ = bookResidualsHistogram(AbsLongpTRes, nPtBins_, PVValHelper::dz, PVValHelper::pT);
1486
1487 TFileDirectory NormTranspTRes = fs->mkdir("Norm_Transv_pT_Residuals");
1488 h_norm_dxy_pT_ = bookResidualsHistogram(NormTranspTRes, nPtBins_, PVValHelper::norm_dxy, PVValHelper::pT, true);
1489
1490 TFileDirectory NormLongpTRes = fs->mkdir("Norm_Long_pT_Residuals");
1491 h_norm_dz_pT_ = bookResidualsHistogram(NormLongpTRes, nPtBins_, PVValHelper::norm_dz, PVValHelper::pT, true);
1492
1493
1494
1495 TFileDirectory AbsTranspTCentralRes = fs->mkdir("Abs_Transv_pTCentral_Residuals");
1496 h_dxy_Central_pT_ = bookResidualsHistogram(AbsTranspTCentralRes, nPtBins_, PVValHelper::dxy, PVValHelper::pTCentral);
1497
1498 TFileDirectory AbsLongpTCentralRes = fs->mkdir("Abs_Long_pTCentral_Residuals");
1499 h_dz_Central_pT_ = bookResidualsHistogram(AbsLongpTCentralRes, nPtBins_, PVValHelper::dz, PVValHelper::pTCentral);
1500
1501 TFileDirectory NormTranspTCentralRes = fs->mkdir("Norm_Transv_pTCentral_Residuals");
1502 h_norm_dxy_Central_pT_ =
1503 bookResidualsHistogram(NormTranspTCentralRes, nPtBins_, PVValHelper::norm_dxy, PVValHelper::pTCentral, true);
1504
1505 TFileDirectory NormLongpTCentralRes = fs->mkdir("Norm_Long_pTCentral_Residuals");
1506 h_norm_dz_Central_pT_ =
1507 bookResidualsHistogram(NormLongpTCentralRes, nPtBins_, PVValHelper::norm_dz, PVValHelper::pTCentral, true);
1508
1509
1510
1511 TFileDirectory AbsTransModZRes = fs->mkdir("Abs_Transv_modZ_Residuals");
1512 h_dxy_modZ_ = bookResidualsHistogram(AbsTransModZRes, nModZ_, PVValHelper::dxy, PVValHelper::modZ);
1513
1514 TFileDirectory AbsLongModZRes = fs->mkdir("Abs_Long_modZ_Residuals");
1515 h_dz_modZ_ = bookResidualsHistogram(AbsLongModZRes, nModZ_, PVValHelper::dz, PVValHelper::modZ);
1516
1517
1518
1519
1520
1521
1522
1523 TFileDirectory NormTransModZRes = fs->mkdir("Norm_Transv_modZ_Residuals");
1524 h_norm_dxy_modZ_ = bookResidualsHistogram(NormTransModZRes, nModZ_, PVValHelper::norm_dxy, PVValHelper::modZ, true);
1525
1526 TFileDirectory NormLongModZRes = fs->mkdir("Norm_Long_modZ_Residuals");
1527 h_norm_dz_modZ_ = bookResidualsHistogram(NormLongModZRes, nModZ_, PVValHelper::norm_dz, PVValHelper::modZ, true);
1528
1529 TFileDirectory AbsTransLadderRes = fs->mkdir("Abs_Transv_ladder_Residuals");
1530 h_dxy_ladder_ = bookResidualsHistogram(AbsTransLadderRes, nLadders_, PVValHelper::dxy, PVValHelper::ladder);
1531
1532 TFileDirectory AbsTransLadderResOverlap = fs->mkdir("Abs_Transv_ladderOverlap_Residuals");
1533 h_dxy_ladderOverlap_ =
1534 bookResidualsHistogram(AbsTransLadderResOverlap, nLadders_, PVValHelper::dxy, PVValHelper::ladder);
1535
1536 TFileDirectory AbsTransLadderResNoOverlap = fs->mkdir("Abs_Transv_ladderNoOverlap_Residuals");
1537 h_dxy_ladderNoOverlap_ =
1538 bookResidualsHistogram(AbsTransLadderResNoOverlap, nLadders_, PVValHelper::dxy, PVValHelper::ladder);
1539
1540 TFileDirectory AbsLongLadderRes = fs->mkdir("Abs_Long_ladder_Residuals");
1541 h_dz_ladder_ = bookResidualsHistogram(AbsLongLadderRes, nLadders_, PVValHelper::dz, PVValHelper::ladder);
1542
1543 TFileDirectory NormTransLadderRes = fs->mkdir("Norm_Transv_ladder_Residuals");
1544 h_norm_dxy_ladder_ =
1545 bookResidualsHistogram(NormTransLadderRes, nLadders_, PVValHelper::norm_dxy, PVValHelper::ladder, true);
1546
1547 TFileDirectory NormLongLadderRes = fs->mkdir("Norm_Long_ladder_Residuals");
1548 h_norm_dz_ladder_ =
1549 bookResidualsHistogram(NormLongLadderRes, nLadders_, PVValHelper::norm_dz, PVValHelper::ladder, true);
1550
1551
1552
1553 for (unsigned int iLadder = 0; iLadder < nLadders_; iLadder++) {
1554 for (unsigned int iModule = 0; iModule < nModZ_; iModule++) {
1555 a_dxyL1ResidualsMap[iLadder][iModule] =
1556 AbsL1Map.make<TH1F>(Form("histo_dxy_ladder%i_module%i", iLadder, iModule),
1557 Form("d_{xy} ladder=%i module=%i;d_{xy} [#mum];tracks", iLadder, iModule),
1558 theDetails_.histobins,
1559 -dzmax_eta,
1560 dzmax_eta);
1561
1562 a_dzL1ResidualsMap[iLadder][iModule] =
1563 AbsL1Map.make<TH1F>(Form("histo_dz_ladder%i_module%i", iLadder, iModule),
1564 Form("d_{z} ladder=%i module=%i;d_{z} [#mum];tracks", iLadder, iModule),
1565 theDetails_.histobins,
1566 -dzmax_eta,
1567 dzmax_eta);
1568
1569 n_dxyL1ResidualsMap[iLadder][iModule] =
1570 NormL1Map.make<TH1F>(Form("histo_norm_dxy_ladder%i_module%i", iLadder, iModule),
1571 Form("d_{xy} ladder=%i module=%i;d_{xy}/#sigma_{d_{xy}};tracks", iLadder, iModule),
1572 theDetails_.histobins,
1573 -dzmax_eta / 100,
1574 dzmax_eta / 100);
1575
1576 n_dzL1ResidualsMap[iLadder][iModule] =
1577 NormL1Map.make<TH1F>(Form("histo_norm_dz_ladder%i_module%i", iLadder, iModule),
1578 Form("d_{z} ladder=%i module=%i;d_{z}/#sigma_{d_{z}};tracks", iLadder, iModule),
1579 theDetails_.histobins,
1580 -dzmax_eta / 100,
1581 dzmax_eta / 100);
1582 }
1583 }
1584
1585
1586
1587 for (int i = 0; i < nBins_; ++i) {
1588 float phiF = theDetails_.trendbins[PVValHelper::phi][i];
1589 float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
1590
1591
1592
1593
1594
1595
1596 for (int j = 0; j < nBins_; ++j) {
1597 float etaF = theDetails_.trendbins[PVValHelper::eta][j];
1598 float etaL = theDetails_.trendbins[PVValHelper::eta][j + 1];
1599
1600 a_dxyResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1601 Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
1602 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy};tracks", etaF, etaL, phiF, phiL),
1603 theDetails_.histobins,
1604 -dzmax_eta,
1605 dzmax_eta);
1606
1607 a_dzResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1608 Form("histo_dz_eta_plot%i_phi_plot%i", i, j),
1609 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z};tracks", etaF, etaL, phiF, phiL),
1610 theDetails_.histobins,
1611 -dzmax_eta,
1612 dzmax_eta);
1613
1614 a_d3DResidualsMap[i][j] = AbsDoubleDiffRes.make<TH1F>(
1615 Form("histo_d3D_eta_plot%i_phi_plot%i", i, j),
1616 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D};tracks", etaF, etaL, phiF, phiL),
1617 theDetails_.histobins,
1618 0.,
1619 d3Dmax_eta);
1620
1621 n_dxyResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1622 Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
1623 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",
1624 etaF,
1625 etaL,
1626 phiF,
1627 phiL),
1628 theDetails_.histobins,
1629 -dzmax_eta / 100,
1630 dzmax_eta / 100);
1631
1632 n_dzResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1633 Form("histo_norm_dz_eta_plot%i_phi_plot%i", i, j),
1634 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",
1635 etaF,
1636 etaL,
1637 phiF,
1638 phiL),
1639 theDetails_.histobins,
1640 -dzmax_eta / 100,
1641 dzmax_eta / 100);
1642
1643 n_d3DResidualsMap[i][j] = NormDoubleDiffRes.make<TH1F>(
1644 Form("histo_norm_d3D_eta_plot%i_phi_plot%i", i, j),
1645 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{3D}/#sigma_{d_{3D}};tracks",
1646 etaF,
1647 etaL,
1648 phiF,
1649 phiL),
1650 theDetails_.histobins,
1651 0.,
1652 d3Dmax_eta);
1653 }
1654 }
1655
1656
1657
1658 TFileDirectory BiasVsParameter = fs->mkdir("BiasVsParameter");
1659
1660 a_dxyVsPhi = BiasVsParameter.make<TH2F>("h2_dxy_vs_phi",
1661 "d_{xy} vs track #phi;track #phi [rad];track d_{xy}(PV) [#mum]",
1662 nBins_,
1663 -M_PI,
1664 M_PI,
1665 theDetails_.histobins,
1666 -dxymax_phi,
1667 dxymax_phi);
1668
1669 a_dzVsPhi = BiasVsParameter.make<TH2F>("h2_dz_vs_phi",
1670 "d_{z} vs track #phi;track #phi [rad];track d_{z}(PV) [#mum]",
1671 nBins_,
1672 -M_PI,
1673 M_PI,
1674 theDetails_.histobins,
1675 -dzmax_phi,
1676 dzmax_phi);
1677
1678 n_dxyVsPhi = BiasVsParameter.make<TH2F>(
1679 "h2_n_dxy_vs_phi",
1680 "d_{xy}/#sigma_{d_{xy}} vs track #phi;track #phi [rad];track d_{xy}(PV)/#sigma_{d_{xy}}",
1681 nBins_,
1682 -M_PI,
1683 M_PI,
1684 theDetails_.histobins,
1685 -dxymax_phi / 100.,
1686 dxymax_phi / 100.);
1687
1688 n_dzVsPhi =
1689 BiasVsParameter.make<TH2F>("h2_n_dz_vs_phi",
1690 "d_{z}/#sigma_{d_{z}} vs track #phi;track #phi [rad];track d_{z}(PV)/#sigma_{d_{z}}",
1691 nBins_,
1692 -M_PI,
1693 M_PI,
1694 theDetails_.histobins,
1695 -dzmax_phi / 100.,
1696 dzmax_phi / 100.);
1697
1698 a_dxyVsEta = BiasVsParameter.make<TH2F>("h2_dxy_vs_eta",
1699 "d_{xy} vs track #eta;track #eta;track d_{xy}(PV) [#mum]",
1700 nBins_,
1701 -etaOfProbe_,
1702 etaOfProbe_,
1703 theDetails_.histobins,
1704 -dxymax_eta,
1705 dzmax_eta);
1706
1707 a_dzVsEta = BiasVsParameter.make<TH2F>("h2_dz_vs_eta",
1708 "d_{z} vs track #eta;track #eta;track d_{z}(PV) [#mum]",
1709 nBins_,
1710 -etaOfProbe_,
1711 etaOfProbe_,
1712 theDetails_.histobins,
1713 -dzmax_eta,
1714 dzmax_eta);
1715
1716 n_dxyVsEta =
1717 BiasVsParameter.make<TH2F>("h2_n_dxy_vs_eta",
1718 "d_{xy}/#sigma_{d_{xy}} vs track #eta;track #eta;track d_{xy}(PV)/#sigma_{d_{xy}}",
1719 nBins_,
1720 -etaOfProbe_,
1721 etaOfProbe_,
1722 theDetails_.histobins,
1723 -dxymax_eta / 100.,
1724 dxymax_eta / 100.);
1725
1726 n_dzVsEta = BiasVsParameter.make<TH2F>("h2_n_dz_vs_eta",
1727 "d_{z}/#sigma_{d_{z}} vs track #eta;track #eta;track d_{z}(PV)/#sigma_{d_{z}}",
1728 nBins_,
1729 -etaOfProbe_,
1730 etaOfProbe_,
1731 theDetails_.histobins,
1732 -dzmax_eta / 100.,
1733 dzmax_eta / 100.);
1734
1735 MeanTrendsDir = fs->mkdir("MeanTrends");
1736 WidthTrendsDir = fs->mkdir("WidthTrends");
1737 MedianTrendsDir = fs->mkdir("MedianTrends");
1738 MADTrendsDir = fs->mkdir("MADTrends");
1739
1740 Mean2DMapsDir = fs->mkdir("MeanMaps");
1741 Width2DMapsDir = fs->mkdir("WidthMaps");
1742
1743 double highedge = nBins_ - 0.5;
1744 double lowedge = -0.5;
1745
1746
1747
1748 a_dxyPhiMeanTrend =
1749 MeanTrendsDir.make<TH1F>("means_dxy_phi",
1750 "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
1751 nBins_,
1752 lowedge,
1753 highedge);
1754
1755 a_dxyPhiWidthTrend =
1756 WidthTrendsDir.make<TH1F>("widths_dxy_phi",
1757 "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
1758 nBins_,
1759 lowedge,
1760 highedge);
1761
1762 a_dzPhiMeanTrend =
1763 MeanTrendsDir.make<TH1F>("means_dz_phi",
1764 "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
1765 nBins_,
1766 lowedge,
1767 highedge);
1768
1769 a_dzPhiWidthTrend =
1770 WidthTrendsDir.make<TH1F>("widths_dz_phi",
1771 "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
1772 nBins_,
1773 lowedge,
1774 highedge);
1775
1776 a_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F>(
1777 "means_dxy_eta", "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]", nBins_, lowedge, highedge);
1778
1779 a_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_eta",
1780 "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
1781 nBins_,
1782 lowedge,
1783 highedge);
1784
1785 a_dzEtaMeanTrend = MeanTrendsDir.make<TH1F>(
1786 "means_dz_eta", "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]", nBins_, lowedge, highedge);
1787
1788 a_dzEtaWidthTrend = WidthTrendsDir.make<TH1F>(
1789 "widths_dz_eta", "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]", nBins_, lowedge, highedge);
1790
1791 n_dxyPhiMeanTrend = MeanTrendsDir.make<TH1F>(
1792 "norm_means_dxy_phi",
1793 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1794 nBins_,
1795 lowedge,
1796 highedge);
1797
1798 n_dxyPhiWidthTrend = WidthTrendsDir.make<TH1F>(
1799 "norm_widths_dxy_phi",
1800 "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
1801 nBins_,
1802 lowedge,
1803 highedge);
1804
1805 n_dzPhiMeanTrend = MeanTrendsDir.make<TH1F>(
1806 "norm_means_dz_phi",
1807 "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
1808 nBins_,
1809 lowedge,
1810 highedge);
1811
1812 n_dzPhiWidthTrend = WidthTrendsDir.make<TH1F>(
1813 "norm_widths_dz_phi",
1814 "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
1815 nBins_,
1816 lowedge,
1817 highedge);
1818
1819 n_dxyEtaMeanTrend = MeanTrendsDir.make<TH1F>(
1820 "norm_means_dxy_eta",
1821 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
1822 nBins_,
1823 lowedge,
1824 highedge);
1825
1826 n_dxyEtaWidthTrend = WidthTrendsDir.make<TH1F>(
1827 "norm_widths_dxy_eta",
1828 "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
1829 nBins_,
1830 lowedge,
1831 highedge);
1832
1833 n_dzEtaMeanTrend =
1834 MeanTrendsDir.make<TH1F>("norm_means_dz_eta",
1835 "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
1836 nBins_,
1837 lowedge,
1838 highedge);
1839
1840 n_dzEtaWidthTrend =
1841 WidthTrendsDir.make<TH1F>("norm_widths_dz_eta",
1842 "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
1843 nBins_,
1844 lowedge,
1845 highedge);
1846
1847
1848
1849 a_dxypTMeanTrend = MeanTrendsDir.make<TH1F>("means_dxy_pT",
1850 "#LT d_{xy} #GT vs pT;p_{T} [GeV];#LT d_{xy} #GT [#mum]",
1851 mypT_bins_.size() - 1,
1852 mypT_bins_.data());
1853
1854 a_dxypTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dxy_pT",
1855 "#sigma_{d_{xy}} vs pT;p_{T} [GeV];#sigma_{d_{xy}} [#mum]",
1856 mypT_bins_.size() - 1,
1857 mypT_bins_.data());
1858
1859 a_dzpTMeanTrend = MeanTrendsDir.make<TH1F>(
1860 "means_dz_pT", "#LT d_{z} #GT vs pT;p_{T} [GeV];#LT d_{z} #GT [#mum]", mypT_bins_.size() - 1, mypT_bins_.data());
1861
1862 a_dzpTWidthTrend = WidthTrendsDir.make<TH1F>("widths_dz_pT",
1863 "#sigma_{d_{z}} vs pT;p_{T} [GeV];#sigma_{d_{z}} [#mum]",
1864 mypT_bins_.size() - 1,
1865 mypT_bins_.data());
1866
1867 n_dxypTMeanTrend =
1868 MeanTrendsDir.make<TH1F>("norm_means_dxy_pT",
1869 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs pT;p_{T} [GeV];#LT d_{xy}/#sigma_{d_{xy}} #GT",
1870 mypT_bins_.size() - 1,
1871 mypT_bins_.data());
1872
1873 n_dxypTWidthTrend =
1874 WidthTrendsDir.make<TH1F>("norm_widths_dxy_pT",
1875 "width(d_{xy}/#sigma_{d_{xy}}) vs pT;p_{T} [GeV]; width(d_{xy}/#sigma_{d_{xy}})",
1876 mypT_bins_.size() - 1,
1877 mypT_bins_.data());
1878
1879 n_dzpTMeanTrend =
1880 MeanTrendsDir.make<TH1F>("norm_means_dz_pT",
1881 "#LT d_{z}/#sigma_{d_{z}} #GT vs pT;p_{T} [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1882 mypT_bins_.size() - 1,
1883 mypT_bins_.data());
1884
1885 n_dzpTWidthTrend =
1886 WidthTrendsDir.make<TH1F>("norm_widths_dz_pT",
1887 "width(d_{z}/#sigma_{d_{z}}) vs pT;p_{T} [GeV];width(d_{z}/#sigma_{d_{z}})",
1888 mypT_bins_.size() - 1,
1889 mypT_bins_.data());
1890
1891 a_dxypTCentralMeanTrend =
1892 MeanTrendsDir.make<TH1F>("means_dxy_pTCentral",
1893 "#LT d_{xy} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy} #GT [#mum]",
1894 mypT_bins_.size() - 1,
1895 mypT_bins_.data());
1896
1897 a_dxypTCentralWidthTrend =
1898 WidthTrendsDir.make<TH1F>("widths_dxy_pTCentral",
1899 "#sigma_{d_{xy}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{xy}} [#mum]",
1900 mypT_bins_.size() - 1,
1901 mypT_bins_.data());
1902
1903 a_dzpTCentralMeanTrend =
1904 MeanTrendsDir.make<TH1F>("means_dz_pTCentral",
1905 "#LT d_{z} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z} #GT [#mum]",
1906 mypT_bins_.size() - 1,
1907 mypT_bins_.data());
1908
1909 a_dzpTCentralWidthTrend =
1910 WidthTrendsDir.make<TH1F>("widths_dz_pTCentral",
1911 "#sigma_{d_{z}} vs p_{T};p_{T}(|#eta|<1.) [GeV];#sigma_{d_{z}} [#mum]",
1912 mypT_bins_.size() - 1,
1913 mypT_bins_.data());
1914
1915 n_dxypTCentralMeanTrend = MeanTrendsDir.make<TH1F>(
1916 "norm_means_dxy_pTCentral",
1917 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{xy}/#sigma_{d_{z}} #GT",
1918 mypT_bins_.size() - 1,
1919 mypT_bins_.data());
1920
1921 n_dxypTCentralWidthTrend = WidthTrendsDir.make<TH1F>(
1922 "norm_widths_dxy_pTCentral",
1923 "width(d_{xy}/#sigma_{d_{xy}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{xy}/#sigma_{d_{z}})",
1924 mypT_bins_.size() - 1,
1925 mypT_bins_.data());
1926
1927 n_dzpTCentralMeanTrend = MeanTrendsDir.make<TH1F>(
1928 "norm_means_dz_pTCentral",
1929 "#LT d_{z}/#sigma_{d_{z}} #GT vs p_{T};p_{T}(|#eta|<1.) [GeV];#LT d_{z}/#sigma_{d_{z}} #GT",
1930 mypT_bins_.size() - 1,
1931 mypT_bins_.data());
1932
1933 n_dzpTCentralWidthTrend = WidthTrendsDir.make<TH1F>(
1934 "norm_widths_dz_pTCentral",
1935 "width(d_{z}/#sigma_{d_{z}}) vs p_{T};p_{T}(|#eta|<1.) [GeV];width(d_{z}/#sigma_{d_{z}})",
1936 mypT_bins_.size() - 1,
1937 mypT_bins_.data());
1938
1939
1940
1941 a_dxyMeanMap = Mean2DMapsDir.make<TH2F>("means_dxy_map",
1942 "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
1943 nBins_,
1944 lowedge,
1945 highedge,
1946 nBins_,
1947 lowedge,
1948 highedge);
1949
1950 a_dzMeanMap = Mean2DMapsDir.make<TH2F>("means_dz_map",
1951 "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
1952 nBins_,
1953 lowedge,
1954 highedge,
1955 nBins_,
1956 lowedge,
1957 highedge);
1958
1959 n_dxyMeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dxy_map",
1960 "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1961 nBins_,
1962 lowedge,
1963 highedge,
1964 nBins_,
1965 lowedge,
1966 highedge);
1967
1968 n_dzMeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dz_map",
1969 "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
1970 nBins_,
1971 lowedge,
1972 highedge,
1973 nBins_,
1974 lowedge,
1975 highedge);
1976
1977 a_dxyWidthMap = Width2DMapsDir.make<TH2F>("widths_dxy_map",
1978 "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
1979 nBins_,
1980 lowedge,
1981 highedge,
1982 nBins_,
1983 lowedge,
1984 highedge);
1985
1986 a_dzWidthMap = Width2DMapsDir.make<TH2F>("widths_dz_map",
1987 "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
1988 nBins_,
1989 lowedge,
1990 highedge,
1991 nBins_,
1992 lowedge,
1993 highedge);
1994
1995 n_dxyWidthMap =
1996 Width2DMapsDir.make<TH2F>("norm_widths_dxy_map",
1997 "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
1998 nBins_,
1999 lowedge,
2000 highedge,
2001 nBins_,
2002 lowedge,
2003 highedge);
2004
2005 n_dzWidthMap = Width2DMapsDir.make<TH2F>("norm_widths_dz_map",
2006 "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
2007 nBins_,
2008 lowedge,
2009 highedge,
2010 nBins_,
2011 lowedge,
2012 highedge);
2013
2014
2015
2016 a_dxyPhiMedianTrend =
2017 MedianTrendsDir.make<TH1F>("medians_dxy_phi",
2018 "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2019 nBins_,
2020 lowedge,
2021 highedge);
2022
2023 a_dxyPhiMADTrend = MADTrendsDir.make<TH1F>(
2024 "MADs_dxy_phi",
2025 "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2026 nBins_,
2027 lowedge,
2028 highedge);
2029
2030 a_dzPhiMedianTrend =
2031 MedianTrendsDir.make<TH1F>("medians_dz_phi",
2032 "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2033 nBins_,
2034 lowedge,
2035 highedge);
2036
2037 a_dzPhiMADTrend = MADTrendsDir.make<TH1F>(
2038 "MADs_dz_phi",
2039 "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2040 nBins_,
2041 lowedge,
2042 highedge);
2043
2044 a_dxyEtaMedianTrend =
2045 MedianTrendsDir.make<TH1F>("medians_dxy_eta",
2046 "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2047 nBins_,
2048 lowedge,
2049 highedge);
2050
2051 a_dxyEtaMADTrend =
2052 MADTrendsDir.make<TH1F>("MADs_dxy_eta",
2053 "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2054 nBins_,
2055 lowedge,
2056 highedge);
2057
2058 a_dzEtaMedianTrend =
2059 MedianTrendsDir.make<TH1F>("medians_dz_eta",
2060 "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2061 nBins_,
2062 lowedge,
2063 highedge);
2064
2065 a_dzEtaMADTrend =
2066 MADTrendsDir.make<TH1F>("MADs_dz_eta",
2067 "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2068 nBins_,
2069 lowedge,
2070 highedge);
2071
2072 n_dxyPhiMedianTrend = MedianTrendsDir.make<TH1F>(
2073 "norm_medians_dxy_phi",
2074 "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2075 nBins_,
2076 lowedge,
2077 highedge);
2078
2079 n_dxyPhiMADTrend = MADTrendsDir.make<TH1F>("norm_MADs_dxy_phi",
2080 "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi "
2081 "sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2082 nBins_,
2083 lowedge,
2084 highedge);
2085
2086 n_dzPhiMedianTrend = MedianTrendsDir.make<TH1F>(
2087 "norm_medians_dz_phi",
2088 "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2089 nBins_,
2090 lowedge,
2091 highedge);
2092
2093 n_dzPhiMADTrend = MADTrendsDir.make<TH1F>("norm_MADs_dz_phi",
2094 "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi "
2095 "(sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2096 nBins_,
2097 lowedge,
2098 highedge);
2099
2100 n_dxyEtaMedianTrend = MedianTrendsDir.make<TH1F>(
2101 "norm_medians_dxy_eta",
2102 "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2103 nBins_,
2104 lowedge,
2105 highedge);
2106
2107 n_dxyEtaMADTrend = MADTrendsDir.make<TH1F>(
2108 "norm_MADs_dxy_eta",
2109 "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2110 nBins_,
2111 lowedge,
2112 highedge);
2113
2114 n_dzEtaMedianTrend = MedianTrendsDir.make<TH1F>(
2115 "norm_medians_dz_eta",
2116 "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2117 nBins_,
2118 lowedge,
2119 highedge);
2120
2121 n_dzEtaMADTrend = MADTrendsDir.make<TH1F>(
2122 "norm_MADs_dz_eta",
2123 "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2124 nBins_,
2125 lowedge,
2126 highedge);
2127
2128
2129
2130
2131
2132
2133
2134
2135 if (useTracksFromRecoVtx_) {
2136 TFileDirectory AbsTransPhiBiasRes = fs->mkdir("Abs_Transv_Phi_BiasResiduals");
2137 a_dxyPhiBiasResiduals = bookResidualsHistogram(AbsTransPhiBiasRes, nBins_, PVValHelper::dxy, PVValHelper::phi);
2138
2139 TFileDirectory AbsTransEtaBiasRes = fs->mkdir("Abs_Transv_Eta_BiasResiduals");
2140 a_dxyEtaBiasResiduals = bookResidualsHistogram(AbsTransEtaBiasRes, nBins_, PVValHelper::dxy, PVValHelper::eta);
2141
2142 TFileDirectory AbsLongPhiBiasRes = fs->mkdir("Abs_Long_Phi_BiasResiduals");
2143 a_dzPhiBiasResiduals = bookResidualsHistogram(AbsLongPhiBiasRes, nBins_, PVValHelper::dz, PVValHelper::phi);
2144
2145 TFileDirectory AbsLongEtaBiasRes = fs->mkdir("Abs_Long_Eta_BiasResiduals");
2146 a_dzEtaBiasResiduals = bookResidualsHistogram(AbsLongEtaBiasRes, nBins_, PVValHelper::dz, PVValHelper::eta);
2147
2148 TFileDirectory NormTransPhiBiasRes = fs->mkdir("Norm_Transv_Phi_BiasResiduals");
2149 n_dxyPhiBiasResiduals = bookResidualsHistogram(NormTransPhiBiasRes, nBins_, PVValHelper::dxy, PVValHelper::phi);
2150
2151 TFileDirectory NormTransEtaBiasRes = fs->mkdir("Norm_Transv_Eta_BiasResiduals");
2152 n_dxyEtaBiasResiduals = bookResidualsHistogram(NormTransEtaBiasRes, nBins_, PVValHelper::dxy, PVValHelper::eta);
2153
2154 TFileDirectory NormLongPhiBiasRes = fs->mkdir("Norm_Long_Phi_BiasResiduals");
2155 n_dzPhiBiasResiduals = bookResidualsHistogram(NormLongPhiBiasRes, nBins_, PVValHelper::dz, PVValHelper::phi);
2156
2157 TFileDirectory NormLongEtaBiasRes = fs->mkdir("Norm_Long_Eta_BiasResiduals");
2158 n_dzEtaBiasResiduals = bookResidualsHistogram(NormLongEtaBiasRes, nBins_, PVValHelper::dz, PVValHelper::eta);
2159
2160 TFileDirectory AbsDoubleDiffBiasRes = fs->mkdir("Abs_DoubleDiffBiasResiduals");
2161 TFileDirectory NormDoubleDiffBiasRes = fs->mkdir("Norm_DoubleDiffBiasResiduals");
2162
2163 for (int i = 0; i < nBins_; ++i) {
2164 float phiF = theDetails_.trendbins[PVValHelper::phi][i];
2165 float phiL = theDetails_.trendbins[PVValHelper::phi][i + 1];
2166
2167 for (int j = 0; j < nBins_; ++j) {
2168 float etaF = theDetails_.trendbins[PVValHelper::eta][j];
2169 float etaL = theDetails_.trendbins[PVValHelper::eta][j + 1];
2170
2171 a_dxyBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(
2172 Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
2173 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy} [#mum];tracks", etaF, etaL, phiF, phiL),
2174 theDetails_.histobins,
2175 -dzmax_eta,
2176 dzmax_eta);
2177
2178 a_dzBiasResidualsMap[i][j] = AbsDoubleDiffBiasRes.make<TH1F>(
2179 Form("histo_dxy_eta_plot%i_phi_plot%i", i, j),
2180 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z} [#mum];tracks", etaF, etaL, phiF, phiL),
2181 theDetails_.histobins,
2182 -dzmax_eta,
2183 dzmax_eta);
2184
2185 n_dxyBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(
2186 Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
2187 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{xy}/#sigma_{d_{xy}};tracks",
2188 etaF,
2189 etaL,
2190 phiF,
2191 phiL),
2192 theDetails_.histobins,
2193 -dzmax_eta / 100,
2194 dzmax_eta / 100);
2195
2196 n_dzBiasResidualsMap[i][j] = NormDoubleDiffBiasRes.make<TH1F>(
2197 Form("histo_norm_dxy_eta_plot%i_phi_plot%i", i, j),
2198 Form("%.2f<#eta_{tk}<%.2f %.2f#circ<#varphi_{tk}<%.2f#circ;d_{z}/#sigma_{d_{z}};tracks",
2199 etaF,
2200 etaL,
2201 phiF,
2202 phiL),
2203 theDetails_.histobins,
2204 -dzmax_eta / 100,
2205 dzmax_eta / 100);
2206 }
2207 }
2208
2209
2210
2211 TFileDirectory MeanBiasTrendsDir = fs->mkdir("MeanBiasTrends");
2212 TFileDirectory WidthBiasTrendsDir = fs->mkdir("WidthBiasTrends");
2213 TFileDirectory MedianBiasTrendsDir = fs->mkdir("MedianBiasTrends");
2214 TFileDirectory MADBiasTrendsDir = fs->mkdir("MADBiasTrends");
2215
2216 TFileDirectory Mean2DBiasMapsDir = fs->mkdir("MeanBiasMaps");
2217 TFileDirectory Width2DBiasMapsDir = fs->mkdir("WidthBiasMaps");
2218
2219
2220
2221 a_dxyPhiMeanBiasTrend =
2222 MeanBiasTrendsDir.make<TH1F>("means_dxy_phi",
2223 "#LT d_{xy} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy} #GT [#mum]",
2224 nBins_,
2225 lowedge,
2226 highedge);
2227
2228 a_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2229 "widths_dxy_phi",
2230 "#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{xy}} [#mum]",
2231 nBins_,
2232 lowedge,
2233 highedge);
2234
2235 a_dzPhiMeanBiasTrend =
2236 MeanBiasTrendsDir.make<TH1F>("means_dz_phi",
2237 "#LT d_{z} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z} #GT [#mum]",
2238 nBins_,
2239 lowedge,
2240 highedge);
2241
2242 a_dzPhiWidthBiasTrend =
2243 WidthBiasTrendsDir.make<TH1F>("widths_dz_phi",
2244 "#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#sigma_{d_{z}} [#mum]",
2245 nBins_,
2246 lowedge,
2247 highedge);
2248
2249 a_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2250 "means_dxy_eta", "#LT d_{xy} #GT vs #eta sector;#eta (sector);#LT d_{xy} #GT [#mum]", nBins_, lowedge, highedge);
2251
2252 a_dxyEtaWidthBiasTrend =
2253 WidthBiasTrendsDir.make<TH1F>("widths_dxy_eta",
2254 "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{xy}} [#mum]",
2255 nBins_,
2256 lowedge,
2257 highedge);
2258
2259 a_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2260 "means_dz_eta", "#LT d_{z} #GT vs #eta sector;#eta (sector);#LT d_{z} #GT [#mum]", nBins_, lowedge, highedge);
2261
2262 a_dzEtaWidthBiasTrend =
2263 WidthBiasTrendsDir.make<TH1F>("widths_dz_eta",
2264 "#sigma_{d_{xy}} vs #eta sector;#eta (sector);#sigma_{d_{z}} [#mum]",
2265 nBins_,
2266 lowedge,
2267 highedge);
2268
2269 n_dxyPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2270 "norm_means_dxy_phi",
2271 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{xy}/#sigma_{d_{xy}} #GT",
2272 nBins_,
2273 lowedge,
2274 highedge);
2275
2276 n_dxyPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2277 "norm_widths_dxy_phi",
2278 "width(d_{xy}/#sigma_{d_{xy}}) vs #phi sector;#varphi (sector) [degrees]; width(d_{xy}/#sigma_{d_{xy}})",
2279 nBins_,
2280 lowedge,
2281 highedge);
2282
2283 n_dzPhiMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2284 "norm_means_dz_phi",
2285 "#LT d_{z}/#sigma_{d_{z}} #GT vs #phi sector;#varphi (sector) [degrees];#LT d_{z}/#sigma_{d_{z}} #GT",
2286 nBins_,
2287 lowedge,
2288 highedge);
2289
2290 n_dzPhiWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2291 "norm_widths_dz_phi",
2292 "width(d_{z}/#sigma_{d_{z}}) vs #phi sector;#varphi (sector) [degrees];width(d_{z}/#sigma_{d_{z}})",
2293 nBins_,
2294 lowedge,
2295 highedge);
2296
2297 n_dxyEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2298 "norm_means_dxy_eta",
2299 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs #eta sector;#eta (sector);#LT d_{xy}/#sigma_{d_{z}} #GT",
2300 nBins_,
2301 lowedge,
2302 highedge);
2303
2304 n_dxyEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2305 "norm_widths_dxy_eta",
2306 "width(d_{xy}/#sigma_{d_{xy}}) vs #eta sector;#eta (sector);width(d_{xy}/#sigma_{d_{z}})",
2307 nBins_,
2308 lowedge,
2309 highedge);
2310
2311 n_dzEtaMeanBiasTrend = MeanBiasTrendsDir.make<TH1F>(
2312 "norm_means_dz_eta",
2313 "#LT d_{z}/#sigma_{d_{z}} #GT vs #eta sector;#eta (sector);#LT d_{z}/#sigma_{d_{z}} #GT",
2314 nBins_,
2315 lowedge,
2316 highedge);
2317
2318 n_dzEtaWidthBiasTrend = WidthBiasTrendsDir.make<TH1F>(
2319 "norm_widths_dz_eta",
2320 "width(d_{z}/#sigma_{d_{z}}) vs #eta sector;#eta (sector);width(d_{z}/#sigma_{d_{z}})",
2321 nBins_,
2322 lowedge,
2323 highedge);
2324
2325
2326
2327 a_dxyMeanBiasMap = Mean2DBiasMapsDir.make<TH2F>("means_dxy_map",
2328 "#LT d_{xy} #GT map;#eta (sector);#varphi (sector) [degrees]",
2329 nBins_,
2330 lowedge,
2331 highedge,
2332 nBins_,
2333 lowedge,
2334 highedge);
2335
2336 a_dzMeanBiasMap = Mean2DBiasMapsDir.make<TH2F>("means_dz_map",
2337 "#LT d_{z} #GT map;#eta (sector);#varphi (sector) [degrees]",
2338 nBins_,
2339 lowedge,
2340 highedge,
2341 nBins_,
2342 lowedge,
2343 highedge);
2344
2345 n_dxyMeanBiasMap =
2346 Mean2DBiasMapsDir.make<TH2F>("norm_means_dxy_map",
2347 "#LT d_{xy}/#sigma_{d_{xy}} #GT map;#eta (sector);#varphi (sector) [degrees]",
2348 nBins_,
2349 lowedge,
2350 highedge,
2351 nBins_,
2352 lowedge,
2353 highedge);
2354
2355 n_dzMeanBiasMap =
2356 Mean2DBiasMapsDir.make<TH2F>("norm_means_dz_map",
2357 "#LT d_{z}/#sigma_{d_{z}} #GT map;#eta (sector);#varphi (sector) [degrees]",
2358 nBins_,
2359 lowedge,
2360 highedge,
2361 nBins_,
2362 lowedge,
2363 highedge);
2364
2365 a_dxyWidthBiasMap = Width2DBiasMapsDir.make<TH2F>("widths_dxy_map",
2366 "#sigma_{d_{xy}} map;#eta (sector);#varphi (sector) [degrees]",
2367 nBins_,
2368 lowedge,
2369 highedge,
2370 nBins_,
2371 lowedge,
2372 highedge);
2373
2374 a_dzWidthBiasMap = Width2DBiasMapsDir.make<TH2F>("widths_dz_map",
2375 "#sigma_{d_{z}} map;#eta (sector);#varphi (sector) [degrees]",
2376 nBins_,
2377 lowedge,
2378 highedge,
2379 nBins_,
2380 lowedge,
2381 highedge);
2382
2383 n_dxyWidthBiasMap =
2384 Width2DBiasMapsDir.make<TH2F>("norm_widths_dxy_map",
2385 "width(d_{xy}/#sigma_{d_{xy}}) map;#eta (sector);#varphi (sector) [degrees]",
2386 nBins_,
2387 lowedge,
2388 highedge,
2389 nBins_,
2390 lowedge,
2391 highedge);
2392
2393 n_dzWidthBiasMap =
2394 Width2DBiasMapsDir.make<TH2F>("norm_widths_dz_map",
2395 "width(d_{z}/#sigma_{d_{z}}) map;#eta (sector);#varphi (sector) [degrees]",
2396 nBins_,
2397 lowedge,
2398 highedge,
2399 nBins_,
2400 lowedge,
2401 highedge);
2402
2403
2404
2405 a_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2406 "medians_dxy_phi",
2407 "Median of d_{xy} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}) [#mum]",
2408 nBins_,
2409 lowedge,
2410 highedge);
2411
2412 a_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2413 "MADs_dxy_phi",
2414 "Median absolute deviation of d_{xy} vs #phi sector;#varphi (sector) [degrees];MAD(d_{xy}) [#mum]",
2415 nBins_,
2416 lowedge,
2417 highedge);
2418
2419 a_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2420 "medians_dz_phi",
2421 "Median of d_{z} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}) [#mum]",
2422 nBins_,
2423 lowedge,
2424 highedge);
2425
2426 a_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2427 "MADs_dz_phi",
2428 "Median absolute deviation of d_{z} vs #phi sector;#varphi (sector) [degrees];MAD(d_{z}) [#mum]",
2429 nBins_,
2430 lowedge,
2431 highedge);
2432
2433 a_dxyEtaMedianBiasTrend =
2434 MedianBiasTrendsDir.make<TH1F>("medians_dxy_eta",
2435 "Median of d_{xy} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}) [#mum]",
2436 nBins_,
2437 lowedge,
2438 highedge);
2439
2440 a_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2441 "MADs_dxy_eta",
2442 "Median absolute deviation of d_{xy} vs #eta sector;#eta (sector);MAD(d_{xy}) [#mum]",
2443 nBins_,
2444 lowedge,
2445 highedge);
2446
2447 a_dzEtaMedianBiasTrend =
2448 MedianBiasTrendsDir.make<TH1F>("medians_dz_eta",
2449 "Median of d_{z} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}) [#mum]",
2450 nBins_,
2451 lowedge,
2452 highedge);
2453
2454 a_dzEtaMADBiasTrend =
2455 MADBiasTrendsDir.make<TH1F>("MADs_dz_eta",
2456 "Median absolute deviation of d_{z} vs #eta sector;#eta (sector);MAD(d_{z}) [#mum]",
2457 nBins_,
2458 lowedge,
2459 highedge);
2460
2461 n_dxyPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2462 "norm_medians_dxy_phi",
2463 "Median of d_{xy}/#sigma_{d_{xy}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{xy}/#sigma_{d_{xy}})",
2464 nBins_,
2465 lowedge,
2466 highedge);
2467
2468 n_dxyPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>("norm_MADs_dxy_phi",
2469 "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #phi "
2470 "sector;#varphi (sector) [degrees]; MAD(d_{xy}/#sigma_{d_{xy}})",
2471 nBins_,
2472 lowedge,
2473 highedge);
2474
2475 n_dzPhiMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2476 "norm_medians_dz_phi",
2477 "Median of d_{z}/#sigma_{d_{z}} vs #phi sector;#varphi (sector) [degrees];#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2478 nBins_,
2479 lowedge,
2480 highedge);
2481
2482 n_dzPhiMADBiasTrend = MADBiasTrendsDir.make<TH1F>("norm_MADs_dz_phi",
2483 "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #phi "
2484 "sector;#varphi (sector) [degrees];MAD(d_{z}/#sigma_{d_{z}})",
2485 nBins_,
2486 lowedge,
2487 highedge);
2488
2489 n_dxyEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2490 "norm_medians_dxy_eta",
2491 "Median of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{xy}/#sigma_{d_{z}})",
2492 nBins_,
2493 lowedge,
2494 highedge);
2495
2496 n_dxyEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2497 "norm_MADs_dxy_eta",
2498 "Median absolute deviation of d_{xy}/#sigma_{d_{xy}} vs #eta sector;#eta (sector);MAD(d_{xy}/#sigma_{d_{z}})",
2499 nBins_,
2500 lowedge,
2501 highedge);
2502
2503 n_dzEtaMedianBiasTrend = MedianBiasTrendsDir.make<TH1F>(
2504 "norm_medians_dz_eta",
2505 "Median of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);#mu_{1/2}(d_{z}/#sigma_{d_{z}})",
2506 nBins_,
2507 lowedge,
2508 highedge);
2509
2510 n_dzEtaMADBiasTrend = MADBiasTrendsDir.make<TH1F>(
2511 "norm_MADs_dz_eta",
2512 "Median absolute deviation of d_{z}/#sigma_{d_{z}} vs #eta sector;#eta (sector);MAD(d_{z}/#sigma_{d_{z}})",
2513 nBins_,
2514 lowedge,
2515 highedge);
2516 }
2517 }
2518
2519
2520 void PrimaryVertexValidation::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
2521
2522
2523
2524
2525 if (!isBFieldConsistentWithMode(iSetup)) {
2526 edm::LogWarning("PrimaryVertexValidation")
2527 << "*********************************************************************************\n"
2528 << "* The configuration (ptOfProbe > " << ptOfProbe_
2529 << "GeV) is not correctly set for current value of magnetic field \n"
2530 << "* Switching it to 0. !!! \n"
2531 << "*********************************************************************************" << std::endl;
2532 ptOfProbe_ = 0.;
2533 }
2534
2535
2536
2537
2538
2539 const auto runNumber_ = iRun.run();
2540 h_runNumber->Fill(runNumber_);
2541
2542 if (!runNumbersTimesLog_.count(runNumber_)) {
2543 auto times = getRunTime(iSetup);
2544
2545 if (debug_) {
2546 const time_t start_time = times.first / 1000000;
2547 edm::LogInfo("PrimaryVertexValidation")
2548 << runNumber_ << " has start time: " << times.first << " - " << times.second << std::endl;
2549 edm::LogInfo("PrimaryVertexValidation")
2550 << "human readable time: " << std::asctime(std::gmtime(&start_time)) << std::endl;
2551 }
2552
2553 runNumbersTimesLog_[runNumber_] = times;
2554 }
2555
2556 if (h_runFromEvent->GetEntries() == 0) {
2557 h_runFromEvent->SetBinContent(1, RunNumber_);
2558 }
2559
2560
2561
2562
2563
2564 const TrackerTopology* const tTopo = &iSetup.getData(topoTokenBR_);
2565
2566
2567
2568
2569
2570 edm::LogInfo("read tracker geometry...");
2571 edm::ESHandle<TrackerGeometry> pDD = iSetup.getHandle(geomTokenBR_);
2572 edm::LogInfo("tracker geometry read") << "There are: " << pDD->dets().size() << " detectors";
2573
2574 PixelTopologyMap PTMap = PixelTopologyMap(pDD.product(), tTopo);
2575 nLadders_ = PTMap.getPXBLadders(1);
2576 nModZ_ = PTMap.getPXBModules(1);
2577
2578
2579
2580
2581
2582 if (h_dxy_ladderOverlap_.size() != nLadders_) {
2583 PVValHelper::shrinkHistVectorToFit(h_dxy_ladderOverlap_, nLadders_);
2584 PVValHelper::shrinkHistVectorToFit(h_dxy_ladderNoOverlap_, nLadders_);
2585 PVValHelper::shrinkHistVectorToFit(h_dxy_ladder_, nLadders_);
2586 PVValHelper::shrinkHistVectorToFit(h_dz_ladder_, nLadders_);
2587 PVValHelper::shrinkHistVectorToFit(h_norm_dxy_ladder_, nLadders_);
2588 PVValHelper::shrinkHistVectorToFit(h_norm_dz_ladder_, nLadders_);
2589
2590 if (debug_) {
2591 edm::LogInfo("PrimaryVertexValidation") << "checking size:" << h_dxy_ladder_.size() << std::endl;
2592 }
2593 }
2594
2595 if (h_dxy_modZ_.size() != nModZ_) {
2596 PVValHelper::shrinkHistVectorToFit(h_dxy_modZ_, nModZ_);
2597 PVValHelper::shrinkHistVectorToFit(h_dz_modZ_, nModZ_);
2598 PVValHelper::shrinkHistVectorToFit(h_norm_dxy_modZ_, nModZ_);
2599 PVValHelper::shrinkHistVectorToFit(h_norm_dxy_modZ_, nModZ_);
2600
2601 if (debug_) {
2602 edm::LogInfo("PrimaryVertexValidation") << "checking size:" << h_dxy_modZ_.size() << std::endl;
2603 }
2604 }
2605
2606 if ((pDD->isThere(GeomDetEnumerators::P2PXB)) || (pDD->isThere(GeomDetEnumerators::P2PXEC))) {
2607
2608 phase_ = PVValHelper::phase2;
2609 if (debug_) {
2610 edm::LogInfo("PrimaryVertexValidation")
2611 << " pixel phase2 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2612 }
2613
2614 } else if ((pDD->isThere(GeomDetEnumerators::P1PXB)) || (pDD->isThere(GeomDetEnumerators::P1PXEC))) {
2615
2616 if (debug_) {
2617 edm::LogInfo("PrimaryVertexValidation")
2618 << " pixel phase1 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2619 }
2620 } else {
2621
2622 phase_ = PVValHelper::phase0;
2623 if (debug_) {
2624 edm::LogInfo("PrimaryVertexValidation")
2625 << " pixel phase0 setup, nLadders: " << nLadders_ << " nModules:" << nModZ_;
2626 }
2627 }
2628
2629 switch (phase_) {
2630 case PVValHelper::phase0:
2631 etaOfProbe_ = std::min(etaOfProbe_, PVValHelper::max_eta_phase0);
2632 break;
2633 case PVValHelper::phase1:
2634 etaOfProbe_ = std::min(etaOfProbe_, PVValHelper::max_eta_phase1);
2635 break;
2636 case PVValHelper::phase2:
2637 etaOfProbe_ = std::min(etaOfProbe_, PVValHelper::max_eta_phase2);
2638 break;
2639 default:
2640 edm::LogWarning("LogicError") << "Unknown detector phase: " << phase_;
2641 }
2642
2643 if (h_etaMax->GetEntries() == 0.) {
2644 h_etaMax->SetBinContent(1., etaOfProbe_);
2645 h_nbins->SetBinContent(1., nBins_);
2646 h_nLadders->SetBinContent(1., nLadders_);
2647 h_nModZ->SetBinContent(1., nModZ_);
2648 h_pTinfo->SetBinContent(1., mypT_bins_.size());
2649 h_pTinfo->SetBinContent(2., minPt_);
2650 h_pTinfo->SetBinContent(3., maxPt_);
2651 }
2652 }
2653
2654
2655 void PrimaryVertexValidation::endJob() {
2656
2657 h_probeL1Ladder_->GetXaxis()->SetRangeUser(-1.5, nLadders_ + 0.5);
2658 h_probeL1Module_->GetXaxis()->SetRangeUser(-1.5, nModZ_ + 0.5);
2659 h2_probeLayer1Map_->GetXaxis()->SetRangeUser(0.5, nModZ_ + 0.5);
2660 h2_probeLayer1Map_->GetYaxis()->SetRangeUser(0.5, nLadders_ + 0.5);
2661 h2_probePassingLayer1Map_->GetXaxis()->SetRangeUser(0.5, nModZ_ + 0.5);
2662 h2_probePassingLayer1Map_->GetYaxis()->SetRangeUser(0.5, nLadders_ + 0.5);
2663
2664 TFileDirectory RunFeatures = fs->mkdir("RunFeatures");
2665 h_runStartTimes = RunFeatures.make<TH1I>(
2666 "runStartTimes", "run start times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
2667 h_runEndTimes =
2668 RunFeatures.make<TH1I>("runEndTimes", "run end times", runNumbersTimesLog_.size(), 0, runNumbersTimesLog_.size());
2669
2670 unsigned int count = 1;
2671 for (const auto& run : runNumbersTimesLog_) {
2672
2673 h_runStartTimes->SetBinContent(count, run.second.first / 10e6);
2674 h_runStartTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
2675
2676 h_runEndTimes->SetBinContent(count, run.second.second / 10e6);
2677 h_runEndTimes->GetXaxis()->SetBinLabel(count, (std::to_string(run.first)).c_str());
2678
2679 count++;
2680 }
2681
2682 edm::LogInfo("PrimaryVertexValidation") << "######################################\n"
2683 << "# PrimaryVertexValidation::endJob()\n"
2684 << "# Number of analyzed events: " << Nevt_ << "\n"
2685 << "######################################";
2686
2687
2688
2689 a_dxymodZMeanTrend = MeanTrendsDir.make<TH1F>(
2690 "means_dxy_modZ", "#LT d_{xy} #GT vs modZ;module number (Z);#LT d_{xy} #GT [#mum]", nModZ_, 0., nModZ_);
2691
2692 a_dxymodZWidthTrend = WidthTrendsDir.make<TH1F>(
2693 "widths_dxy_modZ", "#sigma_{d_{xy}} vs modZ;module number (Z);#sigma_{d_{xy}} [#mum]", nModZ_, 0., nModZ_);
2694
2695 a_dzmodZMeanTrend = MeanTrendsDir.make<TH1F>(
2696 "means_dz_modZ", "#LT d_{z} #GT vs modZ;module number (Z);#LT d_{z} #GT [#mum]", nModZ_, 0., nModZ_);
2697
2698 a_dzmodZWidthTrend = WidthTrendsDir.make<TH1F>(
2699 "widths_dz_modZ", "#sigma_{d_{z}} vs modZ;module number (Z);#sigma_{d_{z}} [#mum]", nModZ_, 0., nModZ_);
2700
2701 a_dxyladderMeanTrend = MeanTrendsDir.make<TH1F>("means_dxy_ladder",
2702 "#LT d_{xy} #GT vs ladder;ladder number (#phi);#LT d_{xy} #GT [#mum]",
2703 nLadders_,
2704 0.,
2705 nLadders_);
2706
2707 a_dxyladderWidthTrend =
2708 WidthTrendsDir.make<TH1F>("widths_dxy_ladder",
2709 "#sigma_{d_{xy}} vs ladder;ladder number (#phi);#sigma_{d_{xy}} [#mum]",
2710 nLadders_,
2711 0.,
2712 nLadders_);
2713
2714 a_dzladderMeanTrend = MeanTrendsDir.make<TH1F>(
2715 "means_dz_ladder", "#LT d_{z} #GT vs ladder;ladder number (#phi);#LT d_{z} #GT [#mum]", nLadders_, 0., nLadders_);
2716
2717 a_dzladderWidthTrend =
2718 WidthTrendsDir.make<TH1F>("widths_dz_ladder",
2719 "#sigma_{d_{z}} vs ladder;ladder number (#phi);#sigma_{d_{z}} [#mum]",
2720 nLadders_,
2721 0.,
2722 nLadders_);
2723
2724 n_dxymodZMeanTrend = MeanTrendsDir.make<TH1F>(
2725 "norm_means_dxy_modZ",
2726 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs modZ;module number (Z);#LT d_{xy}/#sigma_{d_{xy}} #GT",
2727 nModZ_,
2728 0.,
2729 nModZ_);
2730
2731 n_dxymodZWidthTrend = WidthTrendsDir.make<TH1F>(
2732 "norm_widths_dxy_modZ",
2733 "width(d_{xy}/#sigma_{d_{xy}}) vs modZ;module number (Z); width(d_{xy}/#sigma_{d_{xy}})",
2734 nModZ_,
2735 0.,
2736 nModZ_);
2737
2738 n_dzmodZMeanTrend =
2739 MeanTrendsDir.make<TH1F>("norm_means_dz_modZ",
2740 "#LT d_{z}/#sigma_{d_{z}} #GT vs modZ;module number (Z);#LT d_{z}/#sigma_{d_{z}} #GT",
2741 nModZ_,
2742 0.,
2743 nModZ_);
2744
2745 n_dzmodZWidthTrend =
2746 WidthTrendsDir.make<TH1F>("norm_widths_dz_modZ",
2747 "width(d_{z}/#sigma_{d_{z}}) vs pT;module number (Z);width(d_{z}/#sigma_{d_{z}})",
2748 nModZ_,
2749 0.,
2750 nModZ_);
2751
2752 n_dxyladderMeanTrend = MeanTrendsDir.make<TH1F>(
2753 "norm_means_dxy_ladder",
2754 "#LT d_{xy}/#sigma_{d_{xy}} #GT vs ladder;ladder number (#phi);#LT d_{xy}/#sigma_{d_{z}} #GT",
2755 nLadders_,
2756 0.,
2757 nLadders_);
2758
2759 n_dxyladderWidthTrend = WidthTrendsDir.make<TH1F>(
2760 "norm_widths_dxy_ladder",
2761 "width(d_{xy}/#sigma_{d_{xy}}) vs ladder;ladder number (#phi);width(d_{xy}/#sigma_{d_{z}})",
2762 nLadders_,
2763 0.,
2764 nLadders_);
2765
2766 n_dzladderMeanTrend = MeanTrendsDir.make<TH1F>(
2767 "norm_means_dz_ladder",
2768 "#LT d_{z}/#sigma_{d_{z}} #GT vs ladder;ladder number (#phi);#LT d_{z}/#sigma_{d_{z}} #GT",
2769 nLadders_,
2770 0.,
2771 nLadders_);
2772
2773 n_dzladderWidthTrend = WidthTrendsDir.make<TH1F>(
2774 "norm_widths_dz_ladder",
2775 "width(d_{z}/#sigma_{d_{z}}) vs ladder;ladder number (#phi);width(d_{z}/#sigma_{d_{z}})",
2776 nLadders_,
2777 0.,
2778 nLadders_);
2779
2780
2781
2782 a_dxyL1MeanMap = Mean2DMapsDir.make<TH2F>("means_dxy_l1map",
2783 "#LT d_{xy} #GT map;module number [z];ladder number [#varphi]",
2784 nModZ_,
2785 0.,
2786 nModZ_,
2787 nLadders_,
2788 0.,
2789 nLadders_);
2790
2791 a_dzL1MeanMap = Mean2DMapsDir.make<TH2F>("means_dz_l1map",
2792 "#LT d_{z} #GT map;module number [z];ladder number [#varphi]",
2793 nModZ_,
2794 0.,
2795 nModZ_,
2796 nLadders_,
2797 0.,
2798 nLadders_);
2799
2800 n_dxyL1MeanMap =
2801 Mean2DMapsDir.make<TH2F>("norm_means_dxy_l1map",
2802 "#LT d_{xy}/#sigma_{d_{xy}} #GT map;module number [z];ladder number [#varphi]",
2803 nModZ_,
2804 0.,
2805 nModZ_,
2806 nLadders_,
2807 0.,
2808 nLadders_);
2809
2810 n_dzL1MeanMap = Mean2DMapsDir.make<TH2F>("norm_means_dz_l1map",
2811 "#LT d_{z}/#sigma_{d_{z}} #GT map;module number [z];ladder number [#varphi]",
2812 nModZ_,
2813 0.,
2814 nModZ_,
2815 nLadders_,
2816 0.,
2817 nLadders_);
2818
2819 a_dxyL1WidthMap = Width2DMapsDir.make<TH2F>("widths_dxy_l1map",
2820 "#sigma_{d_{xy}} map;module number [z];ladder number [#varphi]",
2821 nModZ_,
2822 0.,
2823 nModZ_,
2824 nLadders_,
2825 0.,
2826 nLadders_);
2827
2828 a_dzL1WidthMap = Width2DMapsDir.make<TH2F>("widths_dz_l1map",
2829 "#sigma_{d_{z}} map;module number [z];ladder number [#varphi]",
2830 nModZ_,
2831 0.,
2832 nModZ_,
2833 nLadders_,
2834 0.,
2835 nLadders_);
2836
2837 n_dxyL1WidthMap =
2838 Width2DMapsDir.make<TH2F>("norm_widths_dxy_l1map",
2839 "width(d_{xy}/#sigma_{d_{xy}}) map;module number [z];ladder number [#varphi]",
2840 nModZ_,
2841 0.,
2842 nModZ_,
2843 nLadders_,
2844 0.,
2845 nLadders_);
2846
2847 n_dzL1WidthMap =
2848 Width2DMapsDir.make<TH2F>("norm_widths_dz_l1map",
2849 "width(d_{z}/#sigma_{d_{z}}) map;module number [z];ladder number [#varphi]",
2850 nModZ_,
2851 0.,
2852 nModZ_,
2853 nLadders_,
2854 0.,
2855 nLadders_);
2856
2857 if (useTracksFromRecoVtx_) {
2858 fillTrendPlotByIndex(a_dxyPhiMeanBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2859 fillTrendPlotByIndex(a_dxyPhiWidthBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2860 fillTrendPlotByIndex(a_dzPhiMeanBiasTrend, a_dzPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2861 fillTrendPlotByIndex(a_dzPhiWidthBiasTrend, a_dzPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2862
2863 fillTrendPlotByIndex(a_dxyEtaMeanBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2864 fillTrendPlotByIndex(a_dxyEtaWidthBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2865 fillTrendPlotByIndex(a_dzEtaMeanBiasTrend, a_dzEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2866 fillTrendPlotByIndex(a_dzEtaWidthBiasTrend, a_dzEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2867
2868 fillTrendPlotByIndex(n_dxyPhiMeanBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2869 fillTrendPlotByIndex(n_dxyPhiWidthBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2870 fillTrendPlotByIndex(n_dzPhiMeanBiasTrend, n_dzPhiBiasResiduals, PVValHelper::MEAN, PVValHelper::phi);
2871 fillTrendPlotByIndex(n_dzPhiWidthBiasTrend, n_dzPhiBiasResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2872
2873 fillTrendPlotByIndex(n_dxyEtaMeanBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2874 fillTrendPlotByIndex(n_dxyEtaWidthBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2875 fillTrendPlotByIndex(n_dzEtaMeanBiasTrend, n_dzEtaBiasResiduals, PVValHelper::MEAN, PVValHelper::eta);
2876 fillTrendPlotByIndex(n_dzEtaWidthBiasTrend, n_dzEtaBiasResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2877
2878
2879
2880 fillTrendPlotByIndex(a_dxyPhiMedianBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2881 fillTrendPlotByIndex(a_dxyPhiMADBiasTrend, a_dxyPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2882 fillTrendPlotByIndex(a_dzPhiMedianBiasTrend, a_dzPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2883 fillTrendPlotByIndex(a_dzPhiMADBiasTrend, a_dzPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2884
2885 fillTrendPlotByIndex(a_dxyEtaMedianBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2886 fillTrendPlotByIndex(a_dxyEtaMADBiasTrend, a_dxyEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2887 fillTrendPlotByIndex(a_dzEtaMedianBiasTrend, a_dzEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2888 fillTrendPlotByIndex(a_dzEtaMADBiasTrend, a_dzEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2889
2890 fillTrendPlotByIndex(n_dxyPhiMedianBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2891 fillTrendPlotByIndex(n_dxyPhiMADBiasTrend, n_dxyPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2892 fillTrendPlotByIndex(n_dzPhiMedianBiasTrend, n_dzPhiBiasResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2893 fillTrendPlotByIndex(n_dzPhiMADBiasTrend, n_dzPhiBiasResiduals, PVValHelper::MAD, PVValHelper::phi);
2894
2895 fillTrendPlotByIndex(n_dxyEtaMedianBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2896 fillTrendPlotByIndex(n_dxyEtaMADBiasTrend, n_dxyEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2897 fillTrendPlotByIndex(n_dzEtaMedianBiasTrend, n_dzEtaBiasResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2898 fillTrendPlotByIndex(n_dzEtaMADBiasTrend, n_dzEtaBiasResiduals, PVValHelper::MAD, PVValHelper::eta);
2899
2900
2901
2902 fillMap(a_dxyMeanBiasMap, a_dxyBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2903 fillMap(a_dxyWidthBiasMap, a_dxyBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2904 fillMap(a_dzMeanBiasMap, a_dzBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2905 fillMap(a_dzWidthBiasMap, a_dzBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2906
2907 fillMap(n_dxyMeanBiasMap, n_dxyBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2908 fillMap(n_dxyWidthBiasMap, n_dxyBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2909 fillMap(n_dzMeanBiasMap, n_dzBiasResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
2910 fillMap(n_dzWidthBiasMap, n_dzBiasResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
2911 }
2912
2913
2914
2915 fillTrendPlotByIndex(a_dxyPhiMeanTrend, a_dxyPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2916 fillTrendPlotByIndex(a_dxyPhiWidthTrend, a_dxyPhiResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2917 fillTrendPlotByIndex(a_dzPhiMeanTrend, a_dzPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2918 fillTrendPlotByIndex(a_dzPhiWidthTrend, a_dzPhiResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2919
2920 fillTrendPlotByIndex(a_dxyEtaMeanTrend, a_dxyEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2921 fillTrendPlotByIndex(a_dxyEtaWidthTrend, a_dxyEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2922 fillTrendPlotByIndex(a_dzEtaMeanTrend, a_dzEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2923 fillTrendPlotByIndex(a_dzEtaWidthTrend, a_dzEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2924
2925 fillTrendPlotByIndex(n_dxyPhiMeanTrend, n_dxyPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2926 fillTrendPlotByIndex(n_dxyPhiWidthTrend, n_dxyPhiResiduals, PVValHelper::WIDTH, PVValHelper::phi);
2927 fillTrendPlotByIndex(n_dzPhiMeanTrend, n_dzPhiResiduals, PVValHelper::MEAN, PVValHelper::phi);
2928 fillTrendPlotByIndex(n_dzPhiWidthTrend, n_dzPhiResiduals, PVValHelper::WIDTH);
2929
2930 fillTrendPlotByIndex(n_dxyEtaMeanTrend, n_dxyEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2931 fillTrendPlotByIndex(n_dxyEtaWidthTrend, n_dxyEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2932 fillTrendPlotByIndex(n_dzEtaMeanTrend, n_dzEtaResiduals, PVValHelper::MEAN, PVValHelper::eta);
2933 fillTrendPlotByIndex(n_dzEtaWidthTrend, n_dzEtaResiduals, PVValHelper::WIDTH, PVValHelper::eta);
2934
2935
2936
2937 fillTrendPlotByIndex(a_dxypTMeanTrend, h_dxy_pT_, PVValHelper::MEAN);
2938 fillTrendPlotByIndex(a_dxypTWidthTrend, h_dxy_pT_, PVValHelper::WIDTH);
2939 fillTrendPlotByIndex(a_dzpTMeanTrend, h_dz_pT_, PVValHelper::MEAN);
2940 fillTrendPlotByIndex(a_dzpTWidthTrend, h_dz_pT_, PVValHelper::WIDTH);
2941
2942 fillTrendPlotByIndex(a_dxypTCentralMeanTrend, h_dxy_Central_pT_, PVValHelper::MEAN);
2943 fillTrendPlotByIndex(a_dxypTCentralWidthTrend, h_dxy_Central_pT_, PVValHelper::WIDTH);
2944 fillTrendPlotByIndex(a_dzpTCentralMeanTrend, h_dz_Central_pT_, PVValHelper::MEAN);
2945 fillTrendPlotByIndex(a_dzpTCentralWidthTrend, h_dz_Central_pT_, PVValHelper::WIDTH);
2946
2947 fillTrendPlotByIndex(n_dxypTMeanTrend, h_norm_dxy_pT_, PVValHelper::MEAN);
2948 fillTrendPlotByIndex(n_dxypTWidthTrend, h_norm_dxy_pT_, PVValHelper::WIDTH);
2949 fillTrendPlotByIndex(n_dzpTMeanTrend, h_norm_dz_pT_, PVValHelper::MEAN);
2950 fillTrendPlotByIndex(n_dzpTWidthTrend, h_norm_dz_pT_, PVValHelper::WIDTH);
2951
2952 fillTrendPlotByIndex(n_dxypTCentralMeanTrend, h_norm_dxy_Central_pT_, PVValHelper::MEAN);
2953 fillTrendPlotByIndex(n_dxypTCentralWidthTrend, h_norm_dxy_Central_pT_, PVValHelper::WIDTH);
2954 fillTrendPlotByIndex(n_dzpTCentralMeanTrend, h_norm_dz_Central_pT_, PVValHelper::MEAN);
2955 fillTrendPlotByIndex(n_dzpTCentralWidthTrend, h_norm_dz_Central_pT_, PVValHelper::WIDTH);
2956
2957
2958
2959 fillTrendPlotByIndex(a_dxymodZMeanTrend, h_dxy_modZ_, PVValHelper::MEAN);
2960 fillTrendPlotByIndex(a_dxymodZWidthTrend, h_dxy_modZ_, PVValHelper::WIDTH);
2961 fillTrendPlotByIndex(a_dzmodZMeanTrend, h_dz_modZ_, PVValHelper::MEAN);
2962 fillTrendPlotByIndex(a_dzmodZWidthTrend, h_dz_modZ_, PVValHelper::WIDTH);
2963
2964 fillTrendPlotByIndex(a_dxyladderMeanTrend, h_dxy_ladder_, PVValHelper::MEAN);
2965 fillTrendPlotByIndex(a_dxyladderWidthTrend, h_dxy_ladder_, PVValHelper::WIDTH);
2966 fillTrendPlotByIndex(a_dzladderMeanTrend, h_dz_ladder_, PVValHelper::MEAN);
2967 fillTrendPlotByIndex(a_dzladderWidthTrend, h_dz_ladder_, PVValHelper::WIDTH);
2968
2969 fillTrendPlotByIndex(n_dxymodZMeanTrend, h_norm_dxy_modZ_, PVValHelper::MEAN);
2970 fillTrendPlotByIndex(n_dxymodZWidthTrend, h_norm_dxy_modZ_, PVValHelper::WIDTH);
2971 fillTrendPlotByIndex(n_dzmodZMeanTrend, h_norm_dz_modZ_, PVValHelper::MEAN);
2972 fillTrendPlotByIndex(n_dzmodZWidthTrend, h_norm_dz_modZ_, PVValHelper::WIDTH);
2973
2974 fillTrendPlotByIndex(n_dxyladderMeanTrend, h_norm_dxy_ladder_, PVValHelper::MEAN);
2975 fillTrendPlotByIndex(n_dxyladderWidthTrend, h_norm_dxy_ladder_, PVValHelper::WIDTH);
2976 fillTrendPlotByIndex(n_dzladderMeanTrend, h_norm_dz_ladder_, PVValHelper::MEAN);
2977 fillTrendPlotByIndex(n_dzladderWidthTrend, h_norm_dz_ladder_, PVValHelper::WIDTH);
2978
2979
2980
2981 fillTrendPlotByIndex(a_dxyPhiMedianTrend, a_dxyPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2982 fillTrendPlotByIndex(a_dxyPhiMADTrend, a_dxyPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2983
2984 fillTrendPlotByIndex(a_dzPhiMedianTrend, a_dzPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2985 fillTrendPlotByIndex(a_dzPhiMADTrend, a_dzPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2986
2987 fillTrendPlotByIndex(a_dxyEtaMedianTrend, a_dxyEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2988 fillTrendPlotByIndex(a_dxyEtaMADTrend, a_dxyEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2989 fillTrendPlotByIndex(a_dzEtaMedianTrend, a_dzEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2990 fillTrendPlotByIndex(a_dzEtaMADTrend, a_dzEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2991
2992 fillTrendPlotByIndex(n_dxyPhiMedianTrend, n_dxyPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2993 fillTrendPlotByIndex(n_dxyPhiMADTrend, n_dxyPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2994 fillTrendPlotByIndex(n_dzPhiMedianTrend, n_dzPhiResiduals, PVValHelper::MEDIAN, PVValHelper::phi);
2995 fillTrendPlotByIndex(n_dzPhiMADTrend, n_dzPhiResiduals, PVValHelper::MAD, PVValHelper::phi);
2996
2997 fillTrendPlotByIndex(n_dxyEtaMedianTrend, n_dxyEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
2998 fillTrendPlotByIndex(n_dxyEtaMADTrend, n_dxyEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
2999 fillTrendPlotByIndex(n_dzEtaMedianTrend, n_dzEtaResiduals, PVValHelper::MEDIAN, PVValHelper::eta);
3000 fillTrendPlotByIndex(n_dzEtaMADTrend, n_dzEtaResiduals, PVValHelper::MAD, PVValHelper::eta);
3001
3002
3003
3004 fillMap(a_dxyMeanMap, a_dxyResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3005 fillMap(a_dxyWidthMap, a_dxyResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3006 fillMap(a_dzMeanMap, a_dzResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3007 fillMap(a_dzWidthMap, a_dzResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3008
3009 fillMap(n_dxyMeanMap, n_dxyResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3010 fillMap(n_dxyWidthMap, n_dxyResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3011 fillMap(n_dzMeanMap, n_dzResidualsMap, PVValHelper::MEAN, nBins_, nBins_);
3012 fillMap(n_dzWidthMap, n_dzResidualsMap, PVValHelper::WIDTH, nBins_, nBins_);
3013
3014
3015
3016 fillMap(a_dxyL1MeanMap, a_dxyL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3017 fillMap(a_dxyL1WidthMap, a_dxyL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3018 fillMap(a_dzL1MeanMap, a_dzL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3019 fillMap(a_dzL1WidthMap, a_dzL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3020
3021 fillMap(n_dxyL1MeanMap, n_dxyL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3022 fillMap(n_dxyL1WidthMap, n_dxyL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3023 fillMap(n_dzL1MeanMap, n_dzL1ResidualsMap, PVValHelper::MEAN, nModZ_, nLadders_);
3024 fillMap(n_dzL1WidthMap, n_dzL1ResidualsMap, PVValHelper::WIDTH, nModZ_, nLadders_);
3025 }
3026
3027
3028 std::pair<long long, long long> PrimaryVertexValidation::getRunTime(const edm::EventSetup& iSetup) const
3029
3030 {
3031 edm::ESHandle<RunInfo> runInfo = iSetup.getHandle(runInfoTokenBR_);
3032 if (debug_) {
3033 edm::LogInfo("PrimaryVertexValidation")
3034 << runInfo.product()->m_start_time_str << " " << runInfo.product()->m_stop_time_str << std::endl;
3035 }
3036 return std::make_pair(runInfo.product()->m_start_time_ll, runInfo.product()->m_stop_time_ll);
3037 }
3038
3039
3040 bool PrimaryVertexValidation::isBFieldConsistentWithMode(const edm::EventSetup& iSetup) const
3041
3042 {
3043 edm::ESHandle<RunInfo> runInfo = iSetup.getHandle(runInfoTokenBR_);
3044 double average_current = runInfo.product()->m_avg_current;
3045 bool isOn = (average_current > 2000.);
3046 bool is0T = (ptOfProbe_ == 0.);
3047
3048 return ((isOn && !is0T) || (!isOn && is0T));
3049 }
3050
3051
3052 void PrimaryVertexValidation::SetVarToZero()
3053
3054 {
3055 nTracks_ = 0;
3056 nClus_ = 0;
3057 nOfflineVertices_ = 0;
3058 RunNumber_ = 0;
3059 LuminosityBlockNumber_ = 0;
3060 xOfflineVertex_ = -999.;
3061 yOfflineVertex_ = -999.;
3062 zOfflineVertex_ = -999.;
3063 xErrOfflineVertex_ = 0.;
3064 yErrOfflineVertex_ = 0.;
3065 zErrOfflineVertex_ = 0.;
3066 BSx0_ = -999.;
3067 BSy0_ = -999.;
3068 BSz0_ = -999.;
3069 Beamsigmaz_ = -999.;
3070 Beamdxdz_ = -999.;
3071 BeamWidthX_ = -999.;
3072 BeamWidthY_ = -999.;
3073 wxy2_ = -999.;
3074
3075 for (int i = 0; i < nMaxtracks_; ++i) {
3076 pt_[i] = 0;
3077 p_[i] = 0;
3078 nhits_[i] = 0;
3079 nhits1D_[i] = 0;
3080 nhits2D_[i] = 0;
3081 nhitsBPIX_[i] = 0;
3082 nhitsFPIX_[i] = 0;
3083 nhitsTIB_[i] = 0;
3084 nhitsTID_[i] = 0;
3085 nhitsTOB_[i] = 0;
3086 nhitsTEC_[i] = 0;
3087 isHighPurity_[i] = 0;
3088 eta_[i] = 0;
3089 theta_[i] = 0;
3090 phi_[i] = 0;
3091 chi2_[i] = 0;
3092 chi2ndof_[i] = 0;
3093 charge_[i] = 0;
3094 qoverp_[i] = 0;
3095 dz_[i] = 0;
3096 dxy_[i] = 0;
3097 dzBs_[i] = 0;
3098 dxyBs_[i] = 0;
3099 xPCA_[i] = 0;
3100 yPCA_[i] = 0;
3101 zPCA_[i] = 0;
3102 xUnbiasedVertex_[i] = 0;
3103 yUnbiasedVertex_[i] = 0;
3104 zUnbiasedVertex_[i] = 0;
3105 chi2normUnbiasedVertex_[i] = 0;
3106 chi2UnbiasedVertex_[i] = 0;
3107 chi2ProbUnbiasedVertex_[i] = 0;
3108 DOFUnbiasedVertex_[i] = 0;
3109 sumOfWeightsUnbiasedVertex_[i] = 0;
3110 tracksUsedForVertexing_[i] = 0;
3111 dxyFromMyVertex_[i] = 0;
3112 dzFromMyVertex_[i] = 0;
3113 d3DFromMyVertex_[i] = 0;
3114 dxyErrorFromMyVertex_[i] = 0;
3115 dzErrorFromMyVertex_[i] = 0;
3116 d3DErrorFromMyVertex_[i] = 0;
3117 IPTsigFromMyVertex_[i] = 0;
3118 IPLsigFromMyVertex_[i] = 0;
3119 IP3DsigFromMyVertex_[i] = 0;
3120 hasRecVertex_[i] = 0;
3121 isGoodTrack_[i] = 0;
3122 }
3123 }
3124
3125
3126 void PrimaryVertexValidation::fillTrendPlot(TH1F* trendPlot,
3127 TH1F* residualsPlot[100],
3128 PVValHelper::estimator fitPar_,
3129 const std::string& var_)
3130
3131 {
3132 for (int i = 0; i < nBins_; ++i) {
3133 char phibincenter[129];
3134 auto phiBins = theDetails_.trendbins[PVValHelper::phi];
3135 sprintf(phibincenter, "%.f", (phiBins[i] + phiBins[i + 1]) / 2.);
3136
3137 char etabincenter[129];
3138 auto etaBins = theDetails_.trendbins[PVValHelper::eta];
3139 sprintf(etabincenter, "%.1f", (etaBins[i] + etaBins[i + 1]) / 2.);
3140
3141 switch (fitPar_) {
3142 case PVValHelper::MEAN: {
3143 float mean_ = PVValHelper::fitResiduals(residualsPlot[i]).first.value();
3144 float meanErr_ = PVValHelper::fitResiduals(residualsPlot[i]).first.error();
3145 trendPlot->SetBinContent(i + 1, mean_);
3146 trendPlot->SetBinError(i + 1, meanErr_);
3147 break;
3148 }
3149 case PVValHelper::WIDTH: {
3150 float width_ = PVValHelper::fitResiduals(residualsPlot[i]).second.value();
3151 float widthErr_ = PVValHelper::fitResiduals(residualsPlot[i]).second.error();
3152 trendPlot->SetBinContent(i + 1, width_);
3153 trendPlot->SetBinError(i + 1, widthErr_);
3154 break;
3155 }
3156 case PVValHelper::MEDIAN: {
3157 float median_ = PVValHelper::getMedian(residualsPlot[i]).value();
3158 float medianErr_ = PVValHelper::getMedian(residualsPlot[i]).error();
3159 trendPlot->SetBinContent(i + 1, median_);
3160 trendPlot->SetBinError(i + 1, medianErr_);
3161 break;
3162 }
3163 case PVValHelper::MAD: {
3164 float mad_ = PVValHelper::getMAD(residualsPlot[i]).value();
3165 float madErr_ = PVValHelper::getMAD(residualsPlot[i]).error();
3166 trendPlot->SetBinContent(i + 1, mad_);
3167 trendPlot->SetBinError(i + 1, madErr_);
3168 break;
3169 }
3170 default:
3171 edm::LogWarning("PrimaryVertexValidation")
3172 << "fillTrendPlot() " << fitPar_ << " unknown estimator!" << std::endl;
3173 break;
3174 }
3175
3176 if (var_.find("eta") != std::string::npos) {
3177 trendPlot->GetXaxis()->SetBinLabel(i + 1, etabincenter);
3178 } else if (var_.find("phi") != std::string::npos) {
3179 trendPlot->GetXaxis()->SetBinLabel(i + 1, phibincenter);
3180 } else {
3181 edm::LogWarning("PrimaryVertexValidation")
3182 << "fillTrendPlot() " << var_ << " unknown track parameter!" << std::endl;
3183 }
3184 }
3185 }
3186
3187
3188 void PrimaryVertexValidation::fillTrendPlotByIndex(TH1F* trendPlot,
3189 std::vector<TH1F*>& h,
3190 PVValHelper::estimator fitPar_,
3191 PVValHelper::plotVariable plotVar)
3192
3193 {
3194 for (auto iterator = h.begin(); iterator != h.end(); iterator++) {
3195 unsigned int bin = std::distance(h.begin(), iterator) + 1;
3196 std::pair<Measurement1D, Measurement1D> myFit = PVValHelper::fitResiduals((*iterator));
3197
3198 switch (fitPar_) {
3199 case PVValHelper::MEAN: {
3200 float mean_ = myFit.first.value();
3201 float meanErr_ = myFit.first.error();
3202 trendPlot->SetBinContent(bin, mean_);
3203 trendPlot->SetBinError(bin, meanErr_);
3204 break;
3205 }
3206 case PVValHelper::WIDTH: {
3207 float width_ = myFit.second.value();
3208 float widthErr_ = myFit.second.error();
3209 trendPlot->SetBinContent(bin, width_);
3210 trendPlot->SetBinError(bin, widthErr_);
3211 break;
3212 }
3213 case PVValHelper::MEDIAN: {
3214 float median_ = PVValHelper::getMedian(*iterator).value();
3215 float medianErr_ = PVValHelper::getMedian(*iterator).error();
3216 trendPlot->SetBinContent(bin, median_);
3217 trendPlot->SetBinError(bin, medianErr_);
3218 break;
3219 }
3220 case PVValHelper::MAD: {
3221 float mad_ = PVValHelper::getMAD(*iterator).value();
3222 float madErr_ = PVValHelper::getMAD(*iterator).error();
3223 trendPlot->SetBinContent(bin, mad_);
3224 trendPlot->SetBinError(bin, madErr_);
3225 break;
3226 }
3227 default:
3228 edm::LogWarning("PrimaryVertexValidation")
3229 << "fillTrendPlotByIndex() " << fitPar_ << " unknown estimator!" << std::endl;
3230 break;
3231 }
3232
3233 char bincenter[129];
3234 if (plotVar == PVValHelper::eta) {
3235 auto etaBins = theDetails_.trendbins[PVValHelper::eta];
3236 sprintf(bincenter, "%.1f", (etaBins[bin - 1] + etaBins[bin]) / 2.);
3237 trendPlot->GetXaxis()->SetBinLabel(bin, bincenter);
3238 } else if (plotVar == PVValHelper::phi) {
3239 auto phiBins = theDetails_.trendbins[PVValHelper::phi];
3240 sprintf(bincenter, "%.f", (phiBins[bin - 1] + phiBins[bin]) / 2.);
3241 trendPlot->GetXaxis()->SetBinLabel(bin, bincenter);
3242 } else {
3243
3244
3245 }
3246 }
3247 }
3248
3249
3250 void PrimaryVertexValidation::fillMap(TH2F* trendMap,
3251 TH1F* residualsMapPlot[100][100],
3252 PVValHelper::estimator fitPar_,
3253 const int nXBins_,
3254 const int nYBins_)
3255
3256 {
3257 for (int i = 0; i < nYBins_; ++i) {
3258 char phibincenter[129];
3259 auto phiBins = theDetails_.trendbins[PVValHelper::phi];
3260 sprintf(phibincenter, "%.f", (phiBins[i] + phiBins[i + 1]) / 2.);
3261
3262 if (nXBins_ == nYBins_) {
3263 trendMap->GetYaxis()->SetBinLabel(i + 1, phibincenter);
3264 }
3265
3266 for (int j = 0; j < nXBins_; ++j) {
3267 char etabincenter[129];
3268 auto etaBins = theDetails_.trendbins[PVValHelper::eta];
3269 sprintf(etabincenter, "%.1f", (etaBins[j] + etaBins[j + 1]) / 2.);
3270
3271 if (i == 0) {
3272 if (nXBins_ == nYBins_) {
3273 trendMap->GetXaxis()->SetBinLabel(j + 1, etabincenter);
3274 }
3275 }
3276
3277 switch (fitPar_) {
3278 case PVValHelper::MEAN: {
3279 float mean_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.value();
3280 float meanErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).first.error();
3281 trendMap->SetBinContent(j + 1, i + 1, mean_);
3282 trendMap->SetBinError(j + 1, i + 1, meanErr_);
3283 break;
3284 }
3285 case PVValHelper::WIDTH: {
3286 float width_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.value();
3287 float widthErr_ = PVValHelper::fitResiduals(residualsMapPlot[i][j]).second.error();
3288 trendMap->SetBinContent(j + 1, i + 1, width_);
3289 trendMap->SetBinError(j + 1, i + 1, widthErr_);
3290 break;
3291 }
3292 case PVValHelper::MEDIAN: {
3293 float median_ = PVValHelper::getMedian(residualsMapPlot[i][j]).value();
3294 float medianErr_ = PVValHelper::getMedian(residualsMapPlot[i][j]).error();
3295 trendMap->SetBinContent(j + 1, i + 1, median_);
3296 trendMap->SetBinError(j + 1, i + 1, medianErr_);
3297 break;
3298 }
3299 case PVValHelper::MAD: {
3300 float mad_ = PVValHelper::getMAD(residualsMapPlot[i][j]).value();
3301 float madErr_ = PVValHelper::getMAD(residualsMapPlot[i][j]).error();
3302 trendMap->SetBinContent(j + 1, i + 1, mad_);
3303 trendMap->SetBinError(j + 1, i + 1, madErr_);
3304 break;
3305 }
3306 default:
3307 edm::LogWarning("PrimaryVertexValidation:") << " fillMap() " << fitPar_ << " unknown estimator!" << std::endl;
3308 }
3309 }
3310 }
3311 }
3312
3313
3314 bool PrimaryVertexValidation::vtxSort(const reco::Vertex& a, const reco::Vertex& b)
3315
3316 {
3317 if (a.tracksSize() != b.tracksSize())
3318 return a.tracksSize() > b.tracksSize() ? true : false;
3319 else
3320 return a.chi2() < b.chi2() ? true : false;
3321 }
3322
3323
3324 bool PrimaryVertexValidation::passesTrackCuts(const reco::Track& track,
3325 const reco::Vertex& vertex,
3326 const std::string& qualityString_,
3327 double dxyErrMax_,
3328 double dzErrMax_,
3329 double ptErrMax_)
3330
3331 {
3332 math::XYZPoint vtxPoint(0.0, 0.0, 0.0);
3333 double vzErr = 0.0, vxErr = 0.0, vyErr = 0.0;
3334 vtxPoint = vertex.position();
3335 vzErr = vertex.zError();
3336 vxErr = vertex.xError();
3337 vyErr = vertex.yError();
3338
3339 double dxy = 0.0, dz = 0.0, dxysigma = 0.0, dzsigma = 0.0;
3340 dxy = track.dxy(vtxPoint);
3341 dz = track.dz(vtxPoint);
3342 dxysigma = sqrt(track.d0Error() * track.d0Error() + vxErr * vyErr);
3343 dzsigma = sqrt(track.dzError() * track.dzError() + vzErr * vzErr);
3344
3345 if (track.quality(reco::TrackBase::qualityByName(qualityString_)) != 1)
3346 return false;
3347 if (std::abs(dxy / dxysigma) > dxyErrMax_)
3348 return false;
3349 if (std::abs(dz / dzsigma) > dzErrMax_)
3350 return false;
3351 if (track.ptError() / track.pt() > ptErrMax_)
3352 return false;
3353
3354 return true;
3355 }
3356
3357
3358 std::map<std::string, TH1*> PrimaryVertexValidation::bookVertexHistograms(const TFileDirectory& dir)
3359
3360 {
3361 TH1F::SetDefaultSumw2(kTRUE);
3362
3363 std::map<std::string, TH1*> h;
3364
3365
3366 std::string types[] = {"all", "sel"};
3367 for (const auto& type : types) {
3368 h["pseudorapidity_" + type] =
3369 dir.make<TH1F>(("rapidity_" + type).c_str(), "track pseudorapidity; track #eta; tracks", 100, -3., 3.);
3370 h["z0_" + type] = dir.make<TH1F>(("z0_" + type).c_str(), "track z_{0};track z_{0} (cm);tracks", 80, -40., 40.);
3371 h["phi_" + type] = dir.make<TH1F>(("phi_" + type).c_str(), "track #phi; track #phi;tracks", 80, -M_PI, M_PI);
3372 h["eta_" + type] = dir.make<TH1F>(("eta_" + type).c_str(), "track #eta; track #eta;tracks", 80, -4., 4.);
3373 h["pt_" + type] = dir.make<TH1F>(("pt_" + type).c_str(), "track p_{T}; track p_{T} [GeV];tracks", 100, 0., 20.);
3374 h["p_" + type] = dir.make<TH1F>(("p_" + type).c_str(), "track p; track p [GeV];tracks", 100, 0., 20.);
3375 h["found_" + type] =
3376 dir.make<TH1F>(("found_" + type).c_str(), "n. found hits;n^{found}_{hits};tracks", 30, 0., 30.);
3377 h["lost_" + type] = dir.make<TH1F>(("lost_" + type).c_str(), "n. lost hits;n^{lost}_{hits};tracks", 20, 0., 20.);
3378 h["nchi2_" + type] =
3379 dir.make<TH1F>(("nchi2_" + type).c_str(), "normalized track #chi^{2};track #chi^{2}/ndf;tracks", 100, 0., 20.);
3380 h["rstart_" + type] = dir.make<TH1F>(
3381 ("rstart_" + type).c_str(), "track start radius; track innermost radius r (cm);tracks", 100, 0., 20.);
3382 h["expectedInner_" + type] = dir.make<TH1F>(
3383 ("expectedInner_" + type).c_str(), "n. expected inner hits;n^{expected}_{inner};tracks", 10, 0., 10.);
3384 h["expectedOuter_" + type] = dir.make<TH1F>(
3385 ("expectedOuter_" + type).c_str(), "n. expected outer hits;n^{expected}_{outer};tracks ", 10, 0., 10.);
3386 h["logtresxy_" + type] =
3387 dir.make<TH1F>(("logtresxy_" + type).c_str(),
3388 "log10(track r-#phi resolution/#mum);log10(track r-#phi resolution/#mum);tracks",
3389 100,
3390 0.,
3391 5.);
3392 h["logtresz_" + type] = dir.make<TH1F>(("logtresz_" + type).c_str(),
3393 "log10(track z resolution/#mum);log10(track z resolution/#mum);tracks",
3394 100,
3395 0.,
3396 5.);
3397 h["tpullxy_" + type] =
3398 dir.make<TH1F>(("tpullxy_" + type).c_str(), "track r-#phi pull;pull_{r-#phi};tracks", 100, -10., 10.);
3399 h["tpullz_" + type] =
3400 dir.make<TH1F>(("tpullz_" + type).c_str(), "track r-z pull;pull_{r-z};tracks", 100, -50., 50.);
3401 h["tlogDCAxy_" + type] = dir.make<TH1F>(
3402 ("tlogDCAxy_" + type).c_str(), "track log_{10}(DCA_{r-#phi});track log_{10}(DCA_{r-#phi});tracks", 200, -5., 3.);
3403 h["tlogDCAz_" + type] = dir.make<TH1F>(
3404 ("tlogDCAz_" + type).c_str(), "track log_{10}(DCA_{r-z});track log_{10}(DCA_{r-z});tracks", 200, -5., 5.);
3405 h["lvseta_" + type] = dir.make<TH2F>(
3406 ("lvseta_" + type).c_str(), "cluster length vs #eta;track #eta;cluster length", 60, -3., 3., 20, 0., 20);
3407 h["lvstanlambda_" + type] = dir.make<TH2F>(("lvstanlambda_" + type).c_str(),
3408 "cluster length vs tan #lambda; tan#lambda;cluster length",
3409 60,
3410 -6.,
3411 6.,
3412 20,
3413 0.,
3414 20);
3415 h["restrkz_" + type] =
3416 dir.make<TH1F>(("restrkz_" + type).c_str(), "z-residuals (track vs vertex);res_{z} (cm);tracks", 200, -5., 5.);
3417 h["restrkzvsphi_" + type] = dir.make<TH2F>(("restrkzvsphi_" + type).c_str(),
3418 "z-residuals (track - vertex) vs track #phi;track #phi;res_{z} (cm)",
3419 12,
3420 -M_PI,
3421 M_PI,
3422 100,
3423 -0.5,
3424 0.5);
3425 h["restrkzvseta_" + type] = dir.make<TH2F>(("restrkzvseta_" + type).c_str(),
3426 "z-residuals (track - vertex) vs track #eta;track #eta;res_{z} (cm)",
3427 12,
3428 -3.,
3429 3.,
3430 200,
3431 -0.5,
3432 0.5);
3433 h["pulltrkzvsphi_" + type] =
3434 dir.make<TH2F>(("pulltrkzvsphi_" + type).c_str(),
3435 "normalized z-residuals (track - vertex) vs track #phi;track #phi;res_{z}/#sigma_{res_{z}}",
3436 12,
3437 -M_PI,
3438 M_PI,
3439 100,
3440 -5.,
3441 5.);
3442 h["pulltrkzvseta_" + type] =
3443 dir.make<TH2F>(("pulltrkzvseta_" + type).c_str(),
3444 "normalized z-residuals (track - vertex) vs track #eta;track #eta;res_{z}/#sigma_{res_{z}}",
3445 12,
3446 -3.,
3447 3.,
3448 100,
3449 -5.,
3450 5.);
3451 h["pulltrkz_" + type] = dir.make<TH1F>(("pulltrkz_" + type).c_str(),
3452 "normalized z-residuals (track vs vertex);res_{z}/#sigma_{res_{z}};tracks",
3453 100,
3454 -5.,
3455 5.);
3456 h["sigmatrkz0_" + type] = dir.make<TH1F>(
3457 ("sigmatrkz0_" + type).c_str(), "z-resolution (excluding beam);#sigma^{trk}_{z_{0}} (cm);tracks", 100, 0., 5.);
3458 h["sigmatrkz_" + type] = dir.make<TH1F>(
3459 ("sigmatrkz_" + type).c_str(), "z-resolution (including beam);#sigma^{trk}_{z} (cm);tracks", 100, 0., 5.);
3460 h["nbarrelhits_" + type] = dir.make<TH1F>(
3461 ("nbarrelhits_" + type).c_str(), "number of pixel barrel hits;n. hits Barrel Pixel;tracks", 10, 0., 10.);
3462 h["nbarrelLayers_" + type] = dir.make<TH1F>(
3463 ("nbarrelLayers_" + type).c_str(), "number of pixel barrel layers;n. layers Barrel Pixel;tracks", 10, 0., 10.);
3464 h["nPxLayers_" + type] = dir.make<TH1F>(
3465 ("nPxLayers_" + type).c_str(), "number of pixel layers (barrel+endcap);n. Pixel layers;tracks", 10, 0., 10.);
3466 h["nSiLayers_" + type] =
3467 dir.make<TH1F>(("nSiLayers_" + type).c_str(), "number of Tracker layers;n. Tracker layers;tracks", 20, 0., 20.);
3468 h["trackAlgo_" + type] =
3469 dir.make<TH1F>(("trackAlgo_" + type).c_str(), "track algorithm;track algo;tracks", 30, 0., 30.);
3470 h["trackQuality_" + type] =
3471 dir.make<TH1F>(("trackQuality_" + type).c_str(), "track quality;track quality;tracks", 7, -1., 6.);
3472 }
3473
3474 return h;
3475 }
3476
3477
3478
3479
3480 std::vector<TH1F*> PrimaryVertexValidation::bookResidualsHistogram(const TFileDirectory& dir,
3481 unsigned int theNOfBins,
3482 PVValHelper::residualType resType,
3483 PVValHelper::plotVariable varType,
3484 bool isNormalized) {
3485 TH1F::SetDefaultSumw2(kTRUE);
3486
3487 auto hash = std::make_pair(resType, varType);
3488
3489 double down = theDetails_.range[hash].first;
3490 double up = theDetails_.range[hash].second;
3491
3492 if (isNormalized) {
3493 up = up / 100.;
3494 down = down / 100.;
3495 }
3496
3497 std::vector<TH1F*> h;
3498 h.reserve(theNOfBins);
3499
3500 if (theNOfBins == 0) {
3501 edm::LogError("PrimaryVertexValidation")
3502 << "bookResidualsHistogram() The number of bins cannot be identically 0" << std::endl;
3503 assert(false);
3504 }
3505
3506 std::string s_resType = std::get<0>(PVValHelper::getTypeString(resType));
3507 std::string s_varType = std::get<0>(PVValHelper::getVarString(varType));
3508
3509 std::string t_resType = std::get<1>(PVValHelper::getTypeString(resType));
3510 std::string t_varType = std::get<1>(PVValHelper::getVarString(varType));
3511 std::string units = std::get<2>(PVValHelper::getTypeString(resType));
3512
3513 for (unsigned int i = 0; i < theNOfBins; i++) {
3514 TString title = (varType == PVValHelper::phi || varType == PVValHelper::eta)
3515 ? Form("%s vs %s - bin %i (%f < %s < %f);%s %s;tracks",
3516 t_resType.c_str(),
3517 t_varType.c_str(),
3518 i,
3519 theDetails_.trendbins[varType][i],
3520 t_varType.c_str(),
3521 theDetails_.trendbins[varType][i + 1],
3522 t_resType.c_str(),
3523 units.c_str())
3524 : Form("%s vs %s - bin %i;%s %s;tracks",
3525 t_resType.c_str(),
3526 t_varType.c_str(),
3527 i,
3528 t_resType.c_str(),
3529 units.c_str());
3530
3531 TH1F* htemp = dir.make<TH1F>(
3532 Form("histo_%s_%s_plot%i", s_resType.c_str(), s_varType.c_str(), i),
3533
3534 title.Data(),
3535 theDetails_.histobins,
3536 down,
3537 up);
3538 h.push_back(htemp);
3539 }
3540
3541 return h;
3542 }
3543
3544
3545 void PrimaryVertexValidation::fillTrackHistos(std::map<std::string, TH1*>& h,
3546 const std::string& ttype,
3547 const reco::TransientTrack* tt,
3548 const reco::Vertex& v,
3549 const reco::BeamSpot& beamSpot,
3550 double fBfield_)
3551
3552 {
3553 using namespace reco;
3554
3555 PVValHelper::fill(h, "pseudorapidity_" + ttype, tt->track().eta());
3556 PVValHelper::fill(h, "z0_" + ttype, tt->track().vz());
3557 PVValHelper::fill(h, "phi_" + ttype, tt->track().phi());
3558 PVValHelper::fill(h, "eta_" + ttype, tt->track().eta());
3559 PVValHelper::fill(h, "pt_" + ttype, tt->track().pt());
3560 PVValHelper::fill(h, "p_" + ttype, tt->track().p());
3561 PVValHelper::fill(h, "found_" + ttype, tt->track().found());
3562 PVValHelper::fill(h, "lost_" + ttype, tt->track().lost());
3563 PVValHelper::fill(h, "nchi2_" + ttype, tt->track().normalizedChi2());
3564 PVValHelper::fill(h, "rstart_" + ttype, (tt->track().innerPosition()).Rho());
3565
3566 double d0Error = tt->track().d0Error();
3567 double d0 = tt->track().dxy(beamSpot.position());
3568 double dz = tt->track().dz(beamSpot.position());
3569 if (d0Error > 0) {
3570 PVValHelper::fill(h, "logtresxy_" + ttype, log(d0Error / 0.0001) / log(10.));
3571 PVValHelper::fill(h, "tpullxy_" + ttype, d0 / d0Error);
3572 PVValHelper::fill(h, "tlogDCAxy_" + ttype, log(std::abs(d0 / d0Error)));
3573 }
3574
3575 double dzError = tt->track().dzError();
3576 if (dzError > 0) {
3577 PVValHelper::fill(h, "logtresz_" + ttype, log(dzError / 0.0001) / log(10.));
3578 PVValHelper::fill(h, "tpullz_" + ttype, dz / dzError);
3579 PVValHelper::fill(h, "tlogDCAz_" + ttype, log(std::abs(dz / dzError)));
3580 }
3581
3582
3583 double wxy2_ = pow(beamSpot.BeamWidthX(), 2) + pow(beamSpot.BeamWidthY(), 2);
3584
3585 PVValHelper::fill(
3586 h, "sigmatrkz_" + ttype, sqrt(pow(tt->track().dzError(), 2) + wxy2_ / pow(tan(tt->track().theta()), 2)));
3587 PVValHelper::fill(h, "sigmatrkz0_" + ttype, tt->track().dzError());
3588
3589
3590 if (v.isValid()) {
3591
3592
3593
3594
3595 double z = (tt->stateAtBeamLine().trackStateAtPCA()).position().z();
3596 double tantheta = tan((tt->stateAtBeamLine().trackStateAtPCA()).momentum().theta());
3597 double dz2 = pow(tt->track().dzError(), 2) + wxy2_ / pow(tantheta, 2);
3598
3599 PVValHelper::fill(h, "restrkz_" + ttype, z - v.position().z());
3600 PVValHelper::fill(h, "restrkzvsphi_" + ttype, tt->track().phi(), z - v.position().z());
3601 PVValHelper::fill(h, "restrkzvseta_" + ttype, tt->track().eta(), z - v.position().z());
3602 PVValHelper::fill(h, "pulltrkzvsphi_" + ttype, tt->track().phi(), (z - v.position().z()) / sqrt(dz2));
3603 PVValHelper::fill(h, "pulltrkzvseta_" + ttype, tt->track().eta(), (z - v.position().z()) / sqrt(dz2));
3604
3605 PVValHelper::fill(h, "pulltrkz_" + ttype, (z - v.position().z()) / sqrt(dz2));
3606
3607 double x1 = tt->track().vx() - beamSpot.x0();
3608 double y1 = tt->track().vy() - beamSpot.y0();
3609
3610 double kappa = -0.002998 * fBfield_ * tt->track().qoverp() / cos(tt->track().theta());
3611 double D0 = x1 * sin(tt->track().phi()) - y1 * cos(tt->track().phi()) - 0.5 * kappa * (x1 * x1 + y1 * y1);
3612 double q = sqrt(1. - 2. * kappa * D0);
3613 double s0 = (x1 * cos(tt->track().phi()) + y1 * sin(tt->track().phi())) / q;
3614
3615 if (std::abs(kappa * s0) > 0.001) {
3616
3617 } else {
3618
3619
3620 }
3621
3622
3623 }
3624
3625
3626
3627 PVValHelper::fill(h, "nbarrelLayers_" + ttype, tt->track().hitPattern().pixelBarrelLayersWithMeasurement());
3628 PVValHelper::fill(h, "nPxLayers_" + ttype, tt->track().hitPattern().pixelLayersWithMeasurement());
3629 PVValHelper::fill(h, "nSiLayers_" + ttype, tt->track().hitPattern().trackerLayersWithMeasurement());
3630 PVValHelper::fill(
3631 h, "expectedInner_" + ttype, tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_INNER_HITS));
3632 PVValHelper::fill(
3633 h, "expectedOuter_" + ttype, tt->track().hitPattern().numberOfLostHits(HitPattern::MISSING_OUTER_HITS));
3634 PVValHelper::fill(h, "trackAlgo_" + ttype, tt->track().algo());
3635 PVValHelper::fill(h, "trackQuality_" + ttype, tt->track().qualityMask());
3636
3637
3638 int longesthit = 0, nbarrel = 0;
3639 for (auto const& hit : tt->track().recHits()) {
3640 if (hit->isValid() && hit->geographicalId().det() == DetId::Tracker) {
3641 bool barrel = DetId(hit->geographicalId()).subdetId() == static_cast<int>(PixelSubdetector::PixelBarrel);
3642
3643 if (barrel) {
3644 const SiPixelRecHit* pixhit = dynamic_cast<const SiPixelRecHit*>(&(*hit));
3645 edm::Ref<edmNew::DetSetVector<SiPixelCluster>, SiPixelCluster> const& clust = (*pixhit).cluster();
3646 if (clust.isNonnull()) {
3647 nbarrel++;
3648 if (clust->sizeY() - longesthit > 0)
3649 longesthit = clust->sizeY();
3650 if (clust->sizeY() > 20.) {
3651 PVValHelper::fill(h, "lvseta_" + ttype, tt->track().eta(), 19.9);
3652 PVValHelper::fill(h, "lvstanlambda_" + ttype, tan(tt->track().lambda()), 19.9);
3653 } else {
3654 PVValHelper::fill(h, "lvseta_" + ttype, tt->track().eta(), float(clust->sizeY()));
3655 PVValHelper::fill(h, "lvstanlambda_" + ttype, tan(tt->track().lambda()), float(clust->sizeY()));
3656 }
3657 }
3658 }
3659 }
3660 }
3661 PVValHelper::fill(h, "nbarrelhits_" + ttype, float(nbarrel));
3662
3663 }
3664
3665 void PrimaryVertexValidation::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
3666 edm::ParameterSetDescription desc;
3667 desc.setComment("Validates alignment payloads by evaluating unbiased track paramter resisuals to vertices");
3668
3669
3670
3671 desc.addUntracked<int>("compressionSettings", -1);
3672 desc.add<bool>("storeNtuple", false);
3673 desc.add<bool>("isLightNtuple", true);
3674 desc.add<bool>("useTracksFromRecoVtx", false);
3675 desc.addUntracked<double>("vertexZMax", 99);
3676 desc.addUntracked<double>("intLumi", 0.);
3677 desc.add<bool>("askFirstLayerHit", false);
3678 desc.addUntracked<bool>("doBPix", true);
3679 desc.addUntracked<bool>("doFPix", true);
3680 desc.addUntracked<double>("probePt", 0.);
3681 desc.addUntracked<double>("probeP", 0.);
3682 desc.addUntracked<double>("probeEta", 2.4);
3683 desc.addUntracked<double>("probeNHits", 0.);
3684 desc.addUntracked<int>("numberOfBins", 24);
3685 desc.addUntracked<double>("minPt", 1.);
3686 desc.addUntracked<double>("maxPt", 20.);
3687 desc.add<bool>("Debug", false);
3688 desc.addUntracked<bool>("runControl", false);
3689 desc.addUntracked<bool>("forceBeamSpot", false);
3690
3691 std::vector<unsigned int> defaultRuns;
3692 defaultRuns.push_back(0);
3693 desc.addUntracked<std::vector<unsigned int>>("runControlNumber", defaultRuns);
3694
3695
3696
3697 desc.add<edm::InputTag>("TrackCollectionTag", edm::InputTag("ALCARECOTkAlMinBias"));
3698 desc.add<edm::InputTag>("VertexCollectionTag", edm::InputTag("offlinePrimaryVertices"));
3699 desc.add<edm::InputTag>("BeamSpotTag", edm::InputTag("offlineBeamSpot"));
3700
3701
3702 edm::ParameterSetDescription psd0;
3703 TrackFilterForPVFinding::fillPSetDescription(psd0);
3704 psd0.add<int>("numTracksThreshold", 0);
3705 desc.add<edm::ParameterSetDescription>("TkFilterParameters", psd0);
3706
3707
3708 {
3709 edm::ParameterSetDescription psd0;
3710 {
3711 edm::ParameterSetDescription psd1;
3712 DAClusterizerInZ_vect::fillPSetDescription(psd1);
3713 psd0.add<edm::ParameterSetDescription>("TkDAClusParameters", psd1);
3714
3715 edm::ParameterSetDescription psd2;
3716 GapClusterizerInZ::fillPSetDescription(psd2);
3717 psd0.add<edm::ParameterSetDescription>("TkGapClusParameters", psd2);
3718 }
3719 psd0.add<std::string>("algorithm", "DA_vect");
3720 desc.add<edm::ParameterSetDescription>("TkClusParameters", psd0);
3721 }
3722
3723 descriptions.add("primaryVertexValidation", desc);
3724 }
3725
3726
3727 DEFINE_FWK_MODULE(PrimaryVertexValidation);
