34 #include "Minuit2/FCNBase.h"
35 #include "Minuit2/FunctionMinimum.h"
36 #include "Minuit2/MnMigrad.h"
37 #include "Minuit2/MnPrint.h"
39 #include "TDirectory.h"
40 #include "TMinuitMinimizer.h"
57 TMinuitMinimizer::UseStaticMinuit(
false);
61 .getUntrackedParameter<edm::InputTag>(
"VertexCollection",
edm::InputTag(
"offlinePrimaryVertices")));
87 hPVx = std::unique_ptr<TH2F>(
89 hPVy = std::unique_ptr<TH2F>(
91 hPVx->SetDirectory(
nullptr);
92 hPVy->SetDirectory(
nullptr);
115 if (
pv->isFake() ||
pv->tracksSize() == 0)
124 const auto& testTrack =
pv->trackRefAt(0);
125 if (testTrack.isNull() || !testTrack.isAvailable()) {
126 edm::LogInfo(
"") <<
"Track collection not found. Skipping cut on sumPt.";
129 for (
auto iTrack =
pv->tracks_begin(); iTrack !=
pv->tracks_end(); ++iTrack) {
130 const auto pt = (*iTrack)->pt();
164 pvData.
posCorr[0] =
pv->covariance(0, 1) /
pv->xError() /
pv->yError();
165 pvData.
posCorr[1] =
pv->covariance(0, 2) /
pv->xError() /
pv->zError();
166 pvData.
posCorr[2] =
pv->covariance(1, 2) /
pv->yError() /
pv->zError();
190 edm::LogInfo(
"PVFitter") <<
" Number of bunch crossings: " <<
bxMap_.size() << std::endl;
194 for (
std::map<
int, std::vector<BeamSpotFitPVData> >::const_iterator pvStore =
bxMap_.begin(); pvStore !=
bxMap_.end();
200 edm::LogInfo(
"PVFitter") <<
" Number of PVs collected for PVFitter: " << (pvStore->second).
size()
201 <<
" in bx: " << pvStore->first << std::endl;
204 edm::LogWarning(
"PVFitter") <<
" not enough PVs, continue" << std::endl;
209 edm::LogInfo(
"PVFitter") <<
"Calculating beam spot with PVs ..." << std::endl;
218 MnUserParameters upar;
219 upar.Add(
"x", 0., 0.02, -10., 10.);
220 upar.Add(
"y", 0., 0.02, -10., 10.);
221 upar.Add(
"z", 0., 0.20, -30., 30.);
222 upar.Add(
"ex", 0.015, 0.01, 0., 10.);
223 upar.Add(
"corrxy", 0., 0.02, -1., 1.);
224 upar.Add(
"ey", 0.015, 0.01, 0., 10.);
225 upar.Add(
"dxdz", 0., 0.0002, -0.1, 0.1);
226 upar.Add(
"dydz", 0., 0.0002, -0.1, 0.1);
227 upar.Add(
"ez", 1., 0.1, 0., 30.);
229 MnMigrad migrad(
fcn, upar);
238 FunctionMinimum ierr = migrad(0, 1.);
239 if (!ierr.IsValid()) {
240 edm::LogInfo(
"PVFitter") <<
"3D beam spot fit failed in 1st iteration" << std::endl;
247 fcn.setLimits(upar.Value(0) -
sigmaCut_ * upar.Value(3),
248 upar.Value(0) +
sigmaCut_ * upar.Value(3),
249 upar.Value(1) -
sigmaCut_ * upar.Value(5),
250 upar.Value(1) +
sigmaCut_ * upar.Value(5),
251 upar.Value(2) -
sigmaCut_ * upar.Value(8),
252 upar.Value(2) +
sigmaCut_ * upar.Value(8));
253 ierr = migrad(0, 1.);
254 if (!ierr.IsValid()) {
255 edm::LogInfo(
"PVFitter") <<
"3D beam spot fit failed in 2nd iteration" << std::endl;
265 ierr = migrad(0, 1.);
266 if (!ierr.IsValid()) {
267 edm::LogInfo(
"PVFitter") <<
"3D beam spot fit failed in 3rd iteration" << std::endl;
300 edm::LogInfo(
"PVFitter") <<
"3D PV fit done for this bunch crossing." << std::endl;
301 fit_ok = fit_ok &
true;
309 edm::LogInfo(
"PVFitter") <<
" Number of PVs collected for PVFitter: " <<
pvStore_.size() << std::endl;
316 TDirectory::TContext guard{
nullptr};
317 std::ostringstream
str;
319 std::unique_ptr<TH1D> h1PVx{
hPVx->ProjectionX((
std::string(
"h1PVx") +
str.str()).c_str(), 0, -1,
"e")};
320 std::unique_ptr<TH1D> h1PVy{
hPVy->ProjectionX((
std::string(
"h1PVy") +
str.str()).c_str(), 0, -1,
"e")};
321 std::unique_ptr<TH1D> h1PVz{
hPVx->ProjectionY((
std::string(
"h1PVz") +
str.str()).c_str(), 0, -1,
"e")};
325 TF1 gausx(
"localGausX",
"gaus", 0., 1., TF1::EAddToList::kNo);
326 TF1 gausy(
"localGausY",
"gaus", 0., 1., TF1::EAddToList::kNo);
327 TF1 gausz(
"localGausZ",
"gaus", 0., 1., TF1::EAddToList::kNo);
329 h1PVx->Fit(&gausx,
"QLMN0 SERIAL");
330 h1PVy->Fit(&gausy,
"QLMN0 SERIAL");
331 h1PVz->Fit(&gausz,
"QLMN0 SERIAL");
333 fwidthX = gausx.GetParameter(2);
334 fwidthY = gausy.GetParameter(2);
335 fwidthZ = gausz.GetParameter(2);
340 double estX = gausx.GetParameter(1);
341 double estY = gausy.GetParameter(1);
342 double estZ = gausz.GetParameter(1);
349 matrix(2, 2) = gausz.GetParError(1) * gausz.GetParError(1);
372 MnUserParameters upar;
373 upar.Add(
"x", estX,
errX, -10., 10.);
374 upar.Add(
"y", estY,
errY, -10., 10.);
375 upar.Add(
"z", estZ,
errZ, -30., 30.);
376 upar.Add(
"ex", 0.015, 0.01, 0., 10.);
377 upar.Add(
"corrxy", 0., 0.02, -1., 1.);
378 upar.Add(
"ey", 0.015, 0.01, 0., 10.);
379 upar.Add(
"dxdz", 0., 0.0002, -0.1, 0.1);
380 upar.Add(
"dydz", 0., 0.0002, -0.1, 0.1);
381 upar.Add(
"ez", 1., 0.1, 0., 30.);
383 MnMigrad migrad(
fcn, upar);
391 FunctionMinimum ierr = migrad(0, 1.);
392 if (!ierr.IsValid()) {
393 edm::LogWarning(
"PVFitter") <<
"3D beam spot fit failed in 1st iteration" << std::endl;
402 results = ierr.UserParameters().Params();
403 errors = ierr.UserParameters().Errors();
411 ierr = migrad(0, 1.);
412 if (!ierr.IsValid()) {
413 edm::LogWarning(
"PVFitter") <<
"3D beam spot fit failed in 2nd iteration" << std::endl;
422 ierr = migrad(0, 1.);
423 if (!ierr.IsValid()) {
424 edm::LogWarning(
"PVFitter") <<
"3D beam spot fit failed in 3rd iteration" << std::endl;
428 results = ierr.UserParameters().Params();
429 errors = ierr.UserParameters().Errors();
440 edm::LogWarning(
"PVFitter") <<
"3D beam spot fit returns nan in 3rd iteration" << std::endl;
482 unsigned int iwrite(0);
483 for (
unsigned int i = 0;
i <
pvStore_.size(); ++
i) {
499 return pv.covariance(0, 0) *
pv.covariance(1, 1) -
pv.covariance(0, 1) *
pv.covariance(0, 1);
506 double ex =
pv.posError[0];
507 double ey =
pv.posError[1];
508 return ex * ex * ey * ey * (1 -
pv.posCorr[0] *
pv.posCorr[0]);