diff --git a/Alignment/OfflineValidation/plugins/DiMuonValidation.cc b/Alignment/OfflineValidation/plugins/DiMuonValidation.cc index 9c37310347454..44ccf9d52ca28 100644 --- a/Alignment/OfflineValidation/plugins/DiMuonValidation.cc +++ b/Alignment/OfflineValidation/plugins/DiMuonValidation.cc @@ -155,6 +155,8 @@ class DiMuonValidation : public edm::one::EDAnalyzer int variable_PairPt_nbins_; edm::EDGetTokenT theTrackCollectionToken_; + + TH1F* th1f_mass; TH2D* th2d_mass_variables_[varNumber_]; // actual histograms std::string variables_name_[varNumber_] = { "CosThetaCS", "DeltaEta", "EtaMinus", "EtaPlus", "PhiCS", "PhiMinus", "PhiPlus", "Pt"}; @@ -175,37 +177,35 @@ void DiMuonValidation::analyze(const edm::Event& iEvent, const edm::EventSetup& const reco::TrackCollection& tC = iEvent.get(theTrackCollectionToken_); DiMuonValid::LV LV_mother(0., 0., 0., 0.); - //for (reco::TrackCollection::const_iterator track1 = tC.begin(); track1 != tC.end(); track1++) { - for (const auto& track1 : tC) { - DiMuonValid::LV LV_track1(track1.px(), - track1.py(), - track1.pz(), - sqrt((track1.p() * track1.p()) + mu_mass2_)); //old 106 - - for (const auto& track2 : tC) { - if (&track1 == &track2) { - continue; - } // discard the same track - - if (track1.charge() == track2.charge()) { + for (reco::TrackCollection::const_iterator track1 = tC.begin(); track1 != tC.end(); track1++) { + DiMuonValid::LV LV_track1(track1->px(), + track1->py(), + track1->pz(), + sqrt((track1->p() * track1->p()) + mu_mass2_)); //old 106 + + for (reco::TrackCollection::const_iterator track2 = track1 + 1; track2 != tC.end(); track2++) { + if (track1->charge() == track2->charge()) { continue; } // only reconstruct opposite charge pair - DiMuonValid::LV LV_track2(track2.px(), track2.py(), track2.pz(), sqrt((track2.p() * track2.p()) + mu_mass2_)); + DiMuonValid::LV LV_track2( + track2->px(), track2->py(), track2->pz(), sqrt((track2->p() * track2->p()) + mu_mass2_)); LV_mother = LV_track1 + LV_track2; double mother_mass = LV_mother.M(); + th1f_mass->Fill(mother_mass); + double mother_pt = LV_mother.Pt(); - int charge1 = track1.charge(); - double etaMu1 = track1.eta(); - double phiMu1 = track1.phi(); - double ptMu1 = track1.pt(); + int charge1 = track1->charge(); + double etaMu1 = track1->eta(); + double phiMu1 = track1->phi(); + double ptMu1 = track1->pt(); - int charge2 = track2.charge(); - double etaMu2 = track2.eta(); - double phiMu2 = track2.phi(); - double ptMu2 = track2.pt(); + int charge2 = track2->charge(); + double etaMu2 = track2->eta(); + double phiMu2 = track2->phi(); + double ptMu2 = track2->pt(); if (charge1 < 0) { // use Mu+ for charge1, Mu- for charge2 std::swap(charge1, charge2); @@ -263,6 +263,8 @@ void DiMuonValidation::beginJob() { fs->file().SetCompressionSettings(compressionSettings_); } + th1f_mass = fs->make("hMass", "mass;m_{#mu#mu} [GeV];events", 200, 0., 200.); + for (int i = 0; i < varNumber_; i++) { std::string th2d_name = fmt::sprintf("th2d_mass_%s", variables_name_[i].c_str()); th2d_mass_variables_[i] = fs->make(th2d_name.c_str(), @@ -288,10 +290,10 @@ void DiMuonValidation::fillDescriptions(edm::ConfigurationDescriptions& descript desc.add("Pair_mass_min", 60); desc.add("Pair_mass_max", 120); desc.add("Pair_mass_nbins", 120); - desc.add("Pair_etaminpos", 60); - desc.add("Pair_etamaxpos", 60); - desc.add("Pair_etaminneg", 60); - desc.add("Pair_etamaxneg", 60); + desc.add("Pair_etaminpos", -2.4); + desc.add("Pair_etamaxpos", 2.4); + desc.add("Pair_etaminneg", -2.4); + desc.add("Pair_etamaxneg", 2.4); desc.add("Variable_CosThetaCS_xmin", -1.); desc.add("Variable_CosThetaCS_xmax", 1.); diff --git a/Alignment/OfflineValidation/plugins/GeneralPurposeTrackAnalyzer.cc b/Alignment/OfflineValidation/plugins/GeneralPurposeTrackAnalyzer.cc index 195a354c2c45e..550e8d739ae64 100644 --- a/Alignment/OfflineValidation/plugins/GeneralPurposeTrackAnalyzer.cc +++ b/Alignment/OfflineValidation/plugins/GeneralPurposeTrackAnalyzer.cc @@ -274,6 +274,15 @@ class GeneralPurposeTrackAnalyzer : public edm::one::EDAnalyzer vTrackHistos_; std::vector vTrackProfiles_; std::vector vTrack2DHistos_; @@ -646,12 +655,31 @@ class GeneralPurposeTrackAnalyzer : public edm::one::EDAnalyzer beamSpotHandle = event.getHandle(beamspotToken_); if (beamSpotHandle.isValid()) { beamSpot = *beamSpotHandle; - math::XYZPoint point(beamSpot.x0(), beamSpot.y0(), beamSpot.z0()); + + double BSx0 = beamSpot.x0(); + double BSy0 = beamSpot.y0(); + double BSz0 = beamSpot.z0(); + double Beamsigmaz = beamSpot.sigmaZ(); + double Beamdxdz = beamSpot.dxdz(); + double Beamdydz = beamSpot.dydz(); + double BeamWidthX = beamSpot.BeamWidthX(); + double BeamWidthY = beamSpot.BeamWidthY(); + + math::XYZPoint point(BSx0, BSy0, BSz0); double dxy = track->dxy(point); double dz = track->dz(point); hdxyBS->Fill(dxy); hd0BS->Fill(-dxy); hdzBS->Fill(dz); + + h_BSx0->Fill(BSx0); + h_BSy0->Fill(BSy0); + h_BSz0->Fill(BSz0); + h_Beamsigmaz->Fill(Beamsigmaz); + h_BeamWidthX->Fill(BeamWidthX); + h_BeamWidthY->Fill(BeamWidthY); + h_BSdxdz->Fill(Beamdxdz); + h_BSdydz->Fill(Beamdydz); } //dxy with respect to the primary vertex @@ -843,6 +871,15 @@ class GeneralPurposeTrackAnalyzer : public edm::one::EDAnalyzer("h_PhiEndcapPlus", "hPhiEndcapPlus (#eta>1.4);track #phi;track", 100, -M_PI, M_PI); hPhiEndcapMinus = book("h_PhiEndcapMinus", "hPhiEndcapMinus (#eta<1.4);track #phi;tracks", 100, -M_PI, M_PI); + h_BSx0 = book("h_BSx0", "x-coordinate of reco beamspot;x^{BS}_{0};n_{events}", 100, -0.1, 0.1); + h_BSy0 = book("h_BSy0", "y-coordinate of reco beamspot;y^{BS}_{0};n_{events}", 100, -0.1, 0.1); + h_BSz0 = book("h_BSz0", "z-coordinate of reco beamspot;z^{BS}_{0};n_{events}", 100, -1., 1.); + h_Beamsigmaz = book("h_Beamsigmaz", "z-coordinate beam width;#sigma_{Z}^{beam};n_{events}", 100, 0., 1.); + h_BeamWidthX = book("h_BeamWidthX", "x-coordinate beam width;#sigma_{X}^{beam};n_{events}", 100, 0., 0.01); + h_BeamWidthY = book("h_BeamWidthY", "y-coordinate beam width;#sigma_{Y}^{beam};n_{events}", 100, 0., 0.01); + h_BSdxdz = book("h_BSdxdz", "BeamSpot dxdz;beamspot dx/dz;n_{events}", 100, -0.0003, 0.0003); + h_BSdydz = book("h_BSdydz", "BeamSpot dydz;beamspot dy/dz;n_{events}", 100, -0.0003, 0.0003); + if (!isCosmics_) { hPhp = book("h_P_hp", "Momentum (high purity);track momentum [GeV];tracks", 100, 0., 100.); hPthp = book("h_Pt_hp", "Transverse Momentum (high purity);track p_{T} [GeV];tracks", 100, 0., 100.); diff --git a/Alignment/OfflineValidation/plugins/GeneralPurposeVertexAnalyzer.cc b/Alignment/OfflineValidation/plugins/GeneralPurposeVertexAnalyzer.cc new file mode 100644 index 0000000000000..0dfe9fa546114 --- /dev/null +++ b/Alignment/OfflineValidation/plugins/GeneralPurposeVertexAnalyzer.cc @@ -0,0 +1,698 @@ +// -*- C++ -*- +// +// Package: Alignment/GeneralPurposeVertexAnalyzer +// Class: GeneralPurposeVertexAnalyzer +// +/**\class GeneralPurposeVertexAnalyzer GeneralPurposeVertexAnalyzer.cc Alignment/GeneralPurposeVertexAnalyzer/plugins/GeneralPurposeVertexAnalyzer.cc + Description: monitor vertex properties for alignment purposes, largely copied from DQMOffline/RecoB/plugins/PrimaryVertexMonitor.cc + +*/ +// +// Original Author: Marco Musich +// Created: Thu, 13 Apr 2023 14:16:43 GMT +// +// + +// ROOT includes files +#include "TMath.h" +#include "TFile.h" +#include "TH1I.h" +#include "TH1D.h" +#include "TH2D.h" +#include "TProfile.h" +#include "TProfile2D.h" + +// system include files +#include +#include +#include + +// user include files +#include "CommonTools/Statistics/interface/ChiSquaredProbability.h" +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/Common/interface/Association.h" +#include "DataFormats/Common/interface/ValueMap.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "FWCore/ServiceRegistry/interface/Service.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/InputTag.h" +#include "FWCore/Utilities/interface/isFinite.h" + +// +// class declaration +// + +using reco::TrackCollection; + +class GeneralPurposeVertexAnalyzer : public edm::one::EDAnalyzer { +public: + explicit GeneralPurposeVertexAnalyzer(const edm::ParameterSet &); + ~GeneralPurposeVertexAnalyzer() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions &descriptions); + +private: + void pvTracksPlots(const reco::Vertex &v); + void vertexPlots(const reco::Vertex &v, const reco::BeamSpot &beamSpot, int i); + template + T *book(const Args &...args) const; + void beginJob() override; + void analyze(const edm::Event &, const edm::EventSetup &) override; + + // ----------member data --------------------------- + edm::Service fs_; + + const int ndof_; + bool errorPrinted_; + + const edm::InputTag vertexInputTag_, beamSpotInputTag_; + const edm::EDGetTokenT vertexToken_; + using VertexScore = edm::ValueMap; + const edm::EDGetTokenT scoreToken_; + const edm::EDGetTokenT beamspotToken_; + + static constexpr int cmToUm = 10000; + + const double vposx_; + const double vposy_; + const int tkNoBin_; + const double tkNoMin_; + const double tkNoMax_; + + const int dxyBin_; + const double dxyMin_; + const double dxyMax_; + + const int dzBin_; + const double dzMin_; + const double dzMax_; + + const int phiBin_; + const int phiBin2D_; + const double phiMin_; + const double phiMax_; + + const int etaBin_; + const int etaBin2D_; + const double etaMin_; + const double etaMax_; + + // the histos + TH1I *nbvtx, *nbgvtx; + TH1D *nbtksinvtx[2], *trksWeight[2], *score[2]; + TH1D *tt[2]; + TH1D *xrec[2], *yrec[2], *zrec[2], *xDiff[2], *yDiff[2], *xerr[2], *yerr[2], *zerr[2]; + TH2D *xerrVsTrks[2], *yerrVsTrks[2], *zerrVsTrks[2]; + TH1D *ntracksVsZ[2]; + TH1D *vtxchi2[2], *vtxndf[2], *vtxprob[2], *nans[2]; + TH1D *type[2]; + TH1D *bsX, *bsY, *bsZ, *bsSigmaZ, *bsDxdz, *bsDydz, *bsBeamWidthX, *bsBeamWidthY, *bsType; + + TH1D *sumpt, *ntracks, *weight, *chi2ndf, *chi2prob; + TH1D *dxy, *dxy2, *dz, *dxyErr, *dzErr; + TH1D *phi_pt1, *eta_pt1; + TH1D *phi_pt10, *eta_pt10; + TProfile *dxyVsPhi_pt1, *dzVsPhi_pt1; + TProfile *dxyVsEta_pt1, *dzVsEta_pt1; + TProfile2D *dxyVsEtaVsPhi_pt1, *dzVsEtaVsPhi_pt1; + TProfile *dxyVsPhi_pt10, *dzVsPhi_pt10; + TProfile *dxyVsEta_pt10, *dzVsEta_pt10; + TProfile2D *dxyVsEtaVsPhi_pt10, *dzVsEtaVsPhi_pt10; +}; + +// +// constructors and destructor +// +GeneralPurposeVertexAnalyzer::GeneralPurposeVertexAnalyzer(const edm::ParameterSet &iConfig) + : ndof_(iConfig.getParameter("ndof")), + errorPrinted_(false), + vertexInputTag_(iConfig.getParameter("vertexLabel")), + beamSpotInputTag_(iConfig.getParameter("beamSpotLabel")), + vertexToken_(consumes(vertexInputTag_)), + scoreToken_(consumes(vertexInputTag_)), + beamspotToken_(consumes(beamSpotInputTag_)), + // to be configured for each year... + vposx_(iConfig.getParameter("Xpos")), + vposy_(iConfig.getParameter("Ypos")), + tkNoBin_(iConfig.getParameter("TkSizeBin")), + tkNoMin_(iConfig.getParameter("TkSizeMin")), + tkNoMax_(iConfig.getParameter("TkSizeMax")), + dxyBin_(iConfig.getParameter("DxyBin")), + dxyMin_(iConfig.getParameter("DxyMin")), + dxyMax_(iConfig.getParameter("DxyMax")), + dzBin_(iConfig.getParameter("DzBin")), + dzMin_(iConfig.getParameter("DzMin")), + dzMax_(iConfig.getParameter("DzMax")), + phiBin_(iConfig.getParameter("PhiBin")), + phiBin2D_(iConfig.getParameter("PhiBin2D")), + phiMin_(iConfig.getParameter("PhiMin")), + phiMax_(iConfig.getParameter("PhiMax")), + etaBin_(iConfig.getParameter("EtaBin")), + etaBin2D_(iConfig.getParameter("EtaBin2D")), + etaMin_(iConfig.getParameter("EtaMin")), + etaMax_(iConfig.getParameter("EtaMax")), + // histograms + nbvtx(nullptr), + bsX(nullptr), + bsY(nullptr), + bsZ(nullptr), + bsSigmaZ(nullptr), + bsDxdz(nullptr), + bsDydz(nullptr), + bsBeamWidthX(nullptr), + bsBeamWidthY(nullptr), + bsType(nullptr), + sumpt(nullptr), + ntracks(nullptr), + weight(nullptr), + chi2ndf(nullptr), + chi2prob(nullptr), + dxy(nullptr), + dxy2(nullptr), + dz(nullptr), + dxyErr(nullptr), + dzErr(nullptr), + phi_pt1(nullptr), + eta_pt1(nullptr), + phi_pt10(nullptr), + eta_pt10(nullptr), + dxyVsPhi_pt1(nullptr), + dzVsPhi_pt1(nullptr), + dxyVsEta_pt1(nullptr), + dzVsEta_pt1(nullptr), + dxyVsEtaVsPhi_pt1(nullptr), + dzVsEtaVsPhi_pt1(nullptr), + dxyVsPhi_pt10(nullptr), + dzVsPhi_pt10(nullptr), + dxyVsEta_pt10(nullptr), + dzVsEta_pt10(nullptr), + dxyVsEtaVsPhi_pt10(nullptr), + dzVsEtaVsPhi_pt10(nullptr) { + usesResource(TFileService::kSharedResource); +} + +// +// member functions +// + +// ------------ method called for each event ------------ +void GeneralPurposeVertexAnalyzer::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) { + using namespace edm; + + const auto &recVtxs = iEvent.getHandle(vertexToken_); + const auto &scores = iEvent.getHandle(scoreToken_); + const auto &beamSpotHandle = iEvent.getHandle(beamspotToken_); + reco::BeamSpot beamSpot = *beamSpotHandle; + + // check for absent products and simply "return" in that case + if (recVtxs.isValid() == false || beamSpotHandle.isValid() == false) { + edm::LogWarning("GeneralPurposeVertexAnalyzer") + << " Some products not available in the event: VertexCollection " << vertexInputTag_ << " " << recVtxs.isValid() + << " BeamSpot " << beamSpotInputTag_ << " " << beamSpotHandle.isValid() << ". Skipping plots for this event"; + return; + } + + // check upfront that refs to track are (likely) to be valid + bool ok{true}; + for (const auto &v : *recVtxs) { + if (v.tracksSize() > 0) { + const auto &ref = v.trackRefAt(0); + if (ref.isNull() || !ref.isAvailable()) { + if (!errorPrinted_) { + edm::LogWarning("GeneralPurposeVertexAnalyzer") + << "Skipping vertex collection: " << vertexInputTag_ + << " since likely the track collection the vertex has refs pointing to is missing (at least the first " + "TrackBaseRef is null or not available)"; + } else { + errorPrinted_ = true; + } + ok = false; + } + } + } + + if (!ok) { + return; + } + + nbvtx->Fill(recVtxs->size()); + int ng = 0; + for (auto const &vx : (*recVtxs)) { + if (vx.isValid() && !vx.isFake() && vx.ndof() >= ndof_) { + ++ng; + } + } + nbgvtx->Fill(ng); + + if (scores.isValid() && !(*scores).empty()) { + auto pvScore = (*scores).get(0); + score[1]->Fill(std::sqrt(pvScore)); + for (unsigned int i = 1; i < (*scores).size(); ++i) { + score[0]->Fill(std::sqrt((*scores).get(i))); + } + } + + // fill PV tracks MEs (as now, for alignment) + if (!recVtxs->empty()) { + vertexPlots(recVtxs->front(), beamSpot, 1); + pvTracksPlots(recVtxs->front()); + + for (reco::VertexCollection::const_iterator v = recVtxs->begin() + 1; v != recVtxs->end(); ++v) { + vertexPlots(*v, beamSpot, 0); + } + } + + // Beamline plots: + bsX->Fill(beamSpot.x0()); + bsY->Fill(beamSpot.y0()); + bsZ->Fill(beamSpot.z0()); + bsSigmaZ->Fill(beamSpot.sigmaZ()); + bsDxdz->Fill(beamSpot.dxdz()); + bsDydz->Fill(beamSpot.dydz()); + bsBeamWidthX->Fill(beamSpot.BeamWidthX() * cmToUm); + bsBeamWidthY->Fill(beamSpot.BeamWidthY() * cmToUm); + bsType->Fill(beamSpot.type()); +} + +void GeneralPurposeVertexAnalyzer::pvTracksPlots(const reco::Vertex &v) { + if (!v.isValid()) + return; + if (v.isFake()) + return; + + if (v.tracksSize() == 0) { + ntracks->Fill(0); + return; + } + + const math::XYZPoint myVertex(v.position().x(), v.position().y(), v.position().z()); + + float sumPT = 0.; + for (const auto &t : v.tracks()) { + const bool isHighPurity = t->quality(reco::TrackBase::highPurity); + if (!isHighPurity) { + continue; + } + + const float pt = t->pt(); + if (pt < 1.f) { + continue; + } + + const float pt2 = pt * pt; + const float eta = t->eta(); + const float phi = t->phi(); + const float w = v.trackWeight(t); + const float chi2NDF = t->normalizedChi2(); + const float chi2Prob = TMath::Prob(t->chi2(), static_cast(t->ndof())); + const float Dxy = t->dxy(myVertex) * cmToUm; + const float Dz = t->dz(myVertex) * cmToUm; + const float DxyErr = t->dxyError() * cmToUm; + const float DzErr = t->dzError() * cmToUm; + + sumPT += pt2; + + weight->Fill(w); + chi2ndf->Fill(chi2NDF); + chi2prob->Fill(chi2Prob); + dxy->Fill(Dxy); + dxy2->Fill(Dxy); + dz->Fill(Dz); + dxyErr->Fill(DxyErr); + dzErr->Fill(DzErr); + phi_pt1->Fill(phi); + eta_pt1->Fill(eta); + dxyVsPhi_pt1->Fill(phi, Dxy); + dzVsPhi_pt1->Fill(phi, Dz); + dxyVsEta_pt1->Fill(eta, Dxy); + dzVsEta_pt1->Fill(eta, Dz); + dxyVsEtaVsPhi_pt1->Fill(eta, phi, Dxy); + dzVsEtaVsPhi_pt1->Fill(eta, phi, Dz); + + if (pt >= 10.f) { + phi_pt10->Fill(phi); + eta_pt10->Fill(eta); + dxyVsPhi_pt10->Fill(phi, Dxy); + dzVsPhi_pt10->Fill(phi, Dz); + dxyVsEta_pt10->Fill(eta, Dxy); + dzVsEta_pt10->Fill(eta, Dz); + dxyVsEtaVsPhi_pt10->Fill(eta, phi, Dxy); + dzVsEtaVsPhi_pt10->Fill(eta, phi, Dz); + } + } + + ntracks->Fill(static_cast(v.tracks().size())); + sumpt->Fill(sumPT); +} + +void GeneralPurposeVertexAnalyzer::vertexPlots(const reco::Vertex &v, const reco::BeamSpot &beamSpot, int i) { + if (i < 0 || i > 1) + return; + if (!v.isValid()) + type[i]->Fill(2.); + else if (v.isFake()) + type[i]->Fill(1.); + else + type[i]->Fill(0.); + + if (v.isValid() && !v.isFake()) { + float weight = 0; + for (reco::Vertex::trackRef_iterator t = v.tracks_begin(); t != v.tracks_end(); t++) + weight += v.trackWeight(*t); + trksWeight[i]->Fill(weight); + nbtksinvtx[i]->Fill(v.tracksSize()); + ntracksVsZ[i]->Fill(v.position().z() - beamSpot.z0(), v.tracksSize()); + + vtxchi2[i]->Fill(v.chi2()); + vtxndf[i]->Fill(v.ndof()); + vtxprob[i]->Fill(ChiSquaredProbability(v.chi2(), v.ndof())); + + xrec[i]->Fill(v.position().x()); + yrec[i]->Fill(v.position().y()); + zrec[i]->Fill(v.position().z()); + + float xb = beamSpot.x0() + beamSpot.dxdz() * (v.position().z() - beamSpot.z0()); + float yb = beamSpot.y0() + beamSpot.dydz() * (v.position().z() - beamSpot.z0()); + xDiff[i]->Fill((v.position().x() - xb) * cmToUm); + yDiff[i]->Fill((v.position().y() - yb) * cmToUm); + + xerr[i]->Fill(v.xError() * cmToUm); + yerr[i]->Fill(v.yError() * cmToUm); + zerr[i]->Fill(v.zError() * cmToUm); + xerrVsTrks[i]->Fill(weight, v.xError() * cmToUm); + yerrVsTrks[i]->Fill(weight, v.yError() * cmToUm); + zerrVsTrks[i]->Fill(weight, v.zError() * cmToUm); + + nans[i]->Fill(1., edm::isNotFinite(v.position().x()) * 1.); + nans[i]->Fill(2., edm::isNotFinite(v.position().y()) * 1.); + nans[i]->Fill(3., edm::isNotFinite(v.position().z()) * 1.); + + int index = 3; + for (int k = 0; k != 3; k++) { + for (int j = k; j != 3; j++) { + index++; + nans[i]->Fill(index * 1., edm::isNotFinite(v.covariance(k, j)) * 1.); + // in addition, diagonal element must be positive + if (j == k && v.covariance(k, j) < 0) { + nans[i]->Fill(index * 1., 1.); + } + } + } + } +} + +template +T *GeneralPurposeVertexAnalyzer::book(const Args &...args) const { + T *t = fs_->make(args...); + return t; +} + +// ------------ method called once each job just before starting event loop ------------ +void GeneralPurposeVertexAnalyzer::beginJob() { + nbvtx = book("vtxNbr", "Reconstructed Vertices in Event", 80, -0.5, 79.5); + nbgvtx = book("goodvtxNbr", "Reconstructed Good Vertices in Event", 80, -0.5, 79.5); + + nbtksinvtx[0] = book("otherVtxTrksNbr", "Reconstructed Tracks in Vertex (other Vtx)", 40, -0.5, 99.5); + ntracksVsZ[0] = + book("otherVtxTrksVsZ", "Reconstructed Tracks in Vertex (other Vtx) vs Z", 80, -20., 20., 0., 100., ""); + ntracksVsZ[0]->SetXTitle("z-bs"); + ntracksVsZ[0]->SetYTitle("#tracks"); + + score[0] = book("otherVtxScore", "sqrt(score) (other Vtx)", 100, 0., 400.); + trksWeight[0] = book("otherVtxTrksWeight", "Total weight of Tracks in Vertex (other Vtx)", 40, 0, 100.); + vtxchi2[0] = book("otherVtxChi2", "#chi^{2} (other Vtx)", 100, 0., 200.); + vtxndf[0] = book("otherVtxNdf", "ndof (other Vtx)", 100, 0., 200.); + vtxprob[0] = book("otherVtxProb", "#chi^{2} probability (other Vtx)", 100, 0., 1.); + nans[0] = book("otherVtxNans", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz (other Vtx)", 9, 0.5, 9.5); + + nbtksinvtx[1] = book("tagVtxTrksNbr", "Reconstructed Tracks in Vertex (tagged Vtx)", 100, -0.5, 99.5); + ntracksVsZ[1] = + book("tagVtxTrksVsZ", "Reconstructed Tracks in Vertex (tagged Vtx) vs Z", 80, -20., 20., 0., 100., ""); + ntracksVsZ[1]->SetXTitle("z-bs"); + ntracksVsZ[1]->SetYTitle("#tracks"); + + score[1] = book("tagVtxScore", "sqrt(score) (tagged Vtx)", 100, 0., 400.); + trksWeight[1] = book("tagVtxTrksWeight", "Total weight of Tracks in Vertex (tagged Vtx)", 100, 0, 100.); + vtxchi2[1] = book("tagVtxChi2", "#chi^{2} (tagged Vtx)", 100, 0., 200.); + vtxndf[1] = book("tagVtxNdf", "ndof (tagged Vtx)", 100, 0., 200.); + vtxprob[1] = book("tagVtxProb", "#chi^{2} probability (tagged Vtx)", 100, 0., 1.); + nans[1] = book("tagVtxNans", "Illegal values for x,y,z,xx,xy,xz,yy,yz,zz (tagged Vtx)", 9, 0.5, 9.5); + + xrec[0] = book("otherPosX", "Position x Coordinate (other Vtx)", 100, vposx_ - 0.1, vposx_ + 0.1); + yrec[0] = book("otherPosY", "Position y Coordinate (other Vtx)", 100, vposy_ - 0.1, vposy_ + 0.1); + zrec[0] = book("otherPosZ", "Position z Coordinate (other Vtx)", 100, -20., 20.); + xDiff[0] = book("otherDiffX", "X distance from BeamSpot (other Vtx)", 100, -500, 500); + yDiff[0] = book("otherDiffY", "Y distance from BeamSpot (other Vtx)", 100, -500, 500); + xerr[0] = book("otherErrX", "Uncertainty x Coordinate (other Vtx)", 100, 0., 100); + yerr[0] = book("otherErrY", "Uncertainty y Coordinate (other Vtx)", 100, 0., 100); + zerr[0] = book("otherErrZ", "Uncertainty z Coordinate (other Vtx)", 100, 0., 100); + xerrVsTrks[0] = book( + "otherErrVsWeightX", "Uncertainty x Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100); + yerrVsTrks[0] = book( + "otherErrVsWeightY", "Uncertainty y Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100); + zerrVsTrks[0] = book( + "otherErrVsWeightZ", "Uncertainty z Coordinate vs. track weight (other Vtx)", 100, 0, 100., 100, 0., 100); + + xrec[1] = book("tagPosX", "Position x Coordinate (tagged Vtx)", 100, vposx_ - 0.1, vposx_ + 0.1); + yrec[1] = book("tagPosY", "Position y Coordinate (tagged Vtx)", 100, vposy_ - 0.1, vposy_ + 0.1); + zrec[1] = book("tagPosZ", "Position z Coordinate (tagged Vtx)", 100, -20., 20.); + xDiff[1] = book("tagDiffX", "X distance from BeamSpot (tagged Vtx)", 100, -500, 500); + yDiff[1] = book("tagDiffY", "Y distance from BeamSpot (tagged Vtx)", 100, -500, 500); + xerr[1] = book("tagErrX", "Uncertainty x Coordinate (tagged Vtx)", 100, 0., 100); + yerr[1] = book("tagErrY", "Uncertainty y Coordinate (tagged Vtx)", 100, 0., 100); + zerr[1] = book("tagErrZ", "Uncertainty z Coordinate (tagged Vtx)", 100, 0., 100); + xerrVsTrks[1] = book( + "tagErrVsWeightX", "Uncertainty x Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100); + yerrVsTrks[1] = book( + "tagErrVsWeightY", "Uncertainty y Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100); + zerrVsTrks[1] = book( + "tagErrVsWeightZ", "Uncertainty z Coordinate vs. track weight (tagged Vtx)", 100, 0, 100., 100, 0., 100); + + type[0] = book("otherType", "Vertex type (other Vtx)", 3, -0.5, 2.5); + type[1] = book("tagType", "Vertex type (tagged Vtx)", 3, -0.5, 2.5); + + for (int i = 0; i < 2; ++i) { + type[i]->GetXaxis()->SetBinLabel(1, "Valid, real"); + type[i]->GetXaxis()->SetBinLabel(2, "Valid, fake"); + type[i]->GetXaxis()->SetBinLabel(3, "Invalid"); + } + + bsX = book("bsX", "BeamSpot x0", 100, vposx_ - 0.1, vposx_ + 0.1); + bsY = book("bsY", "BeamSpot y0", 100, vposy_ - 0.1, vposy_ + 0.1); + bsZ = book("bsZ", "BeamSpot z0", 100, -2., 2.); + bsSigmaZ = book("bsSigmaZ", "BeamSpot sigmaZ", 100, 0., 10.); + bsDxdz = book("bsDxdz", "BeamSpot dxdz", 100, -0.0003, 0.0003); + bsDydz = book("bsDydz", "BeamSpot dydz", 100, -0.0003, 0.0003); + bsBeamWidthX = book("bsBeamWidthX", "BeamSpot BeamWidthX", 100, 0., 100.); + bsBeamWidthY = book("bsBeamWidthY", "BeamSpot BeamWidthY", 100, 0., 100.); + bsType = book("bsType", "BeamSpot type", 4, -1.5, 2.5); + bsType->GetXaxis()->SetBinLabel(1, "Unknown"); + bsType->GetXaxis()->SetBinLabel(2, "Fake"); + bsType->GetXaxis()->SetBinLabel(3, "LHC"); + bsType->GetXaxis()->SetBinLabel(4, "Tracker"); + + // repeated strings in titles + std::string s_1 = "PV Tracks (p_{T} > 1 GeV)"; + std::string s_10 = "PV Tracks (p_{T} > 10 GeV)"; + + ntracks = book("ntracks", fmt::sprintf("number of %s", s_1).c_str(), tkNoBin_, tkNoMin_, tkNoMax_); + ntracks->SetXTitle(fmt::sprintf("Number of %s per Event", s_1).c_str()); + ntracks->SetYTitle("Number of Events"); + + weight = book("weight", fmt::sprintf("weight of %s", s_1).c_str(), 100, 0., 1.); + weight->SetXTitle(fmt::sprintf("weight of %s per Event", s_1).c_str()); + weight->SetYTitle("Number of Events"); + + sumpt = book("sumpt", fmt::sprintf("#Sum p_{T} of %s", s_1).c_str(), 100, -0.5, 249.5); + chi2ndf = book("chi2ndf", fmt::sprintf("%s #chi^{2}/ndof", s_1).c_str(), 100, 0., 20.); + chi2prob = book("chi2prob", fmt::sprintf("%s #chi^{2} probability", s_1).c_str(), 100, 0., 1.); + dxy = book("dxy", fmt::sprintf("%s d_{xy} (#mum)", s_1).c_str(), dxyBin_, dxyMin_, dxyMax_); + dxy2 = book("dxyzoom", fmt::sprintf("%s d_{xy} (#mum)", s_1).c_str(), dxyBin_, dxyMin_ / 5., dxyMax_ / 5.); + dxyErr = book("dxyErr", fmt::sprintf("%s d_{xy} error (#mum)", s_1).c_str(), 100, 0., 2000.); + dz = book("dz", fmt::sprintf("%s d_{z} (#mum)", s_1).c_str(), dzBin_, dzMin_, dzMax_); + dzErr = book("dzErr", fmt::sprintf("%s d_{z} error(#mum)", s_1).c_str(), 100, 0., 10000.); + + phi_pt1 = + book("phi_pt1", fmt::sprintf("%s #phi; PV tracks #phi;#tracks", s_1).c_str(), phiBin_, phiMin_, phiMax_); + eta_pt1 = + book("eta_pt1", fmt::sprintf("%s #eta; PV tracks #eta;#tracks", s_1).c_str(), etaBin_, etaMin_, etaMax_); + phi_pt10 = + book("phi_pt10", fmt::sprintf("%s #phi; PV tracks #phi;#tracks", s_10).c_str(), phiBin_, phiMin_, phiMax_); + eta_pt10 = + book("eta_pt10", fmt::sprintf("%s #phi; PV tracks #eta;#tracks", s_10).c_str(), etaBin_, etaMin_, etaMax_); + + dxyVsPhi_pt1 = book("dxyVsPhi_pt1", + fmt::sprintf("%s d_{xy} (#mum) VS track #phi", s_1).c_str(), + phiBin_, + phiMin_, + phiMax_, + dxyMin_, + dxyMax_); + dxyVsPhi_pt1->SetXTitle("PV track (p_{T} > 1 GeV) #phi"); + dxyVsPhi_pt1->SetYTitle("PV track (p_{T} > 1 GeV) d_{xy} (#mum)"); + + dzVsPhi_pt1 = book("dzVsPhi_pt1", + fmt::sprintf("%s d_{z} (#mum) VS track #phi", s_1).c_str(), + phiBin_, + phiMin_, + phiMax_, + dzMin_, + dzMax_); + dzVsPhi_pt1->SetXTitle("PV track (p_{T} > 1 GeV) #phi"); + dzVsPhi_pt1->SetYTitle("PV track (p_{T} > 1 GeV) d_{z} (#mum)"); + + dxyVsEta_pt1 = book("dxyVsEta_pt1", + fmt::sprintf("%s d_{xy} (#mum) VS track #eta", s_1).c_str(), + etaBin_, + etaMin_, + etaMax_, + dxyMin_, + dxyMax_); + dxyVsEta_pt1->SetXTitle("PV track (p_{T} > 1 GeV) #eta"); + dxyVsEta_pt1->SetYTitle("PV track (p_{T} > 1 GeV) d_{xy} (#mum)"); + + dzVsEta_pt1 = book("dzVsEta_pt1", + fmt::sprintf("%s d_{z} (#mum) VS track #eta", s_1).c_str(), + etaBin_, + etaMin_, + etaMax_, + dzMin_, + dzMax_); + dzVsEta_pt1->SetXTitle("PV track (p_{T} > 1 GeV) #eta"); + dzVsEta_pt1->SetYTitle("PV track (p_{T} > 1 GeV) d_{z} (#mum)"); + + dxyVsEtaVsPhi_pt1 = book("dxyVsEtaVsPhi_pt1", + fmt::sprintf("%s d_{xy} (#mum) VS track #eta VS track #phi", s_1).c_str(), + etaBin2D_, + etaMin_, + etaMax_, + phiBin2D_, + phiMin_, + phiMax_, + dxyMin_, + dxyMax_); + dxyVsEtaVsPhi_pt1->SetXTitle("PV track (p_{T} > 1 GeV) #eta"); + dxyVsEtaVsPhi_pt1->SetYTitle("PV track (p_{T} > 1 GeV) #phi"); + dxyVsEtaVsPhi_pt1->SetZTitle("PV track (p_{T} > 1 GeV) d_{xy} (#mum)"); + + dzVsEtaVsPhi_pt1 = book("dzVsEtaVsPhi_pt1", + fmt::sprintf("%s d_{z} (#mum) VS track #eta VS track #phi", s_1).c_str(), + etaBin2D_, + etaMin_, + etaMax_, + phiBin2D_, + phiMin_, + phiMax_, + dzMin_, + dzMax_); + dzVsEtaVsPhi_pt1->SetXTitle("PV track (p_{T} > 1 GeV) #eta"); + dzVsEtaVsPhi_pt1->SetYTitle("PV track (p_{T} > 1 GeV) #phi"); + dzVsEtaVsPhi_pt1->SetZTitle("PV track (p_{T} > 1 GeV) d_{z} (#mum)"); + + dxyVsPhi_pt10 = book("dxyVsPhi_pt10", + fmt::sprintf("%s d_{xy} (#mum) VS track #phi", s_10).c_str(), + phiBin_, + phiMin_, + phiMax_, + dxyMin_, + dxyMax_); + dxyVsPhi_pt10->SetXTitle("PV track (p_{T} > 10 GeV) #phi"); + dxyVsPhi_pt10->SetYTitle("PV track (p_{T} > 10 GeV) d_{xy} (#mum)"); + + dzVsPhi_pt10 = book("dzVsPhi_pt10", + fmt::sprintf("%s d_{z} (#mum) VS track #phi", s_10).c_str(), + phiBin_, + phiMin_, + phiMax_, + dzMin_, + dzMax_); + dzVsPhi_pt10->SetXTitle("PV track (p_{T} > 10 GeV) #phi"); + dzVsPhi_pt10->SetYTitle("PV track (p_{T} > 10 GeV) d_{z} (#mum)"); + + dxyVsEta_pt10 = book("dxyVsEta_pt10", + fmt::sprintf("%s d_{xy} (#mum) VS track #eta", s_10).c_str(), + etaBin_, + etaMin_, + etaMax_, + dxyMin_, + dxyMax_); + dxyVsEta_pt10->SetXTitle("PV track (p_{T} > 10 GeV) #eta"); + dxyVsEta_pt10->SetYTitle("PV track (p_{T} > 10 GeV) d_{xy} (#mum)"); + + dzVsEta_pt10 = book("dzVsEta_pt10", + fmt::sprintf("%s d_{z} (#mum) VS track #eta", s_10).c_str(), + etaBin_, + etaMin_, + etaMax_, + dzMin_, + dzMax_); + dzVsEta_pt10->SetXTitle("PV track (p_{T} > 10 GeV) #eta"); + dzVsEta_pt10->SetYTitle("PV track (p_{T} > 10 GeV) d_{z} (#mum)"); + + dxyVsEtaVsPhi_pt10 = book("dxyVsEtaVsPhi_pt10", + fmt::sprintf("%s d_{xy} (#mum) VS track #eta VS track #phi", s_10).c_str(), + etaBin2D_, + etaMin_, + etaMax_, + phiBin2D_, + phiMin_, + phiMax_, + dxyMin_, + dxyMax_); + dxyVsEtaVsPhi_pt10->SetXTitle("PV track (p_{T} > 10 GeV) #eta"); + dxyVsEtaVsPhi_pt10->SetYTitle("PV track (p_{T} > 10 GeV) #phi"); + dxyVsEtaVsPhi_pt10->SetZTitle("PV track (p_{T} > 10 GeV) d_{xy} (#mum)"); + + dzVsEtaVsPhi_pt10 = book("dzVsEtaVsPhi_pt10", + fmt::sprintf("%s d_{z} (#mum) VS track #eta VS track #phi", s_10).c_str(), + etaBin2D_, + etaMin_, + etaMax_, + phiBin2D_, + phiMin_, + phiMax_, + dzMin_, + dzMax_); + dzVsEtaVsPhi_pt10->SetXTitle("PV track (p_{T} > 10 GeV) #eta"); + dzVsEtaVsPhi_pt10->SetYTitle("PV track (p_{T} > 10 GeV) #phi"); + dzVsEtaVsPhi_pt10->SetZTitle("PV track (p_{T} > 10 GeV) d_{z} (#mum)"); +} + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void GeneralPurposeVertexAnalyzer::fillDescriptions(edm::ConfigurationDescriptions &descriptions) { + edm::ParameterSetDescription desc; + desc.add("ndof", 4); + desc.add("vertexLabel", edm::InputTag("offlinePrimaryVertices")); + desc.add("beamSpotLabel", edm::InputTag("offlineBeamSpot")); + desc.add("Xpos", 0.1); + desc.add("Ypos", 0.0); + desc.add("TkSizeBin", 100); + desc.add("TkSizeMin", 499.5); + desc.add("TkSizeMax", -0.5); + desc.add("DxyBin", 100); + desc.add("DxyMin", 5000.); + desc.add("DxyMax", -5000.); + desc.add("DzBin", 100); + desc.add("DzMin", -2000.0); + desc.add("DzMax", 2000.0); + desc.add("PhiBin", 32); + desc.add("PhiBin2D", 12); + desc.add("PhiMin", -M_PI); + desc.add("PhiMax", M_PI); + desc.add("EtaBin", 26); + desc.add("EtaBin2D", 8); + desc.add("EtaMin", -2.7); + desc.add("EtaMax", 2.7); + descriptions.addWithDefaultLabel(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(GeneralPurposeVertexAnalyzer); diff --git a/Alignment/OfflineValidation/python/TkAlAllInOneTool/Zmumu_cfg.py b/Alignment/OfflineValidation/python/TkAlAllInOneTool/Zmumu_cfg.py index 83bbb200c3796..2d3699d09972f 100644 --- a/Alignment/OfflineValidation/python/TkAlAllInOneTool/Zmumu_cfg.py +++ b/Alignment/OfflineValidation/python/TkAlAllInOneTool/Zmumu_cfg.py @@ -128,17 +128,17 @@ ################################################################### # The Di Muon Mass Validation module ################################################################### -from Alignment.OfflineValidation.diMuonValidation_cfi.py import diMuonValidation as _diMuonValidation +from Alignment.OfflineValidation.diMuonValidation_cfi import diMuonValidation as _diMuonValidation process.DiMuonMassValidation = _diMuonValidation.clone( TkTag = 'TrackRefitter', # mu mu mass Pair_mass_min = 80., Pair_mass_max = 120., Pair_mass_nbins = 80, - Pair_etaminpos = -1, - Pair_etamaxpos = 1, - Pair_etaminneg = -1, - Pair_etamaxneg = 1, + Pair_etaminpos = -2.4, + Pair_etamaxpos = 2.4, + Pair_etaminneg = -2.4, + Pair_etamaxneg = 2.4, # cosTheta CS Variable_CosThetaCS_xmin = -1., Variable_CosThetaCS_xmax = 1., diff --git a/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py b/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py index 9b1c4128db275..471596b250a37 100644 --- a/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py +++ b/Alignment/OfflineValidation/test/testPrimaryVertexRelatedValidations_cfg.py @@ -376,10 +376,33 @@ def switchClusterizerParameters(da): ################################################################### # The analysis module ################################################################### -process.myanalysis = cms.EDAnalyzer("GeneralPurposeTrackAnalyzer", - TkTag = cms.InputTag('FinalTrackRefitter'), - isCosmics = cms.bool(False) - ) +process.trackanalysis = cms.EDAnalyzer("GeneralPurposeTrackAnalyzer", + TkTag = cms.InputTag('FinalTrackRefitter'), + isCosmics = cms.bool(False)) + +process.vertexanalysis = cms.EDAnalyzer('GeneralPurposeVertexAnalyzer', + ndof = cms.int32(4), + vertexLabel = cms.InputTag('offlinePrimaryVerticesFromRefittedTrks'), + beamSpotLabel = cms.InputTag('offlineBeamSpot'), + Xpos = cms.double(0.1), + Ypos = cms.double(0), + TkSizeBin = cms.int32(100), + TkSizeMin = cms.double(499.5), + TkSizeMax = cms.double(-0.5), + DxyBin = cms.int32(100), + DxyMin = cms.double(5000), + DxyMax = cms.double(-5000), + DzBin = cms.int32(100), + DzMin = cms.double(-2000), + DzMax = cms.double(2000), + PhiBin = cms.int32(32), + PhiBin2D = cms.int32(12), + PhiMin = cms.double(-3.1415926535897931), + PhiMax = cms.double(3.1415926535897931), + EtaBin = cms.int32(26), + EtaBin2D = cms.int32(8), + EtaMin = cms.double(-2.7), + EtaMax = cms.double(2.7)) ################################################################### # The PV resolution module @@ -398,5 +421,6 @@ def switchClusterizerParameters(da): process.seqTrackselRefit + process.offlinePrimaryVerticesFromRefittedTrks + process.PrimaryVertexResolution + - process.myanalysis + process.trackanalysis + + process.vertexanalysis )