Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BeamSpot] fit 1D projections of PV distribution to get initial values for 3D fit parameters #10782

Merged
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
62 changes: 34 additions & 28 deletions RecoVertex/BeamSpotProducer/src/PVFitter.cc
Expand Up @@ -90,16 +90,15 @@ void PVFitter::initialize(const edm::ParameterSet& iConfig,

// 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_);
}

PVFitter::~PVFitter() {

}


void PVFitter::readEvent(const edm::Event& iEvent)
{

Expand Down Expand Up @@ -147,7 +146,6 @@ void PVFitter::readEvent(const edm::Event& iEvent)
if ( pv->isFake() || pv->tracksSize()==0 ) continue;
if ( pv->ndof() < minVtxNdf_ || (pv->ndof()+3.)/pv->tracksSize()<2*minVtxWgt_ ) continue;
//---

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

Expand Down Expand Up @@ -341,28 +339,36 @@ bool PVFitter::runFitter() {

//bool fit_ok = false;

if ( ! do3DFit_ ) {
TH1F *h1PVx = (TH1F*) hPVx->ProjectionX("h1PVx", 0, -1, "e");
TH1F *h1PVy = (TH1F*) hPVy->ProjectionX("h1PVy", 0, -1, "e");
TH1F *h1PVz = (TH1F*) hPVx->ProjectionY("h1PVz", 0, -1, "e");
TH1F *h1PVx = (TH1F*) hPVx->ProjectionX("h1PVx", 0, -1, "e");
TH1F *h1PVy = (TH1F*) hPVy->ProjectionX("h1PVy", 0, -1, "e");
TH1F *h1PVz = (TH1F*) hPVx->ProjectionY("h1PVz", 0, -1, "e");

//Use our own copy for thread safety
TF1 gaus("localGaus","gaus");

h1PVx->Fit(&gaus,"QLM0");
h1PVy->Fit(&gaus,"QLM0");
h1PVz->Fit(&gaus,"QLM0");

TF1 *gausx = h1PVx->GetFunction("localGaus");
TF1 *gausy = h1PVy->GetFunction("localGaus");
TF1 *gausz = h1PVz->GetFunction("localGaus");

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.;

//Use our own copy for thread safety
TF1 gaus("localGaus","gaus");

h1PVx->Fit(&gaus,"QLM0");
h1PVy->Fit(&gaus,"QLM0");
h1PVz->Fit(&gaus,"QLM0");

TF1 *gausx = h1PVx->GetFunction("localGaus");
TF1 *gausy = h1PVy->GetFunction("localGaus");
TF1 *gausz = h1PVz->GetFunction("localGaus");

fwidthX = gausx->GetParameter(2);
fwidthY = gausy->GetParameter(2);
fwidthZ = gausz->GetParameter(2);
fwidthXerr = gausx->GetParError(2);
fwidthYerr = gausy->GetParError(2);
fwidthZerr = gausz->GetParError(2);
if ( ! do3DFit_ ) {

reco::BeamSpot::CovarianceMatrix matrix;
matrix(2,2) = gausz->GetParError(1) * gausz->GetParError(1);
Expand Down Expand Up @@ -390,9 +396,9 @@ bool PVFitter::runFitter() {
// 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("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
Expand Down