PVFitter Class Reference

#include <PVFitter.h>

Public Member Functions
void	compressStore ()
	reduce size of primary vertex cache by increasing quality limit More...
void	dumpTxtFile ()
void	FitPerBunchCrossing ()
reco::BeamSpot	getBeamSpot ()
std::map< int, reco::BeamSpot >	getBeamSpotMap ()
int *	getFitLSRange ()
int	getNPVs ()
const std::map< int, int > &	getNPVsperBX ()
std::vector< BeamSpotFitPVData >	getpvStore ()
time_t *	getRefTime ()
double	getWidthX ()
double	getWidthXerr ()
double	getWidthY ()
double	getWidthYerr ()
double	getWidthZ ()
double	getWidthZerr ()
void	initialize (const edm::ParameterSet &iConfig, edm::ConsumesCollector &iColl)
bool	IsFitPerBunchCrossing ()
	PVFitter ()
	PVFitter (const edm::ParameterSet &iConfig, edm::ConsumesCollector &&iColl)
	PVFitter (const edm::ParameterSet &iConfig, edm::ConsumesCollector &iColl)
double	pvQuality (const reco::Vertex &pv) const
	vertex quality measure More...
double	pvQuality (const BeamSpotFitPVData &pv) const
	vertex quality measure More...
void	readEvent (const edm::Event &iEvent)
void	resetAll ()
void	resetLSRange ()
void	resetRefTime ()
void	resizepvStore (unsigned int rmSize)
bool	runBXFitter ()
bool	runFitter ()
void	setFitLSRange (int ls0, int ls1)
void	setRefTime (std::time_t t0, std::time_t t1)
void	setTree (TTree *tree)
virtual	~PVFitter ()

Private Attributes
std::map< int, std::vector< BeamSpotFitPVData > >	bxMap_
bool	debug_
bool	do3DFit_
double	dynamicQualityCut_
double	errorScale_
std::ofstream	fasciiFile
reco::BeamSpot	fbeamspot
int	fbeginLumiOfFit
std::map< int, reco::BeamSpot >	fbspotMap
int	fendLumiOfFit
bool	fFitPerBunchCrossing
std::time_t	freftime [2]
TTree *	ftree_
double	fwidthX
double	fwidthXerr
double	fwidthY
double	fwidthYerr
double	fwidthZ
double	fwidthZerr
TH2F *	hPVx
TH2F *	hPVy
unsigned int	maxNrVertices_
double	maxVtxNormChi2_
double	maxVtxR_
double	maxVtxZ_
unsigned int	minNrVertices_
double	minSumPt_
double	minVtxNdf_
unsigned int	minVtxTracks_
double	minVtxWgt_
std::map< int, int >	npvsmap_
std::string	outputTxt_
std::vector< double >	pvQualities_
std::vector< BeamSpotFitPVData >	pvStore_
double	sigmaCut_
BeamSpotTreeData	theBeamSpotTreeData_
bool	useOnlyFirstPV_
edm::EDGetTokenT< reco::VertexCollection >	vertexToken_

Detailed Description

Definition at line 44 of file PVFitter.h.

Constructor & Destructor Documentation

PVFitter::PVFitter ( )

inline

Definition at line 46 of file PVFitter.h.

References iEvent, initialize(), and AlcaBeamMonitor_cfi::PVFitter.

{}

PVFitter::PVFitter	(	const edm::ParameterSet &	iConfig,
		edm::ConsumesCollector &&	iColl
	)

Definition at line 46 of file PVFitter.cc.

References initialize().

                                                                               : ftree_(nullptr) {
  initialize(iConfig, iColl);
}

PVFitter::PVFitter	(	const edm::ParameterSet &	iConfig,
		edm::ConsumesCollector &	iColl
	)

Definition at line 50 of file PVFitter.cc.

References initialize().

                                                                              : ftree_(nullptr) {
  initialize(iConfig, iColl);
}

PVFitter::~PVFitter ( )

virtual

Definition at line 91 of file PVFitter.cc.

{}

Member Function Documentation

void PVFitter::compressStore

(

)

reduce size of primary vertex cache by increasing quality limit

Definition at line 461 of file PVFitter.cc.

References dynamicQualityCut_, mps_fire::i, pvQualities_, pvQuality(), and pvStore_.

Referenced by readEvent().

                             {
  //
  // fill vertex qualities
  //
  pvQualities_.resize(pvStore_.size());
  for (unsigned int i = 0; i < pvStore_.size(); ++i)
    pvQualities_[i] = pvQuality(pvStore_[i]);
  sort(pvQualities_.begin(), pvQualities_.end());
  //
  // Set new quality cut to median. This cut will be used to reduce the
  // number of vertices in the store and also apply to all new vertices
  // until the next reset
  //
  dynamicQualityCut_ = pvQualities_[pvQualities_.size() / 2];
  //
  // remove all vertices failing the cut from the store
  //   (to be moved to a more efficient memory management!)
  //
  unsigned int iwrite(0);
  for (unsigned int i = 0; i < pvStore_.size(); ++i) {
    if (pvQuality(pvStore_[i]) > dynamicQualityCut_)
      continue;
    if (i != iwrite)
      pvStore_[iwrite] = pvStore_[i];
    ++iwrite;
  }
  pvStore_.resize(iwrite);
  edm::LogInfo("PVFitter") << "Reduced primary vertex store size to " << pvStore_.size()
                           << " ; new dynamic quality cut = " << dynamicQualityCut_ << std::endl;
}

void PVFitter::dumpTxtFile

(

)

Definition at line 459 of file PVFitter.cc.

{}

void PVFitter::FitPerBunchCrossing ( )

inline

Definition at line 64 of file PVFitter.h.

{ fFitPerBunchCrossing = true; }

reco::BeamSpot PVFitter::getBeamSpot ( )

inline

Definition at line 91 of file PVFitter.h.

Referenced by AlcaBeamMonitor::dqmEndLuminosityBlock(), and BeamFitter::runPVandTrkFitter().

{ return fbeamspot; }

std::map<int, reco::BeamSpot> PVFitter::getBeamSpotMap ( )

inline

Definition at line 92 of file PVFitter.h.

Referenced by BeamFitter::runPVandTrkFitter().

{ return fbspotMap; }

int* PVFitter::getFitLSRange ( )

inline

Definition at line 94 of file PVFitter.h.

References createJobs::tmp.

                       {
    int* tmp = new int[2];
    tmp[0] = fbeginLumiOfFit;
    tmp[1] = fendLumiOfFit;
    return tmp;
  }

int PVFitter::getNPVs ( )

inline

Definition at line 117 of file PVFitter.h.

{ return pvStore_.size(); }

const std::map<int, int>& PVFitter::getNPVsperBX ( )

inline

Definition at line 119 of file PVFitter.h.

References genParticles_cff::map, and findQualityFiles::size.

                                         {
    npvsmap_.clear();

    for (std::map<int, std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin();
         pvStore != bxMap_.end();
         ++pvStore) {
      npvsmap_[pvStore->first] = (pvStore->second).size();
    }
    return npvsmap_;
  }

std::vector<BeamSpotFitPVData> PVFitter::getpvStore ( )

inline

Definition at line 62 of file PVFitter.h.

{ return pvStore_; }

time_t* PVFitter::getRefTime ( )

inline

Definition at line 101 of file PVFitter.h.

                       {
    time_t* tmptime = new time_t[2];
    tmptime[0] = freftime[0];
    tmptime[1] = freftime[1];
    return tmptime;
  }

double PVFitter::getWidthX ( )

inline

Definition at line 55 of file PVFitter.h.

{ return fwidthX; }

double PVFitter::getWidthXerr ( )

inline

Definition at line 58 of file PVFitter.h.

Referenced by BeamFitter::runPVandTrkFitter().

{ return fwidthXerr; }

double PVFitter::getWidthY ( )

inline

Definition at line 56 of file PVFitter.h.

{ return fwidthY; }

double PVFitter::getWidthYerr ( )

inline

Definition at line 59 of file PVFitter.h.

{ return fwidthYerr; }

double PVFitter::getWidthZ ( )

inline

Definition at line 57 of file PVFitter.h.

{ return fwidthZ; }

double PVFitter::getWidthZerr ( )

inline

Definition at line 60 of file PVFitter.h.

{ return fwidthZerr; }

void PVFitter::initialize	(	const edm::ParameterSet &	iConfig,
		edm::ConsumesCollector &	iColl
	)

Definition at line 54 of file PVFitter.cc.

References edm::ConsumesCollector::consumes(), debug_, do3DFit_, dynamicQualityCut_, errorScale_, fFitPerBunchCrossing, edm::ParameterSet::getParameter(), hPVx, hPVy, maxNrVertices_, maxVtxNormChi2_, maxVtxR_, maxVtxZ_, minNrVertices_, minSumPt_, minVtxNdf_, minVtxTracks_, minVtxWgt_, sigmaCut_, useOnlyFirstPV_, and vertexToken_.

Referenced by PVFitter().

                                                                                     {
  //In order to make fitting ROOT histograms thread safe
  // one must call this undocumented function
  TMinuitMinimizer::UseStaticMinuit(false);
  debug_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Debug");
  vertexToken_ = iColl.consumes<reco::VertexCollection>(
      iConfig.getParameter<edm::ParameterSet>("PVFitter")
          .getUntrackedParameter<edm::InputTag>("VertexCollection", edm::InputTag("offlinePrimaryVertices")));
  do3DFit_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("Apply3DFit");
  //writeTxt_          = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("WriteAscii");
  //outputTxt_         = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<std::string>("AsciiFileName");
  maxNrVertices_ =
      iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("maxNrStoredVertices");
  minNrVertices_ =
      iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minNrVerticesForFit");
  minVtxNdf_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexNdf");
  maxVtxNormChi2_ =
      iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexNormChi2");
  minVtxTracks_ =
      iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<unsigned int>("minVertexNTracks");
  minVtxWgt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minVertexMeanWeight");
  maxVtxR_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexR");
  maxVtxZ_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("maxVertexZ");
  errorScale_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("errorScale");
  sigmaCut_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("nSigmaCut");
  fFitPerBunchCrossing =
      iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("FitPerBunchCrossing");
  useOnlyFirstPV_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<bool>("useOnlyFirstPV");
  minSumPt_ = iConfig.getParameter<edm::ParameterSet>("PVFitter").getUntrackedParameter<double>("minSumPt");

  // preset quality cut to "infinite"
  dynamicQualityCut_ = 1.e30;

  hPVx = new TH2F("hPVx", "PVx vs PVz distribution", 200, -maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
  hPVy = new TH2F("hPVy", "PVy vs PVz distribution", 200, -maxVtxR_, maxVtxR_, 200, -maxVtxZ_, maxVtxZ_);
}

bool PVFitter::IsFitPerBunchCrossing ( )

inline

Definition at line 93 of file PVFitter.h.

Referenced by BeamFitter::dumpTxtFile(), and BeamFitter::runPVandTrkFitter().

{ return fFitPerBunchCrossing; }

double PVFitter::pvQuality

(

const reco::Vertex &

pv

)

const

vertex quality measure

Definition at line 492 of file PVFitter.cc.

References reco::Vertex::covariance().

Referenced by compressStore(), and readEvent().

                                                     {
  //
  // determinant of the transverse part of the PV covariance matrix
  //
  return pv.covariance(0, 0) * pv.covariance(1, 1) - pv.covariance(0, 1) * pv.covariance(0, 1);
}

double PVFitter::pvQuality

(

const BeamSpotFitPVData &

pv

)

const

vertex quality measure

Definition at line 499 of file PVFitter.cc.

References BeamSpotFitPVData::posCorr, and BeamSpotFitPVData::posError.

                                                            {
  //
  // determinant of the transverse part of the PV covariance matrix
  //
  double ex = pv.posError[0];
  double ey = pv.posError[1];
  return ex * ex * ey * ey * (1 - pv.posCorr[0] * pv.posCorr[0]);
}

void PVFitter::readEvent

(

const edm::Event &

iEvent

)

Definition at line 93 of file PVFitter.cc.

References BeamSpotFitPVData::bunchCrossing, BeamSpotTreeData::bunchCrossing(), edm::EventBase::bunchCrossing(), l1GtPatternGenerator_cfi::bx, bxMap_, compressStore(), dynamicQualityCut_, fFitPerBunchCrossing, ftree_, edm::Event::getByToken(), hPVx, hPVy, edm::EventBase::id(), BeamSpotTreeData::lumi(), edm::EventBase::luminosityBlock(), maxNrVertices_, minSumPt_, minVtxNdf_, minVtxTracks_, minVtxWgt_, BeamSpotFitPVData::posCorr, BeamSpotFitPVData::posError, BeamSpotFitPVData::position, DiDispStaMuonMonitor_cfi::pt, MetAnalyzer::pv(), L1TEGammaOffline_cfi::PVCollection, BeamSpotTreeData::pvData(), pvQuality(), pvStore_, BeamSpotTreeData::run(), edm::EventID::run(), TtFullHadEvtBuilder_cfi::sumPt, theBeamSpotTreeData_, useOnlyFirstPV_, and vertexToken_.

Referenced by AlcaBeamMonitor::analyze(), and BeamFitter::readEvent().

                                               {
  //------ Primary Vertices
  edm::Handle<reco::VertexCollection> PVCollection;
  bool hasPVs = false;

  if (iEvent.getByToken(vertexToken_, PVCollection)) {
    hasPVs = true;
  }
  //------

  if (hasPVs) {
    for (reco::VertexCollection::const_iterator pv = PVCollection->begin(); pv != PVCollection->end(); ++pv) {
      if (useOnlyFirstPV_) {
        if (pv != PVCollection->begin())
          break;
      }

      //--- vertex selection
      if (pv->isFake() || pv->tracksSize() == 0)
        continue;
      if (pv->ndof() < minVtxNdf_ || (pv->ndof() + 3.) / pv->tracksSize() < 2 * minVtxWgt_)
        continue;
      //---

      if (pv->tracksSize() < minVtxTracks_)
        continue;

      const auto& testTrack = pv->trackRefAt(0);
      if (testTrack.isNull() || !testTrack.isAvailable()) {
        edm::LogInfo("") << "Track collection not found. Skipping cut on sumPt.";
      } else {
        double sumPt = 0;
        for (auto iTrack = pv->tracks_begin(); iTrack != pv->tracks_end(); ++iTrack) {
          const auto pt = (*iTrack)->pt();
          sumPt += pt;
        }
        if (sumPt < minSumPt_)
          continue;
      }

      hPVx->Fill(pv->x(), pv->z());
      hPVy->Fill(pv->y(), pv->z());

      //
      // 3D fit section
      //
      // apply additional quality cut
      if (pvQuality(*pv) > dynamicQualityCut_)
        continue;
      // if store exceeds max. size: reduce size and apply new quality cut
      if (pvStore_.size() >= maxNrVertices_) {
        compressStore();
        if (pvQuality(*pv) > dynamicQualityCut_)
          continue;
      }
      //
      // copy PV to store
      //
      int bx = iEvent.bunchCrossing();
      BeamSpotFitPVData pvData;
      pvData.bunchCrossing = bx;
      pvData.position[0] = pv->x();
      pvData.position[1] = pv->y();
      pvData.position[2] = pv->z();
      pvData.posError[0] = pv->xError();
      pvData.posError[1] = pv->yError();
      pvData.posError[2] = pv->zError();
      pvData.posCorr[0] = pv->covariance(0, 1) / pv->xError() / pv->yError();
      pvData.posCorr[1] = pv->covariance(0, 2) / pv->xError() / pv->zError();
      pvData.posCorr[2] = pv->covariance(1, 2) / pv->yError() / pv->zError();
      pvStore_.push_back(pvData);

      if (ftree_ != nullptr) {
        theBeamSpotTreeData_.run(iEvent.id().run());
        theBeamSpotTreeData_.lumi(iEvent.luminosityBlock());
        theBeamSpotTreeData_.bunchCrossing(bx);
        theBeamSpotTreeData_.pvData(pvData);
        ftree_->Fill();
      }

      if (fFitPerBunchCrossing)
        bxMap_[bx].push_back(pvData);
    }
  }
}

void PVFitter::resetAll ( )

inline

Definition at line 80 of file PVFitter.h.

References align::BeamSpot.

Referenced by BeamFitter::BeamFitter(), and AlcaBeamMonitor::dqmEndLuminosityBlock().

                  {
    resetLSRange();
    resetRefTime();
    pvStore_.clear();
    bxMap_.clear();
    dynamicQualityCut_ = 1.e30;
    hPVx->Reset();
    hPVy->Reset();
    fbeamspot = reco::BeamSpot();
    fbspotMap.clear();
  };

void PVFitter::resetLSRange ( )

inline

Definition at line 67 of file PVFitter.h.

{ fbeginLumiOfFit = fendLumiOfFit = -1; }

void PVFitter::resetRefTime ( )

inline

Definition at line 68 of file PVFitter.h.

{ freftime[0] = freftime[1] = 0; }

void PVFitter::resizepvStore ( unsigned int  rmSize )

inline

Definition at line 109 of file PVFitter.h.

References MetAnalyzer::pv().

{ pvStore_.erase(pvStore_.begin(), pvStore_.begin() + rmSize); }

bool PVFitter::runBXFitter

(

)

Definition at line 184 of file PVFitter.cc.

References align::BeamSpot, bxMap_, errorScale_, fbeamspot, fbspotMap, fcn(), fwidthX, fwidthXerr, fwidthY, fwidthYerr, fwidthZ, fwidthZerr, genParticles_cff::map, makeMuonMisalignmentScenario::matrix, minNrVertices_, funct::pow(), reco::BeamSpot::setBeamWidthX(), reco::BeamSpot::setBeamWidthY(), FcnBeamSpotFitPV::setLimits(), reco::BeamSpot::setType(), sigmaCut_, findQualityFiles::size, and reco::BeamSpot::Tracker.

Referenced by BeamFitter::runPVandTrkFitter().

                           {
  using namespace ROOT::Minuit2;
  edm::LogInfo("PVFitter") << " Number of bunch crossings: " << bxMap_.size() << std::endl;

  bool fit_ok = true;

  for (std::map<int, std::vector<BeamSpotFitPVData> >::const_iterator pvStore = bxMap_.begin(); pvStore != bxMap_.end();
       ++pvStore) {
    // first set null beam spot in case
    // fit fails
    fbspotMap[pvStore->first] = reco::BeamSpot();

    edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << (pvStore->second).size()
                             << " in bx: " << pvStore->first << std::endl;

    if ((pvStore->second).size() <= minNrVertices_) {
      edm::LogWarning("PVFitter") << " not enough PVs, continue" << std::endl;
      fit_ok = false;
      continue;
    }

    edm::LogInfo("PVFitter") << "Calculating beam spot with PVs ..." << std::endl;

    //
    // LL function and fitter
    //
    FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore->second);
    //
    // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
    //
    MnUserParameters upar;
    upar.Add("x", 0., 0.02, -10., 10.);                                                     // 0
    upar.Add("y", 0., 0.02, -10., 10.);                                                     // 1
    upar.Add("z", 0., 0.20, -30., 30.);                                                     // 2
    upar.Add("ex", 0.015, 0.01, 0., 10.);                                                   // 3
    upar.Add("corrxy", 0., 0.02, -1., 1.);                                                  // 4
    upar.Add("ey", 0.015, 0.01, 0., 10.);                                                   // 5
    upar.Add("dxdz", 0., 0.0002, -0.1, 0.1);                                                // 6
    upar.Add("dydz", 0., 0.0002, -0.1, 0.1);                                                // 7
    upar.Add("ez", 1., 0.1, 0., 30.);                                                       // 8
    upar.Add("scale", errorScale_, errorScale_ / 10., errorScale_ / 2., errorScale_ * 2.);  // 9
    MnMigrad migrad(*fcn, upar);

    //
    // first iteration without correlations
    //
    upar.Fix(4);
    upar.Fix(6);
    upar.Fix(7);
    upar.Fix(9);
    FunctionMinimum ierr = migrad(0, 1.);
    if (!ierr.IsValid()) {
      edm::LogInfo("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
      fit_ok = false;
      continue;
    }
    //
    // refit with harder selection on vertices
    //
    fcn->setLimits(upar.Value(0) - sigmaCut_ * upar.Value(3),
                   upar.Value(0) + sigmaCut_ * upar.Value(3),
                   upar.Value(1) - sigmaCut_ * upar.Value(5),
                   upar.Value(1) + sigmaCut_ * upar.Value(5),
                   upar.Value(2) - sigmaCut_ * upar.Value(8),
                   upar.Value(2) + sigmaCut_ * upar.Value(8));
    ierr = migrad(0, 1.);
    if (!ierr.IsValid()) {
      edm::LogInfo("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
      fit_ok = false;
      continue;
    }
    //
    // refit with correlations
    //
    upar.Release(4);
    upar.Release(6);
    upar.Release(7);
    ierr = migrad(0, 1.);
    if (!ierr.IsValid()) {
      edm::LogInfo("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
      fit_ok = false;
      continue;
    }

    fwidthX = upar.Value(3);
    fwidthY = upar.Value(5);
    fwidthZ = upar.Value(8);
    fwidthXerr = upar.Error(3);
    fwidthYerr = upar.Error(5);
    fwidthZerr = upar.Error(8);

    reco::BeamSpot::CovarianceMatrix matrix;
    // need to get the full cov matrix
    matrix(0, 0) = pow(upar.Error(0), 2);
    matrix(1, 1) = pow(upar.Error(1), 2);
    matrix(2, 2) = pow(upar.Error(2), 2);
    matrix(3, 3) = fwidthZerr * fwidthZerr;
    matrix(4, 4) = pow(upar.Error(6), 2);
    matrix(5, 5) = pow(upar.Error(7), 2);
    matrix(6, 6) = fwidthXerr * fwidthXerr;

    fbeamspot = reco::BeamSpot(reco::BeamSpot::Point(upar.Value(0), upar.Value(1), upar.Value(2)),
                               fwidthZ,
                               upar.Value(6),
                               upar.Value(7),
                               fwidthX,
                               matrix);
    fbeamspot.setBeamWidthX(fwidthX);
    fbeamspot.setBeamWidthY(fwidthY);
    fbeamspot.setType(reco::BeamSpot::Tracker);

    fbspotMap[pvStore->first] = fbeamspot;
    edm::LogInfo("PVFitter") << "3D PV fit done for this bunch crossing." << std::endl;
    //delete fcn;
    fit_ok = fit_ok & true;
  }

  return fit_ok;
}

bool PVFitter::runFitter

(

)

Definition at line 304 of file PVFitter.cc.

References align::BeamSpot, do3DFit_, MessageLogger_cfi::errors, errorScale_, fftjetvertexadder_cfi::errX, fftjetvertexadder_cfi::errY, fftjetvertexadder_cfi::errZ, fbeamspot, fcn(), fwidthX, fwidthXerr, fwidthY, fwidthYerr, fwidthZ, fwidthZerr, hPVx, hPVy, edm::isNotFinite(), makeMuonMisalignmentScenario::matrix, minNrVertices_, funct::pow(), pvStore_, bookConverter::results, reco::BeamSpot::setBeamWidthX(), reco::BeamSpot::setBeamWidthY(), FcnBeamSpotFitPV::setLimits(), reco::BeamSpot::setType(), sigmaCut_, str, AlCaHLTBitMon_QueryRunRegistry::string, and reco::BeamSpot::Tracker.

Referenced by AlcaBeamMonitor::dqmEndLuminosityBlock(), and BeamFitter::runPVandTrkFitter().

                         {
  using namespace ROOT::Minuit2;
  edm::LogInfo("PVFitter") << " Number of PVs collected for PVFitter: " << pvStore_.size() << std::endl;

  if (pvStore_.size() <= minNrVertices_)
    return false;

  //need to create a unique histogram name to avoid ROOT trying to re-use one
  // also tell ROOT to hide its global state
  TDirectory::TContext guard{nullptr};
  std::ostringstream str;
  str << this;
  std::unique_ptr<TH1D> h1PVx{hPVx->ProjectionX((std::string("h1PVx") + str.str()).c_str(), 0, -1, "e")};
  std::unique_ptr<TH1D> h1PVy{hPVy->ProjectionX((std::string("h1PVy") + str.str()).c_str(), 0, -1, "e")};
  std::unique_ptr<TH1D> h1PVz{hPVx->ProjectionY((std::string("h1PVz") + str.str()).c_str(), 0, -1, "e")};

  //Use our own copy for thread safety
  // also do not add to global list of functions
  TF1 gausx("localGausX", "gaus", 0., 1., TF1::EAddToList::kNo);
  TF1 gausy("localGausY", "gaus", 0., 1., TF1::EAddToList::kNo);
  TF1 gausz("localGausZ", "gaus", 0., 1., TF1::EAddToList::kNo);

  h1PVx->Fit(&gausx, "QLMN0 SERIAL");
  h1PVy->Fit(&gausy, "QLMN0 SERIAL");
  h1PVz->Fit(&gausz, "QLMN0 SERIAL");

  fwidthX = gausx.GetParameter(2);
  fwidthY = gausy.GetParameter(2);
  fwidthZ = gausz.GetParameter(2);
  fwidthXerr = gausx.GetParError(2);
  fwidthYerr = gausy.GetParError(2);
  fwidthZerr = gausz.GetParError(2);

  double estX = gausx.GetParameter(1);
  double estY = gausy.GetParameter(1);
  double estZ = gausz.GetParameter(1);
  double errX = fwidthX * 3.;
  double errY = fwidthY * 3.;
  double errZ = fwidthZ * 3.;

  if (!do3DFit_) {
    reco::BeamSpot::CovarianceMatrix matrix;
    matrix(2, 2) = gausz.GetParError(1) * gausz.GetParError(1);
    matrix(3, 3) = fwidthZerr * fwidthZerr;
    matrix(6, 6) = fwidthXerr * fwidthXerr;

    fbeamspot =
        reco::BeamSpot(reco::BeamSpot::Point(gausx.GetParameter(1), gausy.GetParameter(1), gausz.GetParameter(1)),
                       fwidthZ,
                       0.,
                       0.,
                       fwidthX,
                       matrix);
    fbeamspot.setBeamWidthX(fwidthX);
    fbeamspot.setBeamWidthY(fwidthY);
    fbeamspot.setType(reco::BeamSpot::Tracker);

  } else {  // do 3D fit
    //
    // LL function and fitter
    //
    FcnBeamSpotFitPV* fcn = new FcnBeamSpotFitPV(pvStore_);
    //
    // fit parameters: positions, widths, x-y correlations, tilts in xz and yz
    //
    MnUserParameters upar;
    upar.Add("x", estX, errX, -10., 10.);                                                   // 0
    upar.Add("y", estY, errY, -10., 10.);                                                   // 1
    upar.Add("z", estZ, errZ, -30., 30.);                                                   // 2
    upar.Add("ex", 0.015, 0.01, 0., 10.);                                                   // 3
    upar.Add("corrxy", 0., 0.02, -1., 1.);                                                  // 4
    upar.Add("ey", 0.015, 0.01, 0., 10.);                                                   // 5
    upar.Add("dxdz", 0., 0.0002, -0.1, 0.1);                                                // 6
    upar.Add("dydz", 0., 0.0002, -0.1, 0.1);                                                // 7
    upar.Add("ez", 1., 0.1, 0., 30.);                                                       // 8
    upar.Add("scale", errorScale_, errorScale_ / 10., errorScale_ / 2., errorScale_ * 2.);  // 9
    MnMigrad migrad(*fcn, upar);
    //
    // first iteration without correlations
    //
    migrad.Fix(4);
    migrad.Fix(6);
    migrad.Fix(7);
    migrad.Fix(9);
    FunctionMinimum ierr = migrad(0, 1.);
    if (!ierr.IsValid()) {
      edm::LogWarning("PVFitter") << "3D beam spot fit failed in 1st iteration" << std::endl;
      return false;
    }
    //
    // refit with harder selection on vertices
    //

    vector<double> results;
    vector<double> errors;
    results = ierr.UserParameters().Params();
    errors = ierr.UserParameters().Errors();

    fcn->setLimits(results[0] - sigmaCut_ * results[3],
                   results[0] + sigmaCut_ * results[3],
                   results[1] - sigmaCut_ * results[5],
                   results[1] + sigmaCut_ * results[5],
                   results[2] - sigmaCut_ * results[8],
                   results[2] + sigmaCut_ * results[8]);
    ierr = migrad(0, 1.);
    if (!ierr.IsValid()) {
      edm::LogWarning("PVFitter") << "3D beam spot fit failed in 2nd iteration" << std::endl;
      return false;
    }
    //
    // refit with correlations
    //
    migrad.Release(4);
    migrad.Release(6);
    migrad.Release(7);
    ierr = migrad(0, 1.);
    if (!ierr.IsValid()) {
      edm::LogWarning("PVFitter") << "3D beam spot fit failed in 3rd iteration" << std::endl;
      return false;
    }

    results = ierr.UserParameters().Params();
    errors = ierr.UserParameters().Errors();

    fwidthX = results[3];
    fwidthY = results[5];
    fwidthZ = results[8];
    fwidthXerr = errors[3];
    fwidthYerr = errors[5];
    fwidthZerr = errors[8];

    // check errors on widths and sigmaZ for nan
    if (edm::isNotFinite(fwidthXerr) || edm::isNotFinite(fwidthYerr) || edm::isNotFinite(fwidthZerr)) {
      edm::LogWarning("PVFitter") << "3D beam spot fit returns nan in 3rd iteration" << std::endl;
      return false;
    }

    reco::BeamSpot::CovarianceMatrix matrix;
    // need to get the full cov matrix
    matrix(0, 0) = pow(errors[0], 2);
    matrix(1, 1) = pow(errors[1], 2);
    matrix(2, 2) = pow(errors[2], 2);
    matrix(3, 3) = fwidthZerr * fwidthZerr;
    matrix(6, 6) = fwidthXerr * fwidthXerr;

    fbeamspot = reco::BeamSpot(
        reco::BeamSpot::Point(results[0], results[1], results[2]), fwidthZ, results[6], results[7], fwidthX, matrix);
    fbeamspot.setBeamWidthX(fwidthX);
    fbeamspot.setBeamWidthY(fwidthY);
    fbeamspot.setType(reco::BeamSpot::Tracker);
  }

  return true;
}

void PVFitter::setFitLSRange ( int  ls0, int  ls1  )

inline

Definition at line 74 of file PVFitter.h.

                                       {
    fbeginLumiOfFit = ls0;
    fendLumiOfFit = ls1;
  }

void PVFitter::setRefTime ( std::time_t  t0, std::time_t  t1  )

inline

Definition at line 70 of file PVFitter.h.

References FrontierCondition_GT_autoExpress_cfi::t0, and RandomServiceHelper::t1.

                                              {
    freftime[0] = t0;
    freftime[1] = t1;
  }

void PVFitter::setTree

(

TTree *

tree

)

Definition at line 179 of file PVFitter.cc.

References BeamSpotTreeData::branch(), ftree_, theBeamSpotTreeData_, and MainPageGenerator::tree.

Referenced by BeamFitter::BeamFitter().

                                  {
  ftree_ = tree;
  theBeamSpotTreeData_.branch(ftree_);
}

Member Data Documentation

std::map<int, std::vector<BeamSpotFitPVData> > PVFitter::bxMap_

private

Definition at line 174 of file PVFitter.h.

Referenced by readEvent(), and runBXFitter().

bool PVFitter::debug_

private

Definition at line 139 of file PVFitter.h.

Referenced by initialize().

bool PVFitter::do3DFit_

private

Definition at line 140 of file PVFitter.h.

Referenced by initialize(), and runFitter().

double PVFitter::dynamicQualityCut_

private

Definition at line 175 of file PVFitter.h.

Referenced by compressStore(), initialize(), and readEvent().

double PVFitter::errorScale_

private

Definition at line 152 of file PVFitter.h.

Referenced by initialize(), runBXFitter(), and runFitter().

std::ofstream PVFitter::fasciiFile

private

Definition at line 137 of file PVFitter.h.

reco::BeamSpot PVFitter::fbeamspot

private

Definition at line 132 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

int PVFitter::fbeginLumiOfFit

private

Definition at line 164 of file PVFitter.h.

std::map<int, reco::BeamSpot> PVFitter::fbspotMap

private

Definition at line 133 of file PVFitter.h.

Referenced by runBXFitter().

int PVFitter::fendLumiOfFit

private

Definition at line 165 of file PVFitter.h.

bool PVFitter::fFitPerBunchCrossing

private

Definition at line 134 of file PVFitter.h.

Referenced by initialize(), and readEvent().

std::time_t PVFitter::freftime[2]

private

Definition at line 156 of file PVFitter.h.

TTree* PVFitter::ftree_

private

Definition at line 161 of file PVFitter.h.

Referenced by readEvent(), and setTree().

double PVFitter::fwidthX

private

Definition at line 166 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

double PVFitter::fwidthXerr

private

Definition at line 169 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

double PVFitter::fwidthY

private

Definition at line 167 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

double PVFitter::fwidthYerr

private

Definition at line 170 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

double PVFitter::fwidthZ

private

Definition at line 168 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

double PVFitter::fwidthZerr

private

Definition at line 171 of file PVFitter.h.

Referenced by runBXFitter(), and runFitter().

TH2F* PVFitter::hPVx

private

Definition at line 158 of file PVFitter.h.

Referenced by initialize(), readEvent(), and runFitter().

TH2F* PVFitter::hPVy

private

Definition at line 159 of file PVFitter.h.

Referenced by initialize(), readEvent(), and runFitter().

unsigned int PVFitter::maxNrVertices_

private

Definition at line 144 of file PVFitter.h.

Referenced by initialize(), and readEvent().

double PVFitter::maxVtxNormChi2_

private

Definition at line 147 of file PVFitter.h.

Referenced by initialize().

double PVFitter::maxVtxR_

private

Definition at line 150 of file PVFitter.h.

Referenced by initialize().

double PVFitter::maxVtxZ_

private

Definition at line 151 of file PVFitter.h.

Referenced by initialize().

unsigned int PVFitter::minNrVertices_

private

Definition at line 145 of file PVFitter.h.

Referenced by initialize(), runBXFitter(), and runFitter().

double PVFitter::minSumPt_

private

Definition at line 154 of file PVFitter.h.

Referenced by initialize(), and readEvent().

double PVFitter::minVtxNdf_

private

Definition at line 146 of file PVFitter.h.

Referenced by initialize(), and readEvent().

unsigned int PVFitter::minVtxTracks_

private

Definition at line 148 of file PVFitter.h.

Referenced by initialize(), and readEvent().

double PVFitter::minVtxWgt_

private

Definition at line 149 of file PVFitter.h.

Referenced by initialize(), and readEvent().

std::map<int, int> PVFitter::npvsmap_

private

Definition at line 131 of file PVFitter.h.

std::string PVFitter::outputTxt_

private

Definition at line 142 of file PVFitter.h.

std::vector<double> PVFitter::pvQualities_

private

Definition at line 176 of file PVFitter.h.

Referenced by compressStore().

std::vector<BeamSpotFitPVData> PVFitter::pvStore_

private

Definition at line 173 of file PVFitter.h.

Referenced by compressStore(), readEvent(), and runFitter().

double PVFitter::sigmaCut_

private

Definition at line 153 of file PVFitter.h.

Referenced by initialize(), runBXFitter(), and runFitter().

BeamSpotTreeData PVFitter::theBeamSpotTreeData_

private

Definition at line 178 of file PVFitter.h.

Referenced by readEvent(), and setTree().

bool PVFitter::useOnlyFirstPV_

private

Definition at line 135 of file PVFitter.h.

Referenced by initialize(), and readEvent().

edm::EDGetTokenT<reco::VertexCollection> PVFitter::vertexToken_

private

Definition at line 141 of file PVFitter.h.

Referenced by initialize(), and readEvent().