|
|
|
#include <BSFitter.h>
_________________________________________________________________ class: BSFitter.h package: RecoVertex/BeamSpotProducer
author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov)
version
________________________________________________________________
Definition at line 32 of file BSFitter.h.
| BSFitter::BSFitter | ( | ) |
Definition at line 42 of file BSFitter.cc.
References Unknown.
{
fbeamtype = reco::BeamSpot::Unknown;
}
| BSFitter::BSFitter | ( | std::vector< BSTrkParameters > | BSvector | ) |
Definition at line 47 of file BSFitter.cc.
References Pi, and mathSSE::sqrt().
{
ffit_type = "default";
ffit_variable = "default";
fBSvector = BSvector;
fsqrt2pi = sqrt(2.* TMath::Pi());
fpar_name[0] = "z0 ";
fpar_name[1] = "SigmaZ0 ";
fpar_name[2] = "X0 ";
fpar_name[3] = "Y0 ";
fpar_name[4] = "dxdz ";
fpar_name[5] = "dydz ";
fpar_name[6] = "SigmaBeam ";
//if (theGausszFcn == 0 ) {
thePDF = new BSpdfsFcn();
//}
//if (theFitter == 0 ) {
theFitter = new VariableMetricMinimizer();
//}
fapplyd0cut = false;
fapplychi2cut = false;
ftmprow = 0;
ftmp.ResizeTo(4,1);
ftmp.Zero();
fnthite=0;
fMaxZ = 50.; //cm
fconvergence = 0.5; // stop fit when 50% of the input collection has been removed.
fminNtrks = 100;
finputBeamWidth = -1; // no input
h1z = new TH1F("h1z","z distribution",200,-fMaxZ, fMaxZ);
}
| BSFitter::~BSFitter | ( | ) | [virtual] |
Definition at line 91 of file BSFitter.cc.
| void BSFitter::d0phi_Init | ( | ) | [inline] |
| reco::BeamSpot BSFitter::Fit | ( | ) |
Definition at line 100 of file BSFitter.cc.
References FitTarget::Fit.
Referenced by BeamFitter::runAllFitter(), BeamFitter::runBeamWidthFitter(), and BeamFitter::runFitterNoTxt().
{
return this->Fit(0);
}
| reco::BeamSpot BSFitter::Fit | ( | double * | inipar = 0 | ) |
Definition at line 107 of file BSFitter.cc.
References reco::BeamSpot::BeamWidthX(), reco::BeamSpot::covariance(), reco::BeamSpot::dxdz(), reco::BeamSpot::dydz(), Exception, edm::detail::isnan(), funct::pow(), asciidump::s, reco::BeamSpot::setType(), reco::BeamSpot::sigmaZ(), Unknown, reco::BeamSpot::Unknown, reco::BeamSpot::x0(), reco::BeamSpot::y0(), and reco::BeamSpot::z0().
{
fbeamtype = reco::BeamSpot::Unknown;
if ( ffit_variable == "z" ) {
if ( ffit_type == "chi2" ) {
return Fit_z_chi2(inipar);
} else if ( ffit_type == "likelihood" ) {
return Fit_z_likelihood(inipar);
} else if ( ffit_type == "combined" ) {
reco::BeamSpot tmp_beamspot = Fit_z_chi2(inipar);
double tmp_par[2] = {tmp_beamspot.z0(), tmp_beamspot.sigmaZ()};
return Fit_z_likelihood(tmp_par);
} else {
throw cms::Exception("LogicError")
<< "Error in BeamSpotProducer/BSFitter: "
<< "Illegal fit type, options are chi2,likelihood or combined(ie. first chi2 then likelihood)";
}
} else if ( ffit_variable == "d" ) {
if ( ffit_type == "d0phi" ) {
this->d0phi_Init();
return Fit_d0phi();
} else if ( ffit_type == "likelihood" ) {
return Fit_d_likelihood(inipar);
} else if ( ffit_type == "combined" ) {
this->d0phi_Init();
reco::BeamSpot tmp_beamspot = Fit_d0phi();
double tmp_par[4] = {tmp_beamspot.x0(), tmp_beamspot.y0(), tmp_beamspot.dxdz(), tmp_beamspot.dydz()};
return Fit_d_likelihood(tmp_par);
} else {
throw cms::Exception("LogicError")
<< "Error in BeamSpotProducer/BSFitter: "
<< "Illegal fit type, options are d0phi, likelihood or combined";
}
} else if ( ffit_variable == "d*z" || ffit_variable == "default" ) {
if ( ffit_type == "likelihood" || ffit_type == "default" ) {
reco::BeamSpot::CovarianceMatrix matrix;
// we are now fitting Z inside d0phi fitter
// first fit z distribution using a chi2 fit
//reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
//for (int j = 2 ; j < 4 ; ++j) {
//for(int k = j ; k < 4 ; ++k) {
// matrix(j,k) = tmp_z.covariance()(j,k);
//}
//}
// use d0-phi algorithm to extract transverse position
this->d0phi_Init();
//reco::BeamSpot tmp_d0phi= Fit_d0phi(); // change to iterative procedure:
this->Setd0Cut_d0phi(4.0);
reco::BeamSpot tmp_d0phi= Fit_ited0phi();
//for (int j = 0 ; j < 2 ; ++j) {
// for(int k = j ; k < 2 ; ++k) {
// matrix(j,k) = tmp_d0phi.covariance()(j,k);
//}
//}
// slopes
//for (int j = 4 ; j < 6 ; ++j) {
// for(int k = j ; k < 6 ; ++k) {
// matrix(j,k) = tmp_d0phi.covariance()(j,k);
// }
//}
// put everything into one object
reco::BeamSpot spot(reco::BeamSpot::Point(tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_d0phi.z0()),
tmp_d0phi.sigmaZ(),
tmp_d0phi.dxdz(),
tmp_d0phi.dydz(),
0.,
tmp_d0phi.covariance(),
fbeamtype );
//reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
//reco::BeamSpot tmp_d0phi = Fit_d0phi();
// log-likelihood fit
if (ffit_type == "likelihood") {
double tmp_par[7] = {tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_d0phi.z0(),
tmp_d0phi.sigmaZ(), tmp_d0phi.dxdz(), tmp_d0phi.dydz(),0.0};
double tmp_error_par[7];
for(int s=0;s<6;s++){ tmp_error_par[s] = pow( tmp_d0phi.covariance()(s,s),0.5);}
tmp_error_par[6]=0.0;
reco::BeamSpot tmp_lh = Fit_d_z_likelihood(tmp_par,tmp_error_par);
if ( isnan(ff_minimum) || std::isinf(ff_minimum) ) {
edm::LogWarning("BSFitter") << "BSFitter: Result is non physical. Log-Likelihood fit to extract beam width did not converge." << std::endl;
tmp_lh.setType(reco::BeamSpot::Unknown);
return tmp_lh;
}
return tmp_lh;
} else {
edm::LogInfo("BSFitter") << "default track-based fit does not extract beam width." << std::endl;
return spot;
}
} else if ( ffit_type == "resolution" ) {
reco::BeamSpot tmp_z = Fit_z_chi2(inipar);
this->d0phi_Init();
reco::BeamSpot tmp_d0phi = Fit_d0phi();
double tmp_par[7] = {tmp_d0phi.x0(), tmp_d0phi.y0(), tmp_z.z0(),
tmp_z.sigmaZ(), tmp_d0phi.dxdz(), tmp_d0phi.dydz(),0.0};
double tmp_error_par[7];
for(int s=0;s<6;s++){ tmp_error_par[s] = pow(tmp_par[s],0.5);}
tmp_error_par[6]=0.0;
reco::BeamSpot tmp_beam = Fit_d_z_likelihood(tmp_par,tmp_error_par);
double tmp_par2[7] = {tmp_beam.x0(), tmp_beam.y0(), tmp_beam.z0(),
tmp_beam.sigmaZ(), tmp_beam.dxdz(), tmp_beam.dydz(),
tmp_beam.BeamWidthX()};
reco::BeamSpot tmp_lh = Fit_dres_z_likelihood(tmp_par2);
if ( isnan(ff_minimum) || std::isinf(ff_minimum) ) {
edm::LogWarning("BSFitter") << "Result is non physical. Log-Likelihood fit did not converge." << std::endl;
tmp_lh.setType(reco::BeamSpot::Unknown);
return tmp_lh;
}
return tmp_lh;
} else {
throw cms::Exception("LogicError")
<< "Error in BeamSpotProducer/BSFitter: "
<< "Illegal fit type, options are likelihood or resolution";
}
} else {
throw cms::Exception("LogicError")
<< "Error in BeamSpotProducer/BSFitter: "
<< "Illegal variable type, options are \"z\", \"d\", or \"d*z\"";
}
}
| reco::BeamSpot BSFitter::Fit_d0phi | ( | ) |
Definition at line 447 of file BSFitter.cc.
References abs, b, align::BeamSpot, funct::cos(), ExpressReco_HICollisions_FallBack::e, g, j, gen::k, funct::sin(), mathSSE::sqrt(), and cond::rpcobtemp::temp.
Referenced by BeamFitter::runAllFitter().
{
//LogDebug ("BSFitter") << " we will use " << fBSvector.size() << " tracks.";
if (fnthite > 0) edm::LogInfo("BSFitter") << " number of tracks used: " << ftmprow << std::endl;
//std::cout << " ftmp = matrix("<<ftmp.GetNrows()<<","<<ftmp.GetNcols()<<")"<<std::endl;
//std::cout << " ftmp(0,0)="<<ftmp(0,0)<<std::endl;
//std::cout << " ftmp(1,0)="<<ftmp(1,0)<<std::endl;
//std::cout << " ftmp(2,0)="<<ftmp(2,0)<<std::endl;
//std::cout << " ftmp(3,0)="<<ftmp(3,0)<<std::endl;
h1z->Reset();
TMatrixD x_result(4,1);
TMatrixDSym V_result(4);
TMatrixDSym Vint(4);
TMatrixD b(4,1);
//Double_t weightsum = 0;
Vint.Zero();
b.Zero();
TMatrixD g(4,1);
TMatrixDSym temp(4);
std::vector<BSTrkParameters>::iterator iparam = fBSvector.begin();
ftmprow=0;
//edm::LogInfo ("BSFitter") << " test";
//std::cout << "BSFitter: fit" << std::endl;
for( iparam = fBSvector.begin() ;
iparam != fBSvector.end() ; ++iparam) {
//if(i->weight2 == 0) continue;
//if (ftmprow==0) {
//std::cout << "d0=" << iparam->d0() << " sigd0=" << iparam->sigd0()
//<< " phi0="<< iparam->phi0() << " z0=" << iparam->z0() << std::endl;
//std::cout << "d0phi_d0=" << iparam->d0phi_d0() << " d0phi_chi2="<<iparam->d0phi_chi2() << std::endl;
//}
g(0,0) = sin(iparam->phi0());
g(1,0) = -cos(iparam->phi0());
g(2,0) = iparam->z0() * g(0,0);
g(3,0) = iparam->z0() * g(1,0);
// average transverse beam width
double sigmabeam2 = 0.006 * 0.006;
if (finputBeamWidth > 0 ) sigmabeam2 = finputBeamWidth * finputBeamWidth;
else { edm::LogWarning("BSFitter") << "using in fit beam width = " << sqrt(sigmabeam2) << std::endl; }
//double sigma2 = sigmabeam2 + (iparam->sigd0())* (iparam->sigd0()) / iparam->weight2;
// this should be 2*sigmabeam2?
double sigma2 = sigmabeam2 + (iparam->sigd0())* (iparam->sigd0());
TMatrixD ftmptrans(1,4);
ftmptrans = ftmptrans.Transpose(ftmp);
TMatrixD dcor = ftmptrans * g;
double chi2tmp = (iparam->d0() - dcor(0,0)) * (iparam->d0() - dcor(0,0))/sigma2;
(*iparam) = BSTrkParameters(iparam->z0(),iparam->sigz0(),iparam->d0(),iparam->sigd0(),
iparam->phi0(), iparam->pt(),dcor(0,0),chi2tmp);
bool pass = true;
if (fapplyd0cut && fnthite>0 ) {
if ( std::abs(iparam->d0() - dcor(0,0)) > fd0cut ) pass = false;
}
if (fapplychi2cut && fnthite>0 ) {
if ( chi2tmp > fchi2cut ) pass = false;
}
if (pass) {
temp.Zero();
for(int j = 0 ; j < 4 ; ++j) {
for(int k = j ; k < 4 ; ++k) {
temp(j,k) += g(j,0) * g(k,0);
}
}
Vint += (temp * (1 / sigma2));
b += (iparam->d0() / sigma2 * g);
//weightsum += sqrt(i->weight2);
ftmprow++;
h1z->Fill( iparam->z0() );
}
}
Double_t determinant;
TDecompBK bk(Vint);
bk.SetTol(1e-11); //FIXME: find a better way to solve x_result
if (!bk.Decompose()) {
goodfit = false;
edm::LogWarning("BSFitter")
<< "Decomposition failed, matrix singular ?" << std::endl
<< "condition number = " << bk.Condition() << std::endl;
}
else {
V_result = Vint.InvertFast(&determinant);
x_result = V_result * b;
}
// for(int j = 0 ; j < 4 ; ++j) {
// for(int k = 0 ; k < 4 ; ++k) {
// std::cout<<"V_result("<<j<<","<<k<<")="<<V_result(j,k)<<std::endl;
// }
// }
// for (int j=0;j<4;++j){
// std::cout<<"x_result("<<j<<",0)="<<x_result(j,0)<<std::endl;
// }
//LogDebug ("BSFitter") << " d0-phi fit done.";
//std::cout<< " d0-phi fit done." << std::endl;
h1z->Fit("gaus","QLM0","",h1z->GetMean() -2.*h1z->GetRMS(),h1z->GetMean() +2.*h1z->GetRMS());
//std::cout << "fitted "<< std::endl;
TF1 *fgaus = h1z->GetFunction("gaus");
//std::cout << "got function" << std::endl;
if (!fgaus){
edm::LogError("NoBeamSpotFit")<<"gaussian fit failed. no BS d0 fit";
return reco::BeamSpot();
}
double fpar[2] = {fgaus->GetParameter(1), fgaus->GetParameter(2) };
reco::BeamSpot::CovarianceMatrix matrix;
// first two parameters
for (int j = 0 ; j < 2 ; ++j) {
for(int k = j ; k < 2 ; ++k) {
matrix(j,k) = V_result(j,k);
}
}
// slope parameters
for (int j = 4 ; j < 6 ; ++j) {
for(int k = j ; k < 6 ; ++k) {
matrix(j,k) = V_result(j-2,k-2);
}
}
// Z0 and sigmaZ
matrix(2,2) = fgaus->GetParError(1) * fgaus->GetParError(1);
matrix(3,3) = fgaus->GetParError(2) * fgaus->GetParError(2);
ftmp = x_result;
// x0 and y0 are *not* x,y at z=0, but actually at z=0
// to correct for this, we need to translate them to z=z0
// using the measured slopes
//
double x0tmp = x_result(0,0);
double y0tmp = x_result(1,0);
x0tmp += x_result(2,0)*fpar[0];
y0tmp += x_result(3,0)*fpar[0];
return reco::BeamSpot( reco::BeamSpot::Point(x0tmp,
y0tmp,
fpar[0]),
fpar[1],
x_result(2,0),
x_result(3,0),
0.,
matrix,
fbeamtype );
}
| reco::BeamSpot BSFitter::Fit_d_likelihood | ( | double * | inipar | ) |
Definition at line 641 of file BSFitter.cc.
References align::BeamSpot, j, and gen::k.
{
thePDF->SetPDFs("PDFGauss_d");
thePDF->SetData(fBSvector);
MnUserParameters upar;
upar.Add("X0", inipar[0],0.001);
upar.Add("Y0", inipar[1],0.001);
upar.Add("Z0", 0.,0.001);
upar.Add("sigmaZ",0.,0.001);
upar.Add("dxdz",inipar[2],0.001);
upar.Add("dydz",inipar[3],0.001);
MnMigrad migrad(*thePDF, upar);
FunctionMinimum fmin = migrad();
ff_minimum = fmin.Fval();
reco::BeamSpot::CovarianceMatrix matrix;
for (int j = 0 ; j < 6 ; ++j) {
for(int k = j ; k < 6 ; ++k) {
matrix(j,k) = fmin.Error().Matrix()(j,k);
}
}
return reco::BeamSpot( reco::BeamSpot::Point(fmin.Parameters().Vec()(0),
fmin.Parameters().Vec()(1),
0.),
0.,
fmin.Parameters().Vec()(4),
fmin.Parameters().Vec()(5),
0.,
matrix,
fbeamtype );
}
| reco::BeamSpot BSFitter::Fit_d_z_likelihood | ( | double * | inipar, |
| double * | error_par | ||
| ) |
Definition at line 770 of file BSFitter.cc.
References align::BeamSpot, j, and gen::k.
{
int tracksFailed=0;
//estimate first guess of beam width and tame 20% extra of it to start
inipar[6]=scanPDF(inipar,tracksFailed,1);
error_par[6]=(inipar[6])*0.20;
//Here remove the tracks which give low pdf and fill into a new vector
//std::cout<<"Size of Old vector = "<<(fBSvector.size())<<std::endl;
/* double junk= */ scanPDF(inipar,tracksFailed,2);
//std::cout<<"Size of New vector = "<<(fBSvectorBW.size())<<std::endl;
//Refill the fBSVector again with new sets of tracks
fBSvector.clear();
std::vector<BSTrkParameters>::const_iterator iparamBW = fBSvectorBW.begin();
for( iparamBW = fBSvectorBW.begin(); iparamBW != fBSvectorBW.end(); ++iparamBW)
{ fBSvector.push_back(*iparamBW);
}
thePDF->SetPDFs("PDFGauss_d*PDFGauss_z");
thePDF->SetData(fBSvector);
MnUserParameters upar;
upar.Add("X0", inipar[0],error_par[0]);
upar.Add("Y0", inipar[1],error_par[1]);
upar.Add("Z0", inipar[2],error_par[2]);
upar.Add("sigmaZ",inipar[3],error_par[3]);
upar.Add("dxdz",inipar[4],error_par[4]);
upar.Add("dydz",inipar[5],error_par[5]);
upar.Add("BeamWidthX",inipar[6],error_par[6]);
MnMigrad migrad(*thePDF, upar);
FunctionMinimum fmin = migrad();
// std::cout<<"-----how the fit evoves------"<<std::endl;
// std::cout<<fmin<<std::endl;
ff_minimum = fmin.Fval();
bool ff_nfcn=fmin.HasReachedCallLimit();
bool ff_cov=fmin.HasCovariance();
bool testing=fmin.IsValid();
//Print WARNINGS if minimum did not converged
if( ! testing )
{
edm::LogWarning("BSFitter") <<"===========>>>>>** WARNING: MINUIT DID NOT CONVERGES PROPERLY !!!!!!"<<std::endl;
if(ff_nfcn) edm::LogWarning("BSFitter") <<"===========>>>>>** WARNING: No. of Calls Exhausted"<<std::endl;
if(!ff_cov) edm::LogWarning("BSFitter") <<"===========>>>>>** WARNING: Covariance did not found"<<std::endl;
}
edm::LogInfo("BSFitter") <<"The Total # Tracks used for beam width fit = "<<(fBSvectorBW.size())<<std::endl;
//Checks after fit is performed
double lastIter_pars[7];
for(int ip=0;ip<7;ip++){ lastIter_pars[ip]=fmin.Parameters().Vec()(ip);
}
tracksFailed=0;
/* double lastIter_scan= */ scanPDF(lastIter_pars,tracksFailed,2);
edm::LogWarning("BSFitter") <<"WARNING: # of tracks which have very low pdf value (pdf_d < 1.0e-05) are = "<<tracksFailed<<std::endl;
//std::cout << " eval= " << ff_minimum
// << "/n params[0]= " << fmin.Parameters().Vec()(0) << std::endl;
reco::BeamSpot::CovarianceMatrix matrix;
for (int j = 0 ; j < 7 ; ++j) {
for(int k = j ; k < 7 ; ++k) {
matrix(j,k) = fmin.Error().Matrix()(j,k);
}
}
return reco::BeamSpot( reco::BeamSpot::Point(fmin.Parameters().Vec()(0),
fmin.Parameters().Vec()(1),
fmin.Parameters().Vec()(2)),
fmin.Parameters().Vec()(3),
fmin.Parameters().Vec()(4),
fmin.Parameters().Vec()(5),
fmin.Parameters().Vec()(6),
matrix,
fbeamtype );
}
| reco::BeamSpot BSFitter::Fit_dres_z_likelihood | ( | double * | inipar | ) |
Definition at line 873 of file BSFitter.cc.
References align::BeamSpot, j, gen::k, and mathSSE::sqrt().
{
thePDF->SetPDFs("PDFGauss_d_resolution*PDFGauss_z");
thePDF->SetData(fBSvector);
MnUserParameters upar;
upar.Add("X0", inipar[0],0.001);
upar.Add("Y0", inipar[1],0.001);
upar.Add("Z0", inipar[2],0.001);
upar.Add("sigmaZ",inipar[3],0.001);
upar.Add("dxdz",inipar[4],0.001);
upar.Add("dydz",inipar[5],0.001);
upar.Add("BeamWidthX",inipar[6],0.0001);
upar.Add("c0",0.0010,0.0001);
upar.Add("c1",0.0090,0.0001);
// fix beam width
upar.Fix("BeamWidthX");
// number of parameters in fit are 9-1 = 8
MnMigrad migrad(*thePDF, upar);
FunctionMinimum fmin = migrad();
ff_minimum = fmin.Fval();
reco::BeamSpot::CovarianceMatrix matrix;
for (int j = 0 ; j < 6 ; ++j) {
for(int k = j ; k < 6 ; ++k) {
matrix(j,k) = fmin.Error().Matrix()(j,k);
}
}
//std::cout << " fill resolution values" << std::endl;
//std::cout << " matrix size= " << fmin.Error().Matrix().size() << std::endl;
//std::cout << " vec(6)="<< fmin.Parameters().Vec()(6) << std::endl;
//std::cout << " vec(7)="<< fmin.Parameters().Vec()(7) << std::endl;
fresolution_c0 = fmin.Parameters().Vec()(6);
fresolution_c1 = fmin.Parameters().Vec()(7);
fres_c0_err = sqrt( fmin.Error().Matrix()(6,6) );
fres_c1_err = sqrt( fmin.Error().Matrix()(7,7) );
for (int j = 6 ; j < 8 ; ++j) {
for(int k = 6 ; k < 8 ; ++k) {
fres_matrix(j-6,k-6) = fmin.Error().Matrix()(j,k);
}
}
return reco::BeamSpot( reco::BeamSpot::Point(fmin.Parameters().Vec()(0),
fmin.Parameters().Vec()(1),
fmin.Parameters().Vec()(2)),
fmin.Parameters().Vec()(3),
fmin.Parameters().Vec()(4),
fmin.Parameters().Vec()(5),
inipar[6],
matrix,
fbeamtype );
}
| reco::BeamSpot BSFitter::Fit_ited0phi | ( | ) |
Definition at line 391 of file BSFitter.cc.
References reco::BeamSpot::Fake, reco::BeamSpot::setType(), align::Tracker, and Unknown.
Referenced by BeamFitter::runAllFitter().
{
this->d0phi_Init();
edm::LogInfo("BSFitter") << "number of total input tracks: " << fBSvector.size() << std::endl;
reco::BeamSpot theanswer;
if ( (int)fBSvector.size() <= fminNtrks ) {
edm::LogWarning("BSFitter") << "need at least " << fminNtrks << " tracks to run beamline fitter." << std::endl;
fbeamtype = reco::BeamSpot::Fake;
theanswer.setType(fbeamtype);
return theanswer;
}
theanswer = Fit_d0phi(); //get initial ftmp and ftmprow
if ( goodfit ) fnthite++;
//std::cout << "Initial tempanswer (iteration 0): " << theanswer << std::endl;
reco::BeamSpot preanswer = theanswer;
while ( goodfit &&
ftmprow > fconvergence * fBSvector.size() &&
ftmprow > fminNtrks ) {
theanswer = Fit_d0phi();
fd0cut /= 1.5;
fchi2cut /= 1.5;
if ( goodfit &&
ftmprow > fconvergence * fBSvector.size() &&
ftmprow > fminNtrks ) {
preanswer = theanswer;
//std::cout << "Iteration " << fnthite << ": " << preanswer << std::endl;
fnthite++;
}
}
// FIXME: return fit results from previous iteration for both bad fit and for >50% tracks thrown away
//std::cout << "The last iteration, theanswer: " << theanswer << std::endl;
theanswer = preanswer;
//std::cout << "Use previous results from iteration #" << ( fnthite > 0 ? fnthite-1 : 0 ) << std::endl;
//if ( fnthite > 1 ) std::cout << theanswer << std::endl;
edm::LogInfo("BSFitter") << "Total number of successful iterations = " << ( goodfit ? (fnthite+1) : fnthite ) << std::endl;
if (goodfit) {
fbeamtype = reco::BeamSpot::Tracker;
theanswer.setType(fbeamtype);
}
else {
edm::LogWarning("BSFitter") << "Fit doesn't converge!!!" << std::endl;
fbeamtype = reco::BeamSpot::Unknown;
theanswer.setType(fbeamtype);
}
return theanswer;
}
| reco::BeamSpot BSFitter::Fit_z | ( | std::string | type, |
| double * | inipar | ||
| ) |
| reco::BeamSpot BSFitter::Fit_z_chi2 | ( | double * | inipar | ) |
Definition at line 341 of file BSFitter.cc.
References align::BeamSpot.
{
// N.B. this fit is not performed anymore but now
// Z is fitted in the same track set used in the d0-phi fit after
// each iteration
//std::cout << "Fit_z_chi2() called" << std::endl;
// FIXME: include whole tracker z length for the time being
// ==> add protection and z0 cut
h1z = new TH1F("h1z","z distribution",200,-fMaxZ, fMaxZ);
std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
// HERE check size of track vector
for( iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
h1z->Fill( iparam->z0() );
//std::cout<<"z0="<<iparam->z0()<<"; sigZ0="<<iparam->sigz0()<<std::endl;
}
h1z->Fit("gaus","QLM0");
//std::cout << "fitted "<< std::endl;
TF1 *fgaus = h1z->GetFunction("gaus");
//std::cout << "got function" << std::endl;
double fpar[2] = {fgaus->GetParameter(1), fgaus->GetParameter(2) };
//std::cout<<"Debug fpar[2] = (" <<fpar[0]<<","<<fpar[1]<<")"<<std::endl;
reco::BeamSpot::CovarianceMatrix matrix;
// add matrix values.
matrix(2,2) = fgaus->GetParError(1) * fgaus->GetParError(1);
matrix(3,3) = fgaus->GetParError(2) * fgaus->GetParError(2);
//delete h1z;
return reco::BeamSpot( reco::BeamSpot::Point(0.,
0.,
fpar[0]),
fpar[1],
0.,
0.,
0.,
matrix,
fbeamtype );
}
| reco::BeamSpot BSFitter::Fit_z_likelihood | ( | double * | inipar | ) |
Definition at line 272 of file BSFitter.cc.
References align::BeamSpot, j, gen::k, and Gflash::par.
{
//std::cout << "Fit_z(double *) called" << std::endl;
//std::cout << "inipar[0]= " << inipar[0] << std::endl;
//std::cout << "inipar[1]= " << inipar[1] << std::endl;
std::vector<double> par(2,0);
std::vector<double> err(2,0);
par.push_back(0.0);
par.push_back(7.0);
err.push_back(0.0001);
err.push_back(0.0001);
//par[0] = 0.0; err[0] = 0.0;
//par[1] = 7.0; err[1] = 0.0;
thePDF->SetPDFs("PDFGauss_z");
thePDF->SetData(fBSvector);
//std::cout << "data loaded"<< std::endl;
//FunctionMinimum fmin = theFitter->Minimize(*theGausszFcn, par, err, 1, 500, 0.1);
MnUserParameters upar;
upar.Add("X0", 0.,0.);
upar.Add("Y0", 0.,0.);
upar.Add("Z0", inipar[0],0.001);
upar.Add("sigmaZ",inipar[1],0.001);
MnMigrad migrad(*thePDF, upar);
FunctionMinimum fmin = migrad();
ff_minimum = fmin.Fval();
//std::cout << " eval= " << ff_minimum
// << "/n params[0]= " << fmin.Parameters().Vec()(0) << std::endl;
/*
TMinuit *gmMinuit = new TMinuit(2);
//gmMinuit->SetFCN(z_fcn);
gmMinuit->SetFCN(myFitz_fcn);
int ierflg = 0;
double step[2] = {0.001,0.001};
for (int i = 0; i<2; i++) {
gmMinuit->mnparm(i,fpar_name[i].c_str(),inipar[i],step[i],0,0,ierflg);
}
gmMinuit->Migrad();
*/
reco::BeamSpot::CovarianceMatrix matrix;
for (int j = 2 ; j < 4 ; ++j) {
for(int k = j ; k < 4 ; ++k) {
matrix(j,k) = fmin.Error().Matrix()(j,k);
}
}
return reco::BeamSpot( reco::BeamSpot::Point(0.,
0.,
fmin.Parameters().Vec()(2)),
fmin.Parameters().Vec()(3),
0.,
0.,
0.,
matrix,
fbeamtype );
}
| int BSFitter::GetAcceptedTrks | ( | ) | [inline] |
| std::vector< BSTrkParameters > BSFitter::GetData | ( | ) | [inline] |
| double BSFitter::GetMinimum | ( | ) | [inline] |
| reco::BeamSpot::ResCovMatrix BSFitter::GetResMatrix | ( | ) | [inline] |
| double BSFitter::GetResPar0 | ( | ) | [inline] |
| double BSFitter::GetResPar0Err | ( | ) | [inline] |
| double BSFitter::GetResPar1 | ( | ) | [inline] |
| double BSFitter::GetResPar1Err | ( | ) | [inline] |
| TH1F* BSFitter::GetVzHisto | ( | ) | [inline] |
Definition at line 106 of file BSFitter.h.
References h1z.
Referenced by BeamFitter::runFitterNoTxt().
{ return h1z; }
| double BSFitter::scanPDF | ( | double * | init_pars, |
| int & | tracksFailed, | ||
| int | option | ||
| ) |
Definition at line 679 of file BSFitter.cc.
References abs, funct::cos(), ExpressReco_HICollisions_FallBack::e, funct::exp(), root::function(), funct::log(), L1TEmulatorMonitor_cff::p, Pi, funct::sin(), and mathSSE::sqrt().
{
if(option==1)init_pars[6]=0.0005; //starting value for any given configuration
//local vairables with initial values
double fsqrt2pi=0.0;
double d_sig=0.0;
double d_dprime=0.0;
double d_result=0.0;
double z_sig=0.0;
double z_result=0.0;
double function=0.0;
double tot_pdf=0.0;
double last_minvalue=1.0e+10;
double init_bw=-99.99;
int iters=0;
//used to remove tracks if far away from bs by this
double DeltadCut=0.1000;
if(init_pars[6]<0.0200){DeltadCut=0.0900; } //worked for high 2.36TeV
if(init_pars[6]<0.0100){DeltadCut=0.0700;} //just a guesss for 7 TeV but one should scan for actual values
std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
if(option==1)iters=500;
if(option==2)iters=1;
for(int p=0;p<iters;p++){
if(iters==500)init_pars[6]+=0.0002;
tracksfixed=0;
for( iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam)
{
fsqrt2pi = sqrt(2.* TMath::Pi());
d_sig = sqrt(init_pars[6]*init_pars[6] + (iparam->sigd0())*(iparam->sigd0()));
d_dprime = iparam->d0() - ( ( (init_pars[0] + iparam->z0()*(init_pars[4]))*sin(iparam->phi0()) )
- ( (init_pars[1] + iparam->z0()*(init_pars[5]))*cos(iparam->phi0()) ) );
//***Remove tracks before the fit which gives low pdf values to blow up the pdf
if(std::abs(d_dprime)<DeltadCut && option==2){ fBSvectorBW.push_back(*iparam);}
d_result = (exp(-(d_dprime*d_dprime)/(2.0*d_sig*d_sig)))/(d_sig*fsqrt2pi);
z_sig = sqrt(iparam->sigz0() * iparam->sigz0() + init_pars[3]*init_pars[3]);
z_result = (exp(-((iparam->z0() - init_pars[2])*(iparam->z0() - init_pars[2]))/(2.0*z_sig*z_sig)))/(z_sig*fsqrt2pi);
tot_pdf=z_result*d_result;
//for those trcks which gives problems due to very tiny pdf_d values.
//Update: This protection will NOT be used with the dprime cut above but still kept here to get
// the intial value of beam width reasonably
//A warning will appear if there were any tracks with < 10^-5 for pdf_d so that (d-dprime) cut can be lowered
if(d_result < 1.0e-05){ tot_pdf=z_result*1.0e-05;
//if(option==2)std::cout<<"last Iter d-d' = "<<(std::abs(d_dprime))<<std::endl;
tracksfixed++; }
function = function + log(tot_pdf);
tot_pdf=0.0;
}//loop over tracks
function= -2.0*function;
if(function<last_minvalue){init_bw=init_pars[6];
last_minvalue=function; }
function=0.0;
}//loop over beam width
if(init_bw>0) {
init_bw=init_bw+(0.20*init_bw); //start with 20 % more
}
else{
if(option==1){
edm::LogWarning("BSFitter")
<<"scanPDF:====>>>> WARNING***: The initial guess value of Beam width is negative!!!!!!"<<std::endl
<<"scanPDF:====>>>> Assigning beam width a starting value of "<<init_bw<<" cm"<<std::endl;
init_bw=0.0200;
}
}
return init_bw;
}
| void BSFitter::SetChi2Cut_d0phi | ( | double | chi2cut | ) |
Definition at line 632 of file BSFitter.cc.
{
fapplychi2cut = true;
//fBSforCuts = BSfitted;
fchi2cut = chi2cut;
}
| void BSFitter::SetConvergence | ( | double | val | ) | [inline] |
Definition at line 63 of file BSFitter.h.
References fconvergence.
Referenced by BeamFitter::runBeamWidthFitter(), and BeamFitter::runFitterNoTxt().
{ fconvergence = val; }
| void BSFitter::Setd0Cut_d0phi | ( | double | d0cut | ) |
Definition at line 623 of file BSFitter.cc.
Referenced by BeamFitter::runAllFitter().
{
fapplyd0cut = true;
//fBSforCuts = BSfitted;
fd0cut = d0cut;
}
| void BSFitter::SetFitType | ( | std::string | type | ) | [inline] |
Definition at line 42 of file BSFitter.h.
References ffit_type.
Referenced by BeamFitter::runAllFitter(), and BeamFitter::runBeamWidthFitter().
| void BSFitter::SetFitVariable | ( | std::string | name | ) | [inline] |
Definition at line 46 of file BSFitter.h.
References ffit_variable, and AlCaRecoCosmics_cfg::name.
Referenced by BeamFitter::runAllFitter(), and BeamFitter::runBeamWidthFitter().
{
ffit_variable = name;
}
| void BSFitter::SetInputBeamWidth | ( | double | val | ) | [inline] |
Definition at line 67 of file BSFitter.h.
References finputBeamWidth.
Referenced by BeamFitter::runBeamWidthFitter(), and BeamFitter::runFitterNoTxt().
{ finputBeamWidth = val; }
| void BSFitter::SetMaximumZ | ( | double | z | ) | [inline] |
Definition at line 62 of file BSFitter.h.
Referenced by BeamFitter::runBeamWidthFitter(), and BeamFitter::runFitterNoTxt().
| void BSFitter::SetMinimumNTrks | ( | int | n | ) | [inline] |
Definition at line 64 of file BSFitter.h.
Referenced by BeamFitter::runBeamWidthFitter(), and BeamFitter::runFitterNoTxt().
bool BSFitter::fapplychi2cut [private] |
Definition at line 137 of file BSFitter.h.
bool BSFitter::fapplyd0cut [private] |
Definition at line 136 of file BSFitter.h.
reco::BeamSpot::BeamType BSFitter::fbeamtype [private] |
Definition at line 114 of file BSFitter.h.
std::vector< BSTrkParameters > BSFitter::fBSvector [private] |
Definition at line 126 of file BSFitter.h.
Referenced by GetData().
std::vector< BSTrkParameters > BSFitter::fBSvectorBW [private] |
Definition at line 127 of file BSFitter.h.
double BSFitter::fchi2cut [private] |
Definition at line 139 of file BSFitter.h.
double BSFitter::fconvergence [private] |
Definition at line 144 of file BSFitter.h.
Referenced by SetConvergence().
double BSFitter::fd0cut [private] |
Definition at line 138 of file BSFitter.h.
const int BSFitter::fdim = 7 [static, private] |
Definition at line 120 of file BSFitter.h.
double BSFitter::ff_minimum [private] |
Definition at line 118 of file BSFitter.h.
Referenced by GetMinimum().
std::string BSFitter::ffit_type [private] |
Definition at line 115 of file BSFitter.h.
Referenced by SetFitType().
std::string BSFitter::ffit_variable [private] |
Definition at line 116 of file BSFitter.h.
Referenced by SetFitVariable().
double BSFitter::finputBeamWidth [private] |
Definition at line 146 of file BSFitter.h.
Referenced by SetInputBeamWidth().
double BSFitter::fMaxZ [private] |
Definition at line 143 of file BSFitter.h.
Referenced by SetMaximumZ().
int BSFitter::fminNtrks [private] |
Definition at line 145 of file BSFitter.h.
Referenced by SetMinimumNTrks().
int BSFitter::fnthite [private] |
Definition at line 141 of file BSFitter.h.
Referenced by d0phi_Init().
std::string BSFitter::fpar_name[fdim] [private] |
Definition at line 122 of file BSFitter.h.
double BSFitter::fres_c0_err [private] |
Definition at line 131 of file BSFitter.h.
Referenced by GetResPar0Err().
double BSFitter::fres_c1_err [private] |
Definition at line 132 of file BSFitter.h.
Referenced by GetResPar1Err().
Definition at line 133 of file BSFitter.h.
Referenced by GetResMatrix().
double BSFitter::fresolution_c0 [private] |
Definition at line 129 of file BSFitter.h.
Referenced by GetResPar0().
double BSFitter::fresolution_c1 [private] |
Definition at line 130 of file BSFitter.h.
Referenced by GetResPar1().
Double_t BSFitter::fsqrt2pi [private] |
Definition at line 124 of file BSFitter.h.
TMatrixD BSFitter::ftmp [private] |
Definition at line 135 of file BSFitter.h.
Referenced by d0phi_Init().
int BSFitter::ftmprow [private] |
Definition at line 140 of file BSFitter.h.
Referenced by d0phi_Init(), and GetAcceptedTrks().
bool BSFitter::goodfit [private] |
Definition at line 142 of file BSFitter.h.
Referenced by d0phi_Init().
TH1F* BSFitter::h1z [private] |
Definition at line 147 of file BSFitter.h.
Referenced by GetVzHisto().
ROOT::Minuit2::ModularFunctionMinimizer* BSFitter::theFitter [private] |
Definition at line 110 of file BSFitter.h.
BSpdfsFcn* BSFitter::thePDF [private] |
Definition at line 112 of file BSFitter.h.
1.7.3