File indexing completed on 2025-06-29 23:01:13
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 #include "RecoVertex/BeamSpotProducer/interface/PVFitter.h"
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "FWCore/ServiceRegistry/interface/Service.h"
0017 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0018 #include "CondFormats/BeamSpotObjects/interface/BeamSpotObjects.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020
0021 #include "FWCore/Utilities/interface/InputTag.h"
0022 #include "DataFormats/Common/interface/View.h"
0023
0024 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0025 #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
0026 #include "DataFormats/TrackReco/interface/Track.h"
0027 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0028 #include "DataFormats/TrackReco/interface/HitPattern.h"
0029 #include "DataFormats/VertexReco/interface/Vertex.h"
0030
0031 #include "RecoVertex/BeamSpotProducer/interface/FcnBeamSpotFitPV.h"
0032 #include "FWCore/Utilities/interface/isFinite.h"
0033
0034 #include "Minuit2/FCNBase.h"
0035 #include "Minuit2/FunctionMinimum.h"
0036 #include "Minuit2/MnMigrad.h"
0037 #include "Minuit2/MnPrint.h"
0038 #include "TF1.h"
0039 #include "TDirectory.h"
0040 #include "TMinuitMinimizer.h"
0041
0042 #include <iostream>
0043 #include <memory>
0044
0045 #include <sstream>
0046 using namespace std;
0047
0048 PVFitter::PVFitter(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iColl) : ftree_(nullptr) {
0049 initialize(iConfig, iColl);
0050 }
0051
0052 PVFitter::PVFitter(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iColl) : ftree_(nullptr) {
0053 initialize(iConfig, iColl);
0054 }
0055
0056 void PVFitter::initialize(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iColl) {
0057
0058
0059 TMinuitMinimizer::UseStaticMinuit(false);
0060 debug_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Debug");
0061 vertexToken_ = iColl.consumes<reco::VertexCollection>(
0062 iConfig.getParameter<edm::ParameterSet>("PVFitter")
0063 .getUntrackedParameter<edm::InputTag>("VertexCollection", edm::InputTag("offlinePrimaryVertices")));
0064 do3DFit_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Apply3DFit");
0065
0066
0067 maxNrVertices_ =
0068 iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("maxNrStoredVertices");
0069 minNrVertices_ =
0070 iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
0071 minVtxNdf_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
0072 maxVtxNormChi2_ =
0073 iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexNormChi2");
0074 minVtxTracks_ =
0075 iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minVertexNTracks");
0076 minVtxWgt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
0077 maxVtxR_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexR");
0078 maxVtxZ_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexZ");
0079 errorScale_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("errorScale");
0080 sigmaCut_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("nSigmaCut");
0081 fFitPerBunchCrossing =
0082 iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("FitPerBunchCrossing");
0083 useOnlyFirstPV_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("useOnlyFirstPV");
0084 minSumPt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minSumPt");
0085
0086
0087 dynamicQualityCut_ = 1.e30;
0088
0089 hPVx = std::make_unique<TH2F>("hPVx", "PVx vs PVz distribution", 200, -maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
0090 hPVy = std::make_unique<TH2F>("hPVy", "PVy vs PVz distribution", 200, -maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
0091 hPVx->SetDirectory(nullptr);
0092 hPVy->SetDirectory(nullptr);
0093 }
0094
0095 PVFitter::~PVFitter() {}
0096
0097 void PVFitter::fillDescription(edm::ParameterSetDescription& iDesc) {
0098 edm::ParameterSetDescription pvFitter;
0099
0100 pvFitter.addUntracked<bool>("Debug");
0101 pvFitter.addUntracked<edm::InputTag>("VertexCollection", edm::InputTag("offlinePrimaryVertices"));
0102 pvFitter.addUntracked<bool>("Apply3DFit");
0103 pvFitter.addUntracked<unsigned int>("maxNrStoredVertices");
0104 pvFitter.addUntracked<unsigned int>("minNrVerticesForFit");
0105 pvFitter.addUntracked<double>("minVertexNdf");
0106 pvFitter.addUntracked<double>("maxVertexNormChi2");
0107 pvFitter.addUntracked<unsigned int>("minVertexNTracks");
0108 pvFitter.addUntracked<double>("minVertexMeanWeight");
0109 pvFitter.addUntracked<double>("maxVertexR");
0110 pvFitter.addUntracked<double>("maxVertexZ");
0111 pvFitter.addUntracked<double>("errorScale");
0112 pvFitter.addUntracked<double>("nSigmaCut");
0113 pvFitter.addUntracked<bool>("FitPerBunchCrossing");
0114 pvFitter.addUntracked<bool>("useOnlyFirstPV");
0115 pvFitter.addUntracked<double>("minSumPt");
0116
0117 iDesc.add<edm::ParameterSetDescription>("PVFitter", pvFitter);
0118 }
0119
0120 void PVFitter::readEvent(const edm::Event& iEvent) {
0121
0122 edm::Handle<reco::VertexCollection> PVCollection;
0123 bool hasPVs = false;
0124
0125 if (iEvent.getByToken(vertexToken_, PVCollection)) {
0126 hasPVs = true;
0127 }
0128
0129
0130 if (hasPVs) {
0131 for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv) {
0132 if (useOnlyFirstPV_) {
0133 if (pv != PVCollection->begin())
0134 break;
0135 }
0136
0137
0138 if (pv->isFake() || pv->tracksSize() == 0)
0139 continue;
0140 if (pv->ndof() < minVtxNdf_ || (pv->ndof() + 3.) / pv->tracksSize() < 2 * minVtxWgt_)
0141 continue;
0142
0143
0144 if (pv->tracksSize() < minVtxTracks_)
0145 continue;
0146
0147 const auto& testTrack = pv->trackRefAt(0);
0148 if (testTrack.isNull() || !testTrack.isAvailable()) {
0149 edm::LogInfo("") << "Track collection not found. Skipping cut on sumPt.";
0150 } else {
0151 double sumPt = 0;
0152 for (auto iTrack = pv->tracks_begin(); iTrack != pv->tracks_end(); ++iTrack) {
0153 const auto pt = (*iTrack)->pt();
0154 sumPt += pt;
0155 }
0156 if (sumPt < minSumPt_)
0157 continue;
0158 }
0159
0160 hPVx->Fill(pv->x(), pv->z());
0161 hPVy->Fill(pv->y(), pv->z());
0162
0163
0164
0165
0166
0167 if (pvQuality(*pv) > dynamicQualityCut_)
0168 continue;
0169
0170 if (pvStore_.size() >= maxNrVertices_) {
0171 compressStore();
0172 if (pvQuality(*pv) > dynamicQualityCut_)
0173 continue;
0174 }
0175
0176
0177
0178 int bx = iEvent.bunchCrossing();
0179 BeamSpotFitPVData pvData;
0180 pvData.bunchCrossing = bx;
0181 pvData.position[0] = pv->x();
0182 pvData.position[1] = pv->y();
0183 pvData.position[2] = pv->z();
0184 pvData.posError[0] = pv->xError();
0185 pvData.posError[1] = pv->yError();
0186 pvData.posError[2] = pv->zError();
0187 pvData.posCorr[0] = pv->covariance(0, 1) / pv->xError() / pv->yError();
0188 pvData.posCorr[1] = pv->covariance(0, 2) / pv->xError() / pv->zError();
0189 pvData.posCorr[2] = pv->covariance(1, 2) / pv->yError() / pv->zError();
0190 pvStore_.push_back(pvData);
0191
0192 if (ftree_ != nullptr) {
0193 theBeamSpotTreeData_.run(iEvent.id().run());
0194 theBeamSpotTreeData_.lumi(iEvent.luminosityBlock());
0195 theBeamSpotTreeData_.bunchCrossing(bx);
0196 theBeamSpotTreeData_.pvData(pvData);
0197 ftree_->Fill();
0198 }
0199
0200 if (fFitPerBunchCrossing)
0201 bxMap_[bx].push_back(pvData);
0202 }
0203 }
0204 }
0205
0206 void PVFitter::setTree(TTree* tree) {
0207 ftree_ = tree;
0208 theBeamSpotTreeData_.branch(ftree_);
0209 }
0210
0211 bool PVFitter::runBXFitter() {
0212 using namespace ROOT::Minuit2;
0213 edm::LogInfo("PVFitter") << " Number of bunch crossings: " << bxMap_.size() << std::endl;
0214
0215 bool fit_ok = true;
0216
0217 for (std::map<int, std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin(); pvStore != bxMap_.end();
0218 ++pvStore) {
0219
0220
0221 fbspotMap[pvStore->first] = reco::BeamSpot();
0222
0223 edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << (pvStore->second).size()
0224 << " in bx: " << pvStore->first << std::endl;
0225
0226 if ((pvStore->second).size() <= minNrVertices_) {
0227 edm::LogWarning("PVFitter") << " not enough PVs, continue" << std::endl;
0228 fit_ok = false;
0229 continue;
0230 }
0231
0232 edm::LogInfo("PVFitter") << "Calculating beam spot with PVs ..." << std::endl;
0233
0234
0235
0236
0237 FcnBeamSpotFitPV fcn(pvStore->second);
0238
0239
0240
0241 MnUserParameters upar;
0242 upar.Add("x", 0., 0.02, -10., 10.);
0243 upar.Add("y", 0., 0.02, -10., 10.);
0244 upar.Add("z", 0., 0.20, -30., 30.);
0245 upar.Add("ex", 0.015, 0.01, 0., 10.);
0246 upar.Add("corrxy", 0., 0.02, -1., 1.);
0247 upar.Add("ey", 0.015, 0.01, 0., 10.);
0248 upar.Add("dxdz", 0., 0.0002, -0.1, 0.1);
0249 upar.Add("dydz", 0., 0.0002, -0.1, 0.1);
0250 upar.Add("ez", 1., 0.1, 0., 30.);
0251 upar.Add("scale", errorScale_, errorScale_ / 10., errorScale_ / 2., errorScale_ * 2.);
0252 MnMigrad migrad(fcn, upar);
0253
0254
0255
0256
0257 upar.Fix(4);
0258 upar.Fix(6);
0259 upar.Fix(7);
0260 upar.Fix(9);
0261 FunctionMinimum ierr = migrad(0, 1.);
0262 if (!ierr.IsValid()) {
0263 edm::LogInfo("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
0264 fit_ok = false;
0265 continue;
0266 }
0267
0268
0269
0270 fcn.setLimits(upar.Value(0) - sigmaCut_ * upar.Value(3),
0271 upar.Value(0) + sigmaCut_ * upar.Value(3),
0272 upar.Value(1) - sigmaCut_ * upar.Value(5),
0273 upar.Value(1) + sigmaCut_ * upar.Value(5),
0274 upar.Value(2) - sigmaCut_ * upar.Value(8),
0275 upar.Value(2) + sigmaCut_ * upar.Value(8));
0276 ierr = migrad(0, 1.);
0277 if (!ierr.IsValid()) {
0278 edm::LogInfo("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
0279 fit_ok = false;
0280 continue;
0281 }
0282
0283
0284
0285 upar.Release(4);
0286 upar.Release(6);
0287 upar.Release(7);
0288 ierr = migrad(0, 1.);
0289 if (!ierr.IsValid()) {
0290 edm::LogInfo("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
0291 fit_ok = false;
0292 continue;
0293 }
0294
0295 fwidthX = upar.Value(3);
0296 fwidthY = upar.Value(5);
0297 fwidthZ = upar.Value(8);
0298 fwidthXerr = upar.Error(3);
0299 fwidthYerr = upar.Error(5);
0300 fwidthZerr = upar.Error(8);
0301
0302 reco::BeamSpot::CovarianceMatrix matrix;
0303
0304 matrix(0, 0) = pow(upar.Error(0), 2);
0305 matrix(1, 1) = pow(upar.Error(1), 2);
0306 matrix(2, 2) = pow(upar.Error(2), 2);
0307 matrix(3, 3) = fwidthZerr * fwidthZerr;
0308 matrix(4, 4) = pow(upar.Error(6), 2);
0309 matrix(5, 5) = pow(upar.Error(7), 2);
0310 matrix(6, 6) = fwidthXerr * fwidthXerr;
0311
0312 fbeamspot = reco::BeamSpot(reco::BeamSpot::Point(upar.Value(0), upar.Value(1), upar.Value(2)),
0313 fwidthZ,
0314 upar.Value(6),
0315 upar.Value(7),
0316 fwidthX,
0317 matrix);
0318 fbeamspot.setBeamWidthX(fwidthX);
0319 fbeamspot.setBeamWidthY(fwidthY);
0320 fbeamspot.setType(reco::BeamSpot::Tracker);
0321
0322 fbspotMap[pvStore->first] = fbeamspot;
0323 edm::LogInfo("PVFitter") << "3D PV fit done for this bunch crossing." << std::endl;
0324 fit_ok = fit_ok & true;
0325 }
0326
0327 return fit_ok;
0328 }
0329
0330 bool PVFitter::runFitter() {
0331 using namespace ROOT::Minuit2;
0332 edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << pvStore_.size() << std::endl;
0333
0334 if (pvStore_.size() <= minNrVertices_)
0335 return false;
0336
0337
0338
0339 TDirectory::TContext guard{nullptr};
0340 std::ostringstream str;
0341 str << this;
0342 std::unique_ptr<TH1D> h1PVx{hPVx->ProjectionX((std::string("h1PVx") + str.str()).c_str(), 0, -1, "e")};
0343 std::unique_ptr<TH1D> h1PVy{hPVy->ProjectionX((std::string("h1PVy") + str.str()).c_str(), 0, -1, "e")};
0344 std::unique_ptr<TH1D> h1PVz{hPVx->ProjectionY((std::string("h1PVz") + str.str()).c_str(), 0, -1, "e")};
0345
0346
0347
0348 TF1 gausx("localGausX", "gaus", 0., 1., TF1::EAddToList::kNo);
0349 TF1 gausy("localGausY", "gaus", 0., 1., TF1::EAddToList::kNo);
0350 TF1 gausz("localGausZ", "gaus", 0., 1., TF1::EAddToList::kNo);
0351
0352 h1PVx->Fit(&gausx, "QLMN0 SERIAL");
0353 h1PVy->Fit(&gausy, "QLMN0 SERIAL");
0354 h1PVz->Fit(&gausz, "QLMN0 SERIAL");
0355
0356 fwidthX = gausx.GetParameter(2);
0357 fwidthY = gausy.GetParameter(2);
0358 fwidthZ = gausz.GetParameter(2);
0359 fwidthXerr = gausx.GetParError(2);
0360 fwidthYerr = gausy.GetParError(2);
0361 fwidthZerr = gausz.GetParError(2);
0362
0363 double estX = gausx.GetParameter(1);
0364 double estY = gausy.GetParameter(1);
0365 double estZ = gausz.GetParameter(1);
0366 double errX = fwidthX * 3.;
0367 double errY = fwidthY * 3.;
0368 double errZ = fwidthZ * 3.;
0369
0370 if (!do3DFit_) {
0371 reco::BeamSpot::CovarianceMatrix matrix;
0372 matrix(2, 2) = gausz.GetParError(1) * gausz.GetParError(1);
0373 matrix(3, 3) = fwidthZerr * fwidthZerr;
0374 matrix(6, 6) = fwidthXerr * fwidthXerr;
0375
0376 fbeamspot =
0377 reco::BeamSpot(reco::BeamSpot::Point(gausx.GetParameter(1), gausy.GetParameter(1), gausz.GetParameter(1)),
0378 fwidthZ,
0379 0.,
0380 0.,
0381 fwidthX,
0382 matrix);
0383 fbeamspot.setBeamWidthX(fwidthX);
0384 fbeamspot.setBeamWidthY(fwidthY);
0385 fbeamspot.setType(reco::BeamSpot::Tracker);
0386
0387 } else {
0388
0389
0390
0391 FcnBeamSpotFitPV fcn(pvStore_);
0392
0393
0394
0395 MnUserParameters upar;
0396 upar.Add("x", estX, errX, -10., 10.);
0397 upar.Add("y", estY, errY, -10., 10.);
0398 upar.Add("z", estZ, errZ, -30., 30.);
0399 upar.Add("ex", 0.015, 0.01, 0., 10.);
0400 upar.Add("corrxy", 0., 0.02, -1., 1.);
0401 upar.Add("ey", 0.015, 0.01, 0., 10.);
0402 upar.Add("dxdz", 0., 0.0002, -0.1, 0.1);
0403 upar.Add("dydz", 0., 0.0002, -0.1, 0.1);
0404 upar.Add("ez", 1., 0.1, 0., 30.);
0405 upar.Add("scale", errorScale_, errorScale_ / 10., errorScale_ / 2., errorScale_ * 2.);
0406 MnMigrad migrad(fcn, upar);
0407
0408 auto& state = const_cast<MnUserParameterState&>(migrad.State());
0409
0410
0411
0412
0413 state.Fix(4);
0414 state.Fix(6);
0415 state.Fix(7);
0416 state.Fix(9);
0417 FunctionMinimum ierr = migrad(0, 1.);
0418 if (!ierr.IsValid()) {
0419 edm::LogWarning("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
0420 return false;
0421 }
0422
0423
0424
0425
0426 vector<double> results;
0427 vector<double> errors;
0428 results = ierr.UserParameters().Params();
0429 errors = ierr.UserParameters().Errors();
0430
0431 fcn.setLimits(results[0] - sigmaCut_ * results[3],
0432 results[0] + sigmaCut_ * results[3],
0433 results[1] - sigmaCut_ * results[5],
0434 results[1] + sigmaCut_ * results[5],
0435 results[2] - sigmaCut_ * results[8],
0436 results[2] + sigmaCut_ * results[8]);
0437 ierr = migrad(0, 1.);
0438 if (!ierr.IsValid()) {
0439 edm::LogWarning("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
0440 return false;
0441 }
0442
0443
0444
0445 state.Release(4);
0446 state.Release(6);
0447 state.Release(7);
0448 ierr = migrad(0, 1.);
0449 if (!ierr.IsValid()) {
0450 edm::LogWarning("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
0451 return false;
0452 }
0453
0454 results = ierr.UserParameters().Params();
0455 errors = ierr.UserParameters().Errors();
0456
0457 fwidthX = results[3];
0458 fwidthY = results[5];
0459 fwidthZ = results[8];
0460 fwidthXerr = errors[3];
0461 fwidthYerr = errors[5];
0462 fwidthZerr = errors[8];
0463
0464
0465 if (edm::isNotFinite(fwidthXerr) || edm::isNotFinite(fwidthYerr) || edm::isNotFinite(fwidthZerr)) {
0466 edm::LogWarning("PVFitter") << "3D beam spot fit returns nan in 3rd iteration" << std::endl;
0467 return false;
0468 }
0469
0470 reco::BeamSpot::CovarianceMatrix matrix;
0471
0472 matrix(0, 0) = pow(errors[0], 2);
0473 matrix(1, 1) = pow(errors[1], 2);
0474 matrix(2, 2) = pow(errors[2], 2);
0475 matrix(3, 3) = fwidthZerr * fwidthZerr;
0476 matrix(6, 6) = fwidthXerr * fwidthXerr;
0477
0478 fbeamspot = reco::BeamSpot(
0479 reco::BeamSpot::Point(results[0], results[1], results[2]), fwidthZ, results[6], results[7], fwidthX, matrix);
0480 fbeamspot.setBeamWidthX(fwidthX);
0481 fbeamspot.setBeamWidthY(fwidthY);
0482 fbeamspot.setType(reco::BeamSpot::Tracker);
0483 }
0484
0485 return true;
0486 }
0487
0488 void PVFitter::dumpTxtFile() {}
0489
0490 void PVFitter::compressStore() {
0491
0492
0493
0494 pvQualities_.resize(pvStore_.size());
0495 for (unsigned int i = 0; i < pvStore_.size(); ++i)
0496 pvQualities_[i] = pvQuality(pvStore_[i]);
0497 sort(pvQualities_.begin(), pvQualities_.end());
0498
0499
0500
0501
0502
0503 dynamicQualityCut_ = pvQualities_[pvQualities_.size() / 2];
0504
0505
0506
0507
0508 unsigned int iwrite(0);
0509 for (unsigned int i = 0; i < pvStore_.size(); ++i) {
0510 if (pvQuality(pvStore_[i]) > dynamicQualityCut_)
0511 continue;
0512 if (i != iwrite)
0513 pvStore_[iwrite] = pvStore_[i];
0514 ++iwrite;
0515 }
0516 pvStore_.resize(iwrite);
0517 edm::LogInfo("PVFitter") << "Reduced primary vertex store size to " << pvStore_.size()
0518 << " ; new dynamic quality cut = " << dynamicQualityCut_ << std::endl;
0519 }
0520
0521 double PVFitter::pvQuality(const reco::Vertex& pv) const {
0522
0523
0524
0525 return pv.covariance(0, 0) * pv.covariance(1, 1) - pv.covariance(0, 1) * pv.covariance(0, 1);
0526 }
0527
0528 double PVFitter::pvQuality(const BeamSpotFitPVData& pv) const {
0529
0530
0531
0532 double ex = pv.posError[0];
0533 double ey = pv.posError[1];
0534 return ex * ex * ey * ey * (1 - pv.posCorr[0] * pv.posCorr[0]);
0535 }