Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:23:15 

0001 /**_________________________________________________________________
0002    class:   PVFitter.cc
0003    package: RecoVertex/BeamSpotProducer
0004 
0005 
0006 
0007    author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov)
0008            Geng-Yuan Jeng, UC Riverside (Geng-Yuan.Jeng@cern.ch)
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   //In order to make fitting ROOT histograms thread safe
0058   // one must call this undocumented function
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   //writeTxt_          = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("WriteAscii");
0066   //outputTxt_         = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<std::string>("AsciiFileName");
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   // preset quality cut to "infinite"
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   //------ Primary Vertices
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       //--- vertex selection
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       // 3D fit section
0165       //
0166       // apply additional quality cut
0167       if (pvQuality(*pv) > dynamicQualityCut_)
0168         continue;
0169       // if store exceeds max. size: reduce size and apply new quality cut
0170       if (pvStore_.size() >= maxNrVertices_) {
0171         compressStore();
0172         if (pvQuality(*pv) > dynamicQualityCut_)
0173           continue;
0174       }
0175       //
0176       // copy PV to store
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     // first set null beam spot in case
0220     // fit fails
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     // LL function and fitter
0236     //
0237     FcnBeamSpotFitPV fcn(pvStore->second);
0238     //
0239     // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
0240     //
0241     MnUserParameters upar;
0242     upar.Add("x", 0., 0.02, -10., 10.);                                                     // 0
0243     upar.Add("y", 0., 0.02, -10., 10.);                                                     // 1
0244     upar.Add("z", 0., 0.20, -30., 30.);                                                     // 2
0245     upar.Add("ex", 0.015, 0.01, 0., 10.);                                                   // 3
0246     upar.Add("corrxy", 0., 0.02, -1., 1.);                                                  // 4
0247     upar.Add("ey", 0.015, 0.01, 0., 10.);                                                   // 5
0248     upar.Add("dxdz", 0., 0.0002, -0.1, 0.1);                                                // 6
0249     upar.Add("dydz", 0., 0.0002, -0.1, 0.1);                                                // 7
0250     upar.Add("ez", 1., 0.1, 0., 30.);                                                       // 8
0251     upar.Add("scale", errorScale_, errorScale_ / 10., errorScale_ / 2., errorScale_ * 2.);  // 9
0252     MnMigrad migrad(fcn, upar);
0253 
0254     //
0255     // first iteration without correlations
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     // refit with harder selection on vertices
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     // refit with correlations
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     // need to get the full cov matrix
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   //need to create a unique histogram name to avoid ROOT trying to re-use one
0338   // also tell ROOT to hide its global state
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   //Use our own copy for thread safety
0347   // also do not add to global list of functions
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 {  // do 3D fit
0388     //
0389     // LL function and fitter
0390     //
0391     FcnBeamSpotFitPV fcn(pvStore_);
0392     //
0393     // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
0394     //
0395     MnUserParameters upar;
0396     upar.Add("x", estX, errX, -10., 10.);                                                   // 0
0397     upar.Add("y", estY, errY, -10., 10.);                                                   // 1
0398     upar.Add("z", estZ, errZ, -30., 30.);                                                   // 2
0399     upar.Add("ex", 0.015, 0.01, 0., 10.);                                                   // 3
0400     upar.Add("corrxy", 0., 0.02, -1., 1.);                                                  // 4
0401     upar.Add("ey", 0.015, 0.01, 0., 10.);                                                   // 5
0402     upar.Add("dxdz", 0., 0.0002, -0.1, 0.1);                                                // 6
0403     upar.Add("dydz", 0., 0.0002, -0.1, 0.1);                                                // 7
0404     upar.Add("ez", 1., 0.1, 0., 30.);                                                       // 8
0405     upar.Add("scale", errorScale_, errorScale_ / 10., errorScale_ / 2., errorScale_ * 2.);  // 9
0406     MnMigrad migrad(fcn, upar);
0407     //
0408     // first iteration without correlations
0409     //
0410     migrad.Fix(4);
0411     migrad.Fix(6);
0412     migrad.Fix(7);
0413     migrad.Fix(9);
0414     FunctionMinimum ierr = migrad(0, 1.);
0415     if (!ierr.IsValid()) {
0416       edm::LogWarning("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
0417       return false;
0418     }
0419     //
0420     // refit with harder selection on vertices
0421     //
0422 
0423     vector<double> results;
0424     vector<double> errors;
0425     results = ierr.UserParameters().Params();
0426     errors = ierr.UserParameters().Errors();
0427 
0428     fcn.setLimits(results[0] - sigmaCut_ * results[3],
0429                   results[0] + sigmaCut_ * results[3],
0430                   results[1] - sigmaCut_ * results[5],
0431                   results[1] + sigmaCut_ * results[5],
0432                   results[2] - sigmaCut_ * results[8],
0433                   results[2] + sigmaCut_ * results[8]);
0434     ierr = migrad(0, 1.);
0435     if (!ierr.IsValid()) {
0436       edm::LogWarning("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
0437       return false;
0438     }
0439     //
0440     // refit with correlations
0441     //
0442     migrad.Release(4);
0443     migrad.Release(6);
0444     migrad.Release(7);
0445     ierr = migrad(0, 1.);
0446     if (!ierr.IsValid()) {
0447       edm::LogWarning("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
0448       return false;
0449     }
0450 
0451     results = ierr.UserParameters().Params();
0452     errors = ierr.UserParameters().Errors();
0453 
0454     fwidthX = results[3];
0455     fwidthY = results[5];
0456     fwidthZ = results[8];
0457     fwidthXerr = errors[3];
0458     fwidthYerr = errors[5];
0459     fwidthZerr = errors[8];
0460 
0461     // check errors on widths and sigmaZ for nan
0462     if (edm::isNotFinite(fwidthXerr) || edm::isNotFinite(fwidthYerr) || edm::isNotFinite(fwidthZerr)) {
0463       edm::LogWarning("PVFitter") << "3D beam spot fit returns nan in 3rd iteration" << std::endl;
0464       return false;
0465     }
0466 
0467     reco::BeamSpot::CovarianceMatrix matrix;
0468     // need to get the full cov matrix
0469     matrix(0, 0) = pow(errors[0], 2);
0470     matrix(1, 1) = pow(errors[1], 2);
0471     matrix(2, 2) = pow(errors[2], 2);
0472     matrix(3, 3) = fwidthZerr * fwidthZerr;
0473     matrix(6, 6) = fwidthXerr * fwidthXerr;
0474 
0475     fbeamspot = reco::BeamSpot(
0476         reco::BeamSpot::Point(results[0], results[1], results[2]), fwidthZ, results[6], results[7], fwidthX, matrix);
0477     fbeamspot.setBeamWidthX(fwidthX);
0478     fbeamspot.setBeamWidthY(fwidthY);
0479     fbeamspot.setType(reco::BeamSpot::Tracker);
0480   }
0481 
0482   return true;
0483 }
0484 
0485 void PVFitter::dumpTxtFile() {}
0486 
0487 void PVFitter::compressStore() {
0488   //
0489   // fill vertex qualities
0490   //
0491   pvQualities_.resize(pvStore_.size());
0492   for (unsigned int i = 0; i < pvStore_.size(); ++i)
0493     pvQualities_[i] = pvQuality(pvStore_[i]);
0494   sort(pvQualities_.begin(), pvQualities_.end());
0495   //
0496   // Set new quality cut to median. This cut will be used to reduce the
0497   // number of vertices in the store and also apply to all new vertices
0498   // until the next reset
0499   //
0500   dynamicQualityCut_ = pvQualities_[pvQualities_.size() / 2];
0501   //
0502   // remove all vertices failing the cut from the store
0503   //   (to be moved to a more efficient memory management!)
0504   //
0505   unsigned int iwrite(0);
0506   for (unsigned int i = 0; i < pvStore_.size(); ++i) {
0507     if (pvQuality(pvStore_[i]) > dynamicQualityCut_)
0508       continue;
0509     if (i != iwrite)
0510       pvStore_[iwrite] = pvStore_[i];
0511     ++iwrite;
0512   }
0513   pvStore_.resize(iwrite);
0514   edm::LogInfo("PVFitter") << "Reduced primary vertex store size to " << pvStore_.size()
0515                            << " ; new dynamic quality cut = " << dynamicQualityCut_ << std::endl;
0516 }
0517 
0518 double PVFitter::pvQuality(const reco::Vertex& pv) const {
0519   //
0520   // determinant of the transverse part of the PV covariance matrix
0521   //
0522   return pv.covariance(0, 0) * pv.covariance(1, 1) - pv.covariance(0, 1) * pv.covariance(0, 1);
0523 }
0524 
0525 double PVFitter::pvQuality(const BeamSpotFitPVData& pv) const {
0526   //
0527   // determinant of the transverse part of the PV covariance matrix
0528   //
0529   double ex = pv.posError[0];
0530   double ey = pv.posError[1];
0531   return ex * ex * ey * ey * (1 - pv.posCorr[0] * pv.posCorr[0]);
0532 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!