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