From 75b54ee6fa78ba5d82979b87be36d68e72e63634 Mon Sep 17 00:00:00 2001 From: mmusich <marco.musich@cern.ch> Date: Wed, 10 Jan 2024 18:22:06 +0100 Subject: [PATCH 1/3] add SagittaBiasNtuplizer Co-authored-by: c-alexe <alexe.cristina99@gmail.com> --- .../plugins/SagittaBiasNtuplizer.cc | 527 ++++++++++++++++++ .../test/SagittaBiasNtuplizer_cfg.py | 233 ++++++++ 2 files changed, 760 insertions(+) create mode 100644 Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc create mode 100644 Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py diff --git a/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc b/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc new file mode 100644 index 0000000000000..203d7a09b068b --- /dev/null +++ b/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc @@ -0,0 +1,527 @@ +// -*- C++ -*- +// +// Package: Alignment/OfflineValidation +// Class: SagittaBiasNtuplizer +// +/* + *\class SagittaBiasNtuplizer SagittaBiasNtuplizer.cc Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc + + Description: This module is meant to create a simple ntuple to be able to validate the tracker alignment for the presence of Sagitta Bias. + + Implementation: the implementation is straightforward. + +*/ +// +// Original Author: Marco Musich, C. Alexe +// Created: Fri, 05 Jan 2023 11:41:00 GMT +// +// + +#include "CommonTools/UtilAlgos/interface/TFileService.h" +#include "DataFormats/BeamSpot/interface/BeamSpot.h" +#include "DataFormats/HepMCCandidate/interface/GenParticle.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/MuonReco/interface/Muon.h" +#include "DataFormats/MuonReco/interface/MuonFwd.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/Framework/interface/Event.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/one/EDAnalyzer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/ServiceRegistry/interface/Service.h" + +#include "TH1F.h" +#include "TH2I.h" +#include "TTree.h" +#include "TLorentzVector.h" + +class SagittaBiasNtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResources> { +public: + explicit SagittaBiasNtuplizer(const edm::ParameterSet&); + ~SagittaBiasNtuplizer() override = default; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + +private: + void analyze(const edm::Event& event, const edm::EventSetup& setup) override; + double openingAngle(const reco::Track* track1, const reco::Track* track2); + std::pair<unsigned int, reco::Vertex> findClosestVertex(const reco::Track* track1, + const reco::VertexCollection& vertices); + const bool useReco_; + const bool doGen_; + std::vector<double> pTthresholds_; + + // either on or the other! + //used to select what muon tracks to read from configuration file + edm::EDGetTokenT<reco::MuonCollection> muonsToken_; + //used to select what tracks to read from configuration file + edm::EDGetTokenT<reco::TrackCollection> alcaRecoToken_; + + //used to select what tracks to read from configuration file + edm::EDGetTokenT<reco::TrackCollection> tracksToken_; + + // for associated genParticles + edm::EDGetTokenT<edm::View<reco::Candidate>> genParticlesToken_; + //edm::EDGetTokenT<std::vector<reco::GenParticle>> genParticlesToken_; + + const edm::EDGetTokenT<reco::VertexCollection> vtxToken_; + const edm::EDGetTokenT<reco::BeamSpot> bsToken_; + + static double constexpr muMass = 0.1056583745; //!< Muon mass [GeV] + + TTree* tree_; + float mass_; + float posTrackDz_; + float negTrackDz_; + float posTrackD0_; + float negTrackD0_; + float posTrackEta_; + float negTrackEta_; + float posTrackPhi_; + float negTrackPhi_; + float posTrackPt_; + float negTrackPt_; + + // for the gen-level info + float genPosMuonDz_; + float genNegMuonDz_; + float genPosMuonD0_; + float genNegMuonD0_; + float genPosMuonEta_; + float genNegMuonEta_; + float genPosMuonPhi_; + float genNegMuonPhi_; + float genPosMuonPt_; + float genNegMuonPt_; + + // control histograms + TH2I* h_VertexMatrix; + TH1F* h_cutFlow; + TH1F* h_DeltaD0; + TH1F* h_DeltaDz; + TH1F* h_CosOpeningAngle; +}; + +//_________________________________________________________________________________ +SagittaBiasNtuplizer::SagittaBiasNtuplizer(const edm::ParameterSet& iConfig) + : useReco_(iConfig.getParameter<bool>("useReco")), + doGen_(iConfig.getParameter<bool>("doGen")), + pTthresholds_(iConfig.getParameter<std::vector<double>>("pTThresholds")), + vtxToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))), + bsToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))), + mass_{0}, + posTrackDz_{0}, + negTrackDz_{0}, + posTrackD0_{0}, + negTrackD0_{0}, + posTrackEta_{0}, + negTrackEta_{0}, + posTrackPhi_{0}, + negTrackPhi_{0}, + posTrackPt_{0}, + negTrackPt_{0}, + genPosMuonEta_{-99.}, + genNegMuonEta_{-99.}, + genPosMuonPhi_{-99.}, + genNegMuonPhi_{-99.}, + genPosMuonPt_{-99.}, + genNegMuonPt_{-99.} { + if (useReco_) { + tracksToken_ = mayConsume<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("tracks")); + muonsToken_ = mayConsume<reco::MuonCollection>(iConfig.getParameter<edm::InputTag>("muons")); + } else { + alcaRecoToken_ = mayConsume<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("muonTracks")); + } + + if (doGen_) { + genParticlesToken_ = consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("genParticles")); + } + + usesResource(TFileService::kSharedResource); + + // sort the vector of thresholds + std::sort(pTthresholds_.begin(), pTthresholds_.end(), [](const double& lhs, const double& rhs) { return lhs > rhs; }); + + edm::LogInfo("SagittaBiasNtuplizer") << __FUNCTION__; + for (const auto& thr : pTthresholds_) { + edm::LogInfo("SagittaBiasNtuplizer") << " Threshold: " << thr << " "; + } + + edm::Service<TFileService> fs; + h_cutFlow = fs->make<TH1F>("cutFlow", "cutFlow;cut;remaining events", 9, -0.5, 8.5); + std::string labels[9] = {"all events", + "common vertex", + "d0 cut", + "p_{T} cut", + "#eta cut", + "mass window", + "#delta d_{0}", + "#delta d_{z}", + "opening angle"}; + unsigned int count{0}; + for (const auto& label : labels) { + count++; + h_cutFlow->GetXaxis()->SetBinLabel(count, label.c_str()); + } + + h_VertexMatrix = fs->make<TH2I>("VertexMatrix", + ";index of closest vertex to #mu_{0};index of closest vertex to #mu_{1}", + 100, + 0, + 100, + 100, + 0, + 100); + h_DeltaD0 = fs->make<TH1F>("DeltaD0", "#Deltad_{0};#Deltad_{0} [cm];events", 100, -0.5, 0.5); + h_DeltaDz = fs->make<TH1F>("DeltaDz", "#Deltad_{z};#Deltad_{z} [cm];events", 100, -1, 1); + h_CosOpeningAngle = fs->make<TH1F>("OpeningAngle", "cos(#gamma(#mu^{+},#mu^{-}));events", 100, -1., 1.); + + tree_ = fs->make<TTree>("tree", "My TTree"); + tree_->Branch("mass", &mass_, "mass/F"); + tree_->Branch("posTrackDz", &posTrackDz_, "posTrackDz/F"); + tree_->Branch("negTrackDz", &negTrackDz_, "negTrackDz/F"); + tree_->Branch("posTrackD0", &posTrackD0_, "posTrackD0/F"); + tree_->Branch("negTrackD0", &negTrackD0_, "negTrackD0/F"); + tree_->Branch("posTrackEta", &posTrackEta_, "posTrackEta/F"); + tree_->Branch("negTrackEta", &negTrackEta_, "negTrackEta/F"); + tree_->Branch("posTrackPhi", &posTrackPhi_, "posTrackPhi/F"); + tree_->Branch("negTrackPhi", &negTrackPhi_, "negTrackPhi/F"); + tree_->Branch("posTrackPt", &posTrackPt_, "posTrackPt/F"); + tree_->Branch("negTrackPt", &negTrackPt_, "negTrackPt/F"); + + if (doGen_) { + tree_->Branch("genPosMuonDz", &genPosMuonDz_, "genPosMuonDz/F"); + tree_->Branch("genNegMuonDz", &genNegMuonDz_, "genNegMuonDz/F"); + tree_->Branch("genPosMuonD0", &genPosMuonD0_, "genPosMuonD0/F"); + tree_->Branch("genNegMuonD0", &genNegMuonD0_, "genNegMuonD0/F"); + tree_->Branch("genPosMuonEta", &genPosMuonEta_, "genPosMuonEta/F"); + tree_->Branch("genNegMuonEta", &genNegMuonEta_, "genNegMuonEta/F"); + tree_->Branch("genPosMuonPhi", &genPosMuonPhi_, "genPosMuonPhi/F"); + tree_->Branch("genNegMuonPhi", &genNegMuonPhi_, "genNegMuonPhi/F"); + tree_->Branch("genPosMuonPt", &genPosMuonPt_, "genPosMuonPt/F"); + tree_->Branch("genNegMuonPt", &genNegMuonPt_, "genNegMuonPt/F"); + } +} + +//_________________________________________________________________________________ +void SagittaBiasNtuplizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.ifValue( + edm::ParameterDescription<bool>("useReco", true, true), + true >> edm::ParameterDescription<edm::InputTag>("muons", edm::InputTag("muons"), true) or + false >> edm::ParameterDescription<edm::InputTag>("muonTracks", edm::InputTag("ALCARECOTkAlZMuMu"), true)) + ->setComment("If useReco is true need to specify the muon tracks, otherwise take the ALCARECO Inner tracks"); + desc.add<bool>("doGen", false); + desc.add<edm::InputTag>("genParticles", edm::InputTag("genParticles")); + desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks")); + desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices")); + desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot")); + desc.add<std::vector<double>>("pTThresholds", {30., 10.}); + descriptions.addWithDefaultLabel(desc); +} + +//_________________________________________________________________________________ +void SagittaBiasNtuplizer::analyze(const edm::Event& event, const edm::EventSetup& setup) { + h_cutFlow->Fill(0); + + // the di-muon tracks + std::vector<const reco::Track*> myTracks; + + // if we have to start from scratch from RECO data-tier + if (useReco_) { + // select the good muons + std::vector<const reco::Muon*> myGoodMuonVector; + for (const auto& muon : event.get(muonsToken_)) { + const reco::TrackRef t = muon.innerTrack(); + if (!t.isNull()) { + if (t->quality(reco::TrackBase::highPurity)) { + if (t->chi2() / t->ndof() <= 2.5 && t->numberOfValidHits() >= 5 && + t->hitPattern().numberOfValidPixelHits() >= 2 && t->quality(reco::TrackBase::highPurity)) + myGoodMuonVector.emplace_back(&muon); + } + } + } + + LogDebug("SagittaBiasNtuplizer") << "myGoodMuonVector size: " << myGoodMuonVector.size() << std::endl; + std::sort(myGoodMuonVector.begin(), myGoodMuonVector.end(), [](const reco::Muon*& lhs, const reco::Muon*& rhs) { + return lhs->pt() > rhs->pt(); + }); + + // just check the ordering + for (const auto& muon : myGoodMuonVector) { + LogDebug("SagittaBiasNtuplizer") << "pT: " << muon->pt() << " "; + } + LogDebug("SagittaBiasNtuplizer") << std::endl; + + // reject if there's no Z + if (myGoodMuonVector.size() < 2) + return; + + if ((myGoodMuonVector[0]->pt()) < pTthresholds_[0] || (myGoodMuonVector[1]->pt() < pTthresholds_[1])) + return; + + if (myGoodMuonVector[0]->charge() * myGoodMuonVector[1]->charge() > 0) + return; + + //const auto& m1 = myGoodMuonVector[1]->p4(); + //const auto& m0 = myGoodMuonVector[0]->p4(); + //const auto& mother = m0 + m1; + + // just copy the top two muons + std::vector<const reco::Muon*> theZMuonVector; + theZMuonVector.reserve(2); + theZMuonVector.emplace_back(myGoodMuonVector[1]); + theZMuonVector.emplace_back(myGoodMuonVector[0]); + + // do the matching of Z muons with inner tracks + unsigned int i = 0; + for (const auto& muon : theZMuonVector) { + i++; + float minD = 1e6; + const reco::Track* theMatch = nullptr; + for (const auto& track : event.get(tracksToken_)) { + float D = ::deltaR2(muon->eta(), muon->phi(), track.eta(), track.phi()); + if (D < minD) { + minD = D; + theMatch = &track; + } + } + LogDebug("SagittaBiasNtuplizer") << "pushing new track: " << i << std::endl; + myTracks.emplace_back(theMatch); + } + } else { + // we start directly with the pre-selected ALCARECO tracks + for (const auto& muon : event.get(alcaRecoToken_)) { + myTracks.emplace_back(&muon); + } + } + + const reco::VertexCollection& vertices = event.get(vtxToken_); + const reco::BeamSpot& beamSpot = event.get(bsToken_); + math::XYZPoint bs(beamSpot.x0(), beamSpot.y0(), beamSpot.z0()); + + //edm::LogPrint("SagittaBiasNtuplizer") << " track size:" << myTracks.size() << " vertex size:" << vertices.size() << std::endl; + + if ((myTracks.size() != 2)) { + LogTrace("SagittaBiasNtuplizer") << "Found " << myTracks.size() << " muons in the event. Skipping"; + return; + } + + bool passSameVertex{true}; + bool passD0sigCut{true}; + bool passPtCut{true}; + bool passEtaCut{true}; + bool passMassWindow{true}; + bool passDeltaD0{true}; + bool passDeltaDz{true}; + bool passOpeningAngle{true}; + double d0[2] = {0., 0.}; + double dz[2] = {0., 0.}; + unsigned int vtxIndex[2] = {999, 999}; + + unsigned int i = 0; + for (const auto& track : myTracks) { + if (track->pt() < 12) { + passPtCut = false; + continue; + } + + if (std::abs(track->eta()) > 2.5) { + passEtaCut = false; + continue; + } + + const auto& closestVertex = this->findClosestVertex(track, vertices); + vtxIndex[i] = closestVertex.first; + d0[i] = track->dxy(closestVertex.second.position()); + dz[i] = track->dz(closestVertex.second.position()); + + if (d0[i] / track->dxyError() > 4) { + passD0sigCut = false; + continue; + } + + if (track->charge() > 0) { + posTrackDz_ = dz[i]; + posTrackD0_ = d0[i]; + posTrackEta_ = track->eta(); + posTrackPhi_ = track->phi(); + posTrackPt_ = track->pt(); + } else { + negTrackDz_ = dz[i]; + negTrackD0_ = d0[i]; + negTrackEta_ = track->eta(); + negTrackPhi_ = track->phi(); + negTrackPt_ = track->pt(); + } + i++; + } + + h_VertexMatrix->Fill(vtxIndex[0], vtxIndex[1]); + // check if the two muons have the same vertex + passSameVertex = (vtxIndex[0] == vtxIndex[1]); + if (!passSameVertex) + return; + h_cutFlow->Fill(1); + + // checks if both muons pass the IP signficance cut (w.r.t BeamSpot) + if (!passD0sigCut) + return; + h_cutFlow->Fill(2); + + // checks if both muons pass the pT cut + if (!passPtCut) + return; + h_cutFlow->Fill(3); + + // check if both muons pass the eta cut + if (!passEtaCut) + return; + h_cutFlow->Fill(4); + + // compute the invariant mass of the system + TLorentzVector posTrack, negTrack, mother; + posTrack.SetPtEtaPhiM(posTrackPt_, posTrackEta_, posTrackPhi_, muMass); // assume muon mass for tracks + negTrack.SetPtEtaPhiM(negTrackPt_, negTrackEta_, negTrackPhi_, muMass); + mother = posTrack + negTrack; + mass_ = mother.M(); + + // checks if invariant mass of the system lies in the fiducial window + passMassWindow = (mass_ > 70 && mass_ < 110); + if (!passMassWindow) + return; + h_cutFlow->Fill(5); + + // checks if the di-muon system passes the d0 compatibility cut + passDeltaD0 = (std::abs(d0[0] - d0[1]) < 0.01); + h_DeltaD0->Fill(d0[0] - d0[1]); + h_DeltaDz->Fill(dz[0] - dz[1]); + if (!passDeltaD0) + return; + h_cutFlow->Fill(6); + + // checks if the di-muon system passes the z0 compatibility cut + passDeltaDz = (std::abs(dz[0] - dz[1]) < 0.06); + if (!passDeltaDz) + return; + h_cutFlow->Fill(7); + + // checks if the di-muon system passes the opening angle cut + double openingAngle = this->openingAngle(myTracks[0], myTracks[1]); + h_CosOpeningAngle->Fill(openingAngle); + passOpeningAngle = true; //(openingAngle > M_PI/4.); + + if (!passOpeningAngle) + return; + h_cutFlow->Fill(8); + + if (doGen_) { + const edm::View<reco::Candidate>* genPartCollection = &event.get(genParticlesToken_); + + // loop on the reconstructed tracks + for (const auto& track : myTracks) { + float drmin = 0.01; + // loop on the gen particles + for (auto g = genPartCollection->begin(); g != genPartCollection->end(); ++g) { + if (g->status() != 1) + continue; + + if (std::abs(g->pdgId()) != 13) + continue; + + if (g->charge() != track->charge()) + continue; + + float dR = reco::deltaR2(*g, *track); + + auto const& vtx = g->vertex(); + auto const& myBeamSpot = beamSpot.position(vtx.z()); + const auto& theptinv2 = 1 / g->pt() * g->pt(); + + if (dR < drmin) { + drmin = dR; + + if (g->charge() > 0) { + genPosMuonPt_ = g->pt(); + genPosMuonEta_ = g->eta(); + genPosMuonPhi_ = g->phi(); + //d0 + genPosMuonD0_ = -(-(vtx.x() - myBeamSpot.x()) * g->py() + (vtx.y() - myBeamSpot.y()) * g->px()) / g->pt(); + //dz + genPosMuonDz_ = + (vtx.z() - myBeamSpot.z()) - + ((vtx.x() - myBeamSpot.x()) * g->px() + (vtx.y() - myBeamSpot.y()) * g->py()) * g->pz() * theptinv2; + } else { + genNegMuonPt_ = g->pt(); + genNegMuonEta_ = g->eta(); + genNegMuonPhi_ = g->phi(); + //d0 + genNegMuonD0_ = (-(vtx.x() - myBeamSpot.x()) * g->py() + (vtx.y() - myBeamSpot.y()) * g->px()) / g->pt(); + //dz + genNegMuonDz_ = + (vtx.z() - myBeamSpot.z()) - + ((vtx.x() - myBeamSpot.x()) * g->px() + (vtx.y() - myBeamSpot.y()) * g->py()) * g->pz() * theptinv2; + } + } else { + continue; + } + } + } + } + + // if all cuts are passed fill the TTree + tree_->Fill(); +} + +//_________________________________________________________________________________ +double SagittaBiasNtuplizer::openingAngle(const reco::Track* trk1, const reco::Track* trk2) { + math::XYZVector vec1(trk1->px(), trk1->py(), trk1->pz()); + math::XYZVector vec2(trk2->px(), trk2->py(), trk2->pz()); + return vec1.Dot(vec2) / (vec1.R() * vec2.R()); +} + +/* +double SagittaBiasNtuplizer::openingAngle(const reco::Track& track1, const reco::Track& track2) { + double dPhi = std::abs(track1.phi() - track2.phi()); + double dEta = std::abs(track1.eta() - track2.eta()); + double deltaR2 = reco::deltaR2(track1, track2); + double openingAngle = 2 * std::atan2(std::sqrt(std::sin(dEta / 2) * std::sin(dEta / 2) + std::sinh(dPhi / 2) * std::sinh(dPhi / 2)), std::sqrt(1 - deltaR2)); + return openingAngle; +} +*/ + +//_________________________________________________________________________________ +std::pair<unsigned int, reco::Vertex> SagittaBiasNtuplizer::findClosestVertex(const reco::Track* track, + const reco::VertexCollection& vertices) { + // Initialize variables for minimum distance and closest vertex + double minDistance = std::numeric_limits<double>::max(); + reco::Vertex closestVertex; + + unsigned int index{0}, theIndex{999}; + // Loop over the primary vertices and calculate the distance to the track + for (const auto& vertex : vertices) { + const math::XYZPoint& vertexPosition = vertex.position(); + + // Calculate the distance to the track + const auto& trackMomentum = track->momentum(); + const auto& vertexToPoint = vertexPosition - track->referencePoint(); + double distance = vertexToPoint.Cross(trackMomentum).R() / trackMomentum.R(); + + // Check if this is the closest vertex so far + if (distance < minDistance) { + minDistance = distance; + closestVertex = vertex; + theIndex = index; + } + index++; + } + return std::make_pair(theIndex, closestVertex); +} + +DEFINE_FWK_MODULE(SagittaBiasNtuplizer); diff --git a/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py b/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py new file mode 100644 index 0000000000000..c8256aeb46195 --- /dev/null +++ b/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py @@ -0,0 +1,233 @@ +import glob +import math +import FWCore.ParameterSet.Config as cms +import FWCore.Utilities.FileUtils as FileUtils +from FWCore.ParameterSet.VarParsing import VarParsing + +options = VarParsing('analysis') +options.register('scenario', + '0', + VarParsing.multiplicity.singleton, + VarParsing.varType.string, + "Name of input misalignment scenario") +options.parseArguments() + +valid_scenarios = ['-10e-6','-8e-6','-6e-6','-4e-6','-2e-6','0','2e-6','4e-6','6e-6','8e-6','10e-6'] + +if options.scenario not in valid_scenarios: + print("Error: Invalid scenario specified. Please choose from the following list: ") + print(valid_scenarios) + exit(1) + +process = cms.Process("DiMuonAnalysis") + +################################################################### +# Set the process to run multi-threaded +################################################################### +process.options.numberOfThreads = 8 + +################################################################### +# Message logger service +################################################################### +process.load("FWCore.MessageLogger.MessageLogger_cfi") +process.MessageLogger.cerr.enable = False +process.MessageLogger.ZtoMMNtupler=dict() +process.MessageLogger.cout = cms.untracked.PSet( + enable = cms.untracked.bool(True), + threshold = cms.untracked.string("INFO"), + default = cms.untracked.PSet(limit = cms.untracked.int32(0)), + FwkReport = cms.untracked.PSet(limit = cms.untracked.int32(-1), + reportEvery = cms.untracked.int32(10000) + ), + ZtoMMNtupler = cms.untracked.PSet( limit = cms.untracked.int32(-1)), + enableStatistics = cms.untracked.bool(True) + ) + +################################################################### +# Geometry producer and standard includes +################################################################### +process.load("RecoVertex.BeamSpotProducer.BeamSpot_cff") +process.load("Configuration.StandardSequences.Services_cff") +process.load("Configuration.StandardSequences.GeometryRecoDB_cff") +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load("CondCore.CondDB.CondDB_cfi") + +################################################################### +# TransientTrack from https://twiki.cern.ch/twiki/bin/view/CMSPublic/SWGuideTransientTracks +################################################################### +process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") +process.load('TrackPropagation.SteppingHelixPropagator.SteppingHelixPropagatorOpposite_cfi') +process.load('TrackPropagation.SteppingHelixPropagator.SteppingHelixPropagatorAlong_cfi') +process.load('TrackingTools.TrackAssociator.DetIdAssociatorESProducer_cff') + +#################################################################### +# Get the GlogalTag +#################################################################### +process.load("Configuration.StandardSequences.FrontierConditions_GlobalTag_cff") +from Configuration.AlCa.GlobalTag import GlobalTag +process.GlobalTag = GlobalTag(process.GlobalTag,"125X_mcRun3_2022_design_v6", '') +if (options.scenario=='null'): + print("null scenario, do nothing") + pass +elif (options.scenario=='ideal'): + print("ideal scenario, use ideal tags") + process.GlobalTag.toGet = cms.VPSet(cms.PSet(record = cms.string('TrackerAlignmentRcd'), + tag = cms.string("TrackerAlignment_Upgrade2017_design_v4")), + cms.PSet(record = cms.string('TrackerAlignmentErrorExtendedRcd'), + tag = cms.string("TrackerAlignmentErrorsExtended_Upgrade2017_design_v0")), + cms.PSet(record = cms.string('TrackerSurfaceDeformationRcd'), + tag = cms.string("TrackerSurfaceDeformations_zero"))) +else : + print("using {} scenario".format(options.scenario)) + process.GlobalTag.toGet = cms.VPSet(cms.PSet(connect = cms.string("sqlite_file:/afs/cern.ch/user/m/musich/public/layer_rotation_studies/outputfile_"+options.scenario+".db"), + record = cms.string('TrackerAlignmentRcd'), + tag = cms.string("Alignments"))) + +################################################################### +# Source +################################################################### +#filelist = FileUtils.loadListFromFile("listOfFiles_idealMC_TkAlDiMuonAndVertex.txt") +filelist = FileUtils.loadListFromFile("listOfFiles_idealMC_GEN-SIM-RECO.txt") +readFiles = cms.untracked.vstring( *filelist) +process.source = cms.Source("PoolSource",fileNames = readFiles) + +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(options.maxEvents)) + +################################################################### +# Alignment Track Selector +################################################################### +import Alignment.CommonAlignmentProducer.AlignmentTrackSelector_cfi +process.MuSkimSelector = Alignment.CommonAlignmentProducer.AlignmentTrackSelector_cfi.AlignmentTrackSelector.clone( + applyBasicCuts = True, + filter = True, + src = "ALCARECOTkAlDiMuon", + ptMin = 17., + pMin = 17., + etaMin = -2.5, + etaMax = 2.5, + d0Min = -2., + d0Max = 2., + dzMin = -25., + dzMax = 25., + nHitMin = 6, + nHitMin2D = 0) + +################################################################### +# refitting the muon tracks +################################################################### +process.load("RecoTracker.TrackProducer.TrackRefitters_cff") +import RecoTracker.TrackProducer.TrackRefitters_cff +process.refittedMuons = RecoTracker.TrackProducer.TrackRefitter_cfi.TrackRefitter.clone( + src = "ALCARECOTkAlDiMuon", + TrajectoryInEvent = True, + NavigationSchool = '', + TTRHBuilder = "WithAngleAndTemplate") + +################################################################### +# refitting the vertex tracks +################################################################### +process.refittedVtxTracks = RecoTracker.TrackProducer.TrackRefitter_cfi.TrackRefitter.clone( + src = "ALCARECOTkAlDiMuonVertexTracks", + TrajectoryInEvent = True, + NavigationSchool = '', + TTRHBuilder = "WithAngleAndTemplate") + +################################################################### +# refitting all tracks +################################################################### +process.refittedTracks = RecoTracker.TrackProducer.TrackRefitter_cfi.TrackRefitter.clone( + src = "generalTracks", + TrajectoryInEvent = True, + NavigationSchool = '', + TTRHBuilder = "WithAngleAndTemplate") + +#################################################################### +# Re-do vertices +#################################################################### +from RecoVertex.PrimaryVertexProducer.OfflinePrimaryVertices_cfi import offlinePrimaryVertices +process.offlinePrimaryVerticesFromRefittedTrks = offlinePrimaryVertices.clone() +#process.offlinePrimaryVerticesFromRefittedTrks.TrackLabel = cms.InputTag("refittedVtxTracks") +process.offlinePrimaryVerticesFromRefittedTrks.TrackLabel = cms.InputTag("refittedTracks") + +################################################################### +# The analysis modules +################################################################### +process.ZtoMMNtuple = cms.EDAnalyzer("SagittaBiasNtuplizer", + #tracks = cms.InputTag('refittedMuons'), + useReco = cms.bool(True), + muons = cms.InputTag('muons'), + doGen = cms.bool(True), + tracks = cms.InputTag('refittedTracks'), + vertices = cms.InputTag('offlinePrimaryVerticesFromRefittedTrks')) + +process.DiMuonVertexValidation = cms.EDAnalyzer("DiMuonVertexValidation", + useReco = cms.bool(False), + ## the two parameters below are mutually exclusive, + ## depending if RECO or ALCARECO is used + #muons = cms.InputTag(''), + muonTracks = cms.InputTag('refittedMuons'), + tracks = cms.InputTag(''), + vertices = cms.InputTag('offlinePrimaryVerticesFromRefittedTrks')) + +from Alignment.OfflineValidation.diMuonValidation_cfi import diMuonValidation as _diMuonValidation +process.DiMuonMassValidation = _diMuonValidation.clone( + TkTag = 'refittedMuons', + #TkTag = 'TrackRefitter1', + # mu mu mass + Pair_mass_min = 80., + Pair_mass_max = 120., + Pair_mass_nbins = 80, + 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., + Variable_CosThetaCS_nbins = 20, + # DeltaEta + Variable_DeltaEta_xmin = -4.8, + Variable_DeltaEta_xmax = 4.8, + Variable_DeltaEta_nbins = 20, + # EtaMinus + Variable_EtaMinus_xmin = -2.4, + Variable_EtaMinus_xmax = 2.4, + Variable_EtaMinus_nbins = 12, + # EtaPlus + Variable_EtaPlus_xmin = -2.4, + Variable_EtaPlus_xmax = 2.4, + Variable_EtaPlus_nbins = 12, + # Phi CS + Variable_PhiCS_xmin = -math.pi/2., + Variable_PhiCS_xmax = math.pi/2., + Variable_PhiCS_nbins = 20, + # Phi Minus + Variable_PhiMinus_xmin = -math.pi, + Variable_PhiMinus_xmax = math.pi, + Variable_PhiMinus_nbins = 16, + # Phi Plus + Variable_PhiPlus_xmin = -math.pi, + Variable_PhiPlus_xmax = math.pi, + Variable_PhiPlus_nbins = 16, + # mu mu pT + Variable_PairPt_xmin = 0., + Variable_PairPt_xmax = 100., + Variable_PairPt_nbins = 100) + +################################################################### +# Output name +################################################################### +process.TFileService = cms.Service("TFileService", + fileName = cms.string("ZmmNtuple_MC_GEN-SIM_"+options.scenario+".root")) + +################################################################### +# Path +################################################################### +process.p1 = cms.Path(process.offlineBeamSpot + #* process.refittedMuons + #* process.refittedVtxTracks + * process.refittedTracks + * process.offlinePrimaryVerticesFromRefittedTrks + * process.ZtoMMNtuple) + #* process.DiMuonVertexValidation + #* process.DiMuonMassValidation) From 5038751b02cc943e305b1518e63d083d41a75465 Mon Sep 17 00:00:00 2001 From: mmusich <marco.musich@cern.ch> Date: Fri, 1 Mar 2024 12:13:36 +0100 Subject: [PATCH 2/3] Add SagittaBiasNtuplizer to unit testing infrastructure --- .../test/SagittaBiasNtuplizer_cfg.py | 53 +++++++++++++++---- .../test/testTrackAnalyzers.cc | 12 +++++ .../testingScripts/test_unitMiscellanea.sh | 5 +- 3 files changed, 59 insertions(+), 11 deletions(-) diff --git a/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py b/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py index c8256aeb46195..e5f08e57ebf02 100644 --- a/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py +++ b/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py @@ -3,23 +3,51 @@ import FWCore.ParameterSet.Config as cms import FWCore.Utilities.FileUtils as FileUtils from FWCore.ParameterSet.VarParsing import VarParsing +from Alignment.OfflineValidation.TkAlAllInOneTool.defaultInputFiles_cff import filesDefaultMC_DoubleMuon_string options = VarParsing('analysis') options.register('scenario', - '0', + 'null', VarParsing.multiplicity.singleton, VarParsing.varType.string, "Name of input misalignment scenario") + +options.register('globalTag', + "'auto:phase1_2022_realistic", # default value + VarParsing.multiplicity.singleton, # singleton or list + VarParsing.varType.string, # string, int, or float + "name of the input Global Tag") + +options.register ('myfile', + filesDefaultMC_DoubleMuon_string, # default value + VarParsing.multiplicity.singleton, + VarParsing.varType.string, + "file name") + +options.register ('FileList', + '', # default value + VarParsing.multiplicity.singleton, + VarParsing.varType.string, + "FileList in DAS format") + options.parseArguments() -valid_scenarios = ['-10e-6','-8e-6','-6e-6','-4e-6','-2e-6','0','2e-6','4e-6','6e-6','8e-6','10e-6'] +if(options.FileList): + print("FileList: ", options.FileList) +else: + print("inputFile: ", options.myfile) +print("outputFile: ", "ZmmNtuple_MC_GEN-SIM_{fscenario}.root".format(fscenario=options.scenario)) +print("conditionGT: ", options.globalTag) +print("max events: ", options.maxEvents) + +valid_scenarios = ['-10e-6','-8e-6','-6e-6','-4e-6','-2e-6','0','2e-6','4e-6','6e-6','8e-6','10e-6','null'] if options.scenario not in valid_scenarios: print("Error: Invalid scenario specified. Please choose from the following list: ") print(valid_scenarios) exit(1) -process = cms.Process("DiMuonAnalysis") +process = cms.Process("SagittaBiasNtuplizer") ################################################################### # Set the process to run multi-threaded @@ -31,15 +59,15 @@ ################################################################### process.load("FWCore.MessageLogger.MessageLogger_cfi") process.MessageLogger.cerr.enable = False -process.MessageLogger.ZtoMMNtupler=dict() +process.MessageLogger.SagittaBiasNtuplizer=dict() process.MessageLogger.cout = cms.untracked.PSet( enable = cms.untracked.bool(True), threshold = cms.untracked.string("INFO"), default = cms.untracked.PSet(limit = cms.untracked.int32(0)), FwkReport = cms.untracked.PSet(limit = cms.untracked.int32(-1), - reportEvery = cms.untracked.int32(10000) + reportEvery = cms.untracked.int32(1000) ), - ZtoMMNtupler = cms.untracked.PSet( limit = cms.untracked.int32(-1)), + SagittaBiasNtuplizer = cms.untracked.PSet( limit = cms.untracked.int32(-1)), enableStatistics = cms.untracked.bool(True) ) @@ -86,10 +114,15 @@ ################################################################### # Source ################################################################### -#filelist = FileUtils.loadListFromFile("listOfFiles_idealMC_TkAlDiMuonAndVertex.txt") -filelist = FileUtils.loadListFromFile("listOfFiles_idealMC_GEN-SIM-RECO.txt") -readFiles = cms.untracked.vstring( *filelist) -process.source = cms.Source("PoolSource",fileNames = readFiles) +if(options.FileList): + print('Loading file list from ASCII file') + filelist = FileUtils.loadListFromFile (options.FileList) + readFiles = cms.untracked.vstring( *filelist) +else: + readFiles = cms.untracked.vstring([options.myfile]) + +process.source = cms.Source("PoolSource", + fileNames = readFiles) process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(options.maxEvents)) diff --git a/Alignment/OfflineValidation/test/testTrackAnalyzers.cc b/Alignment/OfflineValidation/test/testTrackAnalyzers.cc index a90a28a2e6bb6..3a14d911accf9 100644 --- a/Alignment/OfflineValidation/test/testTrackAnalyzers.cc +++ b/Alignment/OfflineValidation/test/testTrackAnalyzers.cc @@ -161,3 +161,15 @@ TEST_CASE("TrackerGeometryCompare tests", "[TrackerGeometryCompare]") { const std::string baseConfig = generateBaseConfig("trackerGeometryCompare", "tesTrackAnalyzer16.root"); runTestForAnalyzer(baseConfig, "trackerGeometryCompare"); } + +//___________________________________________________________________________________________ +TEST_CASE("ShortenedTrackValidation tests", "[ShortenedTrackValidation]") { + const std::string baseConfig = generateBaseConfig("shortenedTrackValidation", "tesTrackAnalyzer17.root"); + runTestForAnalyzer(baseConfig, "shortenedTrackValidation"); +} + +//___________________________________________________________________________________________ +TEST_CASE("SagittaBiasNtuplizer tests", "[SagittaBiasNtuplizer]") { + const std::string baseConfig = generateBaseConfig("sagittaBiasNtuplizer", "tesTrackAnalyzer18.root"); + runTestForAnalyzer(baseConfig, "sagittaBiasNtuplizer"); +} diff --git a/Alignment/OfflineValidation/test/testingScripts/test_unitMiscellanea.sh b/Alignment/OfflineValidation/test/testingScripts/test_unitMiscellanea.sh index 1f339a3bb3779..9ace1d6594d03 100755 --- a/Alignment/OfflineValidation/test/testingScripts/test_unitMiscellanea.sh +++ b/Alignment/OfflineValidation/test/testingScripts/test_unitMiscellanea.sh @@ -14,4 +14,7 @@ echo "TESTING Pixel BaryCenter Analyser ..." cmsRun ${CMSSW_BASE}/src/Alignment/OfflineValidation/test/PixelBaryCentreAnalyzer_cfg.py unitTest=True || die "Failure running PixelBaryCentreAnalyzer_cfg.py" $? echo "TESTING CosmicTrackSplitting Analyser ..." -cmsRun ${CMSSW_BASE}/src/Alignment/OfflineValidation/test/testSplitterValidation_cfg.py unitTest=True || die "Failure running testSplitterValidation_cfg.py" $? +cmsRun ${CMSSW_BASE}/src/Alignment/OfflineValidation/test/testSplitterValidation_cfg.py unitTest=True || die "Failure running testSplitterValidation_cfg.py" $? + +echo "TESTING SagittaBiasNtuplizer Analyser ..." +cmsRun ${CMSSW_BASE}/src/Alignment/OfflineValidation/test/SagittaBiasNtuplizer_cfg.py || die "Failure running SagittaBiasNtuplizer_cfg.py" $? From 6937f908e47b41b31ca3cf6b1a336fd669b11eb2 Mon Sep 17 00:00:00 2001 From: mmusich <marco.musich@cern.ch> Date: Fri, 1 Mar 2024 14:53:08 +0100 Subject: [PATCH 3/3] make selection in SagittaBiasNtuplizer configurable --- .../plugins/SagittaBiasNtuplizer.cc | 35 +++++++++++++++---- 1 file changed, 29 insertions(+), 6 deletions(-) diff --git a/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc b/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc index 203d7a09b068b..0a5c313de907e 100644 --- a/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc +++ b/Alignment/OfflineValidation/plugins/SagittaBiasNtuplizer.cc @@ -55,6 +55,15 @@ class SagittaBiasNtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResourc const reco::VertexCollection& vertices); const bool useReco_; const bool doGen_; + + const double muonEtaCut_; + const double muonPtCut_; + const double muondxySigCut_; + const double minMassWindowCut_; + const double maxMassWindowCut_; + const double d0CompatibilityCut_; + const double z0CompatibilityCut_; + std::vector<double> pTthresholds_; // either on or the other! @@ -112,6 +121,13 @@ class SagittaBiasNtuplizer : public edm::one::EDAnalyzer<edm::one::SharedResourc SagittaBiasNtuplizer::SagittaBiasNtuplizer(const edm::ParameterSet& iConfig) : useReco_(iConfig.getParameter<bool>("useReco")), doGen_(iConfig.getParameter<bool>("doGen")), + muonEtaCut_(iConfig.getParameter<double>("muonEtaCut")), + muonPtCut_(iConfig.getParameter<double>("muonPtCut")), + muondxySigCut_(iConfig.getParameter<double>("muondxySigCut")), + minMassWindowCut_(iConfig.getParameter<double>("minMassWindowCut")), + maxMassWindowCut_(iConfig.getParameter<double>("maxMassWindowCut")), + d0CompatibilityCut_(iConfig.getParameter<double>("d0CompatibilityCut")), + z0CompatibilityCut_(iConfig.getParameter<double>("z0CompatibilityCut")), pTthresholds_(iConfig.getParameter<std::vector<double>>("pTThresholds")), vtxToken_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))), bsToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpot"))), @@ -222,6 +238,13 @@ void SagittaBiasNtuplizer::fillDescriptions(edm::ConfigurationDescriptions& desc desc.add<edm::InputTag>("tracks", edm::InputTag("generalTracks")); desc.add<edm::InputTag>("vertices", edm::InputTag("offlinePrimaryVertices")); desc.add<edm::InputTag>("beamSpot", edm::InputTag("offlineBeamSpot")); + desc.add<double>("muonEtaCut", 2.5)->setComment("muon system acceptance"); + desc.add<double>("muonPtCut", 12.)->setComment("in GeV"); + desc.add<double>("muondxySigCut", 4.)->setComment("significance of the d0 compatibility with closest vertex"); + desc.add<double>("minMassWindowCut", 70.)->setComment("in GeV"); + desc.add<double>("maxMassWindowCut", 110.)->setComment("in GeV"); + desc.add<double>("d0CompatibilityCut", 0.01)->setComment("d0 compatibility between the two muons"); + desc.add<double>("z0CompatibilityCut", 0.06)->setComment("z0 compatibility between the two muons"); desc.add<std::vector<double>>("pTThresholds", {30., 10.}); descriptions.addWithDefaultLabel(desc); } @@ -327,12 +350,12 @@ void SagittaBiasNtuplizer::analyze(const edm::Event& event, const edm::EventSetu unsigned int i = 0; for (const auto& track : myTracks) { - if (track->pt() < 12) { + if (track->pt() < muonPtCut_) { passPtCut = false; continue; } - if (std::abs(track->eta()) > 2.5) { + if (std::abs(track->eta()) > muonEtaCut_) { passEtaCut = false; continue; } @@ -342,7 +365,7 @@ void SagittaBiasNtuplizer::analyze(const edm::Event& event, const edm::EventSetu d0[i] = track->dxy(closestVertex.second.position()); dz[i] = track->dz(closestVertex.second.position()); - if (d0[i] / track->dxyError() > 4) { + if (d0[i] / track->dxyError() > muondxySigCut_) { passD0sigCut = false; continue; } @@ -393,13 +416,13 @@ void SagittaBiasNtuplizer::analyze(const edm::Event& event, const edm::EventSetu mass_ = mother.M(); // checks if invariant mass of the system lies in the fiducial window - passMassWindow = (mass_ > 70 && mass_ < 110); + passMassWindow = (mass_ > minMassWindowCut_ && mass_ < maxMassWindowCut_); if (!passMassWindow) return; h_cutFlow->Fill(5); // checks if the di-muon system passes the d0 compatibility cut - passDeltaD0 = (std::abs(d0[0] - d0[1]) < 0.01); + passDeltaD0 = (std::abs(d0[0] - d0[1]) < d0CompatibilityCut_); h_DeltaD0->Fill(d0[0] - d0[1]); h_DeltaDz->Fill(dz[0] - dz[1]); if (!passDeltaD0) @@ -407,7 +430,7 @@ void SagittaBiasNtuplizer::analyze(const edm::Event& event, const edm::EventSetu h_cutFlow->Fill(6); // checks if the di-muon system passes the z0 compatibility cut - passDeltaDz = (std::abs(dz[0] - dz[1]) < 0.06); + passDeltaDz = (std::abs(dz[0] - dz[1]) < z0CompatibilityCut_); if (!passDeltaDz) return; h_cutFlow->Fill(7);